Mercurial > hg > camir-aes2014

annotate toolboxes/FullBNT-1.0.7/KPMtools/repmatC.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 /*
wolffd@0 2 mex -c mexutil.c
wolffd@0 3 mex repmat.c mexutil.obj
wolffd@0 4 to check for warnings:
wolffd@0 5 gcc -Wall -I/cygdrive/c/MATLAB6p1/extern/include -c repmat.c
wolffd@0 6 */
wolffd@0 7 #include "mexutil.h"
wolffd@0 8 #include <string.h>
wolffd@0 9
wolffd@0 10 /* repeat a block of memory rep times */
wolffd@0 11 void memrep(char *dest, size_t chunk, int rep)
wolffd@0 12 {
wolffd@0 13 #if 0
wolffd@0 14 /* slow way */
wolffd@0 15 int i;
wolffd@0 16 char *p = dest;
wolffd@0 17 for(i=1;i<rep;i++) {
wolffd@0 18 p += chunk;
wolffd@0 19 memcpy(p, dest, chunk);
wolffd@0 20 }
wolffd@0 21 #else
wolffd@0 22 /* fast way */
wolffd@0 23 if(rep == 1) return;
wolffd@0 24 memcpy(dest + chunk, dest, chunk);
wolffd@0 25 if(rep & 1) {
wolffd@0 26 dest += chunk;
wolffd@0 27 memcpy(dest + chunk, dest, chunk);
wolffd@0 28 }
wolffd@0 29 /* now repeat using a block twice as big */
wolffd@0 30 memrep(dest, chunk<<1, rep>>1);
wolffd@0 31 #endif
wolffd@0 32 }
wolffd@0 33
wolffd@0 34 void repmat(char *dest, const char *src, int ndim, int *destdimsize,
wolffd@0 35 int *dimsize, const int *dims, int *rep)
wolffd@0 36 {
wolffd@0 37 int d = ndim-1;
wolffd@0 38 int i, chunk;
wolffd@0 39 /* copy the first repetition into dest */
wolffd@0 40 if(d == 0) {
wolffd@0 41 chunk = dimsize[0];
wolffd@0 42 memcpy(dest,src,chunk);
wolffd@0 43 }
wolffd@0 44 else {
wolffd@0 45 /* recursively repeat each slice of src */
wolffd@0 46 for(i=0;i<dims[d];i++) {
wolffd@0 47 repmat(dest + i*destdimsize[d-1], src + i*dimsize[d-1],
wolffd@0 48 ndim-1, destdimsize, dimsize, dims, rep);
wolffd@0 49 }
wolffd@0 50 chunk = destdimsize[d-1]*dims[d];
wolffd@0 51 }
wolffd@0 52 /* copy the result rep-1 times */
wolffd@0 53 memrep(dest,chunk,rep[d]);
wolffd@0 54 }
wolffd@0 55
wolffd@0 56 void mexFunction(int nlhs, mxArray *plhs[],
wolffd@0 57 int nrhs, const mxArray *prhs[])
wolffd@0 58 {
wolffd@0 59 const mxArray *srcmat;
wolffd@0 60 int ndim, *dimsize, eltsize;
wolffd@0 61 const int *dims;
wolffd@0 62 int ndimdest, *destdims, *destdimsize;
wolffd@0 63 char *src, *dest;
wolffd@0 64 int *rep;
wolffd@0 65 int i,nrep;
wolffd@0 66 int extra_rep = 1;
wolffd@0 67 int empty;
wolffd@0 68
wolffd@0 69 if(nrhs < 2) mexErrMsgTxt("Usage: xrepmat(A, [N M ...])");
wolffd@0 70 srcmat = prhs[0];
wolffd@0 71 if(mxIsSparse(srcmat)) {
wolffd@0 72 mexErrMsgTxt("Sorry, can't handle sparse matrices yet.");
wolffd@0 73 }
wolffd@0 74 if(mxIsCell(srcmat)) {
wolffd@0 75 mexErrMsgTxt("Sorry, can't handle cell arrays yet.");
wolffd@0 76 }
wolffd@0 77 ndim = mxGetNumberOfDimensions(srcmat);
wolffd@0 78 dims = mxGetDimensions(srcmat);
wolffd@0 79 eltsize = mxGetElementSize(srcmat);
wolffd@0 80
wolffd@0 81 /* compute dimension sizes */
wolffd@0 82 dimsize = mxCalloc(ndim, sizeof(int));
wolffd@0 83 dimsize[0] = eltsize*dims[0];
wolffd@0 84 for(i=1;i<ndim;i++) dimsize[i] = dimsize[i-1]*dims[i];
wolffd@0 85
wolffd@0 86 /* determine repetition vector */
wolffd@0 87 ndimdest = ndim;
wolffd@0 88 if(nrhs == 2) {
wolffd@0 89 nrep = mxGetN(prhs[1]);
wolffd@0 90 if(nrep > ndimdest) ndimdest = nrep;
wolffd@0 91 rep = mxCalloc(ndimdest, sizeof(int));
wolffd@0 92 for(i=0;i<nrep;i++) {
wolffd@0 93 double repv = mxGetPr(prhs[1])[i];
wolffd@0 94 rep[i] = (int)repv;
wolffd@0 95 }
wolffd@0 96 if(nrep == 1) {
wolffd@0 97 /* special behavior */
wolffd@0 98 nrep = 2;
wolffd@0 99 rep[1] = rep[0];
wolffd@0 100 }
wolffd@0 101 }
wolffd@0 102 else {
wolffd@0 103 nrep = nrhs-1;
wolffd@0 104 if(nrep > ndimdest) ndimdest = nrep;
wolffd@0 105 rep = mxCalloc(ndimdest, sizeof(int));
wolffd@0 106 for(i=0;i<nrep;i++) {
wolffd@0 107 rep[i] = (int)*mxGetPr(prhs[i+1]);
wolffd@0 108 }
wolffd@0 109 }
wolffd@0 110 for(i=nrep;i<ndimdest;i++) rep[i] = 1;
wolffd@0 111
wolffd@0 112 /* compute output size */
wolffd@0 113 destdims = mxCalloc(ndimdest, sizeof(int));
wolffd@0 114 for(i=0;i<ndim;i++) destdims[i] = dims[i]*rep[i];
wolffd@0 115 for(;i<ndimdest;i++) {
wolffd@0 116 destdims[i] = rep[i];
wolffd@0 117 extra_rep *= rep[i];
wolffd@0 118 }
wolffd@0 119 destdimsize = mxCalloc(ndim, sizeof(int));
wolffd@0 120 destdimsize[0] = eltsize*destdims[0];
wolffd@0 121 for(i=1;i<ndim;i++) destdimsize[i] = destdimsize[i-1]*destdims[i];
wolffd@0 122
wolffd@0 123
wolffd@0 124 /* for speed, array should be uninitialized */
wolffd@0 125 plhs[0] = mxCreateNumericArray(ndimdest, destdims, mxGetClassID(srcmat),
wolffd@0 126 mxIsComplex(srcmat)?mxCOMPLEX:mxREAL);
wolffd@0 127
wolffd@0 128 /* if any rep[i] == 0, output should be empty array.
wolffd@0 129 Added by KPM 11/13/02.
wolffd@0 130 */
wolffd@0 131 empty = 0;
wolffd@0 132 for (i=0; i < nrep; i++) {
wolffd@0 133 if (rep[i]==0)
wolffd@0 134 empty = 1;
wolffd@0 135 }
wolffd@0 136 if (empty)
wolffd@0 137 return;
wolffd@0 138
wolffd@0 139 src = (char*)mxGetData(srcmat);
wolffd@0 140 dest = (char*)mxGetData(plhs[0]);
wolffd@0 141 repmat(dest,src,ndim,destdimsize,dimsize,dims,rep);
wolffd@0 142 if(ndimdest > ndim) memrep(dest,destdimsize[ndim-1],extra_rep);
wolffd@0 143 if(mxIsComplex(srcmat)) {
wolffd@0 144 src = (char*)mxGetPi(srcmat);
wolffd@0 145 dest = (char*)mxGetPi(plhs[0]);
wolffd@0 146 repmat(dest,src,ndim,destdimsize,dimsize,dims,rep);
wolffd@0 147 if(ndimdest > ndim) memrep(dest,destdimsize[ndim-1],extra_rep);
wolffd@0 148 }
wolffd@0 149 }