1 #include "cs_mex.h"
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. */
4
5 #define INDEX(i,j,lda) ((i)+(j)*(lda))
6 #define ISNAN(x) ((x) != (x))
7 #ifdef DBL_MAX
8 #define BIG_VALUE DBL_MAX
9 #else
10 #define BIG_VALUE 1.7e308
11 #endif
12
mexFunction(int nargout,mxArray * pargout[],int nargin,const mxArray * pargin[])13 void mexFunction
14 (
15 int nargout,
16 mxArray *pargout [ ],
17 int nargin,
18 const mxArray *pargin [ ]
19 )
20 {
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)
24 {
25 mexErrMsgTxt ("Usage: S = cs_thumb(A,k)") ;
26 }
27 cs_mex_check (0, -1, -1, 0, 1, 1, pargin [0]) ;
28 m = mxGetM (pargin [0]) ;
29 n = mxGetN (pargin [0]) ;
30 mn = CS_MAX (m,n) ;
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) ;
36 /* create 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++)
44 {
45 sj = j/s ;
46 for (p = Ap [j] ; p < Ap [j+1] ; p++)
47 {
48 si = Ai [p] / s ;
49 ij = INDEX (si,sj,m2) ;
50 ax = Ax [p] ;
51 az = Az ? Az [p] : 0 ;
52 if (az == 0)
53 {
54 aij = fabs (ax) ;
55 }
56 else
57 {
58 aij = sqrt (ax*ax + az*az) ;
59 }
60 if (ISNAN (aij)) aij = BIG_VALUE ;
61 aij = CS_MIN (BIG_VALUE, aij) ;
62 S [ij] = CS_MAX (S [ij], aij) ;
63 }
64 }
65 }
66