2
/* cs_thumb: convert a sparse matrix to a dense 2D thumbnail matrix of size
3
* at most k-by-k. k defaults to 256. A helper mexFunction for cspy. */
5
#define INDEX(i,j,lda) ((i)+(j)*(lda))
6
#define ISNAN(x) ((x) != (x))
8
#define BIG_VALUE DBL_MAX
10
#define BIG_VALUE 1.7e308
18
const mxArray *pargin [ ]
21
CS_INT m, n, mn, m2, n2, k, s, j, ij, sj, si, p, *Ap, *Ai ;
22
double aij, ax, az, *S, *Ax, *Az ;
23
if (nargout > 1 || nargin < 1 || nargin > 2)
25
mexErrMsgTxt ("Usage: S = cs_thumb(A,k)") ;
27
cs_mex_check (0, -1, -1, 0, 1, 1, pargin [0]) ;
28
m = mxGetM (pargin [0]) ;
29
n = mxGetN (pargin [0]) ;
31
k = (nargin == 1) ? 256 : mxGetScalar (pargin [1]) ; /* get k */
32
/* s = size of each submatrix; A(1:s,1:s) maps to S(1,1) */
33
s = (mn < k) ? 1 : (CS_INT) ceil ((double) mn / (double) k) ;
34
m2 = (CS_INT) ceil ((double) m / (double) s) ;
35
n2 = (CS_INT) ceil ((double) n / (double) s) ;
37
pargout [0] = mxCreateDoubleMatrix (m2, n2, mxREAL) ;
38
S = mxGetPr (pargout [0]) ;
39
Ap = (CS_INT *) mxGetJc (pargin [0]) ;
40
Ai = (CS_INT *) mxGetIr (pargin [0]) ;
41
Ax = mxGetPr (pargin [0]) ;
42
Az = (mxIsComplex (pargin [0])) ? mxGetPi (pargin [0]) : NULL ;
43
for (j = 0 ; j < n ; j++)
46
for (p = Ap [j] ; p < Ap [j+1] ; p++)
49
ij = INDEX (si,sj,m2) ;
51
az = Az ? Az [p] : 0 ;
58
aij = sqrt (ax*ax + az*az) ;
60
if (ISNAN (aij)) aij = BIG_VALUE ;
61
aij = CS_MIN (BIG_VALUE, aij) ;
62
S [ij] = CS_MAX (S [ij], aij) ;