1 /* ========================================================================== */
2 /* === Cholesky/cholmod_rowfac ============================================== */
3 /* ========================================================================== */
4
5 /* -----------------------------------------------------------------------------
6 * CHOLMOD/Cholesky Module. Copyright (C) 2005-2013, Timothy A. Davis
7 * -------------------------------------------------------------------------- */
8
9 /* Full or incremental numerical LDL' or LL' factorization (simplicial, not
10 * supernodal) cholmod_factorize is the "easy" wrapper for this code, but it
11 * does not provide access to incremental factorization.
12 *
13 * cholmod_rowfac computes the full or incremental LDL' or LL' factorization of
14 * A+beta*I (where A is symmetric) or A*F+beta*I (where A and F are unsymmetric
15 * and only the upper triangular part of A*F+beta*I is used). It computes
16 * L (and D, for LDL') one row at a time. beta is real.
17 *
18 * A is nrow-by-ncol or nrow-by-nrow. In "packed" form it is a conventional
19 * column-oriented sparse matrix. Row indices of column j are in
20 * Ai [Ap [j] ... Ap [j+1]-1] and values in the same locations of Ax.
21 * will be faster if A has sorted columns. In "unpacked" form the column
22 * of A ends at Ap [j] + Anz [j] - 1 instead of Ap [j+1] - 1.
23 *
24 * Row indices in each column of A can be sorted or unsorted, but the routine
25 * routine works fastest if A is sorted, or if only triu(A) is provided
26 * for the symmetric case.
27 *
28 * The unit-diagonal nrow-by-nrow output matrix L is returned in "unpacked"
29 * column form, with row indices of column j in Li [Lp [j] ...
30 * Lp [j] + Lnz [j] - 1] and values in the same location in Lx. The row
31 * indices in each column of L are in sorted order. The unit diagonal of L
32 * is not stored.
33 *
34 * L can be a simplicial symbolic or numeric (L->is_super must be FALSE).
35 * A symbolic factor is converted immediately into a numeric factor containing
36 * the identity matrix.
37 *
38 * For a full factorization, kstart = 0 and kend = nrow. The existing nonzero
39 * entries (numerical values in L->x and L->z for the zomplex case, and indices
40 * in L->i), if any, are overwritten.
41 *
42 * To compute an incremental factorization, select kstart and kend as the range
43 * of rows of L you wish to compute. A correct factorization will be computed
44 * only if all descendants of all nodes k = kstart to kend-1 in the etree have
45 * been factorized by a prior call to this routine, and if rows kstart to kend-1
46 * have not been factorized. This condition is NOT checked on input.
47 *
48 * ---------------
49 * Symmetric case:
50 * ---------------
51 *
52 * The factorization (in MATLAB notation) is:
53 *
54 * S = beta*I + A
55 * S = triu (S) + triu (S,1)'
56 * L*D*L' = S, or L*L' = S
57 *
58 * A is a conventional sparse matrix in compressed column form. Only the
59 * diagonal and upper triangular part of A is accessed; the lower
60 * triangular part is ignored and assumed to be equal to the upper
61 * triangular part. For an incremental factorization, only columns kstart
62 * to kend-1 of A are accessed. F is not used.
63 *
64 * ---------------
65 * Unsymmetric case:
66 * ---------------
67 *
68 * The factorization (in MATLAB notation) is:
69 *
70 * S = beta*I + A*F
71 * S = triu (S) + triu (S,1)'
72 * L*D*L' = S, or L*L' = S
73 *
74 * The typical case is F=A'. Alternatively, if F=A(:,f)', then this
75 * routine factorizes S = beta*I + A(:,f)*A(:,f)'.
76 *
77 * All of A and F are accessed, but only the upper triangular part of A*F
78 * is used. F must be of size A->ncol by A->nrow. F is used for the
79 * unsymmetric case only. F can be packed or unpacked and it need not be
80 * sorted.
81 *
82 * For a complete factorization of beta*I + A*A',
83 * this routine performs a number of flops exactly equal to:
84 *
85 * sum (for each column j of A) of (Anz (j)^2 + Anz (j)), to form S
86 * +
87 * sum (for each column j of L) of (Lnz (j)^2 + 3*Lnz (j)), to factorize S
88 *
89 * where Anz (j) is the number of nonzeros in column j of A, and Lnz (j)
90 * is the number of nonzero in column j of L below the diagonal.
91 *
92 *
93 * workspace: Flag (nrow), W (nrow if real, 2*nrow if complex/zomplex),
94 * Iwork (nrow)
95 *
96 * Supports any xtype, except a pattern-only input matrix A cannot be
97 * factorized.
98 */
99
100 #ifndef NCHOLESKY
101
102 #include "cholmod_internal.h"
103 #include "cholmod_cholesky.h"
104
105 /* ========================================================================== */
106 /* === subtree ============================================================== */
107 /* ========================================================================== */
108
109 /* Compute the nonzero pattern of the sparse triangular solve Lx=b, where L in
110 * this case is L(0:k-1,0:k-1), and b is a column of A. This is done by
111 * traversing the kth row-subtree of the elimination tree of L, starting from
112 * each nonzero entry in b. The pattern is returned postordered, and is valid
113 * for a subsequent numerical triangular solve of Lx=b. The elimination tree
114 * can be provided in a Parent array, or extracted from the pattern of L itself.
115 *
116 * The pattern of x = inv(L)*b is returned in Stack [top...].
117 * Also scatters b, or a multiple of b, into the work vector W.
118 *
119 * The SCATTER macro is defines how the numerical values of A or A*A' are to be
120 * scattered.
121 *
122 * PARENT(i) is a macro the defines how the etree is accessed. It is either:
123 * #define PARENT(i) Parent [i]
124 * #define PARENT(i) (Lnz [i] > 1) ? (Li [Lp [i] + 1]) : EMPTY
125 */
126
127 #define SUBTREE \
128 for ( ; p < pend ; p++) \
129 { \
130 i = Ai [p] ; \
131 if (i <= k) \
132 { \
133 /* scatter the column of A, or A*A' into Wx and Wz */ \
134 SCATTER ; \
135 /* start at node i and traverse up the subtree, stop at node k */ \
136 for (len = 0 ; i < k && i != EMPTY && Flag [i] < mark ; i = parent) \
137 { \
138 /* L(k,i) is nonzero, and seen for the first time */ \
139 Stack [len++] = i ; /* place i on the stack */ \
140 Flag [i] = mark ; /* mark i as visited */ \
141 parent = PARENT (i) ; /* traverse up the etree to the parent */ \
142 } \
143 /* move the path down to the bottom of the stack */ \
144 while (len > 0) \
145 { \
146 Stack [--top] = Stack [--len] ; \
147 } \
148 } \
149 else if (sorted) \
150 { \
151 break ; \
152 } \
153 }
154
155
156 /* ========================================================================== */
157 /* === TEMPLATE ============================================================= */
158 /* ========================================================================== */
159
160 #define REAL
161 #include "t_cholmod_rowfac.c"
162 #define COMPLEX
163 #include "t_cholmod_rowfac.c"
164 #define ZOMPLEX
165 #include "t_cholmod_rowfac.c"
166
167 #define MASK
168 #define REAL
169 #include "t_cholmod_rowfac.c"
170 #define COMPLEX
171 #include "t_cholmod_rowfac.c"
172 #define ZOMPLEX
173 #include "t_cholmod_rowfac.c"
174 #undef MASK
175
176
177 /* ========================================================================== */
178 /* === cholmod_row_subtree ================================================== */
179 /* ========================================================================== */
180
181 /* Compute the nonzero pattern of the solution to the lower triangular system
182 * L(0:k-1,0:k-1) * x = A (0:k-1,k) if A is symmetric, or
183 * L(0:k-1,0:k-1) * x = A (0:k-1,:) * A (:,k)' if A is unsymmetric.
184 * This gives the nonzero pattern of row k of L (excluding the diagonal).
185 * The pattern is returned postordered.
186 *
187 * The symmetric case requires A to be in symmetric-upper form.
188 *
189 * The result is returned in R, a pre-allocated sparse matrix of size nrow-by-1,
190 * with R->nzmax >= nrow. R is assumed to be packed (Rnz [0] is not updated);
191 * the number of entries in R is given by Rp [0].
192 *
193 * FUTURE WORK: a very minor change to this routine could allow it to compute
194 * the nonzero pattern of x for any system Lx=b. The SUBTREE macro would need
195 * to change, to eliminate its dependence on k.
196 *
197 * workspace: Flag (nrow)
198 */
199
CHOLMOD(row_subtree)200 int CHOLMOD(row_subtree)
201 (
202 /* ---- input ---- */
203 cholmod_sparse *A, /* matrix to analyze */
204 cholmod_sparse *F, /* used for A*A' case only. F=A' or A(:,f)' */
205 size_t krow, /* row k of L */
206 Int *Parent, /* elimination tree */
207 /* ---- output --- */
208 cholmod_sparse *R, /* pattern of L(k,:), 1-by-n with R->nzmax >= n */
209 /* --------------- */
210 cholmod_common *Common
211 )
212 {
213 Int *Rp, *Stack, *Flag, *Ap, *Ai, *Anz, *Fp, *Fi, *Fnz ;
214 Int p, pend, parent, t, stype, nrow, k, pf, pfend, Fpacked, packed,
215 sorted, top, len, i, mark ;
216
217 /* ---------------------------------------------------------------------- */
218 /* check inputs */
219 /* ---------------------------------------------------------------------- */
220
221 RETURN_IF_NULL_COMMON (FALSE) ;
222 RETURN_IF_NULL (A, FALSE) ;
223 RETURN_IF_NULL (R, FALSE) ;
224 RETURN_IF_NULL (Parent, FALSE) ;
225 RETURN_IF_XTYPE_INVALID (A, CHOLMOD_PATTERN, CHOLMOD_ZOMPLEX, FALSE) ;
226 RETURN_IF_XTYPE_INVALID (R, CHOLMOD_PATTERN, CHOLMOD_ZOMPLEX, FALSE) ;
227 stype = A->stype ;
228 if (stype == 0)
229 {
230 RETURN_IF_NULL (F, FALSE) ;
231 RETURN_IF_XTYPE_INVALID (F, CHOLMOD_PATTERN, CHOLMOD_ZOMPLEX, FALSE) ;
232 }
233 if (krow >= A->nrow)
234 {
235 ERROR (CHOLMOD_INVALID, "subtree: k invalid") ;
236 return (FALSE) ;
237 }
238 if (R->ncol != 1 || A->nrow != R->nrow || A->nrow > R->nzmax)
239 {
240 ERROR (CHOLMOD_INVALID, "subtree: R invalid") ;
241 return (FALSE) ;
242 }
243 Common->status = CHOLMOD_OK ;
244
245 /* ---------------------------------------------------------------------- */
246 /* allocate workspace */
247 /* ---------------------------------------------------------------------- */
248
249 nrow = A->nrow ;
250 CHOLMOD(allocate_work) (nrow, 0, 0, Common) ;
251 if (Common->status < CHOLMOD_OK)
252 {
253 return (FALSE) ;
254 }
255 ASSERT (CHOLMOD(dump_work) (TRUE, TRUE, 0, Common)) ;
256
257 /* ---------------------------------------------------------------------- */
258 /* get inputs */
259 /* ---------------------------------------------------------------------- */
260
261 if (stype > 0)
262 {
263 /* symmetric upper case: F is not needed. It may be NULL */
264 Fp = NULL ;
265 Fi = NULL ;
266 Fnz = NULL ;
267 Fpacked = TRUE ;
268 }
269 else if (stype == 0)
270 {
271 /* unsymmetric case: F is required. */
272 Fp = F->p ;
273 Fi = F->i ;
274 Fnz = F->nz ;
275 Fpacked = F->packed ;
276 }
277 else
278 {
279 /* symmetric lower triangular form not supported */
280 ERROR (CHOLMOD_INVALID, "symmetric lower not supported") ;
281 return (FALSE) ;
282 }
283
284 Ap = A->p ;
285 Ai = A->i ;
286 Anz = A->nz ;
287 packed = A->packed ;
288 sorted = A->sorted ;
289
290 k = krow ;
291 Stack = R->i ;
292
293 /* ---------------------------------------------------------------------- */
294 /* get workspace */
295 /* ---------------------------------------------------------------------- */
296
297 Flag = Common->Flag ; /* size nrow, Flag [i] < mark must hold */
298 /* mark = CHOLMOD(clear_flag) (Common) ; */
299 CHOLMOD_CLEAR_FLAG (Common) ;
300 mark = Common->mark ;
301
302 /* ---------------------------------------------------------------------- */
303 /* compute the pattern of L(k,:) */
304 /* ---------------------------------------------------------------------- */
305
306 top = nrow ; /* Stack is empty */
307 Flag [k] = mark ; /* do not include diagonal entry in Stack */
308
309 #define SCATTER /* do not scatter numerical values */
310 #define PARENT(i) Parent [i] /* use Parent for etree */
311
312 if (stype != 0)
313 {
314 /* scatter kth col of triu (A), get pattern L(k,:) */
315 p = Ap [k] ;
316 pend = (packed) ? (Ap [k+1]) : (p + Anz [k]) ;
317 SUBTREE ;
318 }
319 else
320 {
321 /* scatter kth col of triu (beta*I+AA'), get pattern L(k,:) */
322 pf = Fp [k] ;
323 pfend = (Fpacked) ? (Fp [k+1]) : (pf + Fnz [k]) ;
324 for ( ; pf < pfend ; pf++)
325 {
326 /* get nonzero entry F (t,k) */
327 t = Fi [pf] ;
328 p = Ap [t] ;
329 pend = (packed) ? (Ap [t+1]) : (p + Anz [t]) ;
330 SUBTREE ;
331 }
332 }
333
334 #undef SCATTER
335 #undef PARENT
336
337 /* shift the stack upwards, to the first part of R */
338 len = nrow - top ;
339 for (i = 0 ; i < len ; i++)
340 {
341 Stack [i] = Stack [top + i] ;
342 }
343
344 Rp = R->p ;
345 Rp [0] = 0 ;
346 Rp [1] = len ;
347 R->sorted = FALSE ;
348
349 CHOLMOD(clear_flag) (Common) ;
350 ASSERT (CHOLMOD(dump_work) (TRUE, TRUE, 0, Common)) ;
351 return (TRUE) ;
352 }
353
354
355 /* ========================================================================== */
356 /* === cholmod_lsolve_pattern =============================================== */
357 /* ========================================================================== */
358
359 /* Compute the nonzero pattern of Y=L\B. L must be simplicial, and B
360 * must be a single sparse column vector with B->stype = 0. The values of
361 * B are not used; it just specifies a nonzero pattern. The pattern of
362 * Y is not sorted, but is in topological order instead (suitable for a
363 * sparse forward/backsolve).
364 */
365
CHOLMOD(lsolve_pattern)366 int CHOLMOD(lsolve_pattern)
367 (
368 /* ---- input ---- */
369 cholmod_sparse *B, /* sparse right-hand-side (a single sparse column) */
370 cholmod_factor *L, /* the factor L from which parent(i) is derived */
371 /* ---- output --- */
372 cholmod_sparse *Yset, /* pattern of Y=L\B, n-by-1 with Y->nzmax >= n */
373 /* --------------- */
374 cholmod_common *Common
375 )
376 {
377 size_t krow ;
378 RETURN_IF_NULL (B, FALSE) ;
379 krow = B->nrow ;
380 return (CHOLMOD(row_lsubtree) (B, NULL, 0, krow, L, Yset, Common)) ;
381 }
382
383
384 /* ========================================================================== */
385 /* === cholmod_row_lsubtree ================================================= */
386 /* ========================================================================== */
387
388 /* Identical to cholmod_row_subtree, except that the elimination tree is
389 * obtained from L itself, as the first off-diagonal entry in each column.
390 * L must be simplicial, not supernodal.
391 *
392 * If krow = A->nrow, then A must be a single sparse column vector, (A->stype
393 * must be zero), and the nonzero pattern of x=L\b is computed, where b=A(:,0)
394 * is the single sparse right-hand-side. The inputs Fi and fnz are ignored.
395 * See CHOLMOD(lsolve_pattern) above for a simpler interface for this case.
396 */
397
CHOLMOD(row_lsubtree)398 int CHOLMOD(row_lsubtree)
399 (
400 /* ---- input ---- */
401 cholmod_sparse *A, /* matrix to analyze */
402 Int *Fi, size_t fnz, /* nonzero pattern of kth row of A', not required
403 * for the symmetric case. Need not be sorted. */
404 size_t krow, /* row k of L */
405 cholmod_factor *L, /* the factor L from which parent(i) is derived */
406 /* ---- output --- */
407 cholmod_sparse *R, /* pattern of L(k,:), n-by-1 with R->nzmax >= n */
408 /* --------------- */
409 cholmod_common *Common
410 )
411 {
412 Int *Rp, *Stack, *Flag, *Ap, *Ai, *Anz, *Lp, *Li, *Lnz ;
413 Int p, pend, parent, t, stype, nrow, k, pf, packed, sorted, top, len, i,
414 mark, ka ;
415
416 /* ---------------------------------------------------------------------- */
417 /* check inputs */
418 /* ---------------------------------------------------------------------- */
419
420 RETURN_IF_NULL_COMMON (FALSE) ;
421 RETURN_IF_NULL (A, FALSE) ;
422 RETURN_IF_NULL (R, FALSE) ;
423 RETURN_IF_NULL (L, FALSE) ;
424 RETURN_IF_XTYPE_INVALID (A, CHOLMOD_PATTERN, CHOLMOD_ZOMPLEX, FALSE) ;
425 RETURN_IF_XTYPE_INVALID (R, CHOLMOD_PATTERN, CHOLMOD_ZOMPLEX, FALSE) ;
426 RETURN_IF_XTYPE_INVALID (L, CHOLMOD_REAL, CHOLMOD_ZOMPLEX, FALSE) ;
427
428 nrow = A->nrow ;
429 stype = A->stype ;
430 if (stype < 0)
431 {
432 /* symmetric lower triangular form not supported */
433 ERROR (CHOLMOD_INVALID, "symmetric lower not supported") ;
434 return (FALSE) ;
435 }
436
437 if (krow > nrow)
438 {
439 ERROR (CHOLMOD_INVALID, "lsubtree: krow invalid") ;
440 return (FALSE) ;
441 }
442 else if (krow == nrow)
443 {
444 /* find pattern of x=L\b where b=A(:,0) */
445 k = nrow ; /* compute all of the result; don't stop in SUBTREE */
446 ka = 0 ; /* use column A(:,0) */
447 if (stype != 0 || A->ncol != 1)
448 {
449 /* A must be unsymmetric (it's a single sparse column vector) */
450 ERROR (CHOLMOD_INVALID, "lsubtree: A invalid") ;
451 return (FALSE) ;
452 }
453 }
454 else
455 {
456 /* find pattern of L(k,:) using A(:,k) and Fi if A unsymmetric */
457 k = krow ; /* which row of L to compute */
458 ka = k ; /* which column of A to use */
459 if (stype == 0)
460 {
461 RETURN_IF_NULL (Fi, FALSE) ;
462 }
463 }
464
465 if (R->ncol != 1 || nrow != R->nrow || nrow > R->nzmax ||
466 ((krow == nrow || stype != 0) && ka >= A->ncol))
467 {
468 ERROR (CHOLMOD_INVALID, "lsubtree: R invalid") ;
469 return (FALSE) ;
470 }
471 if (L->is_super)
472 {
473 ERROR (CHOLMOD_INVALID, "lsubtree: L invalid (cannot be supernodal)") ;
474 return (FALSE) ;
475 }
476 Common->status = CHOLMOD_OK ;
477
478 /* ---------------------------------------------------------------------- */
479 /* allocate workspace */
480 /* ---------------------------------------------------------------------- */
481
482 CHOLMOD(allocate_work) (nrow, 0, 0, Common) ;
483 if (Common->status < CHOLMOD_OK)
484 {
485 return (FALSE) ;
486 }
487 ASSERT (CHOLMOD(dump_work) (TRUE, TRUE, 0, Common)) ;
488
489 /* ---------------------------------------------------------------------- */
490 /* get inputs */
491 /* ---------------------------------------------------------------------- */
492
493 Ap = A->p ;
494 Ai = A->i ;
495 Anz = A->nz ;
496 packed = A->packed ;
497 sorted = A->sorted ;
498
499 Stack = R->i ;
500
501 Lp = L->p ;
502 Li = L->i ;
503 Lnz = L->nz ;
504
505 /* ---------------------------------------------------------------------- */
506 /* get workspace */
507 /* ---------------------------------------------------------------------- */
508
509 Flag = Common->Flag ; /* size nrow, Flag [i] < mark must hold */
510 mark = CHOLMOD(clear_flag) (Common) ;
511
512 /* ---------------------------------------------------------------------- */
513 /* compute the pattern of L(k,:) */
514 /* ---------------------------------------------------------------------- */
515
516 top = nrow ; /* Stack is empty */
517 if (k < nrow)
518 {
519 Flag [k] = mark ; /* do not include diagonal entry in Stack */
520 }
521
522 #define SCATTER /* do not scatter numerical values */
523 #define PARENT(i) (Lnz [i] > 1) ? (Li [Lp [i] + 1]) : EMPTY
524
525 if (krow == nrow || stype != 0)
526 {
527 /* scatter kth col of triu (A), get pattern L(k,:) */
528 p = Ap [ka] ;
529 pend = (packed) ? (Ap [ka+1]) : (p + Anz [ka]) ;
530 SUBTREE ;
531 }
532 else
533 {
534 /* scatter kth col of triu (beta*I+AA'), get pattern L(k,:) */
535 for (pf = 0 ; pf < (Int) fnz ; pf++)
536 {
537 /* get nonzero entry F (t,k) */
538 t = Fi [pf] ;
539 p = Ap [t] ;
540 pend = (packed) ? (Ap [t+1]) : (p + Anz [t]) ;
541 SUBTREE ;
542 }
543 }
544
545 #undef SCATTER
546 #undef PARENT
547
548 /* shift the stack upwards, to the first part of R */
549 len = nrow - top ;
550 for (i = 0 ; i < len ; i++)
551 {
552 Stack [i] = Stack [top + i] ;
553 }
554
555 Rp = R->p ;
556 Rp [0] = 0 ;
557 Rp [1] = len ;
558 R->sorted = FALSE ;
559
560 CHOLMOD(clear_flag) (Common) ;
561 ASSERT (CHOLMOD(dump_work) (TRUE, TRUE, 0, Common)) ;
562 return (TRUE) ;
563 }
564
565
566 /* ========================================================================== */
567 /* === cholmod_rowfac ======================================================= */
568 /* ========================================================================== */
569
570 /* This is the incremental factorization for general purpose usage. */
571
CHOLMOD(rowfac)572 int CHOLMOD(rowfac)
573 (
574 /* ---- input ---- */
575 cholmod_sparse *A, /* matrix to factorize */
576 cholmod_sparse *F, /* used for A*A' case only. F=A' or A(:,f)' */
577 double beta [2], /* factorize beta*I+A or beta*I+AA' */
578 size_t kstart, /* first row to factorize */
579 size_t kend, /* last row to factorize is kend-1 */
580 /* ---- in/out --- */
581 cholmod_factor *L,
582 /* --------------- */
583 cholmod_common *Common
584 )
585 {
586 return (CHOLMOD(rowfac_mask2) (A, F, beta, kstart, kend, NULL, 0, NULL, L,
587 Common)) ;
588 }
589
590
591 /* ========================================================================== */
592 /* === cholmod_rowfac_mask ================================================== */
593 /* ========================================================================== */
594
595 /* This is meant for use in LPDASA only. */
596
CHOLMOD(rowfac_mask)597 int CHOLMOD(rowfac_mask)
598 (
599 /* ---- input ---- */
600 cholmod_sparse *A, /* matrix to factorize */
601 cholmod_sparse *F, /* used for A*A' case only. F=A' or A(:,f)' */
602 double beta [2], /* factorize beta*I+A or beta*I+AA' */
603 size_t kstart, /* first row to factorize */
604 size_t kend, /* last row to factorize is kend-1 */
605 Int *mask, /* size A->nrow. if mask[i] >= 0 row i is set to zero */
606 Int *RLinkUp, /* size A->nrow. link list of rows to compute */
607 /* ---- in/out --- */
608 cholmod_factor *L,
609 /* --------------- */
610 cholmod_common *Common
611 )
612 {
613 Int maskmark = 0 ;
614 return (CHOLMOD(rowfac_mask2) (A, F, beta, kstart, kend, mask, maskmark,
615 RLinkUp, L, Common)) ;
616 }
617
618 /* ========================================================================== */
619 /* === cholmod_rowfac_mask2 ================================================= */
620 /* ========================================================================== */
621
622 /* This is meant for use in LPDASA only. */
623
CHOLMOD(rowfac_mask2)624 int CHOLMOD(rowfac_mask2)
625 (
626 /* ---- input ---- */
627 cholmod_sparse *A, /* matrix to factorize */
628 cholmod_sparse *F, /* used for A*A' case only. F=A' or A(:,f)' */
629 double beta [2], /* factorize beta*I+A or beta*I+AA' */
630 size_t kstart, /* first row to factorize */
631 size_t kend, /* last row to factorize is kend-1 */
632 Int *mask, /* size A->nrow. if mask[i] >= maskmark row i is set
633 to zero */
634 Int maskmark, /* for mask [i] test */
635 Int *RLinkUp, /* size A->nrow. link list of rows to compute */
636 /* ---- in/out --- */
637 cholmod_factor *L,
638 /* --------------- */
639 cholmod_common *Common
640 )
641 {
642 Int n ;
643 size_t s ;
644 int ok = TRUE ;
645
646 /* ---------------------------------------------------------------------- */
647 /* check inputs */
648 /* ---------------------------------------------------------------------- */
649
650 RETURN_IF_NULL_COMMON (FALSE) ;
651 RETURN_IF_NULL (A, FALSE) ;
652 RETURN_IF_NULL (L, FALSE) ;
653 RETURN_IF_XTYPE_INVALID (A, CHOLMOD_REAL, CHOLMOD_ZOMPLEX, FALSE) ;
654 RETURN_IF_XTYPE_INVALID (L, CHOLMOD_PATTERN, CHOLMOD_ZOMPLEX, FALSE) ;
655 if (L->xtype != CHOLMOD_PATTERN && A->xtype != L->xtype)
656 {
657 ERROR (CHOLMOD_INVALID, "xtype of A and L do not match") ;
658 return (FALSE) ;
659 }
660 if (L->is_super)
661 {
662 ERROR (CHOLMOD_INVALID, "can only do simplicial factorization");
663 return (FALSE) ;
664 }
665 if (A->stype == 0)
666 {
667 RETURN_IF_NULL (F, FALSE) ;
668 if (A->xtype != F->xtype)
669 {
670 ERROR (CHOLMOD_INVALID, "xtype of A and F do not match") ;
671 return (FALSE) ;
672 }
673 }
674 if (A->stype < 0)
675 {
676 /* symmetric lower triangular form not supported */
677 ERROR (CHOLMOD_INVALID, "symmetric lower not supported") ;
678 return (FALSE) ;
679 }
680 if (kend > L->n)
681 {
682 ERROR (CHOLMOD_INVALID, "kend invalid") ;
683 return (FALSE) ;
684 }
685 if (A->nrow != L->n)
686 {
687 ERROR (CHOLMOD_INVALID, "dimensions of A and L do not match") ;
688 return (FALSE) ;
689 }
690 Common->status = CHOLMOD_OK ;
691 Common->rowfacfl = 0 ;
692
693 /* ---------------------------------------------------------------------- */
694 /* allocate workspace */
695 /* ---------------------------------------------------------------------- */
696
697 /* Xwork is of size n for the real case, 2*n for complex/zomplex */
698 n = L->n ;
699
700 /* s = ((A->xtype != CHOLMOD_REAL) ? 2:1)*n */
701 s = CHOLMOD(mult_size_t) (n, ((A->xtype != CHOLMOD_REAL) ? 2:1), &ok) ;
702 if (!ok)
703 {
704 ERROR (CHOLMOD_TOO_LARGE, "problem too large") ;
705 return (FALSE) ;
706 }
707
708 CHOLMOD(allocate_work) (n, n, s, Common) ;
709 if (Common->status < CHOLMOD_OK)
710 {
711 return (FALSE) ;
712 }
713 ASSERT (CHOLMOD(dump_work) (TRUE, TRUE, A->nrow, Common)) ;
714
715 /* ---------------------------------------------------------------------- */
716 /* factorize the matrix, using template routine */
717 /* ---------------------------------------------------------------------- */
718
719 if (RLinkUp == NULL)
720 {
721
722 switch (A->xtype)
723 {
724 case CHOLMOD_REAL:
725 ok = r_cholmod_rowfac (A, F, beta, kstart, kend, L, Common) ;
726 break ;
727
728 case CHOLMOD_COMPLEX:
729 ok = c_cholmod_rowfac (A, F, beta, kstart, kend, L, Common) ;
730 break ;
731
732 case CHOLMOD_ZOMPLEX:
733 ok = z_cholmod_rowfac (A, F, beta, kstart, kend, L, Common) ;
734 break ;
735 }
736
737 }
738 else
739 {
740
741 switch (A->xtype)
742 {
743 case CHOLMOD_REAL:
744 ok = r_cholmod_rowfac_mask (A, F, beta, kstart, kend,
745 mask, maskmark, RLinkUp, L, Common) ;
746 break ;
747
748 case CHOLMOD_COMPLEX:
749 ok = c_cholmod_rowfac_mask (A, F, beta, kstart, kend,
750 mask, maskmark, RLinkUp, L, Common) ;
751 break ;
752
753 case CHOLMOD_ZOMPLEX:
754 ok = z_cholmod_rowfac_mask (A, F, beta, kstart, kend,
755 mask, maskmark, RLinkUp, L, Common) ;
756 break ;
757 }
758 }
759
760 return (ok) ;
761 }
762 #endif
763