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