1
/* ========================================================================== */
2
/* === ldlmex.c: LDL mexFunction =========================================== */
3
/* ========================================================================== */
5
/* MATLAB interface for numerical LDL' factorization using the LDL sparse matrix
8
* MATLAB calling syntax is:
10
* [L, D, Parent, flops] = ldl (A)
11
* [L, D, Parent, flops] = ldl (A, P)
12
* [x, flops] = ldl (A, [ ], b)
13
* [x, flops] = ldl (A, P, b)
15
* The factorization is L*D*L' = A or L*D*L' = A(P,P). A must be sparse,
16
* square, and real. L is lower triangular with unit diagonal, but the diagonal
17
* is not returned. D is diagonal sparse matrix. Let n = size (A,1). If P is
18
* not present or empty, the factorization is:
20
* (L + speye (n)) * D * (L + speye (n))' = A
22
* otherwise, the factorization is
24
* (L + speye (n)) * D * (L + speye (n))' = A(P,P)
26
* P is a permutation of 1:n, an output of AMD, SYMAMD, or SYMRCM, for example.
27
* Only the diagonal and upper triangular part of A or A(P,P) is accessed; the
28
* lower triangular part is ignored.
30
* The elimination tree is returned in the Parent array.
32
* In the x = ldl (A, P, b) usage, the LDL' factorization is not returned.
33
* Instead, the system A*x=b is solved for x, where b is a dense n-by-m matrix,
34
* using P as the fill-reducing ordering for the LDL' factorization of A(P,P).
35
* If P is not present or equal to [ ], it is assumed to be the identity
38
* If no zero entry on the diagonal of D is encountered, then the flops argument
39
* is the floating point count.
41
* If any entry on the diagonal of D is zero, then the LDL' factorization is
42
* terminated at that point. If there is no flops output argument, an error
43
* message is printed and no outputs are returned. Otherwise, flops is
44
* negative, d = -flops, and D (d,d) is the first zero entry on the diagonal of
45
* D. A partial factorization is returned. Let B = A if P is not present or
46
* empty, or B = A(P,P) otherwise. Then the factorization is
48
* LDL = (L + speye (n)) * D * (L + speye (n))'
49
* LDL (1:d, 1:d) = B (1:d,1:d)
51
* That is, the LDL' factorization of B (1:d,1:d) is in the first d rows and
52
* columns of L and D. The rest of L and D are zero.
54
* LDL Copyright (c) by Timothy A Davis,
55
* University of Florida. All Rights Reserved. See README for the License.
65
/* ========================================================================== */
66
/* === LDL mexFunction ====================================================== */
67
/* ========================================================================== */
74
const mxArray *pargin[ ]
77
UF_long i, n, *Pattern, *Flag, *Li, *Lp, *Ap, *Ai, *Lnz, *Parent, do_chol,
78
nrhs, lnz, do_solve, *P, *Pinv, nn, k, j, permute, *Dp, *Di, d,
79
do_flops, psrc, pdst, p2 ;
80
double *Y, *D, *Lx, *Ax, flops, *X, *B, *p ;
82
/* ---------------------------------------------------------------------- */
83
/* get inputs and allocate workspace */
84
/* ---------------------------------------------------------------------- */
86
do_chol = (nargin > 0) && (nargin <= 2) && (nargout <= 4) ;
87
do_solve = (nargin == 3) && (nargout <= 2) ;
88
if (!(do_chol || do_solve))
90
mexErrMsgTxt ("Usage:\n"
91
" [L, D, etree, flopcount] = ldl (A) ;\n"
92
" [L, D, etree, flopcount] = ldl (A, P) ;\n"
93
" [x, flopcount] = ldl (A, [ ], b) ;\n"
94
" [x, flopcount] = ldl (A, P, b) ;\n"
95
"The etree and flopcount arguments are optional.") ;
97
n = mxGetM (pargin [0]) ;
98
if (!mxIsSparse (pargin [0]) || n != mxGetN (pargin [0])
99
|| mxIsComplex (pargin [0]))
101
mexErrMsgTxt ("ldl: A must be sparse, square, and real") ;
105
if (mxIsSparse (pargin [2]) || n != mxGetM (pargin [2])
106
|| !mxIsDouble (pargin [2]) || mxIsComplex (pargin [2]))
109
"ldl: b must be dense, real, and with proper dimension") ;
112
nn = (n == 0) ? 1 : n ;
114
/* get sparse matrix A */
115
Ap = (UF_long *) mxGetJc (pargin [0]) ;
116
Ai = (UF_long *) mxGetIr (pargin [0]) ;
117
Ax = mxGetPr (pargin [0]) ;
119
/* get fill-reducing ordering, if present */
120
permute = ((nargin > 1) && !mxIsEmpty (pargin [1])) ;
123
if (mxGetM (pargin [1]) * mxGetN (pargin [1]) != n ||
124
mxIsSparse (pargin [1]))
126
mexErrMsgTxt ("ldl: invalid input permutation\n") ;
128
P = (UF_long *) mxMalloc (nn * sizeof (UF_long)) ;
129
Pinv = (UF_long *) mxMalloc (nn * sizeof (UF_long)) ;
130
p = mxGetPr (pargin [1]) ;
131
for (k = 0 ; k < n ; k++)
133
P [k] = p [k] - 1 ; /* convert to 0-based */
138
P = (UF_long *) NULL ;
139
Pinv = (UF_long *) NULL ;
142
/* allocate first part of L */
143
Lp = (UF_long *) mxMalloc ((n+1) * sizeof (UF_long)) ;
144
Parent = (UF_long *) mxMalloc (nn * sizeof (UF_long)) ;
147
Y = (double *) mxMalloc (nn * sizeof (double)) ;
148
Flag = (UF_long *) mxMalloc (nn * sizeof (UF_long)) ;
149
Pattern = (UF_long *) mxMalloc (nn * sizeof (UF_long)) ;
150
Lnz = (UF_long *) mxMalloc (nn * sizeof (UF_long)) ;
152
/* make sure the input P is valid */
153
if (permute && !ldl_l_valid_perm (n, P, Flag))
155
mexErrMsgTxt ("ldl: invalid input permutation\n") ;
158
/* note that we assume that the input matrix is valid */
160
/* ---------------------------------------------------------------------- */
161
/* symbolic factorization to get Lp, Parent, Lnz, and optionally Pinv */
162
/* ---------------------------------------------------------------------- */
164
ldl_l_symbolic (n, Ap, Ai, Lp, Parent, Lnz, Flag, P, Pinv) ;
167
/* ---------------------------------------------------------------------- */
169
/* ---------------------------------------------------------------------- */
173
/* create the output matrix L, using the Lp array from ldl_l_symbolic */
174
pargout [0] = mxCreateSparse (n, n, lnz+1, mxREAL) ;
175
mxFree (mxGetJc (pargout [0])) ;
176
mxSetJc (pargout [0], (void *) Lp) ; /* Lp is not mxFree'd */
177
Li = (UF_long *) mxGetIr (pargout [0]) ;
178
Lx = mxGetPr (pargout [0]) ;
180
/* create sparse diagonal matrix D */
183
pargout [1] = mxCreateSparse (n, n, nn, mxREAL) ;
184
Dp = (UF_long *) mxGetJc (pargout [1]) ;
185
Di = (UF_long *) mxGetIr (pargout [1]) ;
186
for (j = 0 ; j < n ; j++)
192
D = mxGetPr (pargout [1]) ;
196
D = (double *) mxMalloc (nn * sizeof (double)) ;
199
/* return elimination tree (add 1 to change from 0-based to 1-based) */
202
pargout [2] = mxCreateDoubleMatrix (1, n, mxREAL) ;
203
p = mxGetPr (pargout [2]) ;
204
for (i = 0 ; i < n ; i++)
206
p [i] = Parent [i] + 1 ;
210
do_flops = (nargout == 4) ? (3) : (-1) ;
214
/* create L and D as temporary matrices */
215
Li = (UF_long *) mxMalloc ((lnz+1) * sizeof (UF_long)) ;
216
Lx = (double *) mxMalloc ((lnz+1) * sizeof (double)) ;
217
D = (double *) mxMalloc (nn * sizeof (double)) ;
219
/* create solution x */
220
nrhs = mxGetN (pargin [2]) ;
221
pargout [0] = mxCreateDoubleMatrix (n, nrhs, mxREAL) ;
222
X = mxGetPr (pargout [0]) ;
223
B = mxGetPr (pargin [2]) ;
225
do_flops = (nargout == 2) ? (1) : (-1) ;
230
/* find flop count for ldl_l_numeric */
232
for (k = 0 ; k < n ; k++)
234
flops += ((double) Lnz [k]) * (Lnz [k] + 2) ;
238
/* add flop count for solve */
239
for (k = 0 ; k < n ; k++)
241
flops += 4 * ((double) Lnz [k]) + 1 ;
244
pargout [do_flops] = mxCreateDoubleMatrix (1, 1, mxREAL) ;
245
p = mxGetPr (pargout [do_flops]) ;
249
/* ---------------------------------------------------------------------- */
250
/* numeric factorization to get Li, Lx, and D */
251
/* ---------------------------------------------------------------------- */
253
d = ldl_l_numeric (n, Ap, Ai, Ax, Lp, Parent, Lnz, Li, Lx, D, Y, Flag,
256
/* ---------------------------------------------------------------------- */
257
/* singular case : truncate the factorization */
258
/* ---------------------------------------------------------------------- */
262
/* D [d] is zero: report error, or clean up */
263
if (do_chol && do_flops < 0)
265
mexErrMsgTxt ("ldl: zero pivot encountered\n") ;
269
/* L and D are incomplete, compact them */
272
for (k = d ; k < n ; k++)
278
for (k = d ; k < n ; k++)
283
for (k = 0 ; k < d ; k++)
285
for (psrc = Lp [k] ; psrc < Lp [k] + Lnz [k] ; psrc++)
287
Li [pdst] = Li [psrc] ;
288
Lx [pdst] = Lx [psrc] ;
292
for (k = 0 ; k < d ; k++)
294
Lp [k+1] = Lp [k] + Lnz [k] ;
296
for (k = d ; k <= n ; k++)
302
/* return -d instead of the flop count (convert to 1-based) */
303
p = mxGetPr (pargout [do_flops]) ;
309
/* ---------------------------------------------------------------------- */
310
/* solve Ax=b, if requested */
311
/* ---------------------------------------------------------------------- */
317
for (j = 0 ; j < nrhs ; j++)
319
ldl_l_perm (n, Y, B, P) ; /* y = Pb */
320
ldl_l_lsolve (n, Y, Lp, Li, Lx) ; /* y = L\y */
321
ldl_l_dsolve (n, Y, D) ; /* y = D\y */
322
ldl_l_ltsolve (n, Y, Lp, Li, Lx) ; /* y = L'\y */
323
ldl_l_permt (n, X, Y, P) ; /* x = P'y */
330
for (j = 0 ; j < nrhs ; j++)
332
for (k = 0 ; k < n ; k++) /* x = b */
336
ldl_l_lsolve (n, X, Lp, Li, Lx) ; /* x = L\x */
337
ldl_l_dsolve (n, X, D) ; /* x = D\x */
338
ldl_l_ltsolve (n, X, Lp, Li, Lx) ; /* x = L'\x */
343
/* free the matrix L */
350
/* ---------------------------------------------------------------------- */
352
/* ---------------------------------------------------------------------- */
354
if (do_chol && nargout < 2)
added
removed
LDL/MATLAB/ldlmex.c
