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