comparison toolboxes/FullBNT-1.0.7/KPMtools/repmatC.c @ 0:e9a9cd732c1e tip

first hg version after svn

+/*
+mex -c mexutil.c
+mex repmat.c mexutil.obj
+to check for warnings:
+gcc -Wall -I/cygdrive/c/MATLAB6p1/extern/include -c repmat.c
+*/
+#include "mexutil.h"
+#include <string.h>
+
+/* repeat a block of memory rep times */
+void memrep(char *dest, size_t chunk, int rep)
+{
+#if 0
+  /* slow way */
+  int i;
+  char *p = dest;
+  for(i=1;i<rep;i++) {
+    p += chunk;
+    memcpy(p, dest, chunk);
+  }
+#else
+  /* fast way */
+  if(rep == 1) return;
+  memcpy(dest + chunk, dest, chunk); 
+  if(rep & 1) {
+    dest += chunk;
+    memcpy(dest + chunk, dest, chunk);
+  }
+  /* now repeat using a block twice as big */
+  memrep(dest, chunk<<1, rep>>1);
+#endif
+}
+
+void repmat(char *dest, const char *src, int ndim, int *destdimsize, 
+	    int *dimsize, const int *dims, int *rep) 
+{
+  int d = ndim-1;
+  int i, chunk;
+  /* copy the first repetition into dest */
+  if(d == 0) {
+    chunk = dimsize[0];
+    memcpy(dest,src,chunk);
+  }
+  else {
+    /* recursively repeat each slice of src */
+    for(i=0;i<dims[d];i++) {
+      repmat(dest + i*destdimsize[d-1], src + i*dimsize[d-1], 
+	     ndim-1, destdimsize, dimsize, dims, rep);
+    }
+    chunk = destdimsize[d-1]*dims[d];
+  }
+  /* copy the result rep-1 times */
+  memrep(dest,chunk,rep[d]);
+}
+
+void mexFunction(int nlhs, mxArray *plhs[],
+                 int nrhs, const mxArray *prhs[])
+{
+  const mxArray *srcmat;
+  int ndim, *dimsize, eltsize;
+  const int *dims;
+  int ndimdest, *destdims, *destdimsize;
+  char *src, *dest;
+  int *rep;
+  int i,nrep;
+  int extra_rep = 1;
+  int empty;
+
+  if(nrhs < 2) mexErrMsgTxt("Usage: xrepmat(A, [N M ...])");
+  srcmat = prhs[0];
+  if(mxIsSparse(srcmat)) {
+    mexErrMsgTxt("Sorry, can't handle sparse matrices yet.");
+  }
+  if(mxIsCell(srcmat)) {
+    mexErrMsgTxt("Sorry, can't handle cell arrays yet.");
+  }
+  ndim = mxGetNumberOfDimensions(srcmat);
+  dims = mxGetDimensions(srcmat);
+  eltsize = mxGetElementSize(srcmat);
+
+  /* compute dimension sizes */
+  dimsize = mxCalloc(ndim, sizeof(int));
+  dimsize[0] = eltsize*dims[0];
+  for(i=1;i<ndim;i++) dimsize[i] = dimsize[i-1]*dims[i];
+
+  /* determine repetition vector */
+  ndimdest = ndim;
+  if(nrhs == 2) {
+    nrep = mxGetN(prhs[1]);
+    if(nrep > ndimdest) ndimdest = nrep;
+    rep = mxCalloc(ndimdest, sizeof(int));
+    for(i=0;i<nrep;i++) {
+      double repv = mxGetPr(prhs[1])[i];
+      rep[i] = (int)repv;
+    }
+    if(nrep == 1) {
+      /* special behavior */
+      nrep = 2;
+      rep[1] = rep[0];
+    }
+  }
+  else {
+    nrep = nrhs-1;
+    if(nrep > ndimdest) ndimdest = nrep;
+    rep = mxCalloc(ndimdest, sizeof(int));
+    for(i=0;i<nrep;i++) {
+      rep[i] = (int)*mxGetPr(prhs[i+1]);
+    }
+  }
+  for(i=nrep;i<ndimdest;i++) rep[i] = 1;
+
+  /* compute output size */
+  destdims = mxCalloc(ndimdest, sizeof(int));
+  for(i=0;i<ndim;i++) destdims[i] = dims[i]*rep[i];
+  for(;i<ndimdest;i++) { 
+    destdims[i] = rep[i];
+    extra_rep *= rep[i];
+  }
+  destdimsize = mxCalloc(ndim, sizeof(int));
+  destdimsize[0] = eltsize*destdims[0];
+  for(i=1;i<ndim;i++) destdimsize[i] = destdimsize[i-1]*destdims[i];
+
+    
+  /* for speed, array should be uninitialized */
+  plhs[0] = mxCreateNumericArray(ndimdest, destdims, mxGetClassID(srcmat), 
+				  mxIsComplex(srcmat)?mxCOMPLEX:mxREAL);
+
+  /* if any rep[i] == 0, output should be empty array.
+     Added by KPM 11/13/02.
+  */
+  empty = 0;
+  for (i=0; i < nrep; i++) {
+    if (rep[i]==0) 
+      empty = 1;
+  }
+  if (empty) 
+    return;
+
+  src = (char*)mxGetData(srcmat);
+  dest = (char*)mxGetData(plhs[0]);
+  repmat(dest,src,ndim,destdimsize,dimsize,dims,rep);
+  if(ndimdest > ndim) memrep(dest,destdimsize[ndim-1],extra_rep);
+  if(mxIsComplex(srcmat)) {
+    src = (char*)mxGetPi(srcmat);
+    dest = (char*)mxGetPi(plhs[0]);
+    repmat(dest,src,ndim,destdimsize,dimsize,dims,rep);
+    if(ndimdest > ndim) memrep(dest,destdimsize[ndim-1],extra_rep);
+  }
+}