Amesos Package Browser (Single Doxygen Collection) Development
amesos_klu_kernel.c
Go to the documentation of this file.
00001 /* ========================================================================== */
00002 /* === KLU_kernel =========================================================== */
00003 /* ========================================================================== */
00004 
00005 /* Sparse left-looking LU factorization, with partial pivoting.  Based on
00006  * Gilbert & Peierl's method, with a non-recursive DFS and with Eisenstat &
00007  * Liu's symmetric pruning.  No user-callable routines are in this file.
00008  */
00009 
00010 #include "amesos_klu_internal.h"
00011 
00012 /* ========================================================================== */
00013 /* === dfs ================================================================== */
00014 /* ========================================================================== */
00015 
00016 /* Does a depth-first-search, starting at node j. */
00017 
00018 static Int dfs
00019 (
00020     /* input, not modified on output: */
00021     Int j,    /* node at which to start the DFS */
00022     Int k,    /* mark value, for the Flag array */
00023     Int Pinv [ ], /* Pinv [i] = k if row i is kth pivot row, or EMPTY if
00024        * row i is not yet pivotal.  */
00025     Int Llen [ ], /* size n, Llen [k] = # nonzeros in column k of L */
00026     Int Lip [ ],  /* size n, Lip [k] is position in LU of column k of L */
00027 
00028     /* workspace, not defined on input or output */
00029     Int Stack [ ],  /* size n */
00030 
00031     /* input/output: */
00032     Int Flag [ ], /* Flag [i] == k means i is marked */
00033     Int Lpend [ ],  /* for symmetric pruning */
00034     Int top,    /* top of stack on input*/
00035     Unit LU [],
00036     Int *Lik,   /* Li row index array of the kth column */
00037     Int *plength,
00038 
00039     /* other, not defined on input or output */
00040     Int Ap_pos [ ]  /* keeps track of position in adj list during DFS */
00041 )
00042 {
00043     Int i, pos, jnew, head, l_length ;
00044     Int *Li ;
00045 
00046     l_length = *plength ;
00047 
00048     head = 0 ;
00049     Stack [0] = j ;
00050     ASSERT (Flag [j] != k) ;
00051 
00052     while (head >= 0)
00053     {
00054   j = Stack [head] ;
00055   jnew = Pinv [j] ;
00056   ASSERT (jnew >= 0 && jnew < k) ;  /* j is pivotal */
00057 
00058   if (Flag [j] != k)      /* a node is not yet visited */
00059   {
00060       /* first time that j has been visited */
00061       Flag [j] = k ;
00062       PRINTF (("[ start dfs at %d : new %d\n", j, jnew)) ;
00063       /* set Ap_pos [head] to one past the last entry in col j to scan */
00064       Ap_pos [head] =
00065     (Lpend [jnew] == EMPTY) ?  Llen [jnew] : Lpend [jnew] ;
00066   }
00067 
00068   /* add the adjacent nodes to the recursive stack by iterating through
00069    * until finding another non-visited pivotal node */
00070   Li = (Int *) (LU + Lip [jnew]) ;
00071   for (pos = --Ap_pos [head] ; pos >= 0 ; --pos)
00072   {
00073       i = Li [pos] ;
00074       if (Flag [i] != k)
00075       {
00076     /* node i is not yet visited */
00077     if (Pinv [i] >= 0)
00078     {
00079         /* keep track of where we left off in the scan of the
00080          * adjacency list of node j so we can restart j where we
00081          * left off. */
00082         Ap_pos [head] = pos ;
00083 
00084         /* node i is pivotal; push it onto the recursive stack
00085          * and immediately break so we can recurse on node i. */
00086         Stack [++head] = i ;
00087         break ;
00088     }
00089     else
00090     {
00091         /* node i is not pivotal (no outgoing edges). */
00092         /* Flag as visited and store directly into L,
00093          * and continue with current node j. */
00094             Flag [i] = k ;
00095         Lik [l_length] = i ;
00096           l_length++ ;
00097     }
00098       }
00099   }
00100 
00101   if (pos == -1)
00102   {
00103       /* if all adjacent nodes of j are already visited, pop j from
00104        * recursive stack and push j onto output stack */
00105       head-- ;
00106       Stack[--top] = j ;
00107       PRINTF (("  end   dfs at %d ] head : %d\n", j, head)) ;
00108   }
00109     }
00110 
00111     *plength = l_length ;
00112     return (top) ;
00113 }
00114 
00115 
00116 /* ========================================================================== */
00117 /* === lsolve_symbolic ====================================================== */
00118 /* ========================================================================== */
00119 
00120 /* Finds the pattern of x, for the solution of Lx=b */
00121 
00122 static Int lsolve_symbolic
00123 (
00124     /* input, not modified on output: */
00125     Int n,              /* L is n-by-n, where n >= 0 */
00126     Int k,    /* also used as the mark value, for the Flag array */
00127     Int Ap [ ],
00128     Int Ai [ ],
00129     Int Q [ ],
00130     Int Pinv [ ], /* Pinv [i] = k if i is kth pivot row, or EMPTY if row i
00131        * is not yet pivotal.  */
00132 
00133     /* workspace, not defined on input or output */
00134     Int Stack [ ],  /* size n */
00135 
00136     /* workspace, defined on input and output */
00137     Int Flag [ ], /* size n.  Initially, all of Flag [0..n-1] < k.  After
00138        * lsolve_symbolic is done, Flag [i] == k if i is in
00139        * the pattern of the output, and Flag [0..n-1] <= k. */
00140 
00141     /* other */
00142     Int Lpend [ ],  /* for symmetric pruning */
00143     Int Ap_pos [ ], /* workspace used in dfs */
00144 
00145     Unit LU [ ],  /* LU factors (pattern and values) */
00146     Int lup,    /* pointer to free space in LU */
00147     Int Llen [ ], /* size n, Llen [k] = # nonzeros in column k of L */
00148     Int Lip [ ],  /* size n, Lip [k] is position in LU of column k of L */
00149 
00150     /* ---- the following are only used in the BTF case --- */
00151 
00152     Int k1,   /* the block of A is from k1 to k2-1 */
00153     Int PSinv [ ] /* inverse of P from symbolic factorization */
00154 )
00155 {
00156     Int *Lik ;
00157     Int i, p, pend, oldcol, kglobal, top, l_length ;
00158 
00159     top = n ;
00160     l_length = 0 ;
00161     Lik = (Int *) (LU + lup);
00162 
00163     /* ---------------------------------------------------------------------- */
00164     /* BTF factorization of A (k1:k2-1, k1:k2-1) */
00165     /* ---------------------------------------------------------------------- */
00166 
00167     kglobal = k + k1 ;  /* column k of the block is col kglobal of A */
00168     oldcol = Q [kglobal] ;  /* Q must be present for BTF case */
00169     pend = Ap [oldcol+1] ;
00170     for (p = Ap [oldcol] ; p < pend ; p++)
00171     {
00172   i = PSinv [Ai [p]] - k1 ;
00173   if (i < 0) continue ; /* skip entry outside the block */
00174 
00175   /* (i,k) is an entry in the block.  start a DFS at node i */
00176   PRINTF (("\n ===== DFS at node %d in b, inew: %d\n", i, Pinv [i])) ;
00177   if (Flag [i] != k)
00178   {
00179       if (Pinv [i] >= 0)
00180       {
00181     top = dfs (i, k, Pinv, Llen, Lip, Stack, Flag,
00182          Lpend, top, LU, Lik, &l_length, Ap_pos) ;
00183       }
00184       else
00185       {
00186     /* i is not pivotal, and not flagged. Flag and put in L */
00187     Flag [i] = k ;
00188     Lik [l_length] = i ;
00189     l_length++;
00190       }
00191   }
00192     }
00193 
00194     /* If Llen [k] is zero, the matrix is structurally singular */
00195     Llen [k] = l_length ;
00196     return (top) ;
00197 }
00198 
00199 
00200 /* ========================================================================== */
00201 /* === construct_column ===================================================== */
00202 /* ========================================================================== */
00203 
00204 /* Construct the kth column of A, and the off-diagonal part, if requested.
00205  * Scatter the numerical values into the workspace X, and construct the
00206  * corresponding column of the off-diagonal matrix. */
00207 
00208 static void construct_column
00209 (
00210     /* inputs, not modified on output */
00211     Int k,      /* the column of A (or the column of the block) to get */
00212     Int Ap [ ],
00213     Int Ai [ ],
00214     Entry Ax [ ],
00215     Int Q [ ],      /* column pre-ordering */
00216 
00217     /* zero on input, modified on output */
00218     Entry X [ ],
00219 
00220     /* ---- the following are only used in the BTF case --- */
00221 
00222     /* inputs, not modified on output */
00223     Int k1,     /* the block of A is from k1 to k2-1 */
00224     Int PSinv [ ],  /* inverse of P from symbolic factorization */
00225     double Rs [ ],  /* scale factors for A */
00226     Int scale,      /* 0: no scaling, nonzero: scale the rows with Rs */
00227 
00228     /* inputs, modified on output */
00229     Int Offp [ ],   /* off-diagonal matrix (modified by this routine) */
00230     Int Offi [ ],
00231     Entry Offx [ ]
00232 )
00233 {
00234     Entry aik ;
00235     Int i, p, pend, oldcol, kglobal, poff, oldrow ;
00236 
00237     /* ---------------------------------------------------------------------- */
00238     /* Scale and scatter the column into X. */
00239     /* ---------------------------------------------------------------------- */
00240 
00241     kglobal = k + k1 ;    /* column k of the block is col kglobal of A */
00242     poff = Offp [kglobal] ; /* start of off-diagonal column */
00243     oldcol = Q [kglobal] ;
00244     pend = Ap [oldcol+1] ;
00245 
00246     if (scale <= 0)
00247     {
00248   /* no scaling */
00249   for (p = Ap [oldcol] ; p < pend ; p++)
00250   {
00251       oldrow = Ai [p] ;
00252       i = PSinv [oldrow] - k1 ;
00253       aik = Ax [p] ;
00254       if (i < 0)
00255       {
00256     /* this is an entry in the off-diagonal part */
00257     Offi [poff] = oldrow ;
00258     Offx [poff] = aik ;
00259     poff++ ;
00260       }
00261       else
00262       {
00263     /* (i,k) is an entry in the block.  scatter into X */
00264     X [i] = aik ;
00265       }
00266   }
00267     }
00268     else
00269     {
00270   /* row scaling */
00271   for (p = Ap [oldcol] ; p < pend ; p++)
00272   {
00273       oldrow = Ai [p] ;
00274       i = PSinv [oldrow] - k1 ;
00275       aik = Ax [p] ;
00276       SCALE_DIV (aik, Rs [oldrow]) ;
00277       if (i < 0)
00278       {
00279     /* this is an entry in the off-diagonal part */
00280     Offi [poff] = oldrow ;
00281     Offx [poff] = aik ;
00282     poff++ ;
00283       }
00284       else
00285       {
00286     /* (i,k) is an entry in the block.  scatter into X */
00287     X [i] = aik ;
00288       }
00289   }
00290     }
00291 
00292     Offp [kglobal+1] = poff ;   /* start of the next col of off-diag part */
00293 }
00294 
00295 
00296 /* ========================================================================== */
00297 /* === lsolve_numeric ======================================================= */
00298 /* ========================================================================== */
00299 
00300 /* Computes the numerical values of x, for the solution of Lx=b.  Note that x
00301  * may include explicit zeros if numerical cancelation occurs.  L is assumed
00302  * to be unit-diagonal, with possibly unsorted columns (but the first entry in
00303  * the column must always be the diagonal entry). */
00304 
00305 static void lsolve_numeric
00306 (
00307     /* input, not modified on output: */
00308     Int Pinv [ ], /* Pinv [i] = k if i is kth pivot row, or EMPTY if row i
00309        * is not yet pivotal.  */
00310     Unit *LU,   /* LU factors (pattern and values) */
00311     Int Stack [ ],  /* stack for dfs */
00312     Int Lip [ ],  /* size n, Lip [k] is position in LU of column k of L */
00313     Int top,    /* top of stack on input */
00314     Int n,    /* A is n-by-n */
00315     Int Llen [ ], /* size n, Llen [k] = # nonzeros in column k of L */
00316 
00317     /* output, must be zero on input: */
00318     Entry X [ ] /* size n, initially zero.  On output,
00319      * X [Ui [up1..up-1]] and X [Li [lp1..lp-1]]
00320      * contains the solution. */
00321 
00322 )
00323 {
00324     Entry xj ;
00325     Entry *Lx ;
00326     Int *Li ;
00327     Int p, s, j, jnew, len ;
00328 
00329     /* solve Lx=b */
00330     for (s = top ; s < n ; s++)
00331     {
00332   /* forward solve with column j of L */
00333   j = Stack [s] ;
00334   jnew = Pinv [j] ;
00335   ASSERT (jnew >= 0) ;
00336   xj = X [j] ;
00337   GET_POINTER (LU, Lip, Llen, Li, Lx, jnew, len) ;
00338   ASSERT (Lip [jnew] <= Lip [jnew+1]) ;
00339   for (p = 0 ; p < len ; p++)
00340   {
00341       /*X [Li [p]] -= Lx [p] * xj ; */
00342       MULT_SUB (X [Li [p]], Lx [p], xj) ;
00343   }
00344     }
00345 }
00346 
00347 
00348 /* ========================================================================== */
00349 /* === lpivot =============================================================== */
00350 /* ========================================================================== */
00351 
00352 /* Find a pivot via partial pivoting, and scale the column of L. */
00353 
00354 static Int lpivot
00355 (
00356     Int diagrow,
00357     Int *p_pivrow,
00358     Entry *p_pivot,
00359     double *p_abs_pivot,
00360     double tol,
00361     Entry X [ ],
00362     Unit *LU,   /* LU factors (pattern and values) */
00363     Int Lip [ ],
00364     Int Llen [ ],
00365     Int k,
00366     Int n,
00367 
00368     Int Pinv [ ], /* Pinv [i] = k if row i is kth pivot row, or EMPTY if
00369        * row i is not yet pivotal.  */
00370 
00371     Int *p_firstrow,
00372     KLU_common *Common
00373 )
00374 {
00375     Entry x, pivot, *Lx ;
00376     double abs_pivot, xabs ;
00377     Int p, i, ppivrow, pdiag, pivrow, *Li, last_row_index, firstrow, len ;
00378 
00379     pivrow = EMPTY ;
00380     if (Llen [k] == 0)
00381     {
00382   /* matrix is structurally singular */
00383   if (Common->halt_if_singular)
00384   {
00385       return (FALSE) ;
00386   }
00387   for (firstrow = *p_firstrow ; firstrow < n ; firstrow++)
00388   {
00389       PRINTF (("check %d\n", firstrow)) ;
00390       if (Pinv [firstrow] < 0)
00391       {
00392     /* found the lowest-numbered non-pivotal row.  Pick it. */
00393     pivrow = firstrow ;
00394     PRINTF (("Got pivotal row: %d\n", pivrow)) ;
00395     break ;
00396       }
00397   }
00398   ASSERT (pivrow >= 0 && pivrow < n) ;
00399   CLEAR (pivot) ;
00400   *p_pivrow = pivrow ;
00401   *p_pivot = pivot ;
00402   *p_abs_pivot = 0 ;
00403   *p_firstrow = firstrow ;
00404   return (FALSE) ;
00405     }
00406 
00407     pdiag = EMPTY ;
00408     ppivrow = EMPTY ;
00409     abs_pivot = EMPTY ;
00410     i = Llen [k] - 1 ;
00411     GET_POINTER (LU, Lip, Llen, Li, Lx, k, len) ;
00412     last_row_index = Li [i] ;
00413 
00414     /* decrement the length by 1 */
00415     Llen [k] = i ;
00416     GET_POINTER (LU, Lip, Llen, Li, Lx, k, len) ;
00417 
00418     /* look in Li [0 ..Llen [k] - 1 ] for a pivot row */
00419     for (p = 0 ; p < len ; p++)
00420     {
00421   /* gather the entry from X and store in L */
00422   i = Li [p] ;
00423   x = X [i] ;
00424   CLEAR (X [i]) ;
00425 
00426   Lx [p] = x ;
00427   /* xabs = ABS (x) ; */
00428   ABS (xabs, x) ;
00429 
00430   /* find the diagonal */
00431   if (i == diagrow)
00432   {
00433       pdiag = p ;
00434   }
00435 
00436   /* find the partial-pivoting choice */
00437   if (xabs > abs_pivot)
00438   {
00439       abs_pivot = xabs ;
00440       ppivrow = p ;
00441   }
00442     }
00443 
00444     /* xabs = ABS (X [last_row_index]) ;*/
00445     ABS (xabs, X [last_row_index]) ;
00446     if (xabs > abs_pivot)
00447     {
00448         abs_pivot = xabs ;
00449         ppivrow = EMPTY ;
00450     }
00451 
00452     /* compare the diagonal with the largest entry */
00453     if (last_row_index == diagrow)
00454     {
00455   if (xabs >= tol * abs_pivot)
00456   {
00457           abs_pivot = xabs ;
00458             ppivrow = EMPTY ;
00459         }
00460     }
00461     else if (pdiag != EMPTY)
00462     {
00463   /* xabs = ABS (Lx [pdiag]) ;*/
00464   ABS (xabs, Lx [pdiag]) ;
00465   if (xabs >= tol * abs_pivot)
00466   {
00467       /* the diagonal is large enough */
00468       abs_pivot = xabs ;
00469       ppivrow = pdiag ;
00470   }
00471     }
00472 
00473     if (ppivrow != EMPTY)
00474     {
00475         pivrow = Li [ppivrow] ;
00476         pivot  = Lx [ppivrow] ;
00477   /* overwrite the ppivrow values with last index values */
00478         Li [ppivrow] = last_row_index ;
00479         Lx [ppivrow] = X [last_row_index] ;
00480     }
00481     else
00482     {
00483         pivrow = last_row_index ;
00484         pivot = X [last_row_index] ;
00485     }
00486     CLEAR (X [last_row_index]) ;
00487 
00488     *p_pivrow = pivrow ;
00489     *p_pivot = pivot ;
00490     *p_abs_pivot = abs_pivot ;
00491     ASSERT (pivrow >= 0 && pivrow < n) ;
00492 
00493     if (IS_ZERO (pivot) && Common->halt_if_singular)
00494     {
00495   /* numerically singular case */
00496   return (FALSE) ;
00497     }
00498 
00499     /* divide L by the pivot value */
00500     for (p = 0 ; p < Llen [k] ; p++)
00501     {
00502   /* Lx [p] /= pivot ; */
00503   DIV (Lx [p], Lx [p], pivot) ;
00504     }
00505 
00506     return (TRUE) ;
00507 }
00508 
00509 
00510 /* ========================================================================== */
00511 /* === prune ================================================================ */
00512 /* ========================================================================== */
00513 
00514 /* Prune the columns of L to reduce work in subsequent depth-first searches */
00515 static void prune
00516 (
00517     /* input/output: */
00518     Int Lpend [ ],  /* Lpend [j] marks symmetric pruning point for L(:,j) */
00519 
00520     /* input: */
00521     Int Pinv [ ], /* Pinv [i] = k if row i is kth pivot row, or EMPTY if
00522        * row i is not yet pivotal.  */
00523     Int k,    /* prune using column k of U */
00524     Int pivrow,   /* current pivot row */
00525 
00526     /* input/output: */
00527     Unit *LU,   /* LU factors (pattern and values) */
00528 
00529     /* input */
00530     Int Uip [ ],  /* size n, column pointers for U */
00531     Int Lip [ ],  /* size n, column pointers for L */
00532     Int Ulen [ ], /* size n, column length of U */
00533     Int Llen [ ]  /* size n, column length of L */
00534 )
00535 {
00536     Entry x ;
00537     Entry *Lx, *Ux ;
00538     Int *Li, *Ui ;
00539     Int p, i, j, p2, phead, ptail, llen, ulen ;
00540 
00541     /* check to see if any column of L can be pruned */
00542     GET_POINTER (LU, Uip, Ulen, Ui, Ux, k, ulen) ;
00543     for (p = 0 ; p < ulen ; p++)
00544     {
00545   j = Ui [p] ;
00546   ASSERT (j < k) ;
00547   PRINTF (("%d is pruned: %d. Lpend[j] %d Lip[j+1] %d\n",
00548       j, Lpend [j] != EMPTY, Lpend [j], Lip [j+1])) ;
00549   if (Lpend [j] == EMPTY)
00550   {
00551       /* scan column j of L for the pivot row */
00552             GET_POINTER (LU, Lip, Llen, Li, Lx, j, llen) ;
00553       for (p2 = 0 ; p2 < llen ; p2++)
00554       {
00555     if (pivrow == Li [p2])
00556     {
00557         /* found it!  This column can be pruned */
00558 #ifndef NDEBUG
00559         PRINTF (("==== PRUNE: col j %d of L\n", j)) ;
00560         {
00561       Int p3 ;
00562       for (p3 = 0 ; p3 < Llen [j] ; p3++)
00563       {
00564           PRINTF (("before: %i  pivotal: %d\n", Li [p3],
00565           Pinv [Li [p3]] >= 0)) ;
00566       }
00567         }
00568 #endif
00569 
00570         /* partition column j of L.  The unit diagonal of L
00571          * is not stored in the column of L. */
00572         phead = 0 ;
00573         ptail = Llen [j] ;
00574         while (phead < ptail)
00575         {
00576       i = Li [phead] ;
00577       if (Pinv [i] >= 0)
00578       {
00579           /* leave at the head */
00580           phead++ ;
00581       }
00582       else
00583       {
00584           /* swap with the tail */
00585           ptail-- ;
00586           Li [phead] = Li [ptail] ;
00587           Li [ptail] = i ;
00588           x = Lx [phead] ;
00589           Lx [phead] = Lx [ptail] ;
00590           Lx [ptail] = x ;
00591       }
00592         }
00593 
00594         /* set Lpend to one past the last entry in the
00595          * first part of the column of L.  Entries in
00596          * Li [0 ... Lpend [j]-1] are the only part of
00597          * column j of L that needs to be scanned in the DFS.
00598          * Lpend [j] was EMPTY; setting it >= 0 also flags
00599          * column j as pruned. */
00600         Lpend [j] = ptail ;
00601 
00602 #ifndef NDEBUG
00603         {
00604       Int p3 ;
00605       for (p3 = 0 ; p3 < Llen [j] ; p3++)
00606       {
00607           if (p3 == Lpend [j]) PRINTF (("----\n")) ;
00608           PRINTF (("after: %i  pivotal: %d\n", Li [p3],
00609           Pinv [Li [p3]] >= 0)) ;
00610       }
00611         }
00612 #endif
00613 
00614         break ;
00615     }
00616       }
00617   }
00618     }
00619 }
00620 
00621 
00622 /* ========================================================================== */
00623 /* === KLU_kernel =========================================================== */
00624 /* ========================================================================== */
00625 
00626 size_t KLU_kernel   /* final size of LU on output */
00627 (
00628     /* input, not modified */
00629     Int n,      /* A is n-by-n */
00630     Int Ap [ ],     /* size n+1, column pointers for A */
00631     Int Ai [ ],     /* size nz = Ap [n], row indices for A */
00632     Entry Ax [ ],   /* size nz, values of A */
00633     Int Q [ ],      /* size n, optional input permutation */
00634     size_t lusize,  /* initial size of LU on input */
00635 
00636     /* output, not defined on input */
00637     Int Pinv [ ],   /* size n, inverse row permutation, where Pinv [i] = k if
00638          * row i is the kth pivot row */
00639     Int P [ ],      /* size n, row permutation, where P [k] = i if row i is the
00640          * kth pivot row. */
00641     Unit **p_LU,  /* LU array, size lusize on input */
00642     Entry Udiag [ ],  /* size n, diagonal of U */
00643     Int Llen [ ],       /* size n, column length of L */
00644     Int Ulen [ ], /* size n, column length of U */
00645     Int Lip [ ],  /* size n, column pointers for L */
00646     Int Uip [ ],  /* size n, column pointers for U */
00647     Int *lnz,   /* size of L*/
00648     Int *unz,   /* size of U*/
00649     /* workspace, not defined on input */
00650     Entry X [ ],    /* size n, undefined on input, zero on output */
00651 
00652     /* workspace, not defined on input or output */
00653     Int Stack [ ],  /* size n */
00654     Int Flag [ ],   /* size n */
00655     Int Ap_pos [ ], /* size n */
00656 
00657     /* other workspace: */
00658     Int Lpend [ ],        /* size n workspace, for pruning only */
00659 
00660     /* inputs, not modified on output */
00661     Int k1,       /* the block of A is from k1 to k2-1 */
00662     Int PSinv [ ],    /* inverse of P from symbolic factorization */
00663     double Rs [ ],    /* scale factors for A */
00664 
00665     /* inputs, modified on output */
00666     Int Offp [ ],   /* off-diagonal matrix (modified by this routine) */
00667     Int Offi [ ],
00668     Entry Offx [ ],
00669     /* --------------- */
00670     KLU_common *Common
00671 )
00672 {
00673     Entry pivot ;
00674     double abs_pivot, xsize, nunits, tol, memgrow ;
00675     Entry *Ux ;
00676     Int *Li, *Ui ;
00677     Unit *LU ;    /* LU factors (pattern and values) */
00678     Int k, p, i, j, pivrow, kbar, diagrow, firstrow, lup, top, scale, len ;
00679     size_t newlusize ;
00680 
00681 #ifndef NDEBUG
00682     Entry *Lx ;
00683 #endif
00684 
00685     ASSERT (Common != NULL) ;
00686     scale = Common->scale ;
00687     tol = Common->tol ;
00688     memgrow = Common->memgrow ;
00689     *lnz = 0 ;
00690     *unz = 0 ;
00691 
00692     /* ---------------------------------------------------------------------- */
00693     /* get initial Li, Lx, Ui, and Ux */
00694     /* ---------------------------------------------------------------------- */
00695 
00696     PRINTF (("input: lusize %d \n", lusize)) ;
00697     ASSERT (lusize > 0) ;
00698     LU = *p_LU ;
00699 
00700     /* ---------------------------------------------------------------------- */
00701     /* initializations */
00702     /* ---------------------------------------------------------------------- */
00703 
00704     firstrow = 0 ;
00705     lup = 0 ;
00706 
00707     for (k = 0 ; k < n ; k++)
00708     {
00709   /* X [k] = 0 ; */
00710   CLEAR (X [k]) ;
00711   Flag [k] = EMPTY ;
00712   Lpend [k] = EMPTY ; /* flag k as not pruned */
00713     }
00714 
00715     /* ---------------------------------------------------------------------- */
00716     /* mark all rows as non-pivotal and determine initial diagonal mapping */
00717     /* ---------------------------------------------------------------------- */
00718 
00719     /* PSinv does the symmetric permutation, so don't do it here */
00720     for (k = 0 ; k < n ; k++)
00721     {
00722   P [k] = k ;
00723   Pinv [k] = FLIP (k) ; /* mark all rows as non-pivotal */
00724     }
00725     /* initialize the construction of the off-diagonal matrix */
00726     Offp [0] = 0 ;
00727 
00728     /* P [k] = row means that UNFLIP (Pinv [row]) = k, and visa versa.
00729      * If row is pivotal, then Pinv [row] >= 0.  A row is initially "flipped"
00730      * (Pinv [k] < EMPTY), and then marked "unflipped" when it becomes
00731      * pivotal. */
00732 
00733 #ifndef NDEBUG
00734     for (k = 0 ; k < n ; k++)
00735     {
00736   PRINTF (("Initial P [%d] = %d\n", k, P [k])) ;
00737     }
00738 #endif
00739 
00740     /* ---------------------------------------------------------------------- */
00741     /* factorize */
00742     /* ---------------------------------------------------------------------- */
00743 
00744     for (k = 0 ; k < n ; k++)
00745     {
00746 
00747   PRINTF (("\n\n==================================== k: %d\n", k)) ;
00748 
00749   /* ------------------------------------------------------------------ */
00750   /* determine if LU factors have grown too big */
00751   /* ------------------------------------------------------------------ */
00752 
00753   /* (n - k) entries for L and k entries for U */
00754   nunits = DUNITS (Int, n - k) + DUNITS (Int, k) +
00755      DUNITS (Entry, n - k) + DUNITS (Entry, k) ;
00756 
00757         /* LU can grow by at most 'nunits' entries if the column is dense */
00758         PRINTF (("lup %d lusize %g lup+nunits: %g\n", lup, (double) lusize,
00759       lup+nunits));
00760   xsize = ((double) lup) + nunits ;
00761   if (xsize > (double) lusize)
00762         {
00763             /* check here how much to grow */
00764       xsize = (memgrow * ((double) lusize) + 4*n + 1) ;
00765             if (INT_OVERFLOW (xsize))
00766             {
00767                 PRINTF (("Matrix is too large (Int overflow)\n")) ;
00768     Common->status = KLU_TOO_LARGE ;
00769                 return (lusize) ;
00770             }
00771             newlusize = (size_t) ( memgrow * lusize + 2*n + 1 ) ;
00772       /* Future work: retry mechanism in case of malloc failure */
00773       LU = KLU_realloc (newlusize, lusize, sizeof (Unit), LU, Common) ;
00774       Common->nrealloc++ ;
00775             *p_LU = LU ;
00776             if (Common->status == KLU_OUT_OF_MEMORY)
00777             {
00778                 PRINTF (("Matrix is too large (LU)\n")) ;
00779                 return (lusize) ;
00780             }
00781       lusize = newlusize ;
00782             PRINTF (("inc LU to %d done\n", lusize)) ;
00783         }
00784 
00785   /* ------------------------------------------------------------------ */
00786   /* start the kth column of L and U */
00787   /* ------------------------------------------------------------------ */
00788 
00789   Lip [k] = lup ;
00790 
00791   /* ------------------------------------------------------------------ */
00792   /* compute the nonzero pattern of the kth column of L and U */
00793   /* ------------------------------------------------------------------ */
00794 
00795 #ifndef NDEBUG
00796   for (i = 0 ; i < n ; i++)
00797   {
00798       ASSERT (Flag [i] < k) ;
00799       /* ASSERT (X [i] == 0) ; */
00800       ASSERT (IS_ZERO (X [i])) ;
00801   }
00802 #endif
00803 
00804   top = lsolve_symbolic (n, k, Ap, Ai, Q, Pinv, Stack, Flag,
00805         Lpend, Ap_pos, LU, lup, Llen, Lip, k1, PSinv) ;
00806 
00807 #ifndef NDEBUG
00808   PRINTF (("--- in U:\n")) ;
00809   for (p = top ; p < n ; p++)
00810   {
00811       PRINTF (("pattern of X for U: %d : %d pivot row: %d\n",
00812     p, Stack [p], Pinv [Stack [p]])) ;
00813       ASSERT (Flag [Stack [p]] == k) ;
00814   }
00815   PRINTF (("--- in L:\n")) ;
00816   Li = (Int *) (LU + Lip [k]);
00817   for (p = 0 ; p < Llen [k] ; p++)
00818   {
00819       PRINTF (("pattern of X in L: %d : %d pivot row: %d\n",
00820     p, Li [p], Pinv [Li [p]])) ;
00821       ASSERT (Flag [Li [p]] == k) ;
00822   }
00823   p = 0 ;
00824   for (i = 0 ; i < n ; i++)
00825   {
00826       ASSERT (Flag [i] <= k) ;
00827       if (Flag [i] == k) p++ ;
00828   }
00829 #endif
00830 
00831   /* ------------------------------------------------------------------ */
00832   /* get the column of the matrix to factorize and scatter into X */
00833   /* ------------------------------------------------------------------ */
00834 
00835   construct_column (k, Ap, Ai, Ax, Q, X,
00836       k1, PSinv, Rs, scale, Offp, Offi, Offx) ;
00837 
00838   /* ------------------------------------------------------------------ */
00839   /* compute the numerical values of the kth column (s = L \ A (:,k)) */
00840   /* ------------------------------------------------------------------ */
00841 
00842   lsolve_numeric (Pinv, LU, Stack, Lip, top, n, Llen, X) ;
00843 
00844 #ifndef NDEBUG
00845   for (p = top ; p < n ; p++)
00846   {
00847       PRINTF (("X for U %d : ",  Stack [p])) ;
00848       PRINT_ENTRY (X [Stack [p]]) ;
00849   }
00850   Li = (Int *) (LU + Lip [k]) ;
00851   for (p = 0 ; p < Llen [k] ; p++)
00852   {
00853       PRINTF (("X for L %d : ", Li [p])) ;
00854       PRINT_ENTRY (X [Li [p]]) ;
00855   }
00856 #endif
00857 
00858   /* ------------------------------------------------------------------ */
00859   /* partial pivoting with diagonal preference */
00860   /* ------------------------------------------------------------------ */
00861 
00862   /* determine what the "diagonal" is */
00863   diagrow = P [k] ;   /* might already be pivotal */
00864   PRINTF (("k %d, diagrow = %d, UNFLIP (diagrow) = %d\n",
00865       k, diagrow, UNFLIP (diagrow))) ;
00866 
00867   /* find a pivot and scale the pivot column */
00868   if (!lpivot (diagrow, &pivrow, &pivot, &abs_pivot, tol, X, LU, Lip,
00869         Llen, k, n, Pinv, &firstrow, Common))
00870   {
00871       /* matrix is structurally or numerically singular */
00872       Common->status = KLU_SINGULAR ;
00873       if (Common->numerical_rank == EMPTY)
00874       {
00875     Common->numerical_rank = k+k1 ;
00876     Common->singular_col = Q [k+k1] ;
00877       }
00878       if (Common->halt_if_singular)
00879       {
00880     /* do not continue the factorization */
00881     return (lusize) ;
00882       }
00883   }
00884 
00885   /* we now have a valid pivot row, even if the column has NaN's or
00886    * has no entries on or below the diagonal at all. */
00887   PRINTF (("\nk %d : Pivot row %d : ", k, pivrow)) ;
00888   PRINT_ENTRY (pivot) ;
00889   ASSERT (pivrow >= 0 && pivrow < n) ;
00890   ASSERT (Pinv [pivrow] < 0) ;
00891 
00892   /* set the Uip pointer */
00893   Uip [k] = Lip [k] + UNITS (Int, Llen [k]) + UNITS (Entry, Llen [k]) ;
00894 
00895         /* move the lup pointer to the position where indices of U
00896          * should be stored */
00897         lup += UNITS (Int, Llen [k]) + UNITS (Entry, Llen [k]) ;
00898 
00899         Ulen [k] = n - top ;
00900 
00901         /* extract Stack [top..n-1] to Ui and the values to Ux and clear X */
00902   GET_POINTER (LU, Uip, Ulen, Ui, Ux, k, len) ;
00903         for (p = top, i = 0 ; p < n ; p++, i++)
00904         {
00905       j = Stack [p] ;
00906       Ui [i] = Pinv [j] ;
00907       Ux [i] = X [j] ;
00908       CLEAR (X [j]) ;
00909         }
00910 
00911         /* position the lu index at the starting point for next column */
00912         lup += UNITS (Int, Ulen [k]) + UNITS (Entry, Ulen [k]) ;
00913 
00914   /* U(k,k) = pivot */
00915   Udiag [k] = pivot ;
00916 
00917   /* ------------------------------------------------------------------ */
00918   /* log the pivot permutation */
00919   /* ------------------------------------------------------------------ */
00920 
00921   ASSERT (UNFLIP (Pinv [diagrow]) < n) ;
00922   ASSERT (P [UNFLIP (Pinv [diagrow])] == diagrow) ;
00923 
00924   if (pivrow != diagrow)
00925   {
00926       /* an off-diagonal pivot has been chosen */
00927       Common->noffdiag++ ;
00928       PRINTF ((">>>>>>>>>>>>>>>>> pivrow %d k %d off-diagonal\n",
00929       pivrow, k)) ;
00930       if (Pinv [diagrow] < 0)
00931       {
00932     /* the former diagonal row index, diagrow, has not yet been
00933      * chosen as a pivot row.  Log this diagrow as the "diagonal"
00934      * entry in the column kbar for which the chosen pivot row,
00935      * pivrow, was originally logged as the "diagonal" */
00936     kbar = FLIP (Pinv [pivrow]) ;
00937     P [kbar] = diagrow ;
00938     Pinv [diagrow] = FLIP (kbar) ;
00939       }
00940   }
00941   P [k] = pivrow ;
00942   Pinv [pivrow] = k ;
00943 
00944 #ifndef NDEBUG
00945   for (i = 0 ; i < n ; i++) { ASSERT (IS_ZERO (X [i])) ;}
00946   GET_POINTER (LU, Uip, Ulen, Ui, Ux, k, len) ;
00947   for (p = 0 ; p < len ; p++)
00948   {
00949       PRINTF (("Column %d of U: %d : ", k, Ui [p])) ;
00950       PRINT_ENTRY (Ux [p]) ;
00951   }
00952   GET_POINTER (LU, Lip, Llen, Li, Lx, k, len) ;
00953   for (p = 0 ; p < len ; p++)
00954   {
00955       PRINTF (("Column %d of L: %d : ", k, Li [p])) ;
00956       PRINT_ENTRY (Lx [p]) ;
00957   }
00958 #endif
00959 
00960   /* ------------------------------------------------------------------ */
00961   /* symmetric pruning */
00962   /* ------------------------------------------------------------------ */
00963 
00964   prune (Lpend, Pinv, k, pivrow, LU, Uip, Lip, Ulen, Llen) ;
00965 
00966   *lnz += Llen [k] + 1 ; /* 1 added to lnz for diagonal */
00967   *unz += Ulen [k] + 1 ; /* 1 added to unz for diagonal */
00968     }
00969 
00970     /* ---------------------------------------------------------------------- */
00971     /* finalize column pointers for L and U, and put L in the pivotal order */
00972     /* ---------------------------------------------------------------------- */
00973 
00974     for (p = 0 ; p < n ; p++)
00975     {
00976   Li = (Int *) (LU + Lip [p]) ;
00977   for (i = 0 ; i < Llen [p] ; i++)
00978   {
00979       Li [i] = Pinv [Li [i]] ;
00980   }
00981     }
00982 
00983 #ifndef NDEBUG
00984     for (i = 0 ; i < n ; i++)
00985     {
00986   PRINTF (("P [%d] = %d   Pinv [%d] = %d\n", i, P [i], i, Pinv [i])) ;
00987     }
00988     for (i = 0 ; i < n ; i++)
00989     {
00990   ASSERT (Pinv [i] >= 0 && Pinv [i] < n) ;
00991   ASSERT (P [i] >= 0 && P [i] < n) ;
00992   ASSERT (P [Pinv [i]] == i) ;
00993   ASSERT (IS_ZERO (X [i])) ;
00994     }
00995 #endif
00996 
00997     /* ---------------------------------------------------------------------- */
00998     /* shrink the LU factors to just the required size */
00999     /* ---------------------------------------------------------------------- */
01000 
01001     newlusize = lup ;
01002     ASSERT ((size_t) newlusize <= lusize) ;
01003 
01004     /* this cannot fail, since the block is descreasing in size */
01005     LU = KLU_realloc (newlusize, lusize, sizeof (Unit), LU, Common) ;
01006     *p_LU = LU ;
01007     return (newlusize) ;
01008 }
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Defines