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
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
|
//------------------------------------------------------------------------------
// CHOLMOD/MATLAB/mread: read a matrix in Matrix Market format
//------------------------------------------------------------------------------
// CHOLMOD/MATLAB Module. Copyright (C) 2005-2023, Timothy A. Davis.
// All Rights Reserved.
// SPDX-License-Identifier: GPL-2.0+
//------------------------------------------------------------------------------
// [A Z] = mread (filename, prefer_binary)
//
// Read a sparse or dense matrix from a file in Matrix Market format.
//
// All MatrixMarket formats are supported.
// The Matrix Market "integer" format is converted into real, but the values
// are preserved. The "pattern" format is converted into real. If a pattern
// matrix is unsymmetric, all of its values are equal to one. If a pattern is
// symmetric, the kth diagonal entry is set to one plus the number of
// off-diagonal nonzeros in row/column k, and off-diagonal entries are set to
// -1.
//
// Explicit zero entries are returned as the binary pattern of the matrix Z.
#include "sputil2.h"
// maximum file length
#define MAXLEN 1030
void mexFunction
(
int nargout,
mxArray *pargout [ ],
int nargin,
const mxArray *pargin [ ]
)
{
void *G ;
cholmod_dense *X = NULL ;
cholmod_sparse *A = NULL, *Z = NULL ;
cholmod_common Common, *cm ;
int64_t *Ap = NULL, *Ai ;
double *Ax, *Az = NULL ;
char filename [MAXLEN] ;
int64_t nz, k, is_complex = FALSE, nrow = 0, ncol = 0 ;
int mtype ;
//--------------------------------------------------------------------------
// start CHOLMOD and set parameters
//--------------------------------------------------------------------------
cm = &Common ;
cholmod_l_start (cm) ;
sputil2_config (SPUMONI, cm) ;
//--------------------------------------------------------------------------
// get inputs
//--------------------------------------------------------------------------
if (nargin < 1 || nargin > 2 || nargout > 2)
{
mexErrMsgTxt ("usage: [A Z] = mread (filename, prefer_binary)") ;
}
if (!mxIsChar (pargin [0]))
{
mexErrMsgTxt ("mread requires a filename") ;
}
mxGetString (pargin [0], filename, MAXLEN) ;
sputil2_file = fopen (filename, "r") ;
if (sputil2_file == NULL)
{
mexErrMsgTxt ("cannot open file") ;
}
if (nargin > 1)
{
cm->prefer_binary = (mxGetScalar (pargin [1]) != 0) ;
}
//--------------------------------------------------------------------------
// read the matrix, as either a dense or sparse matrix
//--------------------------------------------------------------------------
G = cholmod_l_read_matrix (sputil2_file, 1, &mtype, cm) ;
fclose (sputil2_file) ;
sputil2_file = NULL ;
if (G == NULL)
{
mexErrMsgTxt ("could not read file") ;
}
// get the specific matrix (A or X), and change to complex if needed
if (mtype == CHOLMOD_SPARSE)
{
A = (cholmod_sparse *) G ;
nrow = A->nrow ;
ncol = A->ncol ;
Ap = A->p ;
Ai = A->i ;
if (A->xtype == CHOLMOD_ZOMPLEX)
{
// if complex, ensure A is complex, not zomplex
cholmod_l_sparse_xtype (CHOLMOD_COMPLEX, A, cm) ;
}
is_complex = (A->xtype == CHOLMOD_COMPLEX) ;
Ax = A->x ;
}
else if (mtype == CHOLMOD_DENSE)
{
X = (cholmod_dense *) G ;
nrow = X->nrow ;
ncol = X->ncol ;
if (X->xtype == CHOLMOD_ZOMPLEX)
{
// if complex, ensure X is complex, not zomplex
cholmod_l_dense_xtype (CHOLMOD_COMPLEX, X, cm) ;
}
is_complex = (X->xtype == CHOLMOD_COMPLEX) ;
Ax = X->x ;
}
else
{
mexErrMsgTxt ("invalid file") ;
}
//--------------------------------------------------------------------------
// if requested, extract the zero entries and place them in Z
//--------------------------------------------------------------------------
if (nargout > 1)
{
if (mtype == CHOLMOD_SPARSE)
{
// A is a sparse real/zomplex double matrix
Z = sputil2_extract_zeros (A, cm) ;
}
else
{
// input is full; just return an empty Z matrix
Z = cholmod_l_spzeros (nrow, ncol, 0, CHOLMOD_REAL, cm) ;
}
}
//--------------------------------------------------------------------------
// change a complex matrix to real if its imaginary part is all zero
//--------------------------------------------------------------------------
if (is_complex)
{
if (mtype == CHOLMOD_SPARSE)
{
nz = Ap [ncol] ;
}
else
{
nz = nrow * ncol ;
}
bool allzero = true ;
for (k = 0 ; k < nz ; k++)
{
if (Ax [2*k+1] != 0)
{
allzero = false ;
break ;
}
}
if (allzero)
{
// discard the all-zero imaginary part
if (mtype == CHOLMOD_SPARSE)
{
cholmod_l_sparse_xtype (CHOLMOD_REAL, A, cm) ;
}
else
{
cholmod_l_dense_xtype (CHOLMOD_REAL, X, cm) ;
}
}
}
//--------------------------------------------------------------------------
// return results to MATLAB
//--------------------------------------------------------------------------
if (mtype == CHOLMOD_SPARSE)
{
// drop explicit zeros from A; their pattern is kept in Z
pargout [0] = sputil2_put_sparse (&A, mxDOUBLE_CLASS,
/* drop explicit zeros */ true, cm) ;
}
else
{
pargout [0] = sputil2_put_dense (&X, mxDOUBLE_CLASS, cm) ;
}
if (nargout > 1)
{
pargout [1] = sputil2_put_sparse (&Z, mxDOUBLE_CLASS,
/* Z is binary so it has no zeros to drop */ false, cm) ;
}
//--------------------------------------------------------------------------
// free workspace
//--------------------------------------------------------------------------
cholmod_l_finish (cm) ;
if (SPUMONI > 0) cholmod_l_print_common (" ", cm) ;
}
|
