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
|
//------------------------------------------------------------------------------
// CHOLMOD/MATLAB/lxbpattern: MATLAB interface for CHOLMOD symbolic x=L\b solve
//------------------------------------------------------------------------------
// CHOLMOD/MATLAB Module. Copyright (C) 2005-2023, Timothy A. Davis.
// All Rights Reserved.
// SPDX-License-Identifier: GPL-2.0+
//------------------------------------------------------------------------------
// Find the nonzero pattern of x=L\b for a sparse vector b. If numerical
// cancellation has caused entries to drop in L, then this function may
// give an incorrect result.
//
// xpattern = lxbpattern (L,b), same as xpattern = find (L\b),
// assuming no numerical cancellation, except that xpattern will not
// appear in sorted ordering (it appears in topological ordering).
//
// The core cholmod_lsolve_pattern function takes O(length (xpattern)) time,
// except that the initialzations in this mexFunction interface add O(n) time.
//
// This function is for testing cholmod_lsolve_pattern only.
#include "sputil2.h"
void mexFunction
(
int nargout,
mxArray *pargout [ ],
int nargin,
const mxArray *pargin [ ]
)
{
double dummy = 0 ;
int64_t *Lp, *Lnz, *Xp, *Xi, xnz ;
cholmod_sparse *B, Bmatrix, *X ;
cholmod_factor *L ;
cholmod_common Common, *cm ;
int64_t j, n ;
//--------------------------------------------------------------------------
// start CHOLMOD and set parameters
//--------------------------------------------------------------------------
cm = &Common ;
cholmod_l_start (cm) ;
sputil2_config (SPUMONI, cm) ;
//--------------------------------------------------------------------------
// check inputs
//--------------------------------------------------------------------------
if (nargin != 2 || nargout > 1)
{
mexErrMsgTxt ("usage: xpattern = lxbpattern (L,b)") ;
}
n = mxGetN (pargin [0]) ;
if (!mxIsSparse (pargin [0]) || n != mxGetM (pargin [0]))
{
mexErrMsgTxt ("lxbpattern: L must be sparse and square") ;
}
if (n != mxGetM (pargin [1]) || mxGetN (pargin [1]) != 1)
{
mexErrMsgTxt ("lxbpattern: b wrong dimension") ;
}
if (!mxIsSparse (pargin [1]))
{
mexErrMsgTxt ("lxbpattern: b must be sparse") ;
}
//--------------------------------------------------------------------------
// get the sparse b
//--------------------------------------------------------------------------
// get sparse matrix B (unsymmetric)
size_t B_xsize = 0 ;
B = sputil2_get_sparse (pargin [1], 0, CHOLMOD_DOUBLE, &Bmatrix,
&B_xsize, cm) ;
//--------------------------------------------------------------------------
// construct a shallow copy of the input sparse matrix L
//--------------------------------------------------------------------------
// the construction of the CHOLMOD takes O(n) time and memory
// allocate the CHOLMOD symbolic L
L = cholmod_l_allocate_factor (n, cm) ;
L->ordering = CHOLMOD_NATURAL ;
// get the MATLAB L
L->p = mxGetJc (pargin [0]) ;
L->i = mxGetIr (pargin [0]) ;
L->x = mxGetData (pargin [0]) ;
L->z = NULL ;
// allocate and initialize the rest of L
L->nz = cholmod_l_malloc (n, sizeof (int64_t), cm) ;
Lp = L->p ;
Lnz = L->nz ;
for (j = 0 ; j < n ; j++)
{
Lnz [j] = Lp [j+1] - Lp [j] ;
}
// L is not truly a valid CHOLMOD sparse factor, since L->prev and next are
// NULL. But these pointers are not accessed in cholmod_lsolve_pattern
L->prev = NULL ;
L->next = NULL ;
L->xtype = (mxIsComplex (pargin [0])) ? CHOLMOD_COMPLEX : CHOLMOD_REAL ;
L->nzmax = Lp [n] ;
//--------------------------------------------------------------------------
// allocate size-n space for the result X
//--------------------------------------------------------------------------
X = cholmod_l_allocate_sparse (n, 1, n, FALSE, TRUE, 0, 0, cm) ;
//--------------------------------------------------------------------------
// find the pattern of X=L\B
//--------------------------------------------------------------------------
cholmod_l_lsolve_pattern (B, L, X, cm) ;
//--------------------------------------------------------------------------
// return result, converting to 1-based
//--------------------------------------------------------------------------
Xp = (int64_t *) X->p ;
Xi = (int64_t *) X->i ;
xnz = Xp [1] ;
pargout [0] = sputil2_put_int (Xi, xnz, 1) ;
//--------------------------------------------------------------------------
// free workspace and the CHOLMOD L, except for what is copied to MATLAB
//--------------------------------------------------------------------------
L->p = NULL ;
L->i = NULL ;
L->x = NULL ;
L->z = NULL ;
sputil2_free_sparse (&B, &Bmatrix, B_xsize, cm) ;
cholmod_l_free_factor (&L, cm) ;
cholmod_l_free_sparse (&X, cm) ;
cholmod_l_finish (cm) ;
if (SPUMONI > 0) cholmod_l_print_common (" ", cm) ;
}
|
