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
|
// =============================================================================
// === spqr_singletons mexFunction =============================================
// =============================================================================
// SPQR, Copyright (c) 2008-2022, Timothy A Davis. All Rights Reserved.
// SPDX-License-Identifier: GPL-2.0+
#include "spqr_mx.hpp"
#include "spqr.hpp"
// Finds the row and column singletons of a sparse matrix. Note that this
// function uses "non-usercallable" functions from SuiteSparseQR.
// See spqr_singletons.m for details.
//
// [p q n1rows n1cols tol] = spqr_singletons (A)
// [p q n1rows n1cols] = spqr_singletons (A, tol)
#define TRUE 1
#define FALSE 0
void mexFunction
(
int nargout,
mxArray *pargout [ ],
int nargin,
const mxArray *pargin [ ]
)
{
int64_t *P, *Q, *Rp, *Pinv ;
double *Ax, dummy, tol ;
int64_t m, n, anz, is_complex, n1rows, n1cols, i, k ;
cholmod_sparse *A, Amatrix, *Y ;
cholmod_common Common, *cc ;
// -------------------------------------------------------------------------
// start CHOLMOD and set parameters
// -------------------------------------------------------------------------
cc = &Common ;
cholmod_l_start (cc) ;
spqr_mx_config (SPUMONI, cc) ;
// -------------------------------------------------------------------------
// check inputs
// -------------------------------------------------------------------------
if (nargout > 5)
{
mexErrMsgIdAndTxt ("MATLAB:maxlhs", "Too many output arguments") ;
}
if (nargin < 1)
{
mexErrMsgIdAndTxt ("MATLAB:minrhs", "Not enough input arguments") ;
}
if (nargin > 2)
{
mexErrMsgIdAndTxt ("MATLAB:maxrhs", "Too many input arguments") ;
}
// -------------------------------------------------------------------------
// get the input matrix A and convert to merged-complex if needed
// -------------------------------------------------------------------------
if (!mxIsSparse (pargin [0]))
{
mexErrMsgIdAndTxt ("QR:invalidInput", "A must be sparse") ;
}
A = spqr_mx_get_sparse (pargin [0], &Amatrix, &dummy) ;
m = A->nrow ;
n = A->ncol ;
is_complex = mxIsComplex (pargin [0]) ;
Ax = spqr_mx_merge_if_complex (pargin [0], is_complex, &anz, cc) ;
if (is_complex)
{
// A has been converted from real or zomplex to complex
A->x = Ax ;
A->z = NULL ;
A->xtype = CHOLMOD_COMPLEX ;
}
// -------------------------------------------------------------------------
// get the tolerance
// -------------------------------------------------------------------------
if (nargin < 2)
{
tol = is_complex ? spqr_tol <Complex> (A,cc) : spqr_tol <double> (A,cc);
}
else
{
tol = mxGetScalar (pargin [1]) ;
}
// -------------------------------------------------------------------------
// find the singletons
// -------------------------------------------------------------------------
if (is_complex)
{
spqr_1colamd <Complex,int64_t> (SPQR_ORDERING_NATURAL, tol, 0, A,
&Q, &Rp, &Pinv, &Y, &n1cols, &n1rows, cc) ;
}
else
{
spqr_1colamd <double,int64_t> (SPQR_ORDERING_NATURAL, tol, 0, A,
&Q, &Rp, &Pinv, &Y, &n1cols, &n1rows, cc) ;
}
// -------------------------------------------------------------------------
// free unused outputs from spqr_1colamd, and the merged-complex copy of A
// -------------------------------------------------------------------------
cholmod_l_free (n1rows+1, sizeof (int64_t), Rp, cc) ;
cholmod_l_free_sparse (&Y, cc) ;
if (is_complex)
{
// this was allocated by merge_if_complex
cholmod_l_free (anz, sizeof (Complex), Ax, cc) ;
}
// -------------------------------------------------------------------------
// find P from Pinv
// -------------------------------------------------------------------------
P = (int64_t *) cholmod_l_malloc (m, sizeof (int64_t), cc) ;
for (i = 0 ; i < m ; i++)
{
k = Pinv ? Pinv [i] : i ;
P [k] = i ;
}
cholmod_l_free (m, sizeof (int64_t), Pinv, cc) ;
// -------------------------------------------------------------------------
// return results
// -------------------------------------------------------------------------
pargout [0] = spqr_mx_put_permutation (P, m, TRUE, cc) ;
cholmod_l_free (m, sizeof (int64_t), P, cc) ;
if (nargout > 1) pargout [1] = spqr_mx_put_permutation (Q, n, TRUE, cc) ;
cholmod_l_free (n, sizeof (int64_t), Q, cc) ;
if (nargout > 2) pargout [2] = mxCreateDoubleScalar ((double) n1rows) ;
if (nargout > 3) pargout [3] = mxCreateDoubleScalar ((double) n1cols) ;
if (nargout > 4) pargout [4] = mxCreateDoubleScalar (tol) ;
cholmod_l_finish (cc) ;
}
|
