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