1 | /* C implementation of repmat. |
---|
2 | * by Tom Minka |
---|
3 | */ |
---|
4 | /* |
---|
5 | mex -c mexutil.c |
---|
6 | mex repmat.c mexutil.obj |
---|
7 | to check for warnings: |
---|
8 | gcc -Wall -I/cygdrive/c/MATLAB6p1/extern/include -c repmat.c |
---|
9 | */ |
---|
10 | #include "mexutil.h" |
---|
11 | #include <string.h> |
---|
12 | |
---|
13 | /* ALWAYS_2D == 1 gives the usual matlab behavior where repmat(0,3) is |
---|
14 | * a 3x3 matrix. |
---|
15 | * ALWAYS_2D == 0 gives a 3x1 matrix instead. |
---|
16 | * repmat(x,sizes,1) is another way to get this behavior. |
---|
17 | */ |
---|
18 | #define ALWAYS_2D 1 |
---|
19 | |
---|
20 | /* repeat a block of memory rep times */ |
---|
21 | void memrep(char *dest, size_t chunk, int rep) |
---|
22 | { |
---|
23 | #if 0 |
---|
24 | /* slow way */ |
---|
25 | int i; |
---|
26 | char *p = dest; |
---|
27 | for(i=1;i<rep;i++) { |
---|
28 | p += chunk; |
---|
29 | memcpy(p, dest, chunk); |
---|
30 | } |
---|
31 | #else |
---|
32 | /* fast way */ |
---|
33 | if(rep == 1) return; |
---|
34 | memcpy(dest + chunk, dest, chunk); |
---|
35 | if(rep & 1) { |
---|
36 | dest += chunk; |
---|
37 | memcpy(dest + chunk, dest, chunk); |
---|
38 | } |
---|
39 | /* now repeat using a block twice as big */ |
---|
40 | memrep(dest, chunk<<1, rep>>1); |
---|
41 | #endif |
---|
42 | } |
---|
43 | |
---|
44 | void repmat(char *dest, const char *src, int ndim, int *destdimsize, |
---|
45 | int *dimsize, const int *dims, int *rep) |
---|
46 | { |
---|
47 | int d = ndim-1; |
---|
48 | int i, chunk; |
---|
49 | /* copy the first repetition into dest */ |
---|
50 | if(d == 0) { |
---|
51 | chunk = dimsize[0]; |
---|
52 | memcpy(dest,src,chunk); |
---|
53 | } |
---|
54 | else { |
---|
55 | /* recursively repeat each slice of src */ |
---|
56 | for(i=0;i<dims[d];i++) { |
---|
57 | repmat(dest + i*destdimsize[d-1], src + i*dimsize[d-1], |
---|
58 | ndim-1, destdimsize, dimsize, dims, rep); |
---|
59 | } |
---|
60 | chunk = destdimsize[d-1]*dims[d]; |
---|
61 | } |
---|
62 | /* copy the result rep-1 times */ |
---|
63 | memrep(dest,chunk,rep[d]); |
---|
64 | } |
---|
65 | |
---|
66 | void mexFunction(int nlhs, mxArray *plhs[], |
---|
67 | int nrhs, const mxArray *prhs[]) |
---|
68 | { |
---|
69 | const mxArray *srcmat; |
---|
70 | int ndim, *dimsize, eltsize; |
---|
71 | const int *dims; |
---|
72 | int ndimdest, *destdims, *destdimsize; |
---|
73 | char *src, *dest; |
---|
74 | int *rep; |
---|
75 | int i,nrep; |
---|
76 | int extra_rep = 1; |
---|
77 | int empty; |
---|
78 | |
---|
79 | if(nrhs < 2) mexErrMsgTxt("Usage: repmat(A, [M N ...])"); |
---|
80 | srcmat = prhs[0]; |
---|
81 | if(mxIsSparse(srcmat) || mxIsCell(srcmat) || mxIsStruct(srcmat)) { |
---|
82 | /* call Matlab's repmat */ |
---|
83 | mexCallMATLAB(nlhs,plhs,nrhs,prhs,"xrepmat");return; |
---|
84 | /* mexErrMsgTxt("Sorry, can't handle sparse matrices yet."); */ |
---|
85 | } |
---|
86 | ndim = mxGetNumberOfDimensions(srcmat); |
---|
87 | dims = mxGetDimensions(srcmat); |
---|
88 | eltsize = mxGetElementSize(srcmat); |
---|
89 | |
---|
90 | /* compute dimension sizes */ |
---|
91 | dimsize = mxCalloc(ndim, sizeof(int)); |
---|
92 | dimsize[0] = eltsize*dims[0]; |
---|
93 | for(i=1;i<ndim;i++) dimsize[i] = dimsize[i-1]*dims[i]; |
---|
94 | |
---|
95 | /* determine repetition vector */ |
---|
96 | ndimdest = ndim; |
---|
97 | if(nrhs == 2) { |
---|
98 | /* prhs[1] is a vector of reps */ |
---|
99 | nrep = mxGetN(prhs[1]); |
---|
100 | if(nrep > ndimdest) ndimdest = nrep; |
---|
101 | rep = mxCalloc(ndimdest, sizeof(int)); |
---|
102 | for(i=0;i<nrep;i++) { |
---|
103 | double repv = mxGetPr(prhs[1])[i]; |
---|
104 | rep[i] = (int)repv; |
---|
105 | } |
---|
106 | #if ALWAYS_2D |
---|
107 | if(nrep == 1) { |
---|
108 | /* special behavior */ |
---|
109 | nrep = 2; |
---|
110 | rep[1] = rep[0]; |
---|
111 | } |
---|
112 | #endif |
---|
113 | } |
---|
114 | else { |
---|
115 | /* concatenate all prhs's */ |
---|
116 | int ri=0; |
---|
117 | nrep = 0; |
---|
118 | for(i=0;i<nrhs-1;i++) { |
---|
119 | nrep += mxGetNumberOfElements(prhs[i+1]); |
---|
120 | } |
---|
121 | if(nrep > ndimdest) ndimdest = nrep; |
---|
122 | rep = mxCalloc(ndimdest, sizeof(int)); |
---|
123 | for(i=0;i<nrhs-1;i++) { |
---|
124 | double *p = mxGetPr(prhs[i+1]); |
---|
125 | int j, sz = mxGetNumberOfElements(prhs[i+1]); |
---|
126 | for(j=0;j<sz;j++) rep[ri++] = (int)p[j]; |
---|
127 | } |
---|
128 | } |
---|
129 | for(i=nrep;i<ndimdest;i++) rep[i] = 1; |
---|
130 | |
---|
131 | /* compute output size */ |
---|
132 | destdims = mxCalloc(ndimdest, sizeof(int)); |
---|
133 | for(i=0;i<ndim;i++) destdims[i] = dims[i]*rep[i]; |
---|
134 | for(;i<ndimdest;i++) { |
---|
135 | destdims[i] = rep[i]; |
---|
136 | extra_rep *= rep[i]; |
---|
137 | } |
---|
138 | destdimsize = mxCalloc(ndim, sizeof(int)); |
---|
139 | destdimsize[0] = eltsize*destdims[0]; |
---|
140 | for(i=1;i<ndim;i++) destdimsize[i] = destdimsize[i-1]*destdims[i]; |
---|
141 | |
---|
142 | /* for speed, array should be uninitialized */ |
---|
143 | plhs[0] = mxCreateNumericArrayE(ndimdest, destdims, mxGetClassID(srcmat), |
---|
144 | mxIsComplex(srcmat)?mxCOMPLEX:mxREAL); |
---|
145 | |
---|
146 | /* if any rep[i] == 0, output should be empty array. |
---|
147 | Added by KPM 11/13/02. |
---|
148 | */ |
---|
149 | empty = 0; |
---|
150 | for (i=0; i < nrep; i++) { |
---|
151 | if (rep[i]==0) |
---|
152 | empty = 1; |
---|
153 | } |
---|
154 | if (empty) |
---|
155 | return; |
---|
156 | |
---|
157 | src = (char*)mxGetData(srcmat); |
---|
158 | dest = (char*)mxGetData(plhs[0]); |
---|
159 | repmat(dest,src,ndim,destdimsize,dimsize,dims,rep); |
---|
160 | if(ndimdest > ndim) memrep(dest,destdimsize[ndim-1],extra_rep); |
---|
161 | if(mxIsComplex(srcmat)) { |
---|
162 | src = (char*)mxGetPi(srcmat); |
---|
163 | dest = (char*)mxGetPi(plhs[0]); |
---|
164 | repmat(dest,src,ndim,destdimsize,dimsize,dims,rep); |
---|
165 | if(ndimdest > ndim) memrep(dest,destdimsize[ndim-1],extra_rep); |
---|
166 | } |
---|
167 | } |
---|