Mercurial > hg > camir-aes2014

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

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
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_sparse_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 <math.h>
wolffd@0 17 #include <stdlib.h>
wolffd@0 18 #include "mex.h"
wolffd@0 19
wolffd@0 20 int compare(const void* src1, const void* src2){
wolffd@0 21 int i1 = *(int*)src1 ;
wolffd@0 22 int i2 = *(int*)src2 ;
wolffd@0 23 return i1-i2 ;
wolffd@0 24 }
wolffd@0 25
wolffd@0 26 void ind_subv(int index, const int *cumprod, int n, int *bsubv){
wolffd@0 27 int i;
wolffd@0 28
wolffd@0 29 for (i = n-1; i >= 0; i--) {
wolffd@0 30 bsubv[i] = ((int)floor(index / cumprod[i]));
wolffd@0 31 index = index % cumprod[i];
wolffd@0 32 }
wolffd@0 33 }
wolffd@0 34
wolffd@0 35 int subv_ind(const int n, const int *cumprod, const int *subv){
wolffd@0 36 int i, index=0;
wolffd@0 37
wolffd@0 38 for(i=0; i<n; i++){
wolffd@0 39 index += subv[i] * cumprod[i];
wolffd@0 40 }
wolffd@0 41 return index;
wolffd@0 42 }
wolffd@0 43
wolffd@0 44 void reset_nzmax(mxArray *spArray, const int old_nzmax, const int new_nzmax){
wolffd@0 45 double *ptr;
wolffd@0 46 void *newptr;
wolffd@0 47 int *ir, *jc;
wolffd@0 48 int nbytes;
wolffd@0 49
wolffd@0 50 if(new_nzmax == old_nzmax) return;
wolffd@0 51 nbytes = new_nzmax * sizeof(*ptr);
wolffd@0 52 ptr = mxGetPr(spArray);
wolffd@0 53 newptr = mxRealloc(ptr, nbytes);
wolffd@0 54 mxSetPr(spArray, newptr);
wolffd@0 55 nbytes = new_nzmax * sizeof(*ir);
wolffd@0 56 ir = mxGetIr(spArray);
wolffd@0 57 newptr = mxRealloc(ir, nbytes);
wolffd@0 58 mxSetIr(spArray, newptr);
wolffd@0 59 jc = mxGetJc(spArray);
wolffd@0 60 jc[0] = 0;
wolffd@0 61 jc[1] = new_nzmax;
wolffd@0 62 mxSetNzmax(spArray, new_nzmax);
wolffd@0 63 }
wolffd@0 64
wolffd@0 65 void mexFunction(int nlhs, mxArray *plhs[], int nrhs, const mxArray *prhs[]){
wolffd@0 66 int i, j, count, bdim, sdim, NB, NZB, NZS, position, bindex, sindex, nzCounts=0;
wolffd@0 67 int *mask, *result, *bir, *sir, *rir, *bjc, *sjc, *rjc, *bCumprod, *sCumprod, *bsubv, *ssubv;
wolffd@0 68 double *pbDomain, *psDomain, *pbSize, *psSize, *bpr, *spr, *rpr;
wolffd@0 69
wolffd@0 70 pbDomain = mxGetPr(prhs[1]);
wolffd@0 71 bdim = mxGetNumberOfElements(prhs[1]);
wolffd@0 72 psDomain = mxGetPr(prhs[4]);
wolffd@0 73 sdim = mxGetNumberOfElements(prhs[4]);
wolffd@0 74
wolffd@0 75 pbSize = mxGetPr(prhs[2]);
wolffd@0 76 psSize = mxGetPr(prhs[5]);
wolffd@0 77
wolffd@0 78 NB = 1;
wolffd@0 79 for(i=0; i<bdim; i++){
wolffd@0 80 NB *= (int)pbSize[i];
wolffd@0 81 }
wolffd@0 82
wolffd@0 83 bpr = mxGetPr(prhs[0]);
wolffd@0 84 bir = mxGetIr(prhs[0]);
wolffd@0 85 bjc = mxGetJc(prhs[0]);
wolffd@0 86 NZB = bjc[1];
wolffd@0 87
wolffd@0 88 spr = mxGetPr(prhs[3]);
wolffd@0 89 sir = mxGetIr(prhs[3]);
wolffd@0 90 sjc = mxGetJc(prhs[3]);
wolffd@0 91 NZS = sjc[1];
wolffd@0 92
wolffd@0 93 plhs[0] = mxDuplicateArray(prhs[0]);
wolffd@0 94 rpr = mxGetPr(plhs[0]);
wolffd@0 95 rir = mxGetIr(plhs[0]);
wolffd@0 96 rjc = mxGetJc(plhs[0]);
wolffd@0 97 rjc[0] = 0;
wolffd@0 98 rjc[1] = NZB;
wolffd@0 99
wolffd@0 100 if(sdim == 0){
wolffd@0 101 for(i=0; i<NZB; i++){
wolffd@0 102 rpr[i] *= *spr;
wolffd@0 103 }
wolffd@0 104 return;
wolffd@0 105 }
wolffd@0 106
wolffd@0 107 mask = malloc(sdim * sizeof(int));
wolffd@0 108 bCumprod = malloc(bdim * sizeof(int));
wolffd@0 109 sCumprod = malloc(sdim * sizeof(int));
wolffd@0 110 bsubv = malloc(bdim * sizeof(int));
wolffd@0 111 ssubv = malloc(sdim * sizeof(int));
wolffd@0 112
wolffd@0 113 count = 0;
wolffd@0 114 for(i=0; i<sdim; i++){
wolffd@0 115 for(j=0; j<bdim; j++){
wolffd@0 116 if(psDomain[i] == pbDomain[j]){
wolffd@0 117 mask[count] = j;
wolffd@0 118 count++;
wolffd@0 119 break;
wolffd@0 120 }
wolffd@0 121 }
wolffd@0 122 }
wolffd@0 123
wolffd@0 124 bCumprod[0] = 1;
wolffd@0 125 for(i=0; i<bdim-1; i++){
wolffd@0 126 bCumprod[i+1] = bCumprod[i] * (int)pbSize[i];
wolffd@0 127 }
wolffd@0 128 sCumprod[0] = 1;
wolffd@0 129 for(i=0; i<sdim-1; i++){
wolffd@0 130 sCumprod[i+1] = sCumprod[i] * (int)psSize[i];
wolffd@0 131 }
wolffd@0 132
wolffd@0 133 for(i=0; i<NZB; i++){
wolffd@0 134 bindex = bir[i];
wolffd@0 135 ind_subv(bindex, bCumprod, bdim, bsubv);
wolffd@0 136 for(j=0; j<sdim; j++){
wolffd@0 137 ssubv[j] = bsubv[mask[j]];
wolffd@0 138 }
wolffd@0 139 sindex = subv_ind(sdim, sCumprod, ssubv);
wolffd@0 140 result = (int *) bsearch(&sindex, sir, NZS, sizeof(int), compare);
wolffd@0 141 if(result){
wolffd@0 142 position = result - sir;
wolffd@0 143 rpr[nzCounts] = bpr[i] * spr[position];
wolffd@0 144 rir[nzCounts] = bindex;
wolffd@0 145 nzCounts++;
wolffd@0 146 }
wolffd@0 147 }
wolffd@0 148
wolffd@0 149 reset_nzmax(plhs[0], NZB, nzCounts);
wolffd@0 150 free(mask);
wolffd@0 151 free(bCumprod);
wolffd@0 152 free(sCumprod);
wolffd@0 153 free(bsubv);
wolffd@0 154 free(ssubv);
wolffd@0 155 }