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
68
|
// CXSparse/MATLAB/CSparse/cs_thumb_mex: 2D dense thumbnail of a sparse matrix
// CXSparse, Copyright (c) 2006-2022, Timothy A. Davis. All Rights Reserved.
// SPDX-License-Identifier: LGPL-2.1+
#include "cs_mex.h"
/* cs_thumb: convert a sparse matrix to a dense 2D thumbnail matrix of size
* at most k-by-k. k defaults to 256. A helper mexFunction for cspy. */
#define INDEX(i,j,lda) ((i)+(j)*(lda))
#define ISNAN(x) ((x) != (x))
#ifdef DBL_MAX
#define BIG_VALUE DBL_MAX
#else
#define BIG_VALUE 1.7e308
#endif
void mexFunction
(
int nargout,
mxArray *pargout [ ],
int nargin,
const mxArray *pargin [ ]
)
{
int64_t m, n, mn, m2, n2, k, s, j, ij, sj, si, p, *Ap, *Ai ;
double aij, ax, az, *S, *Ax, *Az ;
if (nargout > 1 || nargin < 1 || nargin > 2)
{
mexErrMsgTxt ("Usage: S = cs_thumb(A,k)") ;
}
cs_mex_check (0, -1, -1, 0, 1, 1, pargin [0]) ;
m = mxGetM (pargin [0]) ;
n = mxGetN (pargin [0]) ;
mn = CS_MAX (m,n) ;
k = (nargin == 1) ? 256 : mxGetScalar (pargin [1]) ; /* get k */
/* s = size of each submatrix; A(1:s,1:s) maps to S(1,1) */
s = (mn < k) ? 1 : (int64_t) ceil ((double) mn / (double) k) ;
m2 = (int64_t) ceil ((double) m / (double) s) ;
n2 = (int64_t) ceil ((double) n / (double) s) ;
/* create S */
pargout [0] = mxCreateDoubleMatrix (m2, n2, mxREAL) ;
S = mxGetPr (pargout [0]) ;
Ap = (int64_t *) mxGetJc (pargin [0]) ;
Ai = (int64_t *) mxGetIr (pargin [0]) ;
Ax = mxGetPr (pargin [0]) ;
Az = (mxIsComplex (pargin [0])) ? mxGetPi (pargin [0]) : NULL ;
for (j = 0 ; j < n ; j++)
{
sj = j/s ;
for (p = Ap [j] ; p < Ap [j+1] ; p++)
{
si = Ai [p] / s ;
ij = INDEX (si,sj,m2) ;
ax = Ax [p] ;
az = Az ? Az [p] : 0 ;
if (az == 0)
{
aij = fabs (ax) ;
}
else
{
aij = sqrt (ax*ax + az*az) ;
}
if (ISNAN (aij)) aij = BIG_VALUE ;
aij = CS_MIN (BIG_VALUE, aij) ;
S [ij] = CS_MAX (S [ij], aij) ;
}
}
}
|
