Amesos Package Browser (Single Doxygen Collection) Development
amesos_klu_refactor.c
Go to the documentation of this file.
00001 /* ========================================================================== */
00002 /* === KLU_refactor ========================================================= */
00003 /* ========================================================================== */
00004 
00005 /* Factor the matrix, after ordering and analyzing it with KLU_analyze, and
00006  * factoring it once with KLU_factor.  This routine cannot do any numerical
00007  * pivoting.  The pattern of the input matrix (Ap, Ai) must be identical to
00008  * the pattern given to KLU_factor.
00009  */
00010 
00011 #include "amesos_klu_internal.h"
00012 
00013 
00014 /* ========================================================================== */
00015 /* === KLU_refactor ========================================================= */
00016 /* ========================================================================== */
00017 
00018 Int KLU_refactor  /* returns TRUE if successful, FALSE otherwise */
00019 (
00020     /* inputs, not modified */
00021     Int Ap [ ],   /* size n+1, column pointers */
00022     Int Ai [ ],   /* size nz, row indices */
00023     double Ax [ ],
00024     KLU_symbolic *Symbolic,
00025 
00026     /* input/output */
00027     KLU_numeric *Numeric,
00028     KLU_common  *Common
00029 )
00030 {
00031     Entry ukk, ujk, s ;
00032     Entry *Offx, *Lx, *Ux, *X, *Az, *Udiag ;
00033     double *Rs ;
00034     Int *P, *Q, *R, *Pnum, *Offp, *Offi, *Ui, *Li, *Pinv, *Lip, *Uip, *Llen,
00035   *Ulen ;
00036     Unit **LUbx ;
00037     Unit *LU ;
00038     Int k1, k2, nk, k, block, oldcol, pend, oldrow, n, p, newrow, scale,
00039   nblocks, poff, i, j, up, ulen, llen, maxblock, nzoff ;
00040 
00041     /* ---------------------------------------------------------------------- */
00042     /* check inputs */
00043     /* ---------------------------------------------------------------------- */
00044 
00045     if (Common == NULL)
00046     {
00047   return (FALSE) ;
00048     }
00049     Common->status = KLU_OK ;
00050 
00051     if (Numeric == NULL)
00052     {
00053   /* invalid Numeric object */
00054   Common->status = KLU_INVALID ;
00055   return (FALSE) ;
00056     }
00057 
00058     Common->numerical_rank = EMPTY ;
00059     Common->singular_col = EMPTY ;
00060 
00061     Az = (Entry *) Ax ;
00062 
00063     /* ---------------------------------------------------------------------- */
00064     /* get the contents of the Symbolic object */
00065     /* ---------------------------------------------------------------------- */
00066 
00067     n = Symbolic->n ;
00068     P = Symbolic->P ;
00069     Q = Symbolic->Q ;
00070     R = Symbolic->R ;
00071     nblocks = Symbolic->nblocks ;
00072     maxblock = Symbolic->maxblock ;
00073 
00074     /* ---------------------------------------------------------------------- */
00075     /* get the contents of the Numeric object */
00076     /* ---------------------------------------------------------------------- */
00077 
00078     Pnum = Numeric->Pnum ;
00079     Offp = Numeric->Offp ;
00080     Offi = Numeric->Offi ;
00081     Offx = (Entry *) Numeric->Offx ;
00082 
00083     LUbx = (Unit **) Numeric->LUbx ;
00084 
00085     scale = Common->scale ;
00086     if (scale > 0)
00087     {
00088   /* factorization was not scaled, but refactorization is scaled */
00089   if (Numeric->Rs == NULL)
00090   {
00091       Numeric->Rs = KLU_malloc (n, sizeof (double), Common) ;
00092       if (Common->status < KLU_OK)
00093       {
00094     Common->status = KLU_OUT_OF_MEMORY ;
00095     return (FALSE) ;
00096       }
00097   }
00098     }
00099     else
00100     {
00101   /* no scaling for refactorization; ensure Numeric->Rs is freed.  This
00102    * does nothing if Numeric->Rs is already NULL. */
00103   Numeric->Rs = KLU_free (Numeric->Rs, n, sizeof (double), Common) ;
00104     }
00105     Rs = Numeric->Rs ;
00106 
00107     Pinv = Numeric->Pinv ;
00108     X = (Entry *) Numeric->Xwork ;
00109     Common->nrealloc = 0 ;
00110     Udiag = Numeric->Udiag ;
00111     nzoff = Symbolic->nzoff ;
00112 
00113     /* ---------------------------------------------------------------------- */
00114     /* check the input matrix compute the row scale factors, Rs */
00115     /* ---------------------------------------------------------------------- */
00116 
00117     /* do no scale, or check the input matrix, if scale < 0 */
00118     if (scale >= 0)
00119     {
00120   /* check for out-of-range indices, but do not check for duplicates */
00121   if (!KLU_scale (scale, n, Ap, Ai, Ax, Rs, NULL, Common))
00122   {
00123       return (FALSE) ;
00124   }
00125     }
00126 
00127     /* ---------------------------------------------------------------------- */
00128     /* clear workspace X */
00129     /* ---------------------------------------------------------------------- */
00130 
00131     for (k = 0 ; k < maxblock ; k++)
00132     {
00133   /* X [k] = 0 */
00134   CLEAR (X [k]) ;
00135     }
00136 
00137     poff = 0 ;
00138 
00139     /* ---------------------------------------------------------------------- */
00140     /* factor each block */
00141     /* ---------------------------------------------------------------------- */
00142 
00143     if (scale <= 0)
00144     {
00145 
00146   /* ------------------------------------------------------------------ */
00147   /* no scaling */
00148   /* ------------------------------------------------------------------ */
00149 
00150   for (block = 0 ; block < nblocks ; block++)
00151   {
00152 
00153       /* -------------------------------------------------------------- */
00154       /* the block is from rows/columns k1 to k2-1 */
00155       /* -------------------------------------------------------------- */
00156 
00157       k1 = R [block] ;
00158       k2 = R [block+1] ;
00159       nk = k2 - k1 ;
00160 
00161       if (nk == 1)
00162       {
00163 
00164     /* ---------------------------------------------------------- */
00165     /* singleton case */
00166     /* ---------------------------------------------------------- */
00167 
00168     oldcol = Q [k1] ;
00169     pend = Ap [oldcol+1] ;
00170     CLEAR (s) ;
00171     for (p = Ap [oldcol] ; p < pend ; p++)
00172     {
00173         newrow = Pinv [Ai [p]] - k1 ;
00174         if (newrow < 0 && poff < nzoff)
00175         {
00176       /* entry in off-diagonal block */
00177       Offx [poff] = Az [p] ;
00178       poff++ ;
00179         }
00180         else
00181         {
00182       /* singleton */
00183       s = Az [p] ;
00184         }
00185     }
00186     Udiag [k1] = s ;
00187 
00188       }
00189       else
00190       {
00191 
00192     /* ---------------------------------------------------------- */
00193     /* construct and factor the kth block */
00194     /* ---------------------------------------------------------- */
00195 
00196     Lip  = Numeric->Lip  + k1 ;
00197     Llen = Numeric->Llen + k1 ;
00198     Uip  = Numeric->Uip  + k1 ;
00199     Ulen = Numeric->Ulen + k1 ;
00200     LU = LUbx [block] ;
00201 
00202     for (k = 0 ; k < nk ; k++)
00203     {
00204 
00205         /* ------------------------------------------------------ */
00206         /* scatter kth column of the block into workspace X */
00207         /* ------------------------------------------------------ */
00208 
00209         oldcol = Q [k+k1] ;
00210         pend = Ap [oldcol+1] ;
00211         for (p = Ap [oldcol] ; p < pend ; p++)
00212         {
00213       newrow = Pinv [Ai [p]] - k1 ;
00214       if (newrow < 0 && poff < nzoff)
00215       {
00216           /* entry in off-diagonal block */
00217           Offx [poff] = Az [p] ;
00218           poff++ ;
00219       }
00220       else
00221       {
00222           /* (newrow,k) is an entry in the block */
00223           X [newrow] = Az [p] ;
00224       }
00225         }
00226 
00227         /* ------------------------------------------------------ */
00228         /* compute kth column of U, and update kth column of A */
00229         /* ------------------------------------------------------ */
00230 
00231         GET_POINTER (LU, Uip, Ulen, Ui, Ux, k, ulen) ;
00232         for (up = 0 ; up < ulen ; up++)
00233         {
00234       j = Ui [up] ;
00235       ujk = X [j] ;
00236       /* X [j] = 0 */
00237       CLEAR (X [j]) ;
00238       Ux [up] = ujk ;
00239       GET_POINTER (LU, Lip, Llen, Li, Lx, j, llen) ;
00240       for (p = 0 ; p < llen ; p++)
00241       {
00242           /* X [Li [p]] -= Lx [p] * ujk */
00243           MULT_SUB (X [Li [p]], Lx [p], ujk) ;
00244       }
00245         }
00246         /* get the diagonal entry of U */
00247         ukk = X [k] ;
00248         /* X [k] = 0 */
00249         CLEAR (X [k]) ;
00250         /* singular case */
00251         if (IS_ZERO (ukk))
00252         {
00253       /* matrix is numerically singular */
00254       Common->status = KLU_SINGULAR ;
00255       if (Common->numerical_rank == EMPTY)
00256       {
00257           Common->numerical_rank = k+k1 ;
00258           Common->singular_col = Q [k+k1] ;
00259       }
00260       if (Common->halt_if_singular)
00261       {
00262           /* do not continue the factorization */
00263           return (FALSE) ;
00264       }
00265         }
00266         Udiag [k+k1] = ukk ;
00267         /* gather and divide by pivot to get kth column of L */
00268         GET_POINTER (LU, Lip, Llen, Li, Lx, k, llen) ;
00269         for (p = 0 ; p < llen ; p++)
00270         {
00271       i = Li [p] ;
00272       DIV (Lx [p], X [i], ukk) ;
00273       CLEAR (X [i]) ;
00274         }
00275 
00276     }
00277       }
00278   }
00279 
00280     }
00281     else
00282     {
00283 
00284   /* ------------------------------------------------------------------ */
00285   /* scaling */
00286   /* ------------------------------------------------------------------ */
00287 
00288   for (block = 0 ; block < nblocks ; block++)
00289   {
00290 
00291       /* -------------------------------------------------------------- */
00292       /* the block is from rows/columns k1 to k2-1 */
00293       /* -------------------------------------------------------------- */
00294 
00295       k1 = R [block] ;
00296       k2 = R [block+1] ;
00297       nk = k2 - k1 ;
00298 
00299       if (nk == 1)
00300       {
00301 
00302     /* ---------------------------------------------------------- */
00303     /* singleton case */
00304     /* ---------------------------------------------------------- */
00305 
00306     oldcol = Q [k1] ;
00307     pend = Ap [oldcol+1] ;
00308     CLEAR (s) ;
00309     for (p = Ap [oldcol] ; p < pend ; p++)
00310     {
00311         oldrow = Ai [p] ;
00312         newrow = Pinv [oldrow] - k1 ;
00313         if (newrow < 0 && poff < nzoff)
00314         {
00315       /* entry in off-diagonal block */
00316       /* Offx [poff] = Az [p] / Rs [oldrow] */
00317       SCALE_DIV_ASSIGN (Offx [poff], Az [p], Rs [oldrow]) ;
00318       poff++ ;
00319         }
00320         else
00321         {
00322       /* singleton */
00323       /* s = Az [p] / Rs [oldrow] */
00324       SCALE_DIV_ASSIGN (s, Az [p], Rs [oldrow]) ;
00325         }
00326     }
00327     Udiag [k1] = s ;
00328 
00329       }
00330       else
00331       {
00332 
00333     /* ---------------------------------------------------------- */
00334     /* construct and factor the kth block */
00335     /* ---------------------------------------------------------- */
00336 
00337     Lip  = Numeric->Lip  + k1 ;
00338     Llen = Numeric->Llen + k1 ;
00339     Uip  = Numeric->Uip  + k1 ;
00340     Ulen = Numeric->Ulen + k1 ;
00341     LU = LUbx [block] ;
00342 
00343     for (k = 0 ; k < nk ; k++)
00344     {
00345 
00346         /* ------------------------------------------------------ */
00347         /* scatter kth column of the block into workspace X */
00348         /* ------------------------------------------------------ */
00349 
00350         oldcol = Q [k+k1] ;
00351         pend = Ap [oldcol+1] ;
00352         for (p = Ap [oldcol] ; p < pend ; p++)
00353         {
00354       oldrow = Ai [p] ;
00355       newrow = Pinv [oldrow] - k1 ;
00356       if (newrow < 0 && poff < nzoff)
00357       {
00358           /* entry in off-diagonal part */
00359           /* Offx [poff] = Az [p] / Rs [oldrow] */
00360           SCALE_DIV_ASSIGN (Offx [poff], Az [p], Rs [oldrow]);
00361           poff++ ;
00362       }
00363       else
00364       {
00365           /* (newrow,k) is an entry in the block */
00366           /* X [newrow] = Az [p] / Rs [oldrow] */
00367           SCALE_DIV_ASSIGN (X [newrow], Az [p], Rs [oldrow]) ;
00368       }
00369         }
00370 
00371         /* ------------------------------------------------------ */
00372         /* compute kth column of U, and update kth column of A */
00373         /* ------------------------------------------------------ */
00374 
00375         GET_POINTER (LU, Uip, Ulen, Ui, Ux, k, ulen) ;
00376         for (up = 0 ; up < ulen ; up++)
00377         {
00378       j = Ui [up] ;
00379       ujk = X [j] ;
00380       /* X [j] = 0 */
00381       CLEAR (X [j]) ;
00382       Ux [up] = ujk ;
00383       GET_POINTER (LU, Lip, Llen, Li, Lx, j, llen) ;
00384       for (p = 0 ; p < llen ; p++)
00385       {
00386           /* X [Li [p]] -= Lx [p] * ujk */
00387           MULT_SUB (X [Li [p]], Lx [p], ujk) ;
00388       }
00389         }
00390         /* get the diagonal entry of U */
00391         ukk = X [k] ;
00392         /* X [k] = 0 */
00393         CLEAR (X [k]) ;
00394         /* singular case */
00395         if (IS_ZERO (ukk))
00396         {
00397       /* matrix is numerically singular */
00398       Common->status = KLU_SINGULAR ;
00399       if (Common->numerical_rank == EMPTY)
00400       {
00401           Common->numerical_rank = k+k1 ;
00402           Common->singular_col = Q [k+k1] ;
00403       }
00404       if (Common->halt_if_singular)
00405       {
00406           /* do not continue the factorization */
00407           return (FALSE) ;
00408       }
00409         }
00410         Udiag [k+k1] = ukk ;
00411         /* gather and divide by pivot to get kth column of L */
00412         GET_POINTER (LU, Lip, Llen, Li, Lx, k, llen) ;
00413         for (p = 0 ; p < llen ; p++)
00414         {
00415       i = Li [p] ;
00416       DIV (Lx [p], X [i], ukk) ;
00417       CLEAR (X [i]) ;
00418         }
00419     }
00420       }
00421   }
00422     }
00423 
00424     /* ---------------------------------------------------------------------- */
00425     /* permute scale factors Rs according to pivotal row order */
00426     /* ---------------------------------------------------------------------- */
00427 
00428     if (scale > 0)
00429     {
00430   for (k = 0 ; k < n ; k++)
00431   {
00432       REAL (X [k]) = Rs [Pnum [k]] ;
00433   }
00434   for (k = 0 ; k < n ; k++)
00435   {
00436       Rs [k] = REAL (X [k]) ;
00437   }
00438     }
00439 
00440 #ifndef NDEBUG
00441     ASSERT (Offp [n] == poff) ;
00442     ASSERT (Symbolic->nzoff == poff) ;
00443     PRINTF (("\n------------------- Off diagonal entries, new:\n")) ;
00444     ASSERT (KLU_valid (n, Offp, Offi, Offx)) ;
00445     if (Common->status == KLU_OK)
00446     {
00447   PRINTF (("\n ########### KLU_BTF_REFACTOR done, nblocks %d\n",nblocks));
00448   for (block = 0 ; block < nblocks ; block++)
00449   {
00450       k1 = R [block] ;
00451       k2 = R [block+1] ;
00452       nk = k2 - k1 ;
00453       PRINTF ((
00454     "\n================KLU_refactor output: k1 %d k2 %d nk %d\n",
00455     k1, k2, nk)) ;
00456       if (nk == 1)
00457       {
00458     PRINTF (("singleton  ")) ;
00459     PRINT_ENTRY (Udiag [k1]) ;
00460       }
00461       else
00462       {
00463     Lip = Numeric->Lip + k1 ;
00464                 Llen = Numeric->Llen + k1 ;
00465                 LU = (Unit *) Numeric->LUbx [block] ;
00466     PRINTF (("\n---- L block %d\n", block)) ;
00467     ASSERT (KLU_valid_LU (nk, TRUE, Lip, Llen, LU)) ;
00468     Uip = Numeric->Uip + k1 ;
00469                 Ulen = Numeric->Ulen + k1 ;
00470     PRINTF (("\n---- U block %d\n", block)) ;
00471     ASSERT (KLU_valid_LU (nk, FALSE, Uip, Ulen, LU)) ;
00472       }
00473   }
00474     }
00475 #endif
00476 
00477     return (TRUE) ;
00478 }
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Defines