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