1 /* ========================================================================== */
2 /* === Modify/cholmod_updown ================================================ */
3 /* ========================================================================== */
4
5 /* -----------------------------------------------------------------------------
6 * CHOLMOD/Modify Module.
7 * Copyright (C) 2005-2006, Timothy A. Davis and William W. Hager.
8 * The CHOLMOD/Modify Module is licensed under Version 2.0 of the GNU
9 * General Public License. See gpl.txt for a text of the license.
10 * CHOLMOD is also available under other licenses; contact authors for details.
11 * http://www.suitesparse.com
12 * -------------------------------------------------------------------------- */
13
14 /* Updates/downdates the LDL' factorization (symbolic, then numeric), by
15 * computing a new factorization of
16 *
17 * Lnew * Dnew * Lnew' = Lold * Dold * Lold' +/- C*C'
18 *
19 * C must be sorted. It can be either packed or unpacked. As in all CHOLMOD
20 * routines, the columns of L are sorted on input, and also on output.
21 *
22 * If the factor is not an unpacked LDL' or dynamic LDL', it is converted
23 * to an LDL' dynamic factor. An unpacked LDL' factor may be updated, but if
24 * any one column runs out of space, the factor is converted to an LDL'
25 * dynamic one. If the initial conversion fails, the factor is returned
26 * unchanged.
27 *
28 * If memory runs out during the update, the factor is returned as a simplicial
29 * symbolic factor. That is, everything is freed except for the fill-reducing
30 * ordering and its corresponding column counts (typically computed by
31 * cholmod_analyze).
32 *
33 * Note that the fill-reducing permutation L->Perm is NOT used. The row
34 * indices of C refer to the rows of L, not A. If your original system is
35 * LDL' = PAP' (where P = L->Perm), and you want to compute the LDL'
36 * factorization of A+CC', then you must permute C first. That is:
37 *
38 * PAP' = LDL'
39 * P(A+CC')P' = PAP'+PCC'P' = LDL' + (PC)(PC)' = LDL' + Cnew*Cnew'
40 * where Cnew = P*C.
41 *
42 * You can use the cholmod_submatrix routine in the MatrixOps module
43 * to permute C, with:
44 *
45 * Cnew = cholmod_submatrix (C, L->Perm, L->n, NULL, -1, TRUE, TRUE, Common) ;
46 *
47 * Note that the sorted input parameter to cholmod_submatrix must be TRUE,
48 * because cholmod_updown requires C with sorted columns.
49 *
50 * The system Lx=b can also be updated/downdated. The old system was Lold*x=b.
51 * The new system is Lnew*xnew = b + deltab. The old solution x is overwritten
52 * with xnew. Note that as in the update/downdate of L itself, the fill-
53 * reducing permutation L->Perm is not used. x and b are in the permuted
54 * ordering, not your original ordering. x and b are n-by-1; this routine
55 * does not handle multiple right-hand-sides.
56 *
57 * workspace: Flag (nrow), Head (nrow+1), W (maxrank*nrow), Iwork (nrow),
58 * where maxrank is 2, 4, or 8.
59 *
60 * Only real matrices are supported. A symbolic L is converted into a
61 * numeric identity matrix.
62 */
63
64 #ifndef NGPL
65 #ifndef NMODIFY
66
67 #include "cholmod_internal.h"
68 #include "cholmod_modify.h"
69
70
71 /* ========================================================================== */
72 /* === cholmod_updown ======================================================= */
73 /* ========================================================================== */
74
75 /* Compute the new LDL' factorization of LDL'+CC' (an update) or LDL'-CC'
76 * (a downdate). The factor object L need not be an LDL' factorization; it
77 * is converted to one if it isn't. */
78
CHOLMOD(updown)79 int CHOLMOD(updown)
80 (
81 /* ---- input ---- */
82 int update, /* TRUE for update, FALSE for downdate */
83 cholmod_sparse *C, /* the incoming sparse update */
84 /* ---- in/out --- */
85 cholmod_factor *L, /* factor to modify */
86 /* --------------- */
87 cholmod_common *Common
88 )
89 {
90 return (CHOLMOD(updown_mask2) (update, C, NULL, NULL, 0, L, NULL, NULL,
91 Common)) ;
92 }
93
94
95 /* ========================================================================== */
96 /* === cholmod_updown_solve ================================================= */
97 /* ========================================================================== */
98
99 /* Does the same as cholmod_updown, except that it also updates/downdates the
100 * solution to Lx=b+DeltaB. x and b must be n-by-1 dense matrices. b is not
101 * need as input to this routine, but a sparse change to b is (DeltaB). Only
102 * entries in DeltaB corresponding to columns modified in L are accessed; the
103 * rest are ignored.
104 */
105
CHOLMOD(updown_solve)106 int CHOLMOD(updown_solve)
107 (
108 /* ---- input ---- */
109 int update, /* TRUE for update, FALSE for downdate */
110 cholmod_sparse *C, /* the incoming sparse update */
111 /* ---- in/out --- */
112 cholmod_factor *L, /* factor to modify */
113 cholmod_dense *X, /* solution to Lx=b (size n-by-1) */
114 cholmod_dense *DeltaB, /* change in b, zero on output */
115 /* --------------- */
116 cholmod_common *Common
117 )
118 {
119 return (CHOLMOD(updown_mask2) (update, C, NULL, NULL, 0, L, X, DeltaB,
120 Common)) ;
121 }
122
123
124 /* ========================================================================== */
125 /* === Power2 =============================================================== */
126 /* ========================================================================== */
127
128 /* Power2 [i] is smallest power of 2 that is >= i (for i in range 0 to 8) */
129
130 static Int Power2 [ ] =
131 {
132 /* 0 1 2 3 4 5 6 7 8 */
133 0, 1, 2, 4, 4, 8, 8, 8, 8
134 } ;
135
136 /* ========================================================================== */
137 /* === debug routines ======================================================= */
138 /* ========================================================================== */
139
140 #ifndef NDEBUG
141
dump_set(Int s,Int ** Set_ps1,Int ** Set_ps2,Int j,Int n,cholmod_common * Common)142 static void dump_set (Int s, Int **Set_ps1, Int **Set_ps2, Int j, Int n,
143 cholmod_common *Common)
144 {
145 Int *p, len, i, ilast ;
146
147 if (CHOLMOD(dump) < -1)
148 {
149 /* no checks if debug level is -2 or less */
150 return ;
151 }
152
153 len = Set_ps2 [s] - Set_ps1 [s] ;
154 PRINT2 (("Set s: "ID" len: "ID":", s, len)) ;
155 ASSERT (len > 0) ;
156 ilast = j ;
157 for (p = Set_ps1 [s] ; p < Set_ps2 [s] ; p++)
158 {
159 i = *p ;
160 PRINT3 ((" "ID"", i)) ;
161 ASSERT (i > ilast && i < n) ;
162 ilast = i ;
163 }
164 PRINT3 (("\n")) ;
165 }
166
dump_col(char * w,Int j,Int p1,Int p2,Int * Li,double * Lx,Int n,cholmod_common * Common)167 static void dump_col
168 (
169 char *w, Int j, Int p1, Int p2, Int *Li, double *Lx, Int n,
170 cholmod_common *Common
171 )
172 {
173 Int p, row, lastrow ;
174
175 if (CHOLMOD(dump) < -1)
176 {
177 /* no checks if debug level is -2 or less */
178 return ;
179 }
180
181 PRINT3 (("\n\nDUMP COL==== j = "ID" %s: p1="ID" p2="ID" \n", j, w, p1,p2));
182 lastrow = -1 ;
183 for (p = p1 ; p < p2 ; p++)
184 {
185 PRINT3 ((" "ID": ", p)) ;
186 row = Li [p] ;
187 PRINT3 ((""ID" ", Li [p])) ;
188 PRINT3 (("%g ", Lx [p])) ;
189 PRINT3 (("\n")) ;
190 ASSERT (row > lastrow && row < n) ;
191 lastrow = row ;
192 }
193 ASSERT (p1 < p2) ;
194 ASSERT (Li [p1] == j) ;
195 PRINT3 (("\n")) ;
196 }
197 #endif
198
199
200 /* ========================================================================== */
201 /* === a path =============================================================== */
202 /* ========================================================================== */
203
204 /* A path is a set of nodes of the etree which are all affected by the same
205 * columns of C. */
206
207 typedef struct Path_struct
208 {
209 Int start ; /* column at which to start, or EMPTY if initial */
210 Int end ; /* column at which to end, or EMPTY if initial */
211 Int ccol ; /* column of C to which path refers */
212 Int parent ; /* parent path */
213 Int c ; /* child of j along this path */
214 Int next ; /* next path in link list */
215 Int rank ; /* number of rank-1 paths merged onto this path */
216 Int order ; /* dfs order of this path */
217 Int wfirst ; /* first column of W to affect this path */
218 Int pending ; /* column at which the path is pending */
219 Int botrow ; /* for partial update/downdate of solution to Lx=b */
220
221 } Path_type ;
222
223
224 /* ========================================================================== */
225 /* === dfs ================================================================== */
226 /* ========================================================================== */
227
228 /* Compute the DFS order of the set of paths. This can be recursive because
229 * there are at most 23 paths to sort: one for each column of C (8 at most),
230 * and one for each node in a balanced binary tree with 8 leaves (15).
231 * Stack overflow is thus not a problem. */
232
dfs(Path_type * Path,Int k,Int path,Int * path_order,Int * w_order,Int depth,Int npaths)233 static void dfs
234 (
235 Path_type *Path, /* the set of Paths */
236 Int k, /* the rank of the update/downdate */
237 Int path, /* which path to work on */
238 Int *path_order, /* the current path order */
239 Int *w_order, /* the current order of the columns of W */
240 Int depth,
241 Int npaths /* total number of paths */
242 )
243 {
244 Int c ; /* child path */
245
246 ASSERT (path >= 0 && path < npaths) ;
247 if (path < k)
248 {
249 /* this is a leaf node, corresponding to column W (:,path) */
250 /* and column C (:, Path [path].ccol) */
251 ASSERT (Path [path].ccol >= 0) ;
252 Path [path].wfirst = *w_order ;
253 Path [path].order = *w_order ;
254 (*w_order)++ ;
255 }
256 else
257 {
258 /* this is a non-leaf path, within the tree */
259 ASSERT (Path [path].c != EMPTY) ;
260 ASSERT (Path [path].ccol == EMPTY) ;
261 /* order each child path */
262 for (c = Path [path].c ; c != EMPTY ; c = Path [c].next)
263 {
264 dfs (Path, k, c, path_order, w_order, depth+1, npaths) ;
265 if (Path [path].wfirst == EMPTY)
266 {
267 Path [path].wfirst = Path [c].wfirst ;
268 }
269 }
270 /* order this path next */
271 Path [path].order = (*path_order)++ ;
272 }
273 }
274
275
276 /* ========================================================================== */
277 /* === numeric update/downdate routines ===================================== */
278 /* ========================================================================== */
279
280 #define WDIM 1
281 #include "t_cholmod_updown.c"
282 #define WDIM 2
283 #include "t_cholmod_updown.c"
284 #define WDIM 4
285 #include "t_cholmod_updown.c"
286 #define WDIM 8
287 #include "t_cholmod_updown.c"
288
289
290 /* ========================================================================== */
291 /* === cholmod_updown_mark ================================================== */
292 /* ========================================================================== */
293
294 /* Update/downdate LDL' +/- C*C', and update/downdate selected portions of the
295 * solution to Lx=b.
296 *
297 * The original system is L*x = b. The new system is Lnew*xnew = b + deltab.
298 * deltab(i) can be nonzero only if column i of L is modified by the update/
299 * downdate. If column i is not modified, the deltab(i) is not accessed.
300 *
301 * The solution to Lx=b is not modified if either X or DeltaB are NULL.
302 *
303 * Rowmark and colmark:
304 * --------------------
305 *
306 * rowmark and colmark affect which portions of L take part in the update/
307 * downdate of the solution to Lx=b. They do not affect how L itself is
308 * updated/downdated. They are both ignored if X or DeltaB are NULL.
309 *
310 * If not NULL, rowmark is an integer array of size n where L is n-by-n.
311 * rowmark [j] defines the part of column j of L that takes part in the update/
312 * downdate of the forward solve, Lx=b. Specifically, if i = rowmark [j],
313 * then L(j:i-1,j) is used, and L(i:end,j) is ignored.
314 *
315 * If not NULL, colmark is an integer array of size C->ncol. colmark [ccol]
316 * for a column C(:,ccol) redefines those parts of L that take part in the
317 * update/downdate of Lx=b. Each column of C affects a set of columns of L.
318 * If column ccol of C affects column j of L, then the new rowmark [j] of
319 * column j of L is defined as colmark [ccol]. In a multiple-rank update/
320 * downdate, if two or more columns of C affect column j, its new rowmark [j]
321 * is the colmark of the least-numbered column of C. colmark is ignored if
322 * it is NULL, in which case rowmark is not modified. If colmark [ccol] is
323 * EMPTY (-1), then rowmark is not modified for that particular column of C.
324 * colmark is ignored if it is NULL, or rowmark, X, or DeltaB are NULL.
325 *
326 * The algorithm for modifying the solution to Lx=b when rowmark and colmark
327 * are NULL is as follows:
328 *
329 * for each column j of L that is modified:
330 * deltab (j:end) += L (j:end,j) * x(j)
331 * modify L
332 * for each column j of L that is modified:
333 * x (j) = deltab (j)
334 * deltab (j) = 0
335 * deltab (j+1:end) -= L (j+1:end,j) * x(j)
336 *
337 * If rowmark is non-NULL but colmark is NULL:
338 *
339 * for each column j of L that is modified:
340 * deltab (j:rowmark(j)-1) += L (j:rowmark(j)-1,j) * x(j)
341 * modify L
342 * for each column j of L that is modified:
343 * x (j) = deltab (j)
344 * deltab (j) = 0
345 * deltab (j+1:rowmark(j)-1) -= L (j+1:rowmark(j)-1,j) * x(j)
346 *
347 * If both rowmark and colmark are non-NULL:
348 *
349 * for each column j of L that is modified:
350 * deltab (j:rowmark(j)-1) += L (j:rowmark(j)-1,j) * x(j)
351 * modify L
352 * for each column j of L that is modified:
353 * modify rowmark (j) according to colmark
354 * for each column j of L that is modified:
355 * x (j) = deltab (j)
356 * deltab (j) = 0
357 * deltab (j+1:rowmark(j)-1) -= L (j+1:rowmark(j)-1,j) * x(j)
358 *
359 * Note that if the rank of C exceeds k = Common->maxrank (which is 2, 4, or 8),
360 * then the update/downdate is done as a series of rank-k updates. In this
361 * case, the above algorithm is repeated for each block of k columns of C.
362 *
363 * Unless it leads to no changes in rowmark, colmark should be used only if
364 * C->ncol <= Common->maxrank, because the update/downdate is done with maxrank
365 * columns at a time. Otherwise, the results are undefined.
366 *
367 * This routine is an "expert" routine. It is meant for use in LPDASA only.
368 */
369
CHOLMOD(updown_mark)370 int CHOLMOD(updown_mark)
371 (
372 /* ---- input ---- */
373 int update, /* TRUE for update, FALSE for downdate */
374 cholmod_sparse *C, /* the incoming sparse update */
375 Int *colmark, /* Int array of size n. */
376 /* ---- in/out --- */
377 cholmod_factor *L, /* factor to modify */
378 cholmod_dense *X, /* solution to Lx=b (size n-by-1) */
379 cholmod_dense *DeltaB, /* change in b, zero on output */
380 /* --------------- */
381 cholmod_common *Common
382 )
383 {
384 return (CHOLMOD(updown_mask2) (update, C, colmark, NULL, 0, L, X, DeltaB,
385 Common)) ;
386 }
387
388
389 /* ========================================================================== */
390 /* === cholmod_updown_mask ================================================== */
391 /* ========================================================================== */
392
CHOLMOD(updown_mask)393 int CHOLMOD(updown_mask)
394 (
395 /* ---- input ---- */
396 int update, /* TRUE for update, FALSE for downdate */
397 cholmod_sparse *C, /* the incoming sparse update */
398 Int *colmark, /* Int array of size n. See cholmod_updown.c */
399 Int *mask, /* size n */
400 /* ---- in/out --- */
401 cholmod_factor *L, /* factor to modify */
402 cholmod_dense *X, /* solution to Lx=b (size n-by-1) */
403 cholmod_dense *DeltaB, /* change in b, zero on output */
404 /* --------------- */
405 cholmod_common *Common
406 )
407 {
408 Int maskmark = 0 ;
409 return (CHOLMOD(updown_mask2) (update, C, colmark, mask, maskmark,
410 L, X, DeltaB, Common)) ;
411 }
412
413 /* ========================================================================== */
414 /* === cholmod_updown_mask2 ================================================= */
415 /* ========================================================================== */
416
CHOLMOD(updown_mask2)417 int CHOLMOD(updown_mask2)
418 (
419 /* ---- input ---- */
420 int update, /* TRUE for update, FALSE for downdate */
421 cholmod_sparse *C, /* the incoming sparse update */
422 Int *colmark, /* Int array of size n. See cholmod_updown.c */
423 Int *mask, /* size n */
424 Int maskmark,
425 /* ---- in/out --- */
426 cholmod_factor *L, /* factor to modify */
427 cholmod_dense *X, /* solution to Lx=b (size n-by-1) */
428 cholmod_dense *DeltaB, /* change in b, zero on output */
429 /* --------------- */
430 cholmod_common *Common
431 )
432 {
433 double xj, fl ;
434 double *Lx, *W, *Xx, *Nx ;
435 Int *Li, *Lp, *Lnz, *Cp, *Ci, *Cnz, *Head, *Flag, *Stack, *Lnext, *Iwork,
436 *Set_ps1 [32], *Set_ps2 [32], *ps1, *ps2 ;
437 size_t maxrank ;
438 Path_type OrderedPath [32], Path [32] ;
439 Int n, wdim, k1, k2, npaths, i, j, row, packed, ccol, p, cncol, do_solve,
440 mark, jj, j2, kk, nextj, p1, p2, c, use_colmark, newlnz,
441 k, newpath, path_order, w_order, scattered, path, newparent, pp1, pp2,
442 smax, maxrow, row1, nsets, s, p3, newlnz1, Set [32], top, len, lnz, m,
443 botrow ;
444 size_t w ;
445 int ok = TRUE ;
446 DEBUG (Int oldparent) ;
447
448 /* ---------------------------------------------------------------------- */
449 /* check inputs */
450 /* ---------------------------------------------------------------------- */
451
452 RETURN_IF_NULL_COMMON (FALSE) ;
453 RETURN_IF_NULL (C, FALSE) ;
454 RETURN_IF_NULL (L, FALSE) ;
455 RETURN_IF_XTYPE_INVALID (L, CHOLMOD_PATTERN, CHOLMOD_REAL, FALSE) ;
456 RETURN_IF_XTYPE_INVALID (C, CHOLMOD_REAL, CHOLMOD_REAL, FALSE) ;
457 n = L->n ;
458 cncol = C->ncol ;
459 if (!(C->sorted))
460 {
461 ERROR (CHOLMOD_INVALID, "C must have sorted columns") ;
462 return (FALSE) ;
463 }
464 if (n != (Int) (C->nrow))
465 {
466 ERROR (CHOLMOD_INVALID, "C and L dimensions do not match") ;
467 return (FALSE) ;
468 }
469 do_solve = (X != NULL) && (DeltaB != NULL) ;
470 if (do_solve)
471 {
472 RETURN_IF_XTYPE_INVALID (X, CHOLMOD_REAL, CHOLMOD_REAL, FALSE) ;
473 RETURN_IF_XTYPE_INVALID (DeltaB, CHOLMOD_REAL, CHOLMOD_REAL, FALSE) ;
474 Xx = X->x ;
475 Nx = DeltaB->x ;
476 if (X->nrow != L->n || X->ncol != 1 || DeltaB->nrow != L->n ||
477 DeltaB->ncol != 1 || Xx == NULL || Nx == NULL)
478 {
479 ERROR (CHOLMOD_INVALID, "X and/or DeltaB invalid") ;
480 return (FALSE) ;
481 }
482 }
483 else
484 {
485 Xx = NULL ;
486 Nx = NULL ;
487 }
488 Common->status = CHOLMOD_OK ;
489 Common->modfl = 0 ;
490
491 fl = 0 ;
492 use_colmark = (colmark != NULL) ;
493
494 /* ---------------------------------------------------------------------- */
495 /* allocate workspace */
496 /* ---------------------------------------------------------------------- */
497
498 /* Note: cholmod_rowadd and cholmod_rowdel use the second n doubles in
499 * Common->Xwork for Cx, and then perform a rank-1 update here, which uses
500 * the first n doubles in Common->Xwork. Both the rowadd and rowdel
501 * routines allocate enough workspace so that Common->Xwork isn't destroyed
502 * below. Also, both cholmod_rowadd and cholmod_rowdel use the second n
503 * ints in Common->Iwork for Ci.
504 */
505
506 /* make sure maxrank is in the proper range */
507 maxrank = CHOLMOD(maxrank) (n, Common) ;
508 k = MIN (cncol, (Int) maxrank) ; /* maximum k is wdim */
509 wdim = Power2 [k] ; /* number of columns needed in W */
510 ASSERT (wdim <= (Int) maxrank) ;
511 PRINT1 (("updown wdim final "ID" k "ID"\n", wdim, k)) ;
512
513 /* w = wdim * n */
514 w = CHOLMOD(mult_size_t) (n, wdim, &ok) ;
515 if (!ok)
516 {
517 ERROR (CHOLMOD_TOO_LARGE, "problem too large") ;
518 return (FALSE) ;
519 }
520
521 CHOLMOD(allocate_work) (n, n, w, Common) ;
522 if (Common->status < CHOLMOD_OK || maxrank == 0)
523 {
524 /* out of memory, L is returned unchanged */
525 return (FALSE) ;
526 }
527
528 /* ---------------------------------------------------------------------- */
529 /* convert to simplicial numeric LDL' factor, if not already */
530 /* ---------------------------------------------------------------------- */
531
532 if (L->xtype == CHOLMOD_PATTERN || L->is_super || L->is_ll)
533 {
534 /* can only update/downdate a simplicial LDL' factorization */
535 CHOLMOD(change_factor) (CHOLMOD_REAL, FALSE, FALSE, FALSE, FALSE, L,
536 Common) ;
537 if (Common->status < CHOLMOD_OK)
538 {
539 /* out of memory, L is returned unchanged */
540 return (FALSE) ;
541 }
542 }
543
544 /* ---------------------------------------------------------------------- */
545 /* get inputs */
546 /* ---------------------------------------------------------------------- */
547
548 /* mark = CHOLMOD(clear_flag) (Common) ; */
549 CHOLMOD_CLEAR_FLAG (Common) ;
550 mark = Common->mark ;
551
552 PRINT1 (("updown, rank %g update %d\n", (double) C->ncol, update)) ;
553 DEBUG (CHOLMOD(dump_factor) (L, "input L for updown", Common)) ;
554 ASSERT (CHOLMOD(dump_sparse) (C, "input C for updown", Common) >= 0) ;
555
556 Ci = C->i ;
557 Cp = C->p ;
558 Cnz = C->nz ;
559 packed = C->packed ;
560 ASSERT (IMPLIES (!packed, Cnz != NULL)) ;
561
562 /* ---------------------------------------------------------------------- */
563 /* quick return */
564 /* ---------------------------------------------------------------------- */
565
566 if (cncol <= 0 || n == 0)
567 {
568 /* nothing to do */
569 return (TRUE) ;
570 }
571
572 /* ---------------------------------------------------------------------- */
573 /* get L */
574 /* ---------------------------------------------------------------------- */
575
576 Li = L->i ;
577 Lx = L->x ;
578 Lp = L->p ;
579 Lnz = L->nz ;
580 Lnext = L->next ;
581 ASSERT (Lnz != NULL) ;
582
583 /* ---------------------------------------------------------------------- */
584 /* get workspace */
585 /* ---------------------------------------------------------------------- */
586
587 Flag = Common->Flag ; /* size n, Flag [i] <= mark must hold */
588 Head = Common->Head ; /* size n, Head [i] == EMPTY must hold */
589 W = Common->Xwork ; /* size n-by-wdim, zero on input and output*/
590
591 /* note that Iwork [n .. 2*n-1] (i/i/l) may be in use in rowadd/rowdel: */
592 Iwork = Common->Iwork ;
593 Stack = Iwork ; /* size n, uninitialized (i/i/l) */
594
595 /* ---------------------------------------------------------------------- */
596 /* entire rank-cncol update, done as a sequence of rank-k updates */
597 /* ---------------------------------------------------------------------- */
598
599 ps1 = NULL ;
600 ps2 = NULL ;
601
602 for (k1 = 0 ; k1 < cncol ; k1 += k)
603 {
604
605 /* ------------------------------------------------------------------ */
606 /* get the next k columns of C for the update/downdate */
607 /* ------------------------------------------------------------------ */
608
609 /* the last update/downdate might be less than rank-k */
610 if (k > cncol - k1)
611 {
612 k = cncol - k1 ;
613 wdim = Power2 [k] ;
614 }
615 k2 = k1 + k - 1 ;
616
617 /* workspaces are in the following state, on input and output */
618 ASSERT (CHOLMOD(dump_work) (TRUE, TRUE, wdim, Common)) ;
619
620 /* ------------------------------------------------------------------ */
621 /* create a zero-length path for each column of W */
622 /* ------------------------------------------------------------------ */
623
624 nextj = n ;
625 path = 0 ;
626 for (ccol = k1 ; ccol <= k2 ; ccol++)
627 {
628 PRINT1 (("Column ["ID"]: "ID"\n", path, ccol)) ;
629 ASSERT (ccol >= 0 && ccol <= cncol) ;
630 pp1 = Cp [ccol] ;
631 pp2 = (packed) ? (Cp [ccol+1]) : (pp1 + Cnz [ccol]) ;
632 /* get the row index j of the first entry in C (:,ccol) */
633 if (pp2 > pp1)
634 {
635 /* Column ccol of C has at least one entry. */
636 j = Ci [pp1] ;
637 }
638 else
639 {
640 /* Column ccol of C is empty. Pretend it has one entry in
641 * the last column with numerical value of zero. */
642 j = n-1 ;
643 }
644 ASSERT (j >= 0 && j < n) ;
645
646 /* find first column to work on */
647 nextj = MIN (nextj, j) ;
648
649 Path [path].ccol = ccol ; /* which column of C this path is for */
650 Path [path].start = EMPTY ; /* paths for C have zero length */
651 Path [path].end = EMPTY ;
652 Path [path].parent = EMPTY ; /* no parent yet */
653 Path [path].rank = 1 ; /* one column of W */
654 Path [path].c = EMPTY ; /* no child of this path (case A) */
655 Path [path].next = Head [j] ; /* this path is pending at col j */
656 Path [path].pending = j ; /* this path is pending at col j */
657 Head [j] = path ; /* this path is pending at col j */
658 PRINT1(("Path "ID" starts: start "ID" end "ID" parent "ID" c "ID""
659 "j "ID" ccol "ID"\n", path, Path [path].start,
660 Path [path].end, Path [path].parent,
661 Path [path].c, j, ccol)) ;
662
663 /* initialize botrow for this path */
664 Path [path].botrow = (use_colmark) ? colmark [ccol] : n ;
665
666 path++ ;
667 }
668
669 /* we start with paths 0 to k-1. Next one (now unused) is npaths */
670 npaths = k ;
671
672 j = nextj ;
673 ASSERT (j < n) ;
674 scattered = FALSE ;
675
676 /* ------------------------------------------------------------------ */
677 /* symbolic update of columns of L */
678 /* ------------------------------------------------------------------ */
679
680 while (j < n)
681 {
682 ASSERT (j >= 0 && j < n && Lnz [j] > 0) ;
683
684 /* the old column, Li [p1..p2-1]. D (j,j) is stored in Lx [p1] */
685 p1 = Lp [j] ;
686 newlnz = Lnz [j] ;
687 p2 = p1 + newlnz ;
688
689 #ifndef NDEBUG
690 PRINT1 (("\n=========Column j="ID" p1 "ID" p2 "ID" lnz "ID" \n",
691 j, p1, p2, newlnz)) ;
692 dump_col ("Old", j, p1, p2, Li, Lx, n, Common) ;
693 oldparent = (Lnz [j] > 1) ? (Li [p1 + 1]) : EMPTY ;
694 ASSERT (CHOLMOD(dump_work) (TRUE, FALSE, 0, Common)) ;
695 ASSERT (!scattered) ;
696 PRINT1 (("Col "ID": Checking paths, npaths: "ID"\n", j, npaths)) ;
697 for (kk = 0 ; kk < npaths ; kk++)
698 {
699 Int kk2, found, j3 = Path [kk].pending ;
700 PRINT2 (("Path "ID" pending at "ID".\n", kk, j3)) ;
701 if (j3 != EMPTY)
702 {
703 /* Path kk must be somewhere in link list for column j3 */
704 ASSERT (Head [j3] != EMPTY) ;
705 PRINT3 ((" List at "ID": ", j3)) ;
706 found = FALSE ;
707 for (kk2 = Head [j3] ; kk2 != EMPTY ; kk2 = Path [kk2].next)
708 {
709 PRINT3 ((""ID" ", kk2)) ;
710 ASSERT (Path [kk2].pending == j3) ;
711 found = found || (kk2 == kk) ;
712 }
713 PRINT3 (("\n")) ;
714 ASSERT (found) ;
715 }
716 }
717 PRINT1 (("\nCol "ID": Paths at this column, head "ID"\n",
718 j, Head [j]));
719 ASSERT (Head [j] != EMPTY) ;
720 for (kk = Head [j] ; kk != EMPTY ; kk = Path [kk].next)
721 {
722 PRINT1 (("path "ID": (c="ID" j="ID") npaths "ID"\n",
723 kk, Path[kk].c, j, npaths)) ;
724 ASSERT (kk >= 0 && kk < npaths) ;
725 ASSERT (Path [kk].pending == j) ;
726 }
727 #endif
728
729 /* -------------------------------------------------------------- */
730 /* determine the path we're on */
731 /* -------------------------------------------------------------- */
732
733 /* get the first old path at column j */
734 path = Head [j] ;
735
736 /* -------------------------------------------------------------- */
737 /* update/downdate of forward solve, Lx=b */
738 /* -------------------------------------------------------------- */
739
740 if (do_solve)
741 {
742 xj = Xx [j] ;
743 if (IS_NONZERO (xj))
744 {
745 xj = Xx [j] ;
746 /* This is first time column j has been seen for entire */
747 /* rank-k update/downdate. */
748
749 /* DeltaB += Lold (j:botrow-1,j) * X (j) */
750 Nx [j] += xj ; /* diagonal of L */
751
752 /* find the botrow for this column */
753 botrow = (use_colmark) ? Path [path].botrow : n ;
754
755 for (p = p1 + 1 ; p < p2 ; p++)
756 {
757 i = Li [p] ;
758 if (i >= botrow)
759 {
760 break ;
761 }
762 Nx [i] += Lx [p] * xj ;
763 }
764
765 /* clear X[j] to flag col j of Lold as having been seen. If
766 * X (j) was initially zero, then the above code is never
767 * executed for column j. This is safe, since if xj=0 the
768 * code above does not do anything anyway. */
769 Xx [j] = 0.0 ;
770 }
771 }
772
773 /* -------------------------------------------------------------- */
774 /* start a new path at this column if two or more paths merge */
775 /* -------------------------------------------------------------- */
776
777 newpath =
778 /* start a new path if paths have merged */
779 (Path [path].next != EMPTY)
780 /* or if j is the first node on a path (case A). */
781 || (Path [path].c == EMPTY) ;
782
783 if (newpath)
784 {
785 /* get the botrow of the first path at column j */
786 botrow = (use_colmark) ? Path [path].botrow : n ;
787
788 path = npaths++ ;
789 ASSERT (npaths <= 3*k) ;
790 Path [path].ccol = EMPTY ; /* no single col of C for this path*/
791 Path [path].start = j ; /* path starts at this column j */
792 Path [path].end = EMPTY ; /* don't know yet where it ends */
793 Path [path].parent = EMPTY ;/* don't know parent path yet */
794 Path [path].rank = 0 ; /* rank is sum of child path ranks */
795 PRINT1 (("Path "ID" starts: start "ID" end "ID" parent "ID"\n",
796 path, Path [path].start, Path [path].end, Path [path].parent)) ;
797
798 /* set the botrow of the new path */
799 Path [path].botrow = (use_colmark) ? botrow : n ;
800 }
801
802 /* -------------------------------------------------------------- */
803 /* for each path kk pending at column j */
804 /* -------------------------------------------------------------- */
805
806 /* make a list of the sets that need to be merged into column j */
807 nsets = 0 ;
808
809 for (kk = Head [j] ; kk != EMPTY ; kk = Path [kk].next)
810 {
811
812 /* ---------------------------------------------------------- */
813 /* path kk is at (c,j) */
814 /* ---------------------------------------------------------- */
815
816 c = Path [kk].c ;
817 ASSERT (c < j) ;
818 PRINT1 (("TUPLE on path "ID" (c="ID" j="ID")\n", kk, c, j)) ;
819 ASSERT (Path [kk].pending == j) ;
820
821 if (newpath)
822 {
823 /* finalize path kk and find rank of this path */
824 Path [kk].end = c ; /* end of old path is previous node c */
825 Path [kk].parent = path ; /* parent is this path */
826 Path [path].rank += Path [kk].rank ; /* sum up ranks */
827 Path [kk].pending = EMPTY ;
828 PRINT1 (("Path "ID" done:start "ID" end "ID" parent "ID"\n",
829 kk, Path [kk].start, Path [kk].end, Path [kk].parent)) ;
830 }
831
832 if (c == EMPTY)
833 {
834
835 /* ------------------------------------------------------ */
836 /* CASE A: first node in path */
837 /* ------------------------------------------------------ */
838
839 /* update: add pattern of incoming column */
840
841 /* Column ccol of C is in Ci [pp1 ... pp2-1] */
842 ccol = Path [kk].ccol ;
843 pp1 = Cp [ccol] ;
844 pp2 = (packed) ? (Cp [ccol+1]) : (pp1 + Cnz [ccol]) ;
845 PRINT1 (("Case A, ccol = "ID" len "ID"\n", ccol, pp2-pp1)) ;
846 ASSERT (IMPLIES (pp2 > pp1, Ci [pp1] == j)) ;
847
848 if (!scattered)
849 {
850 /* scatter the original pattern of column j of L */
851 for (p = p1 ; p < p2 ; p++)
852 {
853 Flag [Li [p]] = mark ;
854 }
855 scattered = TRUE ;
856 }
857
858 /* scatter column ccol of C (skip first entry, j) */
859 newlnz1 = newlnz ;
860 for (p = pp1 + 1 ; p < pp2 ; p++)
861 {
862 row = Ci [p] ;
863 if (Flag [row] < mark)
864 {
865 /* this is a new entry in Lj' */
866 Flag [row] = mark ;
867 newlnz++ ;
868 }
869 }
870 if (newlnz1 != newlnz)
871 {
872 /* column ccol of C adds something to column j of L */
873 Set [nsets++] = FLIP (ccol) ;
874 }
875
876 }
877 else if (Head [c] == 1)
878 {
879
880 /* ------------------------------------------------------ */
881 /* CASE B: c is old, but changed, child of j */
882 /* CASE C: new child of j */
883 /* ------------------------------------------------------ */
884
885 /* Head [c] is 1 if col c of L has new entries,
886 * EMPTY otherwise */
887 Flag [c] = 0 ;
888 Head [c] = EMPTY ;
889
890 /* update: add Lc' */
891
892 /* column c of L is in Li [pp1 .. pp2-1] */
893 pp1 = Lp [c] ;
894 pp2 = pp1 + Lnz [c] ;
895 PRINT1 (("Case B/C: c = "ID"\n", c)) ;
896 DEBUG (dump_col ("Child", c, pp1, pp2, Li, Lx, n, Common)) ;
897 ASSERT (j == Li [pp1 + 1]) ; /* j is new parent of c */
898
899 if (!scattered)
900 {
901 /* scatter the original pattern of column j of L */
902 for (p = p1 ; p < p2 ; p++)
903 {
904 Flag [Li [p]] = mark ;
905 }
906 scattered = TRUE ;
907 }
908
909 /* scatter column c of L (skip first two entries, c and j)*/
910 newlnz1 = newlnz ;
911 for (p = pp1 + 2 ; p < pp2 ; p++)
912 {
913 row = Li [p] ;
914 if (Flag [row] < mark)
915 {
916 /* this is a new entry in Lj' */
917 Flag [row] = mark ;
918 newlnz++ ;
919 }
920 }
921 PRINT2 (("\n")) ;
922
923 if (newlnz1 != newlnz)
924 {
925 /* column c of L adds something to column j of L */
926 Set [nsets++] = c ;
927 }
928 }
929 }
930
931 /* -------------------------------------------------------------- */
932 /* update the pattern of column j of L */
933 /* -------------------------------------------------------------- */
934
935 /* Column j of L will be in Li/Lx [p1 .. p3-1] */
936 p3 = p1 + newlnz ;
937 ASSERT (IMPLIES (nsets == 0, newlnz == Lnz [j])) ;
938 PRINT1 (("p1 "ID" p2 "ID" p3 "ID" nsets "ID"\n", p1, p2, p3,nsets));
939
940 /* -------------------------------------------------------------- */
941 /* ensure we have enough space for the longer column */
942 /* -------------------------------------------------------------- */
943
944 if (nsets > 0 && p3 > Lp [Lnext [j]])
945 {
946 PRINT1 (("Col realloc: j "ID" newlnz "ID"\n", j, newlnz)) ;
947 if (!CHOLMOD(reallocate_column) (j, newlnz, L, Common))
948 {
949 /* out of memory, L is now simplicial symbolic */
950 CHOLMOD(clear_flag) (Common) ;
951 for (j = 0 ; j <= n ; j++)
952 {
953 Head [j] = EMPTY ;
954 }
955 ASSERT (CHOLMOD(dump_work) (TRUE, TRUE, wdim, Common)) ;
956 return (FALSE) ;
957 }
958 /* L->i and L->x may have moved. Column j has moved too */
959 Li = L->i ;
960 Lx = L->x ;
961 p1 = Lp [j] ;
962 p2 = p1 + Lnz [j] ;
963 p3 = p1 + newlnz ;
964 }
965
966 /* -------------------------------------------------------------- */
967 /* create set pointers */
968 /* -------------------------------------------------------------- */
969
970 for (s = 0 ; s < nsets ; s++)
971 {
972 /* Pattern of Set s is *(Set_ps1 [s] ... Set_ps2 [s]-1) */
973 c = Set [s] ;
974 if (c < EMPTY)
975 {
976 /* column ccol of C, skip first entry (j) */
977 ccol = FLIP (c) ;
978 pp1 = Cp [ccol] ;
979 pp2 = (packed) ? (Cp [ccol+1]) : (pp1 + Cnz [ccol]) ;
980 ASSERT (pp2 - pp1 > 1) ;
981 Set_ps1 [s] = &(Ci [pp1 + 1]) ;
982 Set_ps2 [s] = &(Ci [pp2]) ;
983 PRINT1 (("set "ID" is ccol "ID"\n", s, ccol)) ;
984 }
985 else
986 {
987 /* column c of L, skip first two entries (c and j) */
988 pp1 = Lp [c] ;
989 pp2 = pp1 + Lnz [c] ;
990 ASSERT (Lnz [c] > 2) ;
991 Set_ps1 [s] = &(Li [pp1 + 2]) ;
992 Set_ps2 [s] = &(Li [pp2]) ;
993 PRINT1 (("set "ID" is L "ID"\n", s, c)) ;
994 }
995 DEBUG (dump_set (s, Set_ps1, Set_ps2, j, n, Common)) ;
996 }
997
998 /* -------------------------------------------------------------- */
999 /* multiset merge */
1000 /* -------------------------------------------------------------- */
1001
1002 /* Merge the sets into a single sorted set, Lj'. Before the merge
1003 * starts, column j is located in Li/Lx [p1 ... p2-1] and the
1004 * space Li/Lx [p2 ... p3-1] is empty. p1 is Lp [j], p2 is
1005 * Lp [j] + Lnz [j] (the old length of the column), and p3 is
1006 * Lp [j] + newlnz (the new and longer length of the column).
1007 *
1008 * The sets 0 to nsets-1 are defined by the Set_ps1 and Set_ps2
1009 * pointers. Set s is located in *(Set_ps1 [s] ... Set_ps2 [s]-1).
1010 * It may be a column of C, or a column of L. All row indices i in
1011 * the sets are in the range i > j and i < n. All sets are sorted.
1012 *
1013 * The merge into column j of L is done in place.
1014 *
1015 * During the merge, p2 and p3 are updated. Li/Lx [p1..p2-1]
1016 * reflects the indices of the old column j of L that are yet to
1017 * be merged into the new column. Entries in their proper place in
1018 * the new column j of L are located in Li/Lx [p3 ... p1+newlnz-1].
1019 * The merge finishes when p2 == p3.
1020 *
1021 * During the merge, set s consumed as it is merged into column j of
1022 * L. Its unconsumed contents are *(Set_ps1 [s] ... Set_ps2 [s]-1).
1023 * When a set is completely consumed, it is removed from the set of
1024 * sets, and nsets is decremented.
1025 *
1026 * The multiset merge and 2-set merge finishes when p2 == p3.
1027 */
1028
1029 PRINT1 (("Multiset merge p3 "ID" p2 "ID" nsets "ID"\n",
1030 p3, p2, nsets)) ;
1031
1032 while (p3 > p2 && nsets > 1)
1033 {
1034
1035 #ifndef NDEBUG
1036 PRINT2 (("\nMultiset merge. nsets = "ID"\n", nsets)) ;
1037 PRINT2 (("Source col p1 = "ID", p2 = "ID", p3= "ID"\n",
1038 p1, p2, p3)) ;
1039 for (p = p1 + 1 ; p < p2 ; p++)
1040 {
1041 PRINT2 ((" p: "ID" source row "ID" %g\n",
1042 p, Li[p], Lx[p])) ;
1043 ASSERT (Li [p] > j && Li [p] < n) ;
1044 }
1045 PRINT2 (("---\n")) ;
1046 for (p = p3 ; p < p1 + newlnz ; p++)
1047 {
1048 PRINT2 ((" p: "ID" target row "ID" %g\n",
1049 p, Li[p], Lx[p])) ;
1050 ASSERT (Li [p] > j && Li [p] < n) ;
1051 }
1052 for (s = 0 ; s < nsets ; s++)
1053 {
1054 dump_set (s, Set_ps1, Set_ps2, j, n, Common) ;
1055 }
1056 #endif
1057
1058 /* get the entry at the tail end of source column Lj */
1059 row1 = Li [p2 - 1] ;
1060 ASSERT (row1 >= j && p2 >= p1) ;
1061
1062 /* find the largest row in all the sets */
1063 maxrow = row1 ;
1064 smax = EMPTY ;
1065 for (s = nsets-1 ; s >= 0 ; s--)
1066 {
1067 ASSERT (Set_ps1 [s] < Set_ps2 [s]) ;
1068 row = *(Set_ps2 [s] - 1) ;
1069 if (row == maxrow)
1070 {
1071 /* skip past this entry in set s (it is a duplicate) */
1072 Set_ps2 [s]-- ;
1073 if (Set_ps1 [s] == Set_ps2 [s])
1074 {
1075 /* nothing more in this set */
1076 nsets-- ;
1077 Set_ps1 [s] = Set_ps1 [nsets] ;
1078 Set_ps2 [s] = Set_ps2 [nsets] ;
1079 if (smax == nsets)
1080 {
1081 /* Set smax redefined; it is now this set */
1082 smax = s ;
1083 }
1084 }
1085 }
1086 else if (row > maxrow)
1087 {
1088 maxrow = row ;
1089 smax = s ;
1090 }
1091 }
1092 ASSERT (maxrow > j) ;
1093
1094 /* move the row onto the stack of the target column */
1095 if (maxrow == row1)
1096 {
1097 /* next entry is in Lj, move to the bottom of Lj' */
1098 ASSERT (smax == EMPTY) ;
1099 p2-- ;
1100 p3-- ;
1101 Li [p3] = maxrow ;
1102 Lx [p3] = Lx [p2] ;
1103 }
1104 else
1105 {
1106 /* new entry in Lj' */
1107 ASSERT (smax >= 0 && smax < nsets) ;
1108 Set_ps2 [smax]-- ;
1109 p3-- ;
1110 Li [p3] = maxrow ;
1111 Lx [p3] = 0.0 ;
1112 if (Set_ps1 [smax] == Set_ps2 [smax])
1113 {
1114 /* nothing more in this set */
1115 nsets-- ;
1116 Set_ps1 [smax] = Set_ps1 [nsets] ;
1117 Set_ps2 [smax] = Set_ps2 [nsets] ;
1118 PRINT1 (("Set "ID" now empty\n", smax)) ;
1119 }
1120 }
1121 }
1122
1123 /* -------------------------------------------------------------- */
1124 /* 2-set merge: */
1125 /* -------------------------------------------------------------- */
1126
1127 /* This the same as the multi-set merge, except there is only one
1128 * set s = 0 left. The source column j and the set 0 are being
1129 * merged into the target column j. */
1130
1131 if (nsets > 0)
1132 {
1133 ps1 = Set_ps1 [0] ;
1134 ps2 = Set_ps2 [0] ;
1135 }
1136
1137 while (p3 > p2)
1138 {
1139
1140 #ifndef NDEBUG
1141 PRINT2 (("\n2-set merge.\n")) ;
1142 ASSERT (nsets == 1) ;
1143 PRINT2 (("Source col p1 = "ID", p2 = "ID", p3= "ID"\n",
1144 p1, p2, p3)) ;
1145 for (p = p1 + 1 ; p < p2 ; p++)
1146 {
1147 PRINT2 ((" p: "ID" source row "ID" %g\n",
1148 p, Li[p], Lx[p])) ;
1149 ASSERT (Li [p] > j && Li [p] < n) ;
1150 }
1151 PRINT2 (("---\n")) ;
1152 for (p = p3 ; p < p1 + newlnz ; p++)
1153 {
1154 PRINT2 ((" p: "ID" target row "ID" %g\n",
1155 p, Li[p], Lx[p])) ;
1156 ASSERT (Li [p] > j && Li [p] < n) ;
1157 }
1158 dump_set (0, Set_ps1, Set_ps2, j, n, Common) ;
1159 #endif
1160
1161 if (p2 == p1 + 1)
1162 {
1163 /* the top of Lj is empty; copy the set and quit */
1164 while (p3 > p2)
1165 {
1166 /* new entry in Lj' */
1167 row = *(--ps2) ;
1168 p3-- ;
1169 Li [p3] = row ;
1170 Lx [p3] = 0.0 ;
1171 }
1172 }
1173 else
1174 {
1175 /* get the entry at the tail end of Lj */
1176 row1 = Li [p2 - 1] ;
1177 ASSERT (row1 > j && row1 < n) ;
1178 /* get the entry at the tail end of the incoming set */
1179 ASSERT (ps1 < ps2) ;
1180 row = *(ps2-1) ;
1181 ASSERT (row > j && row1 < n) ;
1182 /* move the larger of the two entries to the target set */
1183 if (row1 >= row)
1184 {
1185 /* next entry is in Lj, move to the bottom */
1186 if (row1 == row)
1187 {
1188 /* skip past this entry in the set */
1189 ps2-- ;
1190 }
1191 p2-- ;
1192 p3-- ;
1193 Li [p3] = row1 ;
1194 Lx [p3] = Lx [p2] ;
1195 }
1196 else
1197 {
1198 /* new entry in Lj' */
1199 ps2-- ;
1200 p3-- ;
1201 Li [p3] = row ;
1202 Lx [p3] = 0.0 ;
1203 }
1204 }
1205 }
1206
1207 /* -------------------------------------------------------------- */
1208 /* The new column j of L is now in Li/Lx [p1 ... p2-1] */
1209 /* -------------------------------------------------------------- */
1210
1211 p2 = p1 + newlnz ;
1212 DEBUG (dump_col ("After merge: ", j, p1, p2, Li, Lx, n, Common)) ;
1213
1214 fl += Path [path].rank * (6 + 4 * (double) newlnz) ;
1215
1216 /* -------------------------------------------------------------- */
1217 /* clear Flag; original pattern of column j L no longer marked */
1218 /* -------------------------------------------------------------- */
1219
1220 mark = CHOLMOD(clear_flag) (Common) ;
1221 scattered = FALSE ;
1222
1223 /* -------------------------------------------------------------- */
1224 /* find the new parent */
1225 /* -------------------------------------------------------------- */
1226
1227 newparent = (newlnz > 1) ? (Li [p1 + 1]) : EMPTY ;
1228 PRINT1 (("\nNew parent, Lnz: "ID": "ID" "ID"\n",
1229 j, newparent,newlnz));
1230 ASSERT (oldparent == EMPTY || newparent <= oldparent) ;
1231
1232 /* -------------------------------------------------------------- */
1233 /* go to the next node in the path */
1234 /* -------------------------------------------------------------- */
1235
1236 /* path moves to (j,nextj) unless j is a root */
1237 nextj = (newparent == EMPTY) ? n : newparent ;
1238
1239 /* place path at head of list for nextj, or terminate the path */
1240 PRINT1 (("\n j = "ID" nextj = "ID"\n\n", j, nextj)) ;
1241 Path [path].c = j ;
1242 if (nextj < n)
1243 {
1244 /* put path on link list of pending paths at column nextj */
1245 Path [path].next = Head [nextj] ;
1246 Path [path].pending = nextj ;
1247 Head [nextj] = path ;
1248 PRINT1 (("Path "ID" continues to ("ID","ID"). Rank "ID"\n",
1249 path, Path [path].c, nextj, Path [path].rank)) ;
1250 }
1251 else
1252 {
1253 /* path has ended here, at a root */
1254 Path [path].next = EMPTY ;
1255 Path [path].pending = EMPTY ;
1256 Path [path].end = j ;
1257 PRINT1 (("Path "ID" ends at root ("ID"). Rank "ID"\n",
1258 path, Path [path].end, Path [path].rank)) ;
1259 }
1260
1261 /* The link list Head [j] can now be emptied. Set Head [j] to 1
1262 * if column j has changed (it is no longer used as a link list). */
1263 PRINT1 (("column "ID", oldlnz = "ID"\n", j, Lnz [j])) ;
1264 Head [j] = (Lnz [j] != newlnz) ? 1 : EMPTY ;
1265 Lnz [j] = newlnz ;
1266 PRINT1 (("column "ID", newlnz = "ID"\n", j, newlnz)) ;
1267 DEBUG (dump_col ("New", j, p1, p2, Li, Lx, n, Common)) ;
1268
1269 /* move to the next column */
1270 if (k == Path [path].rank)
1271 {
1272 /* only one path left */
1273 j = nextj ;
1274 }
1275 else
1276 {
1277 /* The current path is moving from column j to column nextj
1278 * (nextj is n if the path has ended). However, there may be
1279 * other paths pending in columns j+1 to nextj-1. There are
1280 * two methods for looking for the next column with a pending
1281 * update. The first one looks at all columns j+1 to nextj-1
1282 * for a non-empty link list. This can be costly if j and
1283 * nextj differ by a large amount (it can be O(n), but this
1284 * entire routine may take Omega(1) time). The second method
1285 * looks at all paths and finds the smallest column at which any
1286 * path is pending. It takes O(# of paths), which is bounded
1287 * by 23: one for each column of C (up to 8), and then 15 for a
1288 * balanced binary tree with 8 leaves. However, if j and
1289 * nextj differ by a tiny amount (nextj is often j+1 near
1290 * the end of the matrix), looking at columns j+1 to nextj
1291 * would be faster. Both methods give the same answer. */
1292
1293 if (nextj - j < npaths)
1294 {
1295 /* there are fewer columns to search than paths */
1296 PRINT1 (("check j="ID" to nextj="ID"\n", j, nextj)) ;
1297 for (j2 = j + 1 ; j2 < nextj ; j2++)
1298 {
1299 PRINT1 (("check j="ID" "ID"\n", j2, Head [j2])) ;
1300 if (Head [j2] != EMPTY)
1301 {
1302 PRINT1 (("found, j="ID"\n", j2)) ;
1303 ASSERT (Path [Head [j2]].pending == j2) ;
1304 break ;
1305 }
1306 }
1307 }
1308 else
1309 {
1310 /* there are fewer paths than columns to search */
1311 j2 = nextj ;
1312 for (kk = 0 ; kk < npaths ; kk++)
1313 {
1314 jj = Path [kk].pending ;
1315 PRINT2 (("Path "ID" pending at "ID"\n", kk, jj)) ;
1316 if (jj != EMPTY) j2 = MIN (j2, jj) ;
1317 }
1318 }
1319 j = j2 ;
1320 }
1321 }
1322
1323 /* ensure workspaces are back to the values required on input */
1324 ASSERT (CHOLMOD(dump_work) (TRUE, TRUE, TRUE, Common)) ;
1325
1326 /* ------------------------------------------------------------------ */
1327 /* depth-first-search of tree to order the paths */
1328 /* ------------------------------------------------------------------ */
1329
1330 /* create lists of child paths */
1331 PRINT1 (("\n\nDFS search:\n\n")) ;
1332 for (path = 0 ; path < npaths ; path++)
1333 {
1334 Path [path].c = EMPTY ; /* first child of path */
1335 Path [path].next = EMPTY ; /* next sibling of path */
1336 Path [path].order = EMPTY ; /* path is not ordered yet */
1337 Path [path].wfirst = EMPTY ; /* 1st column of W not found yet */
1338
1339 #ifndef NDEBUG
1340 j = Path [path].start ;
1341 PRINT1 (("Path "ID" : start "ID" end "ID" parent "ID" ccol "ID"\n",
1342 path, j, Path [path].end, Path [path].parent, Path [path].ccol)) ;
1343 for ( ; ; )
1344 {
1345 PRINT1 ((" column "ID"\n", j)) ;
1346 ASSERT (j == EMPTY || (j >= 0 && j < n)) ;
1347 if (j == Path [path].end)
1348 {
1349 break ;
1350 }
1351 ASSERT (j >= 0 && j < n) ;
1352 j = (Lnz [j] > 1) ? (Li [Lp [j] + 1]) : EMPTY ;
1353 }
1354 #endif
1355 }
1356
1357 for (path = 0 ; path < npaths ; path++)
1358 {
1359 p = Path [path].parent ; /* add path to child list of parent */
1360 if (p != EMPTY)
1361 {
1362 ASSERT (p < npaths) ;
1363 Path [path].next = Path [p].c ;
1364 Path [p].c = path ;
1365 }
1366 }
1367
1368 path_order = k ;
1369 w_order = 0 ;
1370 for (path = npaths-1 ; path >= 0 ; path--)
1371 {
1372 if (Path [path].order == EMPTY)
1373 {
1374 /* this path is the root of a subtree of Tbar */
1375 PRINT1 (("Root path "ID"\n", path)) ;
1376 ASSERT (path >= k) ;
1377 dfs (Path, k, path, &path_order, &w_order, 0, npaths) ;
1378 }
1379 }
1380 ASSERT (path_order == npaths) ;
1381 ASSERT (w_order == k) ;
1382
1383 /* reorder the paths */
1384 for (path = 0 ; path < npaths ; path++)
1385 {
1386 /* old order is path, new order is Path [path].order */
1387 OrderedPath [Path [path].order] = Path [path] ;
1388 }
1389
1390 #ifndef NDEBUG
1391 for (path = 0 ; path < npaths ; path++)
1392 {
1393 PRINT1 (("Ordered Path "ID": start "ID" end "ID" wfirst "ID" rank "
1394 ""ID" ccol "ID"\n", path, OrderedPath [path].start,
1395 OrderedPath [path].end, OrderedPath [path].wfirst,
1396 OrderedPath [path].rank, OrderedPath [path].ccol)) ;
1397 if (path < k)
1398 {
1399 ASSERT (OrderedPath [path].ccol >= 0) ;
1400 }
1401 else
1402 {
1403 ASSERT (OrderedPath [path].ccol == EMPTY) ;
1404 }
1405 }
1406 #endif
1407
1408 /* ------------------------------------------------------------------ */
1409 /* numeric update/downdate for all paths */
1410 /* ------------------------------------------------------------------ */
1411
1412 ASSERT (CHOLMOD(dump_work) (TRUE, TRUE, wdim, Common)) ;
1413
1414 switch (wdim)
1415 {
1416 case 1:
1417 updown_1_r (update, C, k, L, W, OrderedPath, npaths, mask,
1418 maskmark, Common) ;
1419 break ;
1420 case 2:
1421 updown_2_r (update, C, k, L, W, OrderedPath, npaths, mask,
1422 maskmark, Common) ;
1423 break ;
1424 case 4:
1425 updown_4_r (update, C, k, L, W, OrderedPath, npaths, mask,
1426 maskmark, Common) ;
1427 break ;
1428 case 8:
1429 updown_8_r (update, C, k, L, W, OrderedPath, npaths, mask,
1430 maskmark, Common) ;
1431 break ;
1432 }
1433
1434 ASSERT (CHOLMOD(dump_work) (TRUE, TRUE, wdim, Common)) ;
1435 }
1436
1437 /* ---------------------------------------------------------------------- */
1438 /* update/downdate the forward solve */
1439 /* ---------------------------------------------------------------------- */
1440
1441 if (do_solve)
1442 {
1443 /* We now have DeltaB += Lold (:,j) * X (j) for all columns j in union
1444 * of all paths seen during the entire rank-cncol update/downdate. For
1445 * each j in path, do DeltaB -= Lnew (:,j)*DeltaB(j)
1446 * in topological order. */
1447
1448 #ifndef NDEBUG
1449 PRINT1 (("\ndo_solve, DeltaB + Lold(:,Path)*X(Path):\n")) ;
1450 for (i = 0 ; i < n ; i++)
1451 {
1452 PRINT1 (("do_solve: "ID" %30.20e\n", i, Nx [i])) ;
1453 }
1454 #endif
1455
1456 /* Note that the downdate, if it deleted entries, would need to compute
1457 * the Stack prior to doing any downdates. */
1458
1459 /* find the union of all the paths in the new L */
1460 top = n ; /* "top" is stack pointer, not a row or column index */
1461 for (ccol = 0 ; ccol < cncol ; ccol++)
1462 {
1463
1464 /* -------------------------------------------------------------- */
1465 /* j = first row index of C (:,ccol) */
1466 /* -------------------------------------------------------------- */
1467
1468 pp1 = Cp [ccol] ;
1469 pp2 = (packed) ? (Cp [ccol+1]) : (pp1 + Cnz [ccol]) ;
1470 if (pp2 > pp1)
1471 {
1472 /* Column ccol of C has at least one entry. */
1473 j = Ci [pp1] ;
1474 }
1475 else
1476 {
1477 /* Column ccol of C is empty */
1478 j = n-1 ;
1479 }
1480 PRINT1 (("\ndo_solve: ccol= "ID"\n", ccol)) ;
1481 ASSERT (j >= 0 && j < n) ;
1482 len = 0 ;
1483
1484 /* -------------------------------------------------------------- */
1485 /* find the new rowmark */
1486 /* -------------------------------------------------------------- */
1487
1488 /* Each column of C can redefine the region of L that takes part in
1489 * the update/downdate of the triangular solve Lx=b. If
1490 * i = colmark [ccol] for column C(:,ccol), then i = rowmark [j] is
1491 * redefined for all columns along the path modified by C(:,ccol).
1492 * If more than one column modifies any given column j of L, then
1493 * the rowmark of j is determined by the colmark of the least-
1494 * numbered column that affects column j. That is, if both
1495 * C(:,ccol1) and C(:,ccol2) affect column j of L, then
1496 * rowmark [j] = colmark [MIN (ccol1, ccol2)].
1497 *
1498 * rowmark [j] is not modified if rowmark or colmark are NULL,
1499 * or if colmark [ccol] is EMPTY.
1500 */
1501
1502 botrow = (use_colmark) ? (colmark [ccol]) : EMPTY ;
1503
1504 /* -------------------------------------------------------------- */
1505 /* traverse from j towards root, stopping if node already visited */
1506 /* -------------------------------------------------------------- */
1507
1508 while (j != EMPTY && Flag [j] < mark)
1509 {
1510 PRINT1 (("do_solve: subpath j= "ID"\n", j)) ;
1511 ASSERT (j >= 0 && j < n) ;
1512 Stack [len++] = j ; /* place j on the stack */
1513 Flag [j] = mark ; /* flag j as visited */
1514
1515 /* if using colmark, mark column j with botrow */
1516 ASSERT (Li [Lp [j]] == j) ; /* diagonal is always present */
1517 if (use_colmark)
1518 {
1519 Li [Lp [j]] = botrow ; /* use the space for botrow */
1520 }
1521
1522 /* go up the tree, to the parent of j */
1523 j = (Lnz [j] > 1) ? (Li [Lp [j] + 1]) : EMPTY ;
1524 }
1525
1526 /* -------------------------------------------------------------- */
1527 /* move the path down to the bottom of the stack */
1528 /* -------------------------------------------------------------- */
1529
1530 ASSERT (len <= top) ;
1531 while (len > 0)
1532 {
1533 Stack [--top] = Stack [--len] ;
1534 }
1535 }
1536
1537 #ifndef NDEBUG
1538 /* Union of paths now in Stack [top..n-1] in topological order */
1539 PRINT1 (("\nTopological order:\n")) ;
1540 for (i = top ; i < n ; i++)
1541 {
1542 PRINT1 (("column "ID" in full path\n", Stack [i])) ;
1543 }
1544 #endif
1545
1546 /* Do the forward solve for the full path part of L */
1547 for (m = top ; m < n ; m++)
1548 {
1549 j = Stack [m] ;
1550 ASSERT (j >= 0 && j < n) ;
1551 PRINT1 (("do_solve: path j= "ID"\n", j)) ;
1552 p1 = Lp [j] ;
1553 lnz = Lnz [j] ;
1554 p2 = p1 + lnz ;
1555 xj = Nx [j] ;
1556
1557 /* copy new solution onto old one, for all cols in full path */
1558 Xx [j] = xj ;
1559 Nx [j] = 0. ;
1560
1561 /* DeltaB -= Lnew (j+1:botrow-1,j) * deltab(j) */
1562 if (use_colmark)
1563 {
1564 botrow = Li [p1] ; /* get botrow */
1565 Li [p1] = j ; /* restore diagonal entry */
1566 for (p = p1 + 1 ; p < p2 ; p++)
1567 {
1568 i = Li [p] ;
1569 if (i >= botrow) break ;
1570 Nx [i] -= Lx [p] * xj ;
1571 }
1572 }
1573 else
1574 {
1575 for (p = p1 + 1 ; p < p2 ; p++)
1576 {
1577 Nx [Li [p]] -= Lx [p] * xj ;
1578 }
1579 }
1580 }
1581
1582 /* clear the Flag */
1583 mark = CHOLMOD(clear_flag) (Common) ;
1584 }
1585
1586 /* ---------------------------------------------------------------------- */
1587 /* successful update/downdate */
1588 /* ---------------------------------------------------------------------- */
1589
1590 Common->modfl = fl ;
1591 DEBUG (for (j = 0 ; j < n ; j++) ASSERT (IMPLIES (do_solve, Nx[j] == 0.))) ;
1592 ASSERT (CHOLMOD(dump_work) (TRUE, TRUE, TRUE, Common)) ;
1593 DEBUG (CHOLMOD(dump_factor) (L, "output L for updown", Common)) ;
1594 return (TRUE) ;
1595 }
1596 #endif
1597 #endif
1598