|
wolffd@0
|
1
|
|
wolffd@0
|
2 m = mxGetM(prhs[0]);
|
|
wolffd@0
|
3 n = mxGetN(prhs[0]);
|
|
wolffd@0
|
4 pr = mxGetPr(prhs[0]);
|
|
wolffd@0
|
5 pi = mxGetPi(prhs[0]);
|
|
wolffd@0
|
6 cmplx = (pi == NULL ? 0 : 1);
|
|
wolffd@0
|
7
|
|
wolffd@0
|
8 /* Allocate space for sparse matrix.
|
|
wolffd@0
|
9 * NOTE: Assume at most 20% of the data is sparse. Use ceil
|
|
wolffd@0
|
10 * to cause it to round up.
|
|
wolffd@0
|
11 */
|
|
wolffd@0
|
12
|
|
wolffd@0
|
13 percent_sparse = 0.2;
|
|
wolffd@0
|
14 nzmax = (int)ceil((double)m*(double)n*percent_sparse);
|
|
wolffd@0
|
15
|
|
wolffd@0
|
16 plhs[0] = mxCreateSparse(m,n,nzmax,cmplx);
|
|
wolffd@0
|
17 sr = mxGetPr(plhs[0]);
|
|
wolffd@0
|
18 si = mxGetPi(plhs[0]);
|
|
wolffd@0
|
19 irs = mxGetIr(plhs[0]);
|
|
wolffd@0
|
20 jcs = mxGetJc(plhs[0]);
|
|
wolffd@0
|
21
|
|
wolffd@0
|
22 /* Copy nonzeros. */
|
|
wolffd@0
|
23 k = 0;
|
|
wolffd@0
|
24 isfull = 0;
|
|
wolffd@0
|
25 for (j = 0; (j < n); j++) {
|
|
wolffd@0
|
26 int i;
|
|
wolffd@0
|
27 jcs[j] = k;
|
|
wolffd@0
|
28 for (i = 0; (i < m); i++) {
|
|
wolffd@0
|
29 if (IsNonZero(pr[i]) || (cmplx && IsNonZero(pi[i]))) {
|
|
wolffd@0
|
30
|
|
wolffd@0
|
31 /* Check to see if non-zero element will fit in
|
|
wolffd@0
|
32 * allocated output array. If not, increase
|
|
wolffd@0
|
33 * percent_sparse by 10%, recalculate nzmax, and augment
|
|
wolffd@0
|
34 * the sparse array.
|
|
wolffd@0
|
35 */
|
|
wolffd@0
|
36 if (k >= nzmax) {
|
|
wolffd@0
|
37 int oldnzmax = nzmax;
|
|
wolffd@0
|
38 percent_sparse += 0.1;
|
|
wolffd@0
|
39 nzmax = (int)ceil((double)m*(double)n*percent_sparse);
|
|
wolffd@0
|
40
|
|
wolffd@0
|
41 /* Make sure nzmax increases atleast by 1. */
|
|
wolffd@0
|
42 if (oldnzmax == nzmax)
|
|
wolffd@0
|
43 nzmax++;
|
|
wolffd@0
|
44
|
|
wolffd@0
|
45 mxSetNzmax(plhs[0], nzmax);
|
|
wolffd@0
|
46 mxSetPr(plhs[0], mxRealloc(sr, nzmax*sizeof(double)));
|
|
wolffd@0
|
47 if (si != NULL)
|
|
wolffd@0
|
48 mxSetPi(plhs[0], mxRealloc(si, nzmax*sizeof(double)));
|
|
wolffd@0
|
49 mxSetIr(plhs[0], mxRealloc(irs, nzmax*sizeof(int)));
|
|
wolffd@0
|
50
|
|
wolffd@0
|
51 sr = mxGetPr(plhs[0]);
|
|
wolffd@0
|
52 si = mxGetPi(plhs[0]);
|
|
wolffd@0
|
53 irs = mxGetIr(plhs[0]);
|
|
wolffd@0
|
54 }
|
|
wolffd@0
|
55 sr[k] = pr[i];
|
|
wolffd@0
|
56 if (cmplx) {
|
|
wolffd@0
|
57 si[k] = pi[i];
|
|
wolffd@0
|
58 }
|
|
wolffd@0
|
59 irs[k] = i;
|
|
wolffd@0
|
60 k++;
|
|
wolffd@0
|
61 }
|
|
wolffd@0
|
62 }
|
|
wolffd@0
|
63 pr += m;
|
|
wolffd@0
|
64 pi += m;
|
|
wolffd@0
|
65 }
|
|
wolffd@0
|
66 jcs[n] = k;
|
|
wolffd@0
|
67 }
|
