annotate toolboxes/FullBNT-1.0.7/bnt/potentials/Tables/mult_by_table.c @ 0:e9a9cd732c1e tip

first hg version after svn
author wolffd
date Tue, 10 Feb 2015 15:05:51 +0000
parents
children
rev   line source
wolffd@0 1 /* mult_by_table.c ../potential/tables */
wolffd@0 2
wolffd@0 3
wolffd@0 4 /******************************************/
wolffd@0 5 /* 6 input & 1 output */
wolffd@0 6 /* Big table [0] */
wolffd@0 7 /* Big domain [1] */
wolffd@0 8 /* big sizes [2] */
wolffd@0 9 /* Small table [3] */
wolffd@0 10 /* small domain [4] */
wolffd@0 11 /* small sizes [5] */
wolffd@0 12 /* */
wolffd@0 13 /* New big table[0] */
wolffd@0 14 /******************************************/
wolffd@0 15
wolffd@0 16 #include "mex.h"
wolffd@0 17
wolffd@0 18 void mexFunction(int nlhs, mxArray *plhs[], int nrhs, const mxArray *prhs[]){
wolffd@0 19 int i, j, count, NB, NS, siz_b, siz_s, ndim, temp;
wolffd@0 20 int *mask, *sx, *sy, *cpsy, *subs, *s, *cpsy2;
wolffd@0 21 double *pbDomain, *psDomain, *sp, *zp, *bs;
wolffd@0 22
wolffd@0 23 plhs[0] = mxDuplicateArray(prhs[0]);
wolffd@0 24 zp = mxGetPr(plhs[0]);
wolffd@0 25
wolffd@0 26 siz_b = mxGetNumberOfElements(prhs[1]);
wolffd@0 27 siz_s = mxGetNumberOfElements(prhs[4]);
wolffd@0 28 pbDomain = mxGetPr(prhs[1]);
wolffd@0 29 psDomain = mxGetPr(prhs[4]);
wolffd@0 30
wolffd@0 31 NB = mxGetNumberOfElements(prhs[0]);
wolffd@0 32 NS = mxGetNumberOfElements(prhs[3]);
wolffd@0 33 sp = mxGetPr(prhs[3]);
wolffd@0 34
wolffd@0 35 bs = mxGetPr(prhs[2]);
wolffd@0 36
wolffd@0 37 if(NS == 1){
wolffd@0 38 for(i=0; i<NB; i++){
wolffd@0 39 zp[i] *= *sp;
wolffd@0 40 }
wolffd@0 41 return;
wolffd@0 42 }
wolffd@0 43
wolffd@0 44 if(NS == NB){
wolffd@0 45 for(i=0; i<NB; i++){
wolffd@0 46 zp[i] *= sp[i];
wolffd@0 47 }
wolffd@0 48 return;
wolffd@0 49 }
wolffd@0 50
wolffd@0 51 mask = malloc(siz_s * sizeof(int));
wolffd@0 52 count = 0;
wolffd@0 53 for(i=0; i<siz_s; i++){
wolffd@0 54 for(j=0; j<siz_b; j++){
wolffd@0 55 if(psDomain[i] == pbDomain[j]){
wolffd@0 56 mask[count] = j;
wolffd@0 57 count++;
wolffd@0 58 break;
wolffd@0 59 }
wolffd@0 60 }
wolffd@0 61 }
wolffd@0 62
wolffd@0 63 ndim = siz_b;
wolffd@0 64 sx = (int *)malloc(sizeof(int)*ndim);
wolffd@0 65 sy = (int *)malloc(sizeof(int)*ndim);
wolffd@0 66 for(i=0; i<ndim; i++){
wolffd@0 67 sx[i] = (int)bs[i];
wolffd@0 68 sy[i] = 1;
wolffd@0 69 }
wolffd@0 70 for(i=0; i<count; i++){
wolffd@0 71 temp = mask[i];
wolffd@0 72 sy[temp] = sx[temp];
wolffd@0 73 }
wolffd@0 74
wolffd@0 75 s = (int *)malloc(sizeof(int)*ndim);
wolffd@0 76 *(cpsy = (int *)malloc(sizeof(int)*ndim)) = 1;
wolffd@0 77 subs = (int *)malloc(sizeof(int)*ndim);
wolffd@0 78 cpsy2 = (int *)malloc(sizeof(int)*ndim);
wolffd@0 79 for(i = 0; i < ndim; i++){
wolffd@0 80 subs[i] = 0;
wolffd@0 81 s[i] = sx[i] - 1;
wolffd@0 82 }
wolffd@0 83
wolffd@0 84 for(i = 0; i < ndim-1; i++){
wolffd@0 85 cpsy[i+1] = cpsy[i]*sy[i]--;
wolffd@0 86 cpsy2[i] = cpsy[i]*sy[i];
wolffd@0 87 }
wolffd@0 88 cpsy2[ndim-1] = cpsy[ndim-1]*(--sy[ndim-1]);
wolffd@0 89
wolffd@0 90 for(j=0; j<NB; j++){
wolffd@0 91 *zp++ *= *sp;
wolffd@0 92 for(i = 0; i < ndim; i++){
wolffd@0 93 if(subs[i] == s[i]){
wolffd@0 94 subs[i] = 0;
wolffd@0 95 if(sy[i])
wolffd@0 96 sp -= cpsy2[i];
wolffd@0 97 }
wolffd@0 98 else{
wolffd@0 99 subs[i]++;
wolffd@0 100 if(sy[i])
wolffd@0 101 sp += cpsy[i];
wolffd@0 102 break;
wolffd@0 103 }
wolffd@0 104 }
wolffd@0 105 }
wolffd@0 106 free(sx);
wolffd@0 107 free(sy);
wolffd@0 108 free(s);
wolffd@0 109 free(cpsy);
wolffd@0 110 free(subs);
wolffd@0 111 free(cpsy2);
wolffd@0 112 free(mask);
wolffd@0 113 }
wolffd@0 114