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