Mercurial > hg > camir-aes2014

annotate toolboxes/FullBNT-1.0.7/bnt/potentials/Tables/divide_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 /* divide_by_sparse_table.c ../potential/tables*/
wolffd@0 2
wolffd@0 3 /******************************************/
wolffd@0 4 /* 6 input & 1 output */
wolffd@0 5 /* Big table [0] */
wolffd@0 6 /* Big domain [1] */
wolffd@0 7 /* big sizes [2] */
wolffd@0 8 /* Small table [3] */
wolffd@0 9 /* small domain [4] */
wolffd@0 10 /* small sizes [5] */
wolffd@0 11 /* */
wolffd@0 12 /* New big table[0] */
wolffd@0 13 /******************************************/
wolffd@0 14
wolffd@0 15 #include <math.h>
wolffd@0 16 #include <stdlib.h>
wolffd@0 17 #include "mex.h"
wolffd@0 18
wolffd@0 19 int compare(const void* src1, const void* src2){
wolffd@0 20 int i1 = *(int*)src1 ;
wolffd@0 21 int i2 = *(int*)src2 ;
wolffd@0 22 return i1-i2 ;
wolffd@0 23 }
wolffd@0 24
wolffd@0 25 void ind_subv(int index, const int *cumprod, int n, int *bsubv){
wolffd@0 26 int i;
wolffd@0 27
wolffd@0 28 for (i = n-1; i >= 0; i--) {
wolffd@0 29 bsubv[i] = ((int)floor(index / cumprod[i]));
wolffd@0 30 index = index % cumprod[i];
wolffd@0 31 }
wolffd@0 32 }
wolffd@0 33
wolffd@0 34 int subv_ind(const int n, const int *cumprod, const int *subv){
wolffd@0 35 int i, index=0;
wolffd@0 36
wolffd@0 37 for(i=0; i<n; i++){
wolffd@0 38 index += subv[i] * cumprod[i];
wolffd@0 39 }
wolffd@0 40 return index;
wolffd@0 41 }
wolffd@0 42
wolffd@0 43 void mexFunction(int nlhs, mxArray *plhs[], int nrhs, const mxArray *prhs[]){
wolffd@0 44 int i, j, count, bdim, sdim, NB, NZB, NZS, position, bindex, sindex;
wolffd@0 45 int *mask, *result, *bir, *sir, *bjc, *sjc, *bCumprod, *sCumprod, *bsubv, *ssubv;
wolffd@0 46 double *pbDomain, *psDomain, *pbSize, *psSize, *bpr, *spr, value;
wolffd@0 47
wolffd@0 48 plhs[0] = mxDuplicateArray(prhs[0]);
wolffd@0 49 pbDomain = mxGetPr(prhs[1]);
wolffd@0 50 bdim = mxGetNumberOfElements(prhs[1]);
wolffd@0 51 psDomain = mxGetPr(prhs[4]);
wolffd@0 52 sdim = mxGetNumberOfElements(prhs[4]);
wolffd@0 53
wolffd@0 54 pbSize = mxGetPr(prhs[2]);
wolffd@0 55 psSize = mxGetPr(prhs[5]);
wolffd@0 56
wolffd@0 57 NB = 1;
wolffd@0 58 for(i=0; i<bdim; i++){
wolffd@0 59 NB *= (int)pbSize[i];
wolffd@0 60 }
wolffd@0 61
wolffd@0 62 bpr = mxGetPr(plhs[0]);
wolffd@0 63 bir = mxGetIr(plhs[0]);
wolffd@0 64 bjc = mxGetJc(plhs[0]);
wolffd@0 65 NZB = bjc[1];
wolffd@0 66
wolffd@0 67 spr = mxGetPr(prhs[3]);
wolffd@0 68 sir = mxGetIr(prhs[3]);
wolffd@0 69 sjc = mxGetJc(prhs[3]);
wolffd@0 70 NZS = sjc[1];
wolffd@0 71
wolffd@0 72 if(sdim == 0){
wolffd@0 73 value = *spr;
wolffd@0 74 if(value == 0)value = 1;
wolffd@0 75 for(i=0; i<NZB; i++){
wolffd@0 76 bpr[i] /= value;
wolffd@0 77 }
wolffd@0 78 return;
wolffd@0 79 }
wolffd@0 80
wolffd@0 81 mask = malloc(sdim * sizeof(int));
wolffd@0 82 bCumprod = malloc(bdim * sizeof(int));
wolffd@0 83 sCumprod = malloc(sdim * sizeof(int));
wolffd@0 84 bsubv = malloc(bdim * sizeof(int));
wolffd@0 85 ssubv = malloc(sdim * sizeof(int));
wolffd@0 86
wolffd@0 87 count = 0;
wolffd@0 88 for(i=0; i<sdim; i++){
wolffd@0 89 for(j=0; j<bdim; j++){
wolffd@0 90 if(psDomain[i] == pbDomain[j]){
wolffd@0 91 mask[count] = j;
wolffd@0 92 count++;
wolffd@0 93 break;
wolffd@0 94 }
wolffd@0 95 }
wolffd@0 96 }
wolffd@0 97
wolffd@0 98 bCumprod[0] = 1;
wolffd@0 99 for(i=0; i<bdim-1; i++){
wolffd@0 100 bCumprod[i+1] = bCumprod[i] * (int)pbSize[i];
wolffd@0 101 }
wolffd@0 102 sCumprod[0] = 1;
wolffd@0 103 for(i=0; i<sdim-1; i++){
wolffd@0 104 sCumprod[i+1] = sCumprod[i] * (int)psSize[i];
wolffd@0 105 }
wolffd@0 106
wolffd@0 107 for(i=0; i<NZB; i++){
wolffd@0 108 bindex = bir[i];
wolffd@0 109 ind_subv(bindex, bCumprod, bdim, bsubv);
wolffd@0 110 for(j=0; j<sdim; j++){
wolffd@0 111 ssubv[j] = bsubv[mask[j]];
wolffd@0 112 }
wolffd@0 113 sindex = subv_ind(sdim, sCumprod, ssubv);
wolffd@0 114 result = (int *) bsearch(&sindex, sir, NZS, sizeof(int), compare);
wolffd@0 115 if(result){
wolffd@0 116 position = result - sir;
wolffd@0 117 bpr[i] /= spr[position];
wolffd@0 118 }
wolffd@0 119 }
wolffd@0 120
wolffd@0 121 free(mask);
wolffd@0 122 free(bCumprod);
wolffd@0 123 free(sCumprod);
wolffd@0 124 free(bsubv);
wolffd@0 125 free(ssubv);
wolffd@0 126 }