//------------------------------------------------------------------------------ // CHOLMOD/Modify/t_cholmod_updown_worker: sparse Cholesky update/downdate //------------------------------------------------------------------------------ // CHOLMOD/Modify Module. Copyright (C) 2005-2023, Timothy A. Davis, // and William W. Hager. All Rights Reserved. // SPDX-License-Identifier: GPL-2.0+ //------------------------------------------------------------------------------ // This method, TEMPLATE (cholmod_updown_worker), is #included twice in // cholmod_updown.c, to create rs_cholmod_updown_worker for the single case, // and r_cholmod_updown_worker for the double case. // t_cholmod_update_wdim.c is #included four times below, to create the set of // update/downdate methods for each value of WDIM (1, 2, 4, and 8). Then each // of the four t_cholmod_update_wdim.c methods in turn #includes // t_cholmod_updown_numkr.c WDIM times. #include "cholmod_template.h" //------------------------------------------------------------------------------ // templates for each size of W //------------------------------------------------------------------------------ #undef UPDOWN #ifdef DOUBLE // double: d_updown_k_rank #define UPDOWN(k,rank) UPDOWN_METHOD(d_,k,rank) #else // single: s_updown_k_rank #define UPDOWN(k,rank) UPDOWN_METHOD(s_,k,rank) #endif #define WDIM 1 #include "t_cholmod_updown_wdim.c" #define WDIM 2 #include "t_cholmod_updown_wdim.c" #define WDIM 4 #include "t_cholmod_updown_wdim.c" #define WDIM 8 #include "t_cholmod_updown_wdim.c" //------------------------------------------------------------------------------ // debug routine //------------------------------------------------------------------------------ #ifndef NDEBUG static void TEMPLATE (dump_col) ( char *w, Int j, Int p1, Int p2, Int *Li, Real *Lx, Int n, cholmod_common *Common ) { Int p, row, lastrow ; if (CHOLMOD(dump) < -1) { // no checks if debug level is -2 or less return ; } PRINT3 (("\n\nDUMP COL==== j = "ID" %s: p1="ID" p2="ID" \n", j, w, p1,p2)); lastrow = -1 ; for (p = p1 ; p < p2 ; p++) { PRINT3 ((" "ID": ", p)) ; row = Li [p] ; PRINT3 ((""ID" ", Li [p])) ; PRINT3 (("%g ", Lx [p])) ; PRINT3 (("\n")) ; ASSERT (row > lastrow && row < n) ; lastrow = row ; } ASSERT (p1 < p2) ; ASSERT (Li [p1] == j) ; PRINT3 (("\n")) ; } #endif //------------------------------------------------------------------------------ // t_cholmod_updown_worker: for the single and double cases //------------------------------------------------------------------------------ static int TEMPLATE (cholmod_updown_worker) ( // input: Int k, // maximum rank for each update/downdate int update, // TRUE for update, FALSE for downdate cholmod_sparse *C, // the incoming sparse update Int *colmark, // array of size n. See cholmod_updown.c for details Int *mask, // size n Int maskmark, // input/output: cholmod_factor *L, // factor to modify cholmod_dense *X, // solution to Lx=b (size n-by-1) cholmod_dense *DeltaB, // change in b, zero on output cholmod_common *Common ) { //-------------------------------------------------------------------------- // get inputs //-------------------------------------------------------------------------- Int *Set_ps1 [32], *Set_ps2 [32] ; Path_type OrderedPath [32], Path [32] ; Int npaths, i, j, row, packed, ccol, p, jj, j2, kk, nextj, p1, p2, c, newlnz, newpath, path_order, w_order, scattered, path, pp1, pp2, smax, maxrow, row1, nsets, s, p3, newlnz1, Set [32], top, len, lnz, m, botrow ; DEBUG (Int oldparent) ; Int mark = Common->mark ; Int *Ci = C->i ; Int *Cp = C->p ; Int *Cnz = C->nz ; packed = C->packed ; Int cncol = C->ncol ; ASSERT (IMPLIES (!packed, Cnz != NULL)) ; Real *Xx, *Nx ; bool do_solve = (X != NULL) && (DeltaB != NULL) ; if (do_solve) { Xx = X->x ; Nx = DeltaB->x ; } else { Xx = NULL ; Nx = NULL ; } Int n = L->n ; Int *Li = L->i ; Real *Lx = L->x ; Int *Lp = L->p ; Int *Lnz = L->nz ; Int *Lnext = L->next ; ASSERT (Lnz != NULL) ; Int wdim = Power2 [k] ; // number of columns in W double fl = 0 ; bool use_colmark = (colmark != NULL) ; Int *ps1 = NULL ; Int *ps2 = NULL ; //-------------------------------------------------------------------------- // get workspace //-------------------------------------------------------------------------- Int *Flag = Common->Flag ; // size n, Flag [i] <= mark must hold Int *Head = Common->Head ; // size n, Head [i] == EMPTY must hold Real *W = Common->Xwork ; // size n-by-wdim, zero on input and output // note that Iwork [n .. 2*n-1] (i/i/l) may be in use in rowadd/rowdel: Int *Iwork = Common->Iwork ; Int *Stack = Iwork ; // size n, uninitialized //-------------------------------------------------------------------------- // entire rank-cncol update, done as a sequence of rank-k updates //-------------------------------------------------------------------------- for (Int k1 = 0 ; k1 < cncol ; k1 += k) { //---------------------------------------------------------------------- // get the next k columns of C for the update/downdate //---------------------------------------------------------------------- // the last update/downdate might be less than rank-k if (k > cncol - k1) { k = cncol - k1 ; wdim = Power2 [k] ; } Int k2 = k1 + k - 1 ; // workspaces are in the following state, on input and output ASSERT (CHOLMOD(dump_work) (TRUE, TRUE, wdim, L->dtype, Common)) ; //---------------------------------------------------------------------- // create a zero-length path for each column of W //---------------------------------------------------------------------- nextj = n ; path = 0 ; for (ccol = k1 ; ccol <= k2 ; ccol++) { PRINT1 (("Column ["ID"]: "ID"\n", path, ccol)) ; ASSERT (ccol >= 0 && ccol <= cncol) ; pp1 = Cp [ccol] ; pp2 = (packed) ? (Cp [ccol+1]) : (pp1 + Cnz [ccol]) ; // get the row index j of the first entry in C (:,ccol) if (pp2 > pp1) { // Column ccol of C has at least one entry. j = Ci [pp1] ; } else { // Column ccol of C is empty. Pretend it has one entry in // the last column with numerical value of zero. j = n-1 ; } ASSERT (j >= 0 && j < n) ; // find first column to work on nextj = MIN (nextj, j) ; Path [path].ccol = ccol ; // which column of C this path is for Path [path].start = EMPTY ; // paths for C have zero length Path [path].end = EMPTY ; Path [path].parent = EMPTY ; // no parent yet Path [path].rank = 1 ; // one column of W Path [path].c = EMPTY ; // no child of this path (case A) Path [path].next = Head [j] ; // this path is pending at col j Path [path].pending = j ; // this path is pending at col j Head [j] = path ; // this path is pending at col j PRINT1(("Path "ID" starts: start "ID" end "ID" parent "ID" c "ID"" "j "ID" ccol "ID"\n", path, Path [path].start, Path [path].end, Path [path].parent, Path [path].c, j, ccol)) ; // initialize botrow for this path Path [path].botrow = (use_colmark) ? colmark [ccol] : n ; path++ ; } // we start with paths 0 to k-1. Next one (now unused) is npaths npaths = k ; j = nextj ; ASSERT (j < n) ; scattered = FALSE ; //---------------------------------------------------------------------- // symbolic update of columns of L //---------------------------------------------------------------------- while (j < n) { ASSERT (j >= 0 && j < n && Lnz [j] > 0) ; // the old column, Li [p1..p2-1]. D (j,j) is stored in Lx [p1] p1 = Lp [j] ; newlnz = Lnz [j] ; p2 = p1 + newlnz ; #ifndef NDEBUG PRINT1 (("\n=========Column j="ID" p1 "ID" p2 "ID" lnz "ID" \n", j, p1, p2, newlnz)) ; TEMPLATE (dump_col) ("Old", j, p1, p2, Li, Lx, n, Common) ; oldparent = (Lnz [j] > 1) ? (Li [p1 + 1]) : EMPTY ; ASSERT (CHOLMOD(dump_work) (TRUE, FALSE, 0, 0, Common)) ; ASSERT (!scattered) ; PRINT1 (("Col "ID": Checking paths, npaths: "ID"\n", j, npaths)) ; for (kk = 0 ; kk < npaths ; kk++) { Int kk2, found, j3 = Path [kk].pending ; PRINT2 (("Path "ID" pending at "ID".\n", kk, j3)) ; if (j3 != EMPTY) { // Path kk must be somewhere in link list for column j3 ASSERT (Head [j3] != EMPTY) ; PRINT3 ((" List at "ID": ", j3)) ; found = FALSE ; for (kk2 = Head [j3] ; kk2 != EMPTY ; kk2 = Path [kk2].next) { PRINT3 ((""ID" ", kk2)) ; ASSERT (Path [kk2].pending == j3) ; found = found || (kk2 == kk) ; } PRINT3 (("\n")) ; ASSERT (found) ; } } PRINT1 (("\nCol "ID": Paths at this column, head "ID"\n", j, Head [j])); ASSERT (Head [j] != EMPTY) ; for (kk = Head [j] ; kk != EMPTY ; kk = Path [kk].next) { PRINT1 (("path "ID": (c="ID" j="ID") npaths "ID"\n", kk, Path[kk].c, j, npaths)) ; ASSERT (kk >= 0 && kk < npaths) ; ASSERT (Path [kk].pending == j) ; } #endif //------------------------------------------------------------------ // determine the path we're on //------------------------------------------------------------------ // get the first old path at column j path = Head [j] ; //------------------------------------------------------------------ // update/downdate of forward solve, Lx=b //------------------------------------------------------------------ if (do_solve) { Real xj = Xx [j] ; if (xj != 0) { xj = Xx [j] ; // This is first time column j has been seen for entire // rank-k update/downdate. // DeltaB += Lold (j:botrow-1,j) * X (j) Nx [j] += xj ; // diagonal of L // find the botrow for this column botrow = (use_colmark) ? Path [path].botrow : n ; for (p = p1 + 1 ; p < p2 ; p++) { i = Li [p] ; if (i >= botrow) { break ; } Nx [i] += Lx [p] * xj ; } // clear X[j] to flag col j of Lold as having been seen. If // X (j) was initially zero, then the above code is never // executed for column j. This is safe, since if xj=0 the // code above does not do anything anyway. Xx [j] = 0.0 ; } } //------------------------------------------------------------------ // start a new path at this column if two or more paths merge //------------------------------------------------------------------ newpath = // start a new path if paths have merged (Path [path].next != EMPTY) // or if j is the first node on a path (case A). || (Path [path].c == EMPTY) ; if (newpath) { // get the botrow of the first path at column j botrow = (use_colmark) ? Path [path].botrow : n ; path = npaths++ ; ASSERT (npaths <= 3*k) ; Path [path].ccol = EMPTY ; // no single col of C for this path Path [path].start = j ; // path starts at this column j Path [path].end = EMPTY ; // don't know yet where it ends Path [path].parent = EMPTY ;// don't know parent path yet Path [path].rank = 0 ; // rank is sum of child path ranks PRINT1 (("Path "ID" starts: start "ID" end "ID" parent "ID"\n", path, Path [path].start, Path [path].end, Path [path].parent)) ; // set the botrow of the new path Path [path].botrow = (use_colmark) ? botrow : n ; } //------------------------------------------------------------------ // for each path kk pending at column j //------------------------------------------------------------------ // make a list of the sets that need to be merged into column j nsets = 0 ; for (kk = Head [j] ; kk != EMPTY ; kk = Path [kk].next) { //-------------------------------------------------------------- // path kk is at (c,j) //-------------------------------------------------------------- c = Path [kk].c ; ASSERT (c < j) ; PRINT1 (("TUPLE on path "ID" (c="ID" j="ID")\n", kk, c, j)) ; ASSERT (Path [kk].pending == j) ; if (newpath) { // finalize path kk and find rank of this path Path [kk].end = c ; // end of old path is previous node c Path [kk].parent = path ; // parent is this path Path [path].rank += Path [kk].rank ; // sum up ranks Path [kk].pending = EMPTY ; PRINT1 (("Path "ID" done:start "ID" end "ID" parent "ID"\n", kk, Path [kk].start, Path [kk].end, Path [kk].parent)) ; } if (c == EMPTY) { //---------------------------------------------------------- // CASE A: first node in path //---------------------------------------------------------- // update: add pattern of incoming column // Column ccol of C is in Ci [pp1 ... pp2-1] ccol = Path [kk].ccol ; pp1 = Cp [ccol] ; pp2 = (packed) ? (Cp [ccol+1]) : (pp1 + Cnz [ccol]) ; PRINT1 (("Case A, ccol = "ID" len "ID"\n", ccol, pp2-pp1)) ; ASSERT (IMPLIES (pp2 > pp1, Ci [pp1] == j)) ; if (!scattered) { // scatter the original pattern of column j of L for (p = p1 ; p < p2 ; p++) { Flag [Li [p]] = mark ; } scattered = TRUE ; } // scatter column ccol of C (skip first entry, j) newlnz1 = newlnz ; for (p = pp1 + 1 ; p < pp2 ; p++) { row = Ci [p] ; if (Flag [row] < mark) { // this is a new entry in Lj' Flag [row] = mark ; newlnz++ ; } } if (newlnz1 != newlnz) { // column ccol of C adds something to column j of L Set [nsets++] = FLIP (ccol) ; } } else if (Head [c] == 1) { //---------------------------------------------------------- // CASE B: c is old, but changed, child of j // CASE C: new child of j //---------------------------------------------------------- // Head [c] is 1 if col c of L has new entries, // EMPTY otherwise Flag [c] = 0 ; Head [c] = EMPTY ; // update: add Lc' // column c of L is in Li [pp1 .. pp2-1] pp1 = Lp [c] ; pp2 = pp1 + Lnz [c] ; PRINT1 (("Case B/C: c = "ID"\n", c)) ; DEBUG (TEMPLATE (dump_col) ("Child", c, pp1, pp2, Li, Lx, n, Common)) ; ASSERT (j == Li [pp1 + 1]) ; // j is new parent of c if (!scattered) { // scatter the original pattern of column j of L for (p = p1 ; p < p2 ; p++) { Flag [Li [p]] = mark ; } scattered = TRUE ; } // scatter column c of L (skip first two entries, c and j) newlnz1 = newlnz ; for (p = pp1 + 2 ; p < pp2 ; p++) { row = Li [p] ; if (Flag [row] < mark) { // this is a new entry in Lj' Flag [row] = mark ; newlnz++ ; } } PRINT2 (("\n")) ; if (newlnz1 != newlnz) { // column c of L adds something to column j of L Set [nsets++] = c ; } } } //------------------------------------------------------------------ // update the pattern of column j of L //------------------------------------------------------------------ // Column j of L will be in Li/Lx [p1 .. p3-1] p3 = p1 + newlnz ; ASSERT (IMPLIES (nsets == 0, newlnz == Lnz [j])) ; PRINT1 (("p1 "ID" p2 "ID" p3 "ID" nsets "ID"\n", p1, p2, p3,nsets)); //------------------------------------------------------------------ // ensure we have enough space for the longer column //------------------------------------------------------------------ if (nsets > 0 && p3 > Lp [Lnext [j]]) { PRINT1 (("Col realloc: j "ID" newlnz "ID"\n", j, newlnz)) ; if (!CHOLMOD(reallocate_column) (j, newlnz, L, Common)) { // out of memory, L is now simplicial symbolic CHOLMOD(clear_flag) (Common) ; for (j = 0 ; j <= n ; j++) { Head [j] = EMPTY ; } ASSERT (CHOLMOD(dump_work) (TRUE, TRUE, wdim, L->dtype, Common)) ; return (FALSE) ; } // L->i and L->x may have moved. Column j has moved too Li = L->i ; Lx = L->x ; p1 = Lp [j] ; p2 = p1 + Lnz [j] ; p3 = p1 + newlnz ; } //------------------------------------------------------------------ // create set pointers //------------------------------------------------------------------ for (s = 0 ; s < nsets ; s++) { // Pattern of Set s is *(Set_ps1 [s] ... Set_ps2 [s]-1) c = Set [s] ; if (c < EMPTY) { // column ccol of C, skip first entry (j) ccol = FLIP (c) ; pp1 = Cp [ccol] ; pp2 = (packed) ? (Cp [ccol+1]) : (pp1 + Cnz [ccol]) ; ASSERT (pp2 - pp1 > 1) ; Set_ps1 [s] = &(Ci [pp1 + 1]) ; Set_ps2 [s] = &(Ci [pp2]) ; PRINT1 (("set "ID" is ccol "ID"\n", s, ccol)) ; } else { // column c of L, skip first two entries (c and j) pp1 = Lp [c] ; pp2 = pp1 + Lnz [c] ; ASSERT (Lnz [c] > 2) ; Set_ps1 [s] = &(Li [pp1 + 2]) ; Set_ps2 [s] = &(Li [pp2]) ; PRINT1 (("set "ID" is L "ID"\n", s, c)) ; } DEBUG (dump_set (s, Set_ps1, Set_ps2, j, n, Common)) ; } //------------------------------------------------------------------ // multiset merge //------------------------------------------------------------------ // Merge the sets into a single sorted set, Lj'. Before the merge // starts, column j is located in Li/Lx [p1 ... p2-1] and the // space Li/Lx [p2 ... p3-1] is empty. p1 is Lp [j], p2 is // Lp [j] + Lnz [j] (the old length of the column), and p3 is // Lp [j] + newlnz (the new and longer length of the column). // // The sets 0 to nsets-1 are defined by the Set_ps1 and Set_ps2 // pointers. Set s is located in *(Set_ps1 [s] ... Set_ps2 [s]-1). // It may be a column of C, or a column of L. All row indices i in // the sets are in the range i > j and i < n. All sets are sorted. // // The merge into column j of L is done in place. // // During the merge, p2 and p3 are updated. Li/Lx [p1..p2-1] // reflects the indices of the old column j of L that are yet to // be merged into the new column. Entries in their proper place in // the new column j of L are located in Li/Lx [p3 ... p1+newlnz-1]. // The merge finishes when p2 == p3. // // During the merge, set s consumed as it is merged into column j of // L. Its unconsumed contents are *(Set_ps1 [s] ... Set_ps2 [s]-1). // When a set is completely consumed, it is removed from the set of // sets, and nsets is decremented. // // The multiset merge and 2-set merge finishes when p2 == p3. PRINT1 (("Multiset merge p3 "ID" p2 "ID" nsets "ID"\n", p3, p2, nsets)) ; while (p3 > p2 && nsets > 1) { #ifndef NDEBUG PRINT2 (("\nMultiset merge. nsets = "ID"\n", nsets)) ; PRINT2 (("Source col p1 = "ID", p2 = "ID", p3= "ID"\n", p1, p2, p3)) ; for (p = p1 + 1 ; p < p2 ; p++) { PRINT2 ((" p: "ID" source row "ID" %g\n", p, Li[p], Lx[p])) ; ASSERT (Li [p] > j && Li [p] < n) ; } PRINT2 (("---\n")) ; for (p = p3 ; p < p1 + newlnz ; p++) { PRINT2 ((" p: "ID" target row "ID" %g\n", p, Li[p], Lx[p])) ; ASSERT (Li [p] > j && Li [p] < n) ; } for (s = 0 ; s < nsets ; s++) { dump_set (s, Set_ps1, Set_ps2, j, n, Common) ; } #endif // get the entry at the tail end of source column Lj row1 = Li [p2 - 1] ; ASSERT (row1 >= j && p2 >= p1) ; // find the largest row in all the sets maxrow = row1 ; smax = EMPTY ; for (s = nsets-1 ; s >= 0 ; s--) { ASSERT (Set_ps1 [s] < Set_ps2 [s]) ; row = *(Set_ps2 [s] - 1) ; if (row == maxrow) { // skip past this entry in set s (it is a duplicate) Set_ps2 [s]-- ; if (Set_ps1 [s] == Set_ps2 [s]) { // nothing more in this set nsets-- ; Set_ps1 [s] = Set_ps1 [nsets] ; Set_ps2 [s] = Set_ps2 [nsets] ; if (smax == nsets) { // Set smax redefined; it is now this set smax = s ; } } } else if (row > maxrow) { maxrow = row ; smax = s ; } } ASSERT (maxrow > j) ; // move the row onto the stack of the target column if (maxrow == row1) { // next entry is in Lj, move to the bottom of Lj' ASSERT (smax == EMPTY) ; p2-- ; p3-- ; Li [p3] = maxrow ; Lx [p3] = Lx [p2] ; } else { // new entry in Lj' ASSERT (smax >= 0 && smax < nsets) ; Set_ps2 [smax]-- ; p3-- ; Li [p3] = maxrow ; Lx [p3] = 0.0 ; if (Set_ps1 [smax] == Set_ps2 [smax]) { // nothing more in this set nsets-- ; Set_ps1 [smax] = Set_ps1 [nsets] ; Set_ps2 [smax] = Set_ps2 [nsets] ; PRINT1 (("Set "ID" now empty\n", smax)) ; } } } //------------------------------------------------------------------ // 2-set merge: //------------------------------------------------------------------ // This the same as the multi-set merge, except there is only one // set s = 0 left. The source column j and the set 0 are being // merged into the target column j. if (nsets > 0) { ps1 = Set_ps1 [0] ; ps2 = Set_ps2 [0] ; } while (p3 > p2) { #ifndef NDEBUG PRINT2 (("\n2-set merge.\n")) ; ASSERT (nsets == 1) ; PRINT2 (("Source col p1 = "ID", p2 = "ID", p3= "ID"\n", p1, p2, p3)) ; for (p = p1 + 1 ; p < p2 ; p++) { PRINT2 ((" p: "ID" source row "ID" %g\n", p, Li[p], Lx[p])) ; ASSERT (Li [p] > j && Li [p] < n) ; } PRINT2 (("---\n")) ; for (p = p3 ; p < p1 + newlnz ; p++) { PRINT2 ((" p: "ID" target row "ID" %g\n", p, Li[p], Lx[p])) ; ASSERT (Li [p] > j && Li [p] < n) ; } dump_set (0, Set_ps1, Set_ps2, j, n, Common) ; #endif if (p2 == p1 + 1) { // the top of Lj is empty; copy the set and quit while (p3 > p2) { // new entry in Lj' row = *(--ps2) ; p3-- ; Li [p3] = row ; Lx [p3] = 0.0 ; } } else { // get the entry at the tail end of Lj row1 = Li [p2 - 1] ; ASSERT (row1 > j && row1 < n) ; // get the entry at the tail end of the incoming set ASSERT (ps1 < ps2) ; row = *(ps2-1) ; ASSERT (row > j && row1 < n) ; // move the larger of the two entries to the target set if (row1 >= row) { // next entry is in Lj, move to the bottom if (row1 == row) { // skip past this entry in the set ps2-- ; } p2-- ; p3-- ; Li [p3] = row1 ; Lx [p3] = Lx [p2] ; } else { // new entry in Lj' ps2-- ; p3-- ; Li [p3] = row ; Lx [p3] = 0.0 ; } } } //------------------------------------------------------------------ // The new column j of L is now in Li/Lx [p1 ... p2-1] //------------------------------------------------------------------ p2 = p1 + newlnz ; DEBUG (TEMPLATE (dump_col) ("After merge: ", j, p1, p2, Li, Lx, n, Common)) ; fl += Path [path].rank * (6 + 4 * (double) newlnz) ; //------------------------------------------------------------------ // clear Flag; original pattern of column j L no longer marked //------------------------------------------------------------------ mark = CHOLMOD(clear_flag) (Common) ; scattered = FALSE ; //------------------------------------------------------------------ // find the new parent //------------------------------------------------------------------ Int newparent = (newlnz > 1) ? (Li [p1 + 1]) : EMPTY ; PRINT1 (("\nNew parent, Lnz: "ID": "ID" "ID"\n", j, newparent,newlnz)); ASSERT (oldparent == EMPTY || newparent <= oldparent) ; //------------------------------------------------------------------ // go to the next node in the path //------------------------------------------------------------------ // path moves to (j,nextj) unless j is a root nextj = (newparent == EMPTY) ? n : newparent ; // place path at head of list for nextj, or terminate the path PRINT1 (("\n j = "ID" nextj = "ID"\n\n", j, nextj)) ; Path [path].c = j ; if (nextj < n) { // put path on link list of pending paths at column nextj Path [path].next = Head [nextj] ; Path [path].pending = nextj ; Head [nextj] = path ; PRINT1 (("Path "ID" continues to ("ID","ID"). Rank "ID"\n", path, Path [path].c, nextj, Path [path].rank)) ; } else { // path has ended here, at a root Path [path].next = EMPTY ; Path [path].pending = EMPTY ; Path [path].end = j ; PRINT1 (("Path "ID" ends at root ("ID"). Rank "ID"\n", path, Path [path].end, Path [path].rank)) ; } // The link list Head [j] can now be emptied. Set Head [j] to 1 // if column j has changed (it is no longer used as a link list). PRINT1 (("column "ID", oldlnz = "ID"\n", j, Lnz [j])) ; Head [j] = (Lnz [j] != newlnz) ? 1 : EMPTY ; Lnz [j] = newlnz ; PRINT1 (("column "ID", newlnz = "ID"\n", j, newlnz)) ; DEBUG (TEMPLATE (dump_col) ("New", j, p1, p2, Li, Lx, n, Common)) ; // move to the next column if (k == Path [path].rank) { // only one path left j = nextj ; } else { // The current path is moving from column j to column nextj // (nextj is n if the path has ended). However, there may be // other paths pending in columns j+1 to nextj-1. There are // two methods for looking for the next column with a pending // update. The first one looks at all columns j+1 to nextj-1 // for a non-empty link list. This can be costly if j and // nextj differ by a large amount (it can be O(n), but this // entire routine may take Omega(1) time). The second method // looks at all paths and finds the smallest column at which any // path is pending. It takes O(# of paths), which is bounded // by 23: one for each column of C (up to 8), and then 15 for a // balanced binary tree with 8 leaves. However, if j and // nextj differ by a tiny amount (nextj is often j+1 near // the end of the matrix), looking at columns j+1 to nextj // would be faster. Both methods give the same answer. if (nextj - j < npaths) { // there are fewer columns to search than paths PRINT1 (("check j="ID" to nextj="ID"\n", j, nextj)) ; for (j2 = j + 1 ; j2 < nextj ; j2++) { PRINT1 (("check j="ID" "ID"\n", j2, Head [j2])) ; if (Head [j2] != EMPTY) { PRINT1 (("found, j="ID"\n", j2)) ; ASSERT (Path [Head [j2]].pending == j2) ; break ; } } } else { // there are fewer paths than columns to search j2 = nextj ; for (kk = 0 ; kk < npaths ; kk++) { jj = Path [kk].pending ; PRINT2 (("Path "ID" pending at "ID"\n", kk, jj)) ; if (jj != EMPTY) j2 = MIN (j2, jj) ; } } j = j2 ; } } // ensure workspaces are back to the values required on input ASSERT (CHOLMOD(dump_work) (TRUE, TRUE, 0, 0, Common)) ; //---------------------------------------------------------------------- // depth-first-search of tree to order the paths //---------------------------------------------------------------------- // create lists of child paths PRINT1 (("\n\nDFS search:\n\n")) ; for (path = 0 ; path < npaths ; path++) { Path [path].c = EMPTY ; // first child of path Path [path].next = EMPTY ; // next sibling of path Path [path].order = EMPTY ; // path is not ordered yet Path [path].wfirst = EMPTY ; // 1st column of W not found yet #ifndef NDEBUG j = Path [path].start ; PRINT1 (("Path "ID" : start "ID" end "ID" parent "ID" ccol "ID"\n", path, j, Path [path].end, Path [path].parent, Path [path].ccol)) ; for ( ; ; ) { PRINT1 ((" column "ID"\n", j)) ; ASSERT (j == EMPTY || (j >= 0 && j < n)) ; if (j == Path [path].end) { break ; } ASSERT (j >= 0 && j < n) ; j = (Lnz [j] > 1) ? (Li [Lp [j] + 1]) : EMPTY ; } #endif } for (path = 0 ; path < npaths ; path++) { p = Path [path].parent ; // add path to child list of parent if (p != EMPTY) { ASSERT (p < npaths) ; Path [path].next = Path [p].c ; Path [p].c = path ; } } path_order = k ; w_order = 0 ; for (path = npaths-1 ; path >= 0 ; path--) { if (Path [path].order == EMPTY) { // this path is the root of a subtree of Tbar PRINT1 (("Root path "ID"\n", path)) ; ASSERT (path >= k) ; dfs (Path, k, path, &path_order, &w_order, 0, npaths) ; } } ASSERT (path_order == npaths) ; ASSERT (w_order == k) ; // reorder the paths for (path = 0 ; path < npaths ; path++) { // old order is path, new order is Path [path].order OrderedPath [Path [path].order] = Path [path] ; } #ifndef NDEBUG for (path = 0 ; path < npaths ; path++) { PRINT1 (("Ordered Path "ID": start "ID" end "ID" wfirst "ID" rank " ""ID" ccol "ID"\n", path, OrderedPath [path].start, OrderedPath [path].end, OrderedPath [path].wfirst, OrderedPath [path].rank, OrderedPath [path].ccol)) ; if (path < k) { ASSERT (OrderedPath [path].ccol >= 0) ; } else { ASSERT (OrderedPath [path].ccol == EMPTY) ; } } #endif //---------------------------------------------------------------------- // numeric update/downdate for all paths //---------------------------------------------------------------------- ASSERT (CHOLMOD(dump_work) (TRUE, TRUE, wdim, L->dtype, Common)) ; switch (wdim) { case 1: UPDOWN (1,r) (update, C, k, L, W, OrderedPath, npaths, mask, maskmark, Common) ; break ; case 2: UPDOWN (2,r) (update, C, k, L, W, OrderedPath, npaths, mask, maskmark, Common) ; break ; case 4: UPDOWN (4,r) (update, C, k, L, W, OrderedPath, npaths, mask, maskmark, Common) ; break ; case 8: UPDOWN (8,r) (update, C, k, L, W, OrderedPath, npaths, mask, maskmark, Common) ; break ; } ASSERT (CHOLMOD(dump_work) (TRUE, TRUE, wdim, L->dtype, Common)) ; } //-------------------------------------------------------------------------- // update/downdate the forward solve //-------------------------------------------------------------------------- if (do_solve) { // We now have DeltaB += Lold (:,j) * X (j) for all columns j in union // of all paths seen during the entire rank-cncol update/downdate. For // each j in path, do DeltaB -= Lnew (:,j)*DeltaB(j) // in topological order. #ifndef NDEBUG PRINT1 (("\ndo_solve, DeltaB + Lold(:,Path)*X(Path):\n")) ; for (i = 0 ; i < n ; i++) { PRINT1 (("do_solve: "ID" %30.20e\n", i, Nx [i])) ; } #endif // Note that the downdate, if it deleted entries, would need to compute // the Stack prior to doing any downdates. // find the union of all the paths in the new L top = n ; // "top" is stack pointer, not a row or column index for (ccol = 0 ; ccol < cncol ; ccol++) { //------------------------------------------------------------------ // j = first row index of C (:,ccol) //------------------------------------------------------------------ pp1 = Cp [ccol] ; pp2 = (packed) ? (Cp [ccol+1]) : (pp1 + Cnz [ccol]) ; if (pp2 > pp1) { // Column ccol of C has at least one entry. j = Ci [pp1] ; } else { // Column ccol of C is empty j = n-1 ; } PRINT1 (("\ndo_solve: ccol= "ID"\n", ccol)) ; ASSERT (j >= 0 && j < n) ; len = 0 ; //------------------------------------------------------------------ // find the new rowmark //------------------------------------------------------------------ // Each column of C can redefine the region of L that takes part in // the update/downdate of the triangular solve Lx=b. If // i = colmark [ccol] for column C(:,ccol), then i = rowmark [j] is // redefined for all columns along the path modified by C(:,ccol). // If more than one column modifies any given column j of L, then // the rowmark of j is determined by the colmark of the least- // numbered column that affects column j. That is, if both // C(:,ccol1) and C(:,ccol2) affect column j of L, then // rowmark [j] = colmark [MIN (ccol1, ccol2)]. // // rowmark [j] is not modified if rowmark or colmark are NULL, // or if colmark [ccol] is EMPTY. botrow = (use_colmark) ? (colmark [ccol]) : EMPTY ; //------------------------------------------------------------------ // traverse from j towards root, stopping if node already visited //------------------------------------------------------------------ while (j != EMPTY && Flag [j] < mark) { PRINT1 (("do_solve: subpath j= "ID"\n", j)) ; ASSERT (j >= 0 && j < n) ; Stack [len++] = j ; // place j on the stack Flag [j] = mark ; // flag j as visited // if using colmark, mark column j with botrow ASSERT (Li [Lp [j]] == j) ; // diagonal is always present if (use_colmark) { Li [Lp [j]] = botrow ; // use the space for botrow } // go up the tree, to the parent of j j = (Lnz [j] > 1) ? (Li [Lp [j] + 1]) : EMPTY ; } //------------------------------------------------------------------ // move the path down to the bottom of the stack //------------------------------------------------------------------ ASSERT (len <= top) ; while (len > 0) { Stack [--top] = Stack [--len] ; } } #ifndef NDEBUG // Union of paths now in Stack [top..n-1] in topological order PRINT1 (("\nTopological order:\n")) ; for (i = top ; i < n ; i++) { PRINT1 (("column "ID" in full path\n", Stack [i])) ; } #endif // Do the forward solve for the full path part of L for (m = top ; m < n ; m++) { j = Stack [m] ; ASSERT (j >= 0 && j < n) ; PRINT1 (("do_solve: path j= "ID"\n", j)) ; p1 = Lp [j] ; lnz = Lnz [j] ; p2 = p1 + lnz ; Real xj = Nx [j] ; // copy new solution onto old one, for all cols in full path Xx [j] = xj ; Nx [j] = 0. ; // DeltaB -= Lnew (j+1:botrow-1,j) * deltab(j) if (use_colmark) { botrow = Li [p1] ; // get botrow Li [p1] = j ; // restore diagonal entry for (p = p1 + 1 ; p < p2 ; p++) { i = Li [p] ; if (i >= botrow) break ; Nx [i] -= Lx [p] * xj ; } } else { for (p = p1 + 1 ; p < p2 ; p++) { Nx [Li [p]] -= Lx [p] * xj ; } } } // clear the Flag mark = CHOLMOD(clear_flag) (Common) ; } //-------------------------------------------------------------------------- // successful update/downdate //-------------------------------------------------------------------------- Common->modfl = fl ; DEBUG (for (j = 0 ; j < n ; j++) ASSERT (IMPLIES (do_solve, Nx[j] == 0.))) ; return (TRUE) ; } #undef PATTERN #undef REAL #undef COMPLEX #undef ZOMPLEX