1 /* ========================================================================== */
2 /* === Cholesky/cholmod_factorize =========================================== */
3 /* ========================================================================== */
4
5 /* -----------------------------------------------------------------------------
6 * CHOLMOD/Cholesky Module. Copyright (C) 2005-2006, Timothy A. Davis
7 * -------------------------------------------------------------------------- */
8
9 /* Computes the numerical factorization of a symmetric matrix. The primary
10 * inputs to this routine are a sparse matrix A and the symbolic factor L from
11 * cholmod_analyze or a prior numerical factor L. If A is symmetric, this
12 * routine factorizes A(p,p)+beta*I (beta can be zero), where p is the
13 * fill-reducing permutation (L->Perm). If A is unsymmetric, either
14 * A(p,:)*A(p,:)'+beta*I or A(p,f)*A(p,f)'+beta*I is factorized. The set f and
15 * the nonzero pattern of the matrix A must be the same as the matrix passed to
16 * cholmod_analyze for the supernodal case. For the simplicial case, it can
17 * be different, but it should be the same for best performance. beta is real.
18 *
19 * A simplicial factorization or supernodal factorization is chosen, based on
20 * the type of the factor L. If L->is_super is TRUE, a supernodal LL'
21 * factorization is computed. Otherwise, a simplicial numeric factorization
22 * is computed, either LL' or LDL', depending on Common->final_ll.
23 *
24 * Once the factorization is complete, it can be left as is or optionally
25 * converted into any simplicial numeric type, depending on the
26 * Common->final_* parameters. If converted from a supernodal to simplicial
27 * type, and the Common->final_resymbol parameter is true, then numerically
28 * zero entries in L due to relaxed supernodal amalgamation are removed from
29 * the simplicial factor (they are always left in the supernodal form of L).
30 * Entries that are numerically zero but present in the simplicial symbolic
31 * pattern of L are left in place (that is, the graph of L remains chordal).
32 * This is required for the update/downdate/rowadd/rowdel routines to work
33 * properly.
34 *
35 * workspace: Flag (nrow), Head (nrow+1),
36 * if symmetric: Iwork (2*nrow+2*nsuper)
37 * if unsymmetric: Iwork (2*nrow+MAX(2*nsuper,ncol))
38 * where nsuper is 0 if simplicial, or the # of relaxed supernodes in
39 * L otherwise (nsuper <= nrow).
40 * if simplicial: W (nrow).
41 * Allocates up to two temporary copies of its input matrix (including
42 * both pattern and numerical values).
43 *
44 * If the matrix is not positive definite the routine returns TRUE, but
45 * sets Common->status to CHOLMOD_NOT_POSDEF and L->minor is set to the
46 * column at which the failure occurred. Columns L->minor to L->n-1 are
47 * set to zero.
48 *
49 * Supports any xtype (pattern, real, complex, or zomplex), except that the
50 * input matrix A cannot be pattern-only. If L is simplicial, its numeric
51 * xtype matches A on output. If L is supernodal, its xtype is real if A is
52 * real, or complex if A is complex or zomplex.
53 */
54
55 #ifndef NCHOLESKY
56
57 #include "cholmod_internal.h"
58 #include "cholmod_cholesky.h"
59
60 #ifndef NSUPERNODAL
61 #include "cholmod_supernodal.h"
62 #endif
63
64
65 /* ========================================================================== */
66 /* === cholmod_factorize ==================================================== */
67 /* ========================================================================== */
68
69 /* Factorizes PAP' (or PAA'P' if A->stype is 0), using a factor obtained
70 * from cholmod_analyze. The analysis can be re-used simply by calling this
71 * routine a second time with another matrix. A must have the same nonzero
72 * pattern as that passed to cholmod_analyze. */
73
CHOLMOD(factorize)74 int CHOLMOD(factorize)
75 (
76 /* ---- input ---- */
77 cholmod_sparse *A, /* matrix to factorize */
78 /* ---- in/out --- */
79 cholmod_factor *L, /* resulting factorization */
80 /* --------------- */
81 cholmod_common *Common
82 )
83 {
84 double zero [2] ;
85 zero [0] = 0 ;
86 zero [1] = 0 ;
87 return (CHOLMOD(factorize_p) (A, zero, NULL, 0, L, Common)) ;
88 }
89
90
91 /* ========================================================================== */
92 /* === cholmod_factorize_p ================================================== */
93 /* ========================================================================== */
94
95 /* Same as cholmod_factorize, but with more options. */
96
CHOLMOD(factorize_p)97 int CHOLMOD(factorize_p)
98 (
99 /* ---- input ---- */
100 cholmod_sparse *A, /* matrix to factorize */
101 double beta [2], /* factorize beta*I+A or beta*I+A'*A */
102 Int *fset, /* subset of 0:(A->ncol)-1 */
103 size_t fsize, /* size of fset */
104 /* ---- in/out --- */
105 cholmod_factor *L, /* resulting factorization */
106 /* --------------- */
107 cholmod_common *Common
108 )
109 {
110 cholmod_sparse *S, *F, *A1, *A2 ;
111 Int nrow, ncol, stype, convert, n, nsuper, grow2, status ;
112 size_t s, t, uncol ;
113 int ok = TRUE ;
114
115 /* ---------------------------------------------------------------------- */
116 /* check inputs */
117 /* ---------------------------------------------------------------------- */
118
119 RETURN_IF_NULL_COMMON (FALSE) ;
120 RETURN_IF_NULL (A, FALSE) ;
121 RETURN_IF_NULL (L, FALSE) ;
122 RETURN_IF_XTYPE_INVALID (A, CHOLMOD_REAL, CHOLMOD_ZOMPLEX, FALSE) ;
123 RETURN_IF_XTYPE_INVALID (L, CHOLMOD_PATTERN, CHOLMOD_ZOMPLEX, FALSE) ;
124 nrow = A->nrow ;
125 ncol = A->ncol ;
126 n = L->n ;
127 stype = A->stype ;
128 if (L->n != A->nrow)
129 {
130 ERROR (CHOLMOD_INVALID, "A and L dimensions do not match") ;
131 return (FALSE) ;
132 }
133 if (stype != 0 && nrow != ncol)
134 {
135 ERROR (CHOLMOD_INVALID, "matrix invalid") ;
136 return (FALSE) ;
137 }
138 DEBUG (CHOLMOD(dump_sparse) (A, "A for cholmod_factorize", Common)) ;
139 Common->status = CHOLMOD_OK ;
140
141 /* ---------------------------------------------------------------------- */
142 /* allocate workspace */
143 /* ---------------------------------------------------------------------- */
144
145 nsuper = (L->is_super ? L->nsuper : 0) ;
146 uncol = ((stype != 0) ? 0 : ncol) ;
147
148 /* s = 2*nrow + MAX (uncol, 2*nsuper) */
149 s = CHOLMOD(mult_size_t) (nsuper, 2, &ok) ;
150 s = MAX (uncol, s) ;
151 t = CHOLMOD(mult_size_t) (nrow, 2, &ok) ;
152 s = CHOLMOD(add_size_t) (s, t, &ok) ;
153 if (!ok)
154 {
155 ERROR (CHOLMOD_TOO_LARGE, "problem too large") ;
156 return (FALSE) ;
157 }
158
159 CHOLMOD(allocate_work) (nrow, s, 0, Common) ;
160 if (Common->status < CHOLMOD_OK)
161 {
162 return (FALSE) ;
163 }
164
165 S = NULL ;
166 F = NULL ;
167 A1 = NULL ;
168 A2 = NULL ;
169
170 /* convert to another form when done, if requested */
171 convert = !(Common->final_asis) ;
172
173 /* ---------------------------------------------------------------------- */
174 /* perform supernodal LL' or simplicial LDL' factorization */
175 /* ---------------------------------------------------------------------- */
176
177 if (L->is_super)
178 {
179
180 #ifndef NSUPERNODAL
181
182 /* ------------------------------------------------------------------ */
183 /* supernodal factorization */
184 /* ------------------------------------------------------------------ */
185
186 if (L->ordering == CHOLMOD_NATURAL)
187 {
188
189 /* -------------------------------------------------------------- */
190 /* natural ordering */
191 /* -------------------------------------------------------------- */
192
193 if (stype > 0)
194 {
195 /* S = tril (A'), F not needed */
196 /* workspace: Iwork (nrow) */
197 A1 = CHOLMOD(ptranspose) (A, 2, NULL, NULL, 0, Common) ;
198 S = A1 ;
199 }
200 else if (stype < 0)
201 {
202 /* This is the fastest option for the natural ordering */
203 /* S = A; F not needed */
204 S = A ;
205 }
206 else
207 {
208 /* F = A(:,f)' */
209 /* workspace: Iwork (nrow) */
210 /* workspace: Iwork (nrow if no fset; MAX (nrow,ncol) if fset)*/
211 A1 = CHOLMOD(ptranspose) (A, 2, NULL, fset, fsize, Common) ;
212 F = A1 ;
213 /* S = A */
214 S = A ;
215 }
216
217 }
218 else
219 {
220
221 /* -------------------------------------------------------------- */
222 /* permute the input matrix before factorization */
223 /* -------------------------------------------------------------- */
224
225 if (stype > 0)
226 {
227 /* This is the fastest option for factoring a permuted matrix */
228 /* S = tril (PAP'); F not needed */
229 /* workspace: Iwork (2*nrow) */
230 A1 = CHOLMOD(ptranspose) (A, 2, L->Perm, NULL, 0, Common) ;
231 S = A1 ;
232 }
233 else if (stype < 0)
234 {
235 /* A2 = triu (PAP') */
236 /* workspace: Iwork (2*nrow) */
237 A2 = CHOLMOD(ptranspose) (A, 2, L->Perm, NULL, 0, Common) ;
238 /* S = tril (A2'); F not needed */
239 /* workspace: Iwork (nrow) */
240 A1 = CHOLMOD(ptranspose) (A2, 2, NULL, NULL, 0, Common) ;
241 S = A1 ;
242 CHOLMOD(free_sparse) (&A2, Common) ;
243 ASSERT (A2 == NULL) ;
244 }
245 else
246 {
247 /* F = A(p,f)' */
248 /* workspace: Iwork (nrow if no fset; MAX (nrow,ncol) if fset)*/
249 A1 = CHOLMOD(ptranspose) (A, 2, L->Perm, fset, fsize, Common) ;
250 F = A1 ;
251 /* S = F' */
252 /* workspace: Iwork (nrow) */
253 A2 = CHOLMOD(ptranspose) (F, 2, NULL, NULL, 0, Common) ;
254 S = A2 ;
255 }
256 }
257
258 /* ------------------------------------------------------------------ */
259 /* supernodal factorization */
260 /* ------------------------------------------------------------------ */
261
262 /* workspace: Flag (nrow), Head (nrow+1), Iwork (2*nrow+2*nsuper) */
263 if (Common->status == CHOLMOD_OK)
264 {
265 CHOLMOD(super_numeric) (S, F, beta, L, Common) ;
266 }
267 status = Common->status ;
268 ASSERT (IMPLIES (status >= CHOLMOD_OK, L->xtype != CHOLMOD_PATTERN)) ;
269
270 /* ------------------------------------------------------------------ */
271 /* convert to final form, if requested */
272 /* ------------------------------------------------------------------ */
273
274 if (Common->status >= CHOLMOD_OK && convert)
275 {
276 /* workspace: none */
277 ok = CHOLMOD(change_factor) (L->xtype, Common->final_ll,
278 Common->final_super, Common->final_pack,
279 Common->final_monotonic, L, Common) ;
280 if (ok && Common->final_resymbol && !(L->is_super))
281 {
282 /* workspace: Flag (nrow), Head (nrow+1),
283 * if symmetric: Iwork (2*nrow)
284 * if unsymmetric: Iwork (2*nrow+ncol) */
285 CHOLMOD(resymbol_noperm) (S, fset, fsize, Common->final_pack,
286 L, Common) ;
287 }
288 }
289
290 #else
291
292 /* ------------------------------------------------------------------ */
293 /* CHOLMOD Supernodal module not installed */
294 /* ------------------------------------------------------------------ */
295
296 status = CHOLMOD_NOT_INSTALLED ;
297 ERROR (CHOLMOD_NOT_INSTALLED,"Supernodal module not installed") ;
298
299 #endif
300
301 }
302 else
303 {
304
305 /* ------------------------------------------------------------------ */
306 /* simplicial LDL' factorization */
307 /* ------------------------------------------------------------------ */
308
309 /* Permute the input matrix A if necessary. cholmod_rowfac requires
310 * triu(A) in column form for the symmetric case, and A in column form
311 * for the unsymmetric case (the matrix S). The unsymmetric case
312 * requires A in row form, or equivalently A' in column form (the
313 * matrix F).
314 */
315
316 if (L->ordering == CHOLMOD_NATURAL)
317 {
318
319 /* -------------------------------------------------------------- */
320 /* natural ordering */
321 /* -------------------------------------------------------------- */
322
323 if (stype > 0)
324 {
325 /* F is not needed, S = A */
326 S = A ;
327 }
328 else if (stype < 0)
329 {
330 /* F is not needed, S = A' */
331 /* workspace: Iwork (nrow) */
332 A2 = CHOLMOD(ptranspose) (A, 2, NULL, NULL, 0, Common) ;
333 S = A2 ;
334 }
335 else
336 {
337 /* F = A (:,f)' */
338 /* workspace: Iwork (nrow if no fset; MAX (nrow,ncol) if fset)*/
339 A1 = CHOLMOD(ptranspose) (A, 2, NULL, fset, fsize, Common) ;
340 F = A1 ;
341 S = A ;
342 }
343
344 }
345 else
346 {
347
348 /* -------------------------------------------------------------- */
349 /* permute the input matrix before factorization */
350 /* -------------------------------------------------------------- */
351
352 if (stype > 0)
353 {
354 /* F = tril (A (p,p)') */
355 /* workspace: Iwork (2*nrow) */
356 A1 = CHOLMOD(ptranspose) (A, 2, L->Perm, NULL, 0, Common) ;
357 /* A2 = triu (F') */
358 /* workspace: Iwork (nrow) */
359 A2 = CHOLMOD(ptranspose) (A1, 2, NULL, NULL, 0, Common) ;
360 /* the symmetric case does not need F, free it and set to NULL*/
361 CHOLMOD(free_sparse) (&A1, Common) ;
362 }
363 else if (stype < 0)
364 {
365 /* A2 = triu (A (p,p)'), F not needed. This is the fastest
366 * way to factorize a matrix using the simplicial routine
367 * (cholmod_rowfac). */
368 /* workspace: Iwork (2*nrow) */
369 A2 = CHOLMOD(ptranspose) (A, 2, L->Perm, NULL, 0, Common) ;
370 }
371 else
372 {
373 /* F = A (p,f)' */
374 /* workspace: Iwork (nrow if no fset; MAX (nrow,ncol) if fset)*/
375 A1 = CHOLMOD(ptranspose) (A, 2, L->Perm, fset, fsize, Common) ;
376 F = A1 ;
377 /* A2 = F' */
378 /* workspace: Iwork (nrow) */
379 A2 = CHOLMOD(ptranspose) (F, 2, NULL, NULL, 0, Common) ;
380 }
381 S = A2 ;
382 }
383
384 /* ------------------------------------------------------------------ */
385 /* simplicial LDL' or LL' factorization */
386 /* ------------------------------------------------------------------ */
387
388 /* factorize beta*I+S (symmetric) or beta*I+F*F' (unsymmetric) */
389 /* workspace: Flag (nrow), W (nrow), Iwork (2*nrow) */
390 if (Common->status == CHOLMOD_OK)
391 {
392 grow2 = Common->grow2 ;
393 L->is_ll = BOOLEAN (Common->final_ll) ;
394 if (L->xtype == CHOLMOD_PATTERN && Common->final_pack)
395 {
396 /* allocate a factor with exactly the space required */
397 Common->grow2 = 0 ;
398 }
399 CHOLMOD(rowfac) (S, F, beta, 0, nrow, L, Common) ;
400 Common->grow2 = grow2 ;
401 }
402 status = Common->status ;
403
404 /* ------------------------------------------------------------------ */
405 /* convert to final form, if requested */
406 /* ------------------------------------------------------------------ */
407
408 if (Common->status >= CHOLMOD_OK && convert)
409 {
410 /* workspace: none */
411 CHOLMOD(change_factor) (L->xtype, L->is_ll, FALSE,
412 Common->final_pack, Common->final_monotonic, L, Common) ;
413 }
414 }
415
416 /* ---------------------------------------------------------------------- */
417 /* free A1 and A2 if they exist */
418 /* ---------------------------------------------------------------------- */
419
420 CHOLMOD(free_sparse) (&A1, Common) ;
421 CHOLMOD(free_sparse) (&A2, Common) ;
422 Common->status = MAX (Common->status, status) ;
423 return (Common->status >= CHOLMOD_OK) ;
424 }
425 #endif
426