Amesos Package Browser (Single Doxygen Collection) Development
amesos_klu_l_factor.c
Go to the documentation of this file.
00001 /* ========================================================================== */
00002 /* === KLU_factor =========================================================== */
00003 /* ========================================================================== */
00004 
00005 /* Factor the matrix, after ordering and analyzing it with KLU_analyze
00006  * or KLU_analyze_given.
00007  */
00008 
00009 /* This file should make the long int version of KLU */
00010 #define DLONG 1
00011 
00012 #include "amesos_klu_internal.h"
00013 
00014 /* ========================================================================== */
00015 /* === KLU_factor2 ========================================================== */
00016 /* ========================================================================== */
00017 
00018 static void factor2
00019 (
00020     /* inputs, not modified */
00021     Int Ap [ ],   /* size n+1, column pointers */
00022     Int Ai [ ],   /* size nz, row indices */
00023     Entry Ax [ ],
00024     KLU_symbolic *Symbolic,
00025 
00026     /* inputs, modified on output: */
00027     KLU_numeric *Numeric,
00028     KLU_common *Common
00029 )
00030 {
00031     double lsize ;
00032     double *Lnz, *Rs ;
00033     Int *P, *Q, *R, *Pnum, *Offp, *Offi, *Pblock, *Pinv, *Iwork,
00034   *Lip, *Uip, *Llen, *Ulen ;
00035     Entry *Offx, *X, s, *Udiag ;
00036     Unit **LUbx ;
00037     Int k1, k2, nk, k, block, oldcol, pend, oldrow, n, lnz, unz, p, newrow,
00038   nblocks, poff, nzoff, lnz_block, unz_block, scale, max_lnz_block,
00039   max_unz_block ;
00040 
00041     /* ---------------------------------------------------------------------- */
00042     /* initializations */
00043     /* ---------------------------------------------------------------------- */
00044 
00045     /* get the contents of the Symbolic object */
00046     n = Symbolic->n ;
00047     P = Symbolic->P ;
00048     Q = Symbolic->Q ;
00049     R = Symbolic->R ;
00050     Lnz = Symbolic->Lnz ;
00051     nblocks = Symbolic->nblocks ;
00052     nzoff = Symbolic->nzoff ;
00053 
00054     Pnum = Numeric->Pnum ;
00055     Offp = Numeric->Offp ;
00056     Offi = Numeric->Offi ;
00057     Offx = (Entry *) Numeric->Offx ;
00058 
00059     Lip = Numeric->Lip ;
00060     Uip = Numeric->Uip ;
00061     Llen = Numeric->Llen ;
00062     Ulen = Numeric->Ulen ;
00063     LUbx = (Unit **) Numeric->LUbx ;
00064     Udiag = Numeric->Udiag ;
00065 
00066     Rs = Numeric->Rs ;
00067     Pinv = Numeric->Pinv ;
00068     X = (Entry *) Numeric->Xwork ;    /* X is of size n */
00069     Iwork = Numeric->Iwork ;      /* 5*maxblock for KLU_factor */
00070             /* 1*maxblock for Pblock */
00071     Pblock = Iwork + 5*((size_t) Symbolic->maxblock) ;
00072     Common->nrealloc = 0 ;
00073     scale = Common->scale ;
00074     max_lnz_block = 1 ;
00075     max_unz_block = 1 ;
00076 
00077     /* compute the inverse of P from symbolic analysis.  Will be updated to
00078      * become the inverse of the numerical factorization when the factorization
00079      * is done, for use in KLU_refactor */
00080 #ifndef NDEBUG
00081     for (k = 0 ; k < n ; k++)
00082     {
00083   Pinv [k] = EMPTY ;
00084     }
00085 #endif
00086     for (k = 0 ; k < n ; k++)
00087     {
00088   ASSERT (P [k] >= 0 && P [k] < n) ;
00089   Pinv [P [k]] = k ;
00090     }
00091 #ifndef NDEBUG
00092     for (k = 0 ; k < n ; k++) ASSERT (Pinv [k] != EMPTY) ;
00093 #endif
00094 
00095     lnz = 0 ;
00096     unz = 0 ;
00097     Common->noffdiag = 0 ;
00098     Offp [0] = 0 ;
00099 
00100     /* ---------------------------------------------------------------------- */
00101     /* optionally check input matrix and compute scale factors */
00102     /* ---------------------------------------------------------------------- */
00103 
00104     if (scale >= 0)
00105     {
00106   /* use Pnum as workspace. NOTE: scale factors are not yet permuted
00107    * according to the final pivot row ordering, so Rs [oldrow] is the
00108    * scale factor for A (oldrow,:), for the user's matrix A.  Pnum is
00109    * used as workspace in KLU_scale.  When the factorization is done,
00110    * the scale factors are permuted according to the final pivot row
00111    * permutation, so that Rs [k] is the scale factor for the kth row of
00112    * A(p,q) where p and q are the final row and column permutations. */
00113   KLU_scale (scale, n, Ap, Ai, (double *) Ax, Rs, Pnum, Common) ;
00114   if (Common->status < KLU_OK)
00115   {
00116       /* matrix is invalid */
00117       return ;
00118   }
00119     }
00120 
00121 #ifndef NDEBUG
00122     if (scale > 0)
00123     {
00124   for (k = 0 ; k < n ; k++) PRINTF (("Rs [%d] %g\n", k, Rs [k])) ;
00125     }
00126 #endif
00127 
00128     /* ---------------------------------------------------------------------- */
00129     /* factor each block using klu */
00130     /* ---------------------------------------------------------------------- */
00131 
00132     for (block = 0 ; block < nblocks ; block++)
00133     {
00134 
00135   /* ------------------------------------------------------------------ */
00136   /* the block is from rows/columns k1 to k2-1 */
00137   /* ------------------------------------------------------------------ */
00138 
00139   k1 = R [block] ;
00140   k2 = R [block+1] ;
00141   nk = k2 - k1 ;
00142   PRINTF (("FACTOR BLOCK %d, k1 %d k2-1 %d nk %d\n", block, k1,k2-1,nk)) ;
00143 
00144   if (nk == 1)
00145   {
00146 
00147       /* -------------------------------------------------------------- */
00148       /* singleton case */
00149       /* -------------------------------------------------------------- */
00150 
00151       poff = Offp [k1] ;
00152       oldcol = Q [k1] ;
00153       pend = Ap [oldcol+1] ;
00154       CLEAR (s) ;
00155 
00156       if (scale <= 0)
00157       {
00158     /* no scaling */
00159     for (p = Ap [oldcol] ; p < pend ; p++)
00160     {
00161         oldrow = Ai [p] ;
00162         newrow = Pinv [oldrow] ;
00163         if (newrow < k1)
00164         {
00165       Offi [poff] = oldrow ;
00166       Offx [poff] = Ax [p] ;
00167       poff++ ;
00168         }
00169         else
00170         {
00171       ASSERT (newrow == k1) ;
00172       PRINTF (("singleton block %d", block)) ;
00173       PRINT_ENTRY (Ax [p]) ;
00174       s = Ax [p] ;
00175         }
00176     }
00177       }
00178       else
00179       {
00180     /* row scaling.  NOTE: scale factors are not yet permuted
00181      * according to the pivot row permutation, so Rs [oldrow] is
00182      * used below.  When the factorization is done, the scale
00183      * factors are permuted, so that Rs [newrow] will be used in
00184      * klu_solve, klu_tsolve, and klu_rgrowth */
00185     for (p = Ap [oldcol] ; p < pend ; p++)
00186     {
00187         oldrow = Ai [p] ;
00188         newrow = Pinv [oldrow] ;
00189         if (newrow < k1)
00190         {
00191       Offi [poff] = oldrow ;
00192       /* Offx [poff] = Ax [p] / Rs [oldrow] ; */
00193       SCALE_DIV_ASSIGN (Offx [poff], Ax [p], Rs [oldrow]) ;
00194       poff++ ;
00195         }
00196         else
00197         {
00198       ASSERT (newrow == k1) ;
00199       PRINTF (("singleton block %d ", block)) ;
00200       PRINT_ENTRY (Ax[p]) ;
00201       SCALE_DIV_ASSIGN (s, Ax [p], Rs [oldrow]) ;
00202         }
00203     }
00204       }
00205 
00206       Udiag [k1] = s ;
00207 
00208       if (IS_ZERO (s))
00209       {
00210     /* singular singleton */
00211     Common->status = KLU_SINGULAR ;
00212     Common->numerical_rank = k1 ;
00213     Common->singular_col = oldcol ;
00214     if (Common->halt_if_singular)
00215     {
00216         return ;
00217     }
00218       }
00219 
00220       Offp [k1+1] = poff ;
00221       Pnum [k1] = P [k1] ;
00222       lnz++ ;
00223       unz++ ;
00224 
00225   }
00226   else
00227   {
00228 
00229       /* -------------------------------------------------------------- */
00230       /* construct and factorize the kth block */
00231       /* -------------------------------------------------------------- */
00232 
00233       if (Lnz [block] < 0)
00234       {
00235     /* COLAMD was used - no estimate of fill-in */
00236     /* use 10 times the nnz in A, plus n */
00237     lsize = -(Common->initmem) ;
00238       }
00239       else
00240       {
00241     lsize = Common->initmem_amd * Lnz [block] + nk ;
00242       }
00243 
00244       /* allocates 1 arrays: LUbx [block] */
00245       Numeric->LUsize [block] = KLU_kernel_factor (nk, Ap, Ai, Ax, Q,
00246         lsize, &LUbx [block], Udiag + k1, Llen + k1, Ulen + k1,
00247         Lip + k1, Uip + k1, Pblock, &lnz_block, &unz_block,
00248         X, Iwork, k1, Pinv, Rs, Offp, Offi, Offx, Common) ;
00249 
00250       if (Common->status < KLU_OK ||
00251          (Common->status == KLU_SINGULAR && Common->halt_if_singular))
00252       {
00253     /* out of memory, invalid inputs, or singular */
00254     return ;
00255       }
00256 
00257       PRINTF (("\n----------------------- L %d:\n", block)) ;
00258       ASSERT (KLU_valid_LU (nk, TRUE, Lip+k1, Llen+k1, LUbx [block])) ;
00259       PRINTF (("\n----------------------- U %d:\n", block)) ;
00260       ASSERT (KLU_valid_LU (nk, FALSE, Uip+k1, Ulen+k1, LUbx [block])) ;
00261 
00262       /* -------------------------------------------------------------- */
00263       /* get statistics */
00264       /* -------------------------------------------------------------- */
00265 
00266       lnz += lnz_block ;
00267       unz += unz_block ;
00268       max_lnz_block = MAX (max_lnz_block, lnz_block) ;
00269       max_unz_block = MAX (max_unz_block, unz_block) ;
00270 
00271       if (Lnz [block] == EMPTY)
00272       {
00273     /* revise estimate for subsequent factorization */
00274     Lnz [block] = MAX (lnz_block, unz_block) ;
00275       }
00276 
00277       /* -------------------------------------------------------------- */
00278       /* combine the klu row ordering with the symbolic pre-ordering */
00279       /* -------------------------------------------------------------- */
00280 
00281       PRINTF (("Pnum, 1-based:\n")) ;
00282       for (k = 0 ; k < nk ; k++)
00283       {
00284     ASSERT (k + k1 < n) ;
00285     ASSERT (Pblock [k] + k1 < n) ;
00286     Pnum [k + k1] = P [Pblock [k] + k1] ;
00287     PRINTF (("Pnum (%d + %d + 1 = %d) = %d + 1 = %d\n",
00288         k, k1, k+k1+1, Pnum [k+k1], Pnum [k+k1]+1)) ;
00289       }
00290 
00291       /* the local pivot row permutation Pblock is no longer needed */
00292   }
00293     }
00294     ASSERT (nzoff == Offp [n]) ;
00295     PRINTF (("\n------------------- Off diagonal entries:\n")) ;
00296     ASSERT (KLU_valid (n, Offp, Offi, Offx)) ;
00297 
00298     Numeric->lnz = lnz ;
00299     Numeric->unz = unz ;
00300     Numeric->max_lnz_block = max_lnz_block ;
00301     Numeric->max_unz_block = max_unz_block ;
00302 
00303     /* compute the inverse of Pnum */
00304 #ifndef NDEBUG
00305     for (k = 0 ; k < n ; k++)
00306     {
00307   Pinv [k] = EMPTY ;
00308     }
00309 #endif
00310     for (k = 0 ; k < n ; k++)
00311     {
00312   ASSERT (Pnum [k] >= 0 && Pnum [k] < n) ;
00313   Pinv [Pnum [k]] = k ;
00314     }
00315 #ifndef NDEBUG
00316     for (k = 0 ; k < n ; k++) ASSERT (Pinv [k] != EMPTY) ;
00317 #endif
00318 
00319     /* permute scale factors Rs according to pivotal row order */
00320     if (scale > 0)
00321     {
00322   for (k = 0 ; k < n ; k++)
00323   {
00324       REAL (X [k]) = Rs [Pnum [k]] ;
00325   }
00326   for (k = 0 ; k < n ; k++)
00327   {
00328       Rs [k] = REAL (X [k]) ;
00329   }
00330     }
00331 
00332     PRINTF (("\n------------------- Off diagonal entries, old:\n")) ;
00333     ASSERT (KLU_valid (n, Offp, Offi, Offx)) ;
00334 
00335     /* apply the pivot row permutations to the off-diagonal entries */
00336     for (p = 0 ; p < nzoff ; p++)
00337     {
00338   ASSERT (Offi [p] >= 0 && Offi [p] < n) ;
00339   Offi [p] = Pinv [Offi [p]] ;
00340     }
00341 
00342     PRINTF (("\n------------------- Off diagonal entries, new:\n")) ;
00343     ASSERT (KLU_valid (n, Offp, Offi, Offx)) ;
00344 
00345 #ifndef NDEBUG
00346     {
00347   PRINTF (("\n ############# KLU_BTF_FACTOR done, nblocks %d\n",nblocks));
00348   Entry ss, *Udiag = Numeric->Udiag ;
00349   for (block = 0 ; block < nblocks && Common->status == KLU_OK ; block++)
00350   {
00351       k1 = R [block] ;
00352       k2 = R [block+1] ;
00353       nk = k2 - k1 ;
00354       PRINTF (("\n======================KLU_factor output: k1 %d k2 %d nk %d\n",k1,k2,nk)) ;
00355       if (nk == 1)
00356       {
00357     PRINTF (("singleton  ")) ;
00358     /* ENTRY_PRINT (singleton [block]) ; */
00359     ss = Udiag [k1] ;
00360     PRINT_ENTRY (ss) ;
00361       }
00362       else
00363       {
00364     Int *Lip, *Uip, *Llen, *Ulen ;
00365     Unit *LU ;
00366     Lip = Numeric->Lip + k1 ;
00367     Llen = Numeric->Llen + k1 ;
00368     LU = (Unit *) Numeric->LUbx [block] ;
00369     PRINTF (("\n---- L block %d\n", block));
00370     ASSERT (KLU_valid_LU (nk, TRUE, Lip, Llen, LU)) ;
00371     Uip = Numeric->Uip + k1 ;
00372     Ulen = Numeric->Ulen + k1 ;
00373     PRINTF (("\n---- U block %d\n", block)) ;
00374     ASSERT (KLU_valid_LU (nk, FALSE, Uip, Ulen, LU)) ;
00375       }
00376   }
00377     }
00378 #endif
00379 }
00380 
00381 
00382 
00383 /* ========================================================================== */
00384 /* === KLU_factor =========================================================== */
00385 /* ========================================================================== */
00386 
00387 KLU_numeric *KLU_factor   /* returns NULL if error, or a valid
00388            KLU_numeric object if successful */
00389 (
00390     /* --- inputs --- */
00391     Int Ap [ ],   /* size n+1, column pointers */
00392     Int Ai [ ],   /* size nz, row indices */
00393     double Ax [ ],
00394     KLU_symbolic *Symbolic,
00395     /* -------------- */
00396     KLU_common *Common
00397 )
00398 {
00399     Int n, nzoff, nblocks, maxblock, k, ok = TRUE ;
00400     Int *R ;
00401     KLU_numeric *Numeric ;
00402     size_t n1, nzoff1, s, b6, n3 ;
00403 
00404     if (Common == NULL)
00405     {
00406   return (NULL) ;
00407     }
00408     Common->status = KLU_OK ;
00409     Common->numerical_rank = EMPTY ;
00410     Common->singular_col = EMPTY ;
00411 
00412     /* ---------------------------------------------------------------------- */
00413     /* get the contents of the Symbolic object */
00414     /* ---------------------------------------------------------------------- */
00415 
00416     /* check for a valid Symbolic object */
00417     if (Symbolic == NULL)
00418     {
00419   Common->status = KLU_INVALID ;
00420   return (NULL) ;
00421     }
00422 
00423     n = Symbolic->n ;
00424     nzoff = Symbolic->nzoff ;
00425     nblocks = Symbolic->nblocks ;
00426     maxblock = Symbolic->maxblock ;
00427     R = Symbolic->R ;
00428     PRINTF (("KLU_factor:  n %d nzoff %d nblocks %d maxblock %d\n",
00429   n, nzoff, nblocks, maxblock)) ;
00430 
00431     /* ---------------------------------------------------------------------- */
00432     /* get control parameters and make sure they are in the proper range */
00433     /* ---------------------------------------------------------------------- */
00434 
00435     Common->initmem_amd = MAX (1.0, Common->initmem_amd) ;
00436     Common->initmem = MAX (1.0, Common->initmem) ;
00437     Common->tol = MIN (Common->tol, 1.0) ;
00438     Common->tol = MAX (0.0, Common->tol) ;
00439     Common->memgrow = MAX (1.0, Common->memgrow) ;
00440 
00441     /* ---------------------------------------------------------------------- */
00442     /* allocate the Numeric object  */
00443     /* ---------------------------------------------------------------------- */
00444 
00445     /* this will not cause size_t overflow (already checked by KLU_symbolic) */
00446     n1 = ((size_t) n) + 1 ;
00447     nzoff1 = ((size_t) nzoff) + 1 ;
00448 
00449     Numeric = KLU_malloc (sizeof (KLU_numeric), 1, Common) ;
00450     if (Common->status < KLU_OK)
00451     {
00452   /* out of memory */
00453   Common->status = KLU_OUT_OF_MEMORY ;
00454   return (NULL) ;
00455     }
00456     Numeric->n = n ;
00457     Numeric->nblocks = nblocks ;
00458     Numeric->nzoff = nzoff ;
00459     Numeric->Pnum = KLU_malloc (n, sizeof (Int), Common) ;
00460     Numeric->Offp = KLU_malloc (n1, sizeof (Int), Common) ;
00461     Numeric->Offi = KLU_malloc (nzoff1, sizeof (Int), Common) ;
00462     Numeric->Offx = KLU_malloc (nzoff1, sizeof (Entry), Common) ;
00463 
00464     Numeric->Lip  = KLU_malloc (n, sizeof (Int), Common) ;
00465     Numeric->Uip  = KLU_malloc (n, sizeof (Int), Common) ;
00466     Numeric->Llen = KLU_malloc (n, sizeof (Int), Common) ;
00467     Numeric->Ulen = KLU_malloc (n, sizeof (Int), Common) ;
00468 
00469     Numeric->LUsize = KLU_malloc (nblocks, sizeof (size_t), Common) ;
00470 
00471     Numeric->LUbx = KLU_malloc (nblocks, sizeof (Unit *), Common) ;
00472     if (Numeric->LUbx != NULL)
00473     {
00474   for (k = 0 ; k < nblocks ; k++)
00475   {
00476       Numeric->LUbx [k] = NULL ;
00477   }
00478     }
00479 
00480     Numeric->Udiag = KLU_malloc (n, sizeof (Entry), Common) ;
00481 
00482     if (Common->scale > 0)
00483     {
00484   Numeric->Rs = KLU_malloc (n, sizeof (double), Common) ;
00485     }
00486     else
00487     {
00488   /* no scaling */
00489   Numeric->Rs = NULL ;
00490     }
00491 
00492     Numeric->Pinv = KLU_malloc (n, sizeof (Int), Common) ;
00493 
00494     /* allocate permanent workspace for factorization and solve.  Note that the
00495      * solver will use an Xwork of size 4n, whereas the factorization codes use
00496      * an Xwork of size n and integer space (Iwork) of size 6n. KLU_condest
00497      * uses an Xwork of size 2n.  Total size is:
00498      *
00499      *    n*sizeof(Entry) + max (6*maxblock*sizeof(Int), 3*n*sizeof(Entry))
00500      */
00501     s = KLU_mult_size_t (n, sizeof (Entry), &ok) ;
00502     n3 = KLU_mult_size_t (n, 3 * sizeof (Entry), &ok) ;
00503     b6 = KLU_mult_size_t (maxblock, 6 * sizeof (Int), &ok) ;
00504     Numeric->worksize = KLU_add_size_t (s, MAX (n3, b6), &ok) ;
00505     Numeric->Work = KLU_malloc (Numeric->worksize, 1, Common) ;
00506     Numeric->Xwork = Numeric->Work ;
00507     Numeric->Iwork = (Int *) ((Entry *) Numeric->Xwork + n) ;
00508     if (!ok || Common->status < KLU_OK)
00509     {
00510   /* out of memory or problem too large */
00511   Common->status = ok ? KLU_OUT_OF_MEMORY : KLU_TOO_LARGE ;
00512   KLU_free_numeric (&Numeric, Common) ;
00513   return (NULL) ;
00514     }
00515 
00516     /* ---------------------------------------------------------------------- */
00517     /* factorize the blocks */
00518     /* ---------------------------------------------------------------------- */
00519 
00520     factor2 (Ap, Ai, (Entry *) Ax, Symbolic, Numeric, Common) ;
00521 
00522     /* ---------------------------------------------------------------------- */
00523     /* return or free the Numeric object */
00524     /* ---------------------------------------------------------------------- */
00525 
00526     if (Common->status < KLU_OK)
00527     {
00528   /* out of memory or inputs invalid */
00529   KLU_free_numeric (&Numeric, Common) ;
00530     }
00531     else if (Common->status == KLU_SINGULAR)
00532     {
00533   if (Common->halt_if_singular)
00534   {
00535       /* Matrix is singular, and the Numeric object is only partially
00536        * defined because we halted early.  This is the default case for
00537        * a singular matrix. */
00538       KLU_free_numeric (&Numeric, Common) ;
00539   }
00540     }
00541     else if (Common->status == KLU_OK)
00542     {
00543   /* successful non-singular factorization */
00544   Common->numerical_rank = n ;
00545   Common->singular_col = n ;
00546     }
00547     return (Numeric) ;
00548 }
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Defines