www.pudn.com > starter.zip > repmatC.c, change:2011-01-04,size:3965b


/* 
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); 
  } 
}