Amesos Package Browser (Single Doxygen Collection) Development
amesos_cholmod_factorize.c
Go to the documentation of this file.
```00001 /* ========================================================================== */
00002 /* === Cholesky/cholmod_factorize =========================================== */
00003 /* ========================================================================== */
00004
00005 /* -----------------------------------------------------------------------------
00006  * CHOLMOD/Cholesky Module.  Copyright (C) 2005-2006, Timothy A. Davis
00008  * Lesser General Public License.  See lesser.txt for a text of the license.
00009  * CHOLMOD is also available under other licenses; contact authors for details.
00010  * http://www.cise.ufl.edu/research/sparse
00011  * -------------------------------------------------------------------------- */
00012
00013 /* Computes the numerical factorization of a symmetric matrix.  The primary
00014  * inputs to this routine are a sparse matrix A and the symbolic factor L from
00015  * cholmod_analyze or a prior numerical factor L.  If A is symmetric, this
00016  * routine factorizes A(p,p)+beta*I (beta can be zero), where p is the
00017  * fill-reducing permutation (L->Perm).  If A is unsymmetric, either
00018  * A(p,:)*A(p,:)'+beta*I or A(p,f)*A(p,f)'+beta*I is factorized.  The set f and
00019  * the nonzero pattern of the matrix A must be the same as the matrix passed to
00020  * cholmod_analyze for the supernodal case.  For the simplicial case, it can
00021  * be different, but it should be the same for best performance.  beta is real.
00022  *
00023  * A simplicial factorization or supernodal factorization is chosen, based on
00024  * the type of the factor L.  If L->is_super is TRUE, a supernodal LL'
00025  * factorization is computed.  Otherwise, a simplicial numeric factorization
00026  * is computed, either LL' or LDL', depending on Common->final_ll.
00027  *
00028  * Once the factorization is complete, it can be left as is or optionally
00029  * converted into any simplicial numeric type, depending on the
00030  * Common->final_* parameters.  If converted from a supernodal to simplicial
00031  * type, and the Common->final_resymbol parameter is true, then numerically
00032  * zero entries in L due to relaxed supernodal amalgamation are removed from
00033  * the simplicial factor (they are always left in the supernodal form of L).
00034  * Entries that are numerically zero but present in the simplicial symbolic
00035  * pattern of L are left in place (that is, the graph of L remains chordal).
00036  * This is required for the update/downdate/rowadd/rowdel routines to work
00037  * properly.
00038  *
00039  * workspace: Flag (nrow), Head (nrow+1),
00040  *  if symmetric:   Iwork (2*nrow+2*nsuper)
00041  *  if unsymmetric: Iwork (2*nrow+MAX(2*nsuper,ncol))
00042  *      where nsuper is 0 if simplicial, or the # of relaxed supernodes in
00043  *      L otherwise (nsuper <= nrow).
00044  *  if simplicial: W (nrow).
00045  *  Allocates up to two temporary copies of its input matrix (including
00046  *  both pattern and numerical values).
00047  *
00048  * If the matrix is not positive definite the routine returns TRUE, but
00049  * sets Common->status to CHOLMOD_NOT_POSDEF and L->minor is set to the
00050  * column at which the failure occurred.  Columns L->minor to L->n-1 are
00051  * set to zero.
00052  *
00053  * Supports any xtype (pattern, real, complex, or zomplex), except that the
00054  * input matrix A cannot be pattern-only.  If L is simplicial, its numeric
00055  * xtype matches A on output.  If L is supernodal, its xtype is real if A is
00056  * real, or complex if A is complex or zomplex.
00057  */
00058
00059 #ifndef NCHOLESKY
00060
00061 #include "amesos_cholmod_internal.h"
00062 #include "amesos_cholmod_cholesky.h"
00063
00064
00065 /* ========================================================================== */
00066 /* === cholmod_factorize ==================================================== */
00067 /* ========================================================================== */
00068
00069 /* Factorizes PAP' (or PAA'P' if A->stype is 0), using a factor obtained
00070  * from cholmod_analyze.  The analysis can be re-used simply by calling this
00071  * routine a second time with another matrix.  A must have the same nonzero
00072  * pattern as that passed to cholmod_analyze. */
00073
00074 int CHOLMOD(factorize)
00075 (
00076     /* ---- input ---- */
00077     cholmod_sparse *A,  /* matrix to factorize */
00078     /* ---- in/out --- */
00079     cholmod_factor *L,  /* resulting factorization */
00080     /* --------------- */
00081     cholmod_common *Common
00082 )
00083 {
00084     double zero [2] ;
00085     zero [0] = 0 ;
00086     zero [1] = 0 ;
00087     return (CHOLMOD(factorize_p) (A, zero, NULL, 0, L, Common)) ;
00088 }
00089
00090
00091 /* ========================================================================== */
00092 /* === cholmod_factorize_p ================================================== */
00093 /* ========================================================================== */
00094
00095 /* Same as cholmod_factorize, but with more options. */
00096
00097 int CHOLMOD(factorize_p)
00098 (
00099     /* ---- input ---- */
00100     cholmod_sparse *A,  /* matrix to factorize */
00101     double beta [2],  /* factorize beta*I+A or beta*I+A'*A */
00102     Int *fset,    /* subset of 0:(A->ncol)-1 */
00103     size_t fsize, /* size of fset */
00104     /* ---- in/out --- */
00105     cholmod_factor *L,  /* resulting factorization */
00106     /* --------------- */
00107     cholmod_common *Common
00108 )
00109 {
00110     cholmod_sparse *S, *F, *A1, *A2 ;
00111     Int nrow, ncol, stype, convert, n, nsuper, grow2, status ;
00112     size_t s, t, uncol ;
00113     int ok = TRUE ;
00114
00115     /* ---------------------------------------------------------------------- */
00116     /* check inputs */
00117     /* ---------------------------------------------------------------------- */
00118
00119     RETURN_IF_NULL_COMMON (FALSE) ;
00120     RETURN_IF_NULL (A, FALSE) ;
00121     RETURN_IF_NULL (L, FALSE) ;
00122     RETURN_IF_XTYPE_INVALID (A, CHOLMOD_REAL, CHOLMOD_ZOMPLEX, FALSE) ;
00123     RETURN_IF_XTYPE_INVALID (L, CHOLMOD_PATTERN, CHOLMOD_ZOMPLEX, FALSE) ;
00124     nrow = A->nrow ;
00125     ncol = A->ncol ;
00126     n = L->n ;
00127     stype = A->stype ;
00128     if (L->n != A->nrow)
00129     {
00130   ERROR (CHOLMOD_INVALID, "A and L dimensions do not match") ;
00131   return (FALSE) ;
00132     }
00133     if (stype != 0 && nrow != ncol)
00134     {
00135   ERROR (CHOLMOD_INVALID, "matrix invalid") ;
00136   return (FALSE) ;
00137     }
00138     DEBUG (CHOLMOD(dump_sparse) (A, "A for cholmod_factorize", Common)) ;
00139     Common->status = CHOLMOD_OK ;
00140
00141     /* ---------------------------------------------------------------------- */
00142     /* allocate workspace */
00143     /* ---------------------------------------------------------------------- */
00144
00145     nsuper = (L->is_super ? L->nsuper : 0) ;
00146     uncol = ((stype != 0) ? 0 : ncol) ;
00147
00148     /* s = 2*nrow + MAX (uncol, 2*nsuper) */
00149     s = CHOLMOD(mult_size_t) (nsuper, 2, &ok) ;
00150     s = MAX (uncol, s) ;
00151     t = CHOLMOD(mult_size_t) (nrow, 2, &ok) ;
00152     s = CHOLMOD(add_size_t) (s, t, &ok) ;
00153     if (!ok)
00154     {
00155   ERROR (CHOLMOD_TOO_LARGE, "problem too large") ;
00156   return (FALSE) ;
00157     }
00158
00159     CHOLMOD(allocate_work) (nrow, s, 0, Common) ;
00160     if (Common->status < CHOLMOD_OK)
00161     {
00162   return (FALSE) ;
00163     }
00164
00165     S  = NULL ;
00166     F  = NULL ;
00167     A1 = NULL ;
00168     A2 = NULL ;
00169
00170     /* convert to another form when done, if requested */
00171     convert = !(Common->final_asis) ;
00172
00173     /* ---------------------------------------------------------------------- */
00174     /* perform supernodal LL' or simplicial LDL' factorization */
00175     /* ---------------------------------------------------------------------- */
00176
00177     if (L->is_super)
00178     {
00179
00180 #ifndef NSUPERNODAL
00181
00182   /* ------------------------------------------------------------------ */
00183   /* supernodal factorization */
00184   /* ------------------------------------------------------------------ */
00185
00186   if (L->ordering == CHOLMOD_NATURAL)
00187   {
00188
00189       /* -------------------------------------------------------------- */
00190       /* natural ordering */
00191       /* -------------------------------------------------------------- */
00192
00193       if (stype > 0)
00194       {
00195     /* S = tril (A'), F not needed */
00196     /* workspace: Iwork (nrow) */
00197     A1 = CHOLMOD(ptranspose) (A, 2, NULL, NULL, 0, Common) ;
00198     S = A1 ;
00199       }
00200       else if (stype < 0)
00201       {
00202     /* This is the fastest option for the natural ordering */
00203     /* S = A; F not needed */
00204     S = A ;
00205       }
00206       else
00207       {
00208     /* F = A(:,f)' */
00209     /* workspace: Iwork (nrow) */
00210     /* workspace: Iwork (nrow if no fset; MAX (nrow,ncol) if fset)*/
00211     A1 = CHOLMOD(ptranspose) (A, 2, NULL, fset, fsize, Common) ;
00212     F = A1 ;
00213     /* S = A */
00214     S = A ;
00215       }
00216
00217   }
00218   else
00219   {
00220
00221       /* -------------------------------------------------------------- */
00222       /* permute the input matrix before factorization */
00223       /* -------------------------------------------------------------- */
00224
00225       if (stype > 0)
00226       {
00227     /* This is the fastest option for factoring a permuted matrix */
00228     /* S = tril (PAP'); F not needed */
00229     /* workspace: Iwork (2*nrow) */
00230     A1 = CHOLMOD(ptranspose) (A, 2, L->Perm, NULL, 0, Common) ;
00231     S = A1 ;
00232       }
00233       else if (stype < 0)
00234       {
00235     /* A2 = triu (PAP') */
00236     /* workspace: Iwork (2*nrow) */
00237     A2 = CHOLMOD(ptranspose) (A, 2, L->Perm, NULL, 0, Common) ;
00238     /* S = tril (A2'); F not needed */
00239     /* workspace: Iwork (nrow) */
00240     A1 = CHOLMOD(ptranspose) (A2, 2, NULL, NULL, 0, Common) ;
00241     S = A1 ;
00242     CHOLMOD(free_sparse) (&A2, Common) ;
00243     ASSERT (A2 == NULL) ;
00244       }
00245       else
00246       {
00247     /* F = A(p,f)' */
00248     /* workspace: Iwork (nrow if no fset; MAX (nrow,ncol) if fset)*/
00249     A1 = CHOLMOD(ptranspose) (A, 2, L->Perm, fset, fsize, Common) ;
00250     F = A1 ;
00251     /* S = F' */
00252     /* workspace: Iwork (nrow) */
00253     A2 = CHOLMOD(ptranspose) (F, 2, NULL, NULL, 0, Common) ;
00254     S = A2 ;
00255       }
00256   }
00257
00258   /* ------------------------------------------------------------------ */
00259   /* supernodal factorization */
00260   /* ------------------------------------------------------------------ */
00261
00262   /* workspace: Flag (nrow), Head (nrow+1), Iwork (2*nrow+2*nsuper) */
00263   if (Common->status == CHOLMOD_OK)
00264   {
00265       CHOLMOD(super_numeric) (S, F, beta, L, Common) ;
00266   }
00267   status = Common->status ;
00268   ASSERT (IMPLIES (status >= CHOLMOD_OK, L->xtype != CHOLMOD_PATTERN)) ;
00269
00270   /* ------------------------------------------------------------------ */
00271   /* convert to final form, if requested */
00272   /* ------------------------------------------------------------------ */
00273
00274   if (Common->status >= CHOLMOD_OK && convert)
00275   {
00276       /* workspace: none */
00277       ok = CHOLMOD(change_factor) (L->xtype, Common->final_ll,
00278         Common->final_super, Common->final_pack,
00279         Common->final_monotonic, L, Common) ;
00280       if (ok && Common->final_resymbol && !(L->is_super))
00281       {
00282     /* workspace: Flag (nrow), Head (nrow+1),
00283      *  if symmetric:   Iwork (2*nrow)
00284      *  if unsymmetric: Iwork (2*nrow+ncol) */
00285     CHOLMOD(resymbol_noperm) (S, fset, fsize, Common->final_pack,
00286         L, Common) ;
00287       }
00288   }
00289
00290 #else
00291
00292   /* ------------------------------------------------------------------ */
00293   /* CHOLMOD Supernodal module not installed */
00294   /* ------------------------------------------------------------------ */
00295
00296   status = CHOLMOD_NOT_INSTALLED ;
00297   ERROR (CHOLMOD_NOT_INSTALLED,"Supernodal module not installed") ;
00298
00299 #endif
00300
00301     }
00302     else
00303     {
00304
00305   /* ------------------------------------------------------------------ */
00306   /* simplicial LDL' factorization */
00307   /* ------------------------------------------------------------------ */
00308
00309   /* Permute the input matrix A if necessary.  cholmod_rowfac requires
00310    * triu(A) in column form for the symmetric case, and A in column form
00311    * for the unsymmetric case (the matrix S).  The unsymmetric case
00312    * requires A in row form, or equivalently A' in column form (the
00313    * matrix F).
00314    */
00315
00316   if (L->ordering == CHOLMOD_NATURAL)
00317   {
00318
00319       /* -------------------------------------------------------------- */
00320       /* natural ordering */
00321       /* -------------------------------------------------------------- */
00322
00323       if (stype > 0)
00324       {
00325     /* F is not needed, S = A */
00326     S = A ;
00327       }
00328       else if (stype < 0)
00329       {
00330     /* F is not needed, S = A' */
00331     /* workspace: Iwork (nrow) */
00332     A2 = CHOLMOD(ptranspose) (A, 2, NULL, NULL, 0, Common) ;
00333     S = A2 ;
00334       }
00335       else
00336       {
00337     /* F = A (:,f)' */
00338     /* workspace: Iwork (nrow if no fset; MAX (nrow,ncol) if fset)*/
00339     A1 = CHOLMOD(ptranspose) (A, 2, NULL, fset, fsize, Common) ;
00340     F = A1 ;
00341     S = A ;
00342       }
00343
00344   }
00345   else
00346   {
00347
00348       /* -------------------------------------------------------------- */
00349       /* permute the input matrix before factorization */
00350       /* -------------------------------------------------------------- */
00351
00352       if (stype > 0)
00353       {
00354     /* F = tril (A (p,p)') */
00355     /* workspace: Iwork (2*nrow) */
00356     A1 = CHOLMOD(ptranspose) (A, 2, L->Perm, NULL, 0, Common) ;
00357     /* A2 = triu (F') */
00358     /* workspace: Iwork (nrow) */
00359     A2 = CHOLMOD(ptranspose) (A1, 2, NULL, NULL, 0, Common) ;
00360     /* the symmetric case does not need F, free it and set to NULL*/
00361     CHOLMOD(free_sparse) (&A1, Common) ;
00362       }
00363       else if (stype < 0)
00364       {
00365     /* A2 = triu (A (p,p)'), F not needed.  This is the fastest
00366      * way to factorize a matrix using the simplicial routine
00367      * (cholmod_rowfac). */
00368     /* workspace: Iwork (2*nrow) */
00369     A2 = CHOLMOD(ptranspose) (A, 2, L->Perm, NULL, 0, Common) ;
00370       }
00371       else
00372       {
00373     /* F = A (p,f)' */
00374     /* workspace: Iwork (nrow if no fset; MAX (nrow,ncol) if fset)*/
00375     A1 = CHOLMOD(ptranspose) (A, 2, L->Perm, fset, fsize, Common) ;
00376     F = A1 ;
00377     /* A2 = F' */
00378     /* workspace: Iwork (nrow) */
00379     A2 = CHOLMOD(ptranspose) (F, 2, NULL, NULL, 0, Common) ;
00380       }
00381       S = A2 ;
00382   }
00383
00384   /* ------------------------------------------------------------------ */
00385   /* simplicial LDL' or LL' factorization */
00386   /* ------------------------------------------------------------------ */
00387
00388   /* factorize beta*I+S (symmetric) or beta*I+F*F' (unsymmetric) */
00389   /* workspace: Flag (nrow), W (nrow), Iwork (2*nrow) */
00390   if (Common->status == CHOLMOD_OK)
00391   {
00392       grow2 = Common->grow2 ;
00393       L->is_ll = BOOLEAN (Common->final_ll) ;
00394       if (L->xtype == CHOLMOD_PATTERN && Common->final_pack)
00395       {
00396     /* allocate a factor with exactly the space required */
00397     Common->grow2 = 0 ;
00398       }
00399       CHOLMOD(rowfac) (S, F, beta, 0, nrow, L, Common) ;
00400       Common->grow2 = grow2 ;
00401   }
00402   status = Common->status ;
00403
00404   /* ------------------------------------------------------------------ */
00405   /* convert to final form, if requested */
00406   /* ------------------------------------------------------------------ */
00407
00408   if (Common->status >= CHOLMOD_OK && convert)
00409   {
00410       /* workspace: none */
00411       CHOLMOD(change_factor) (L->xtype, L->is_ll, FALSE,
00412         Common->final_pack, Common->final_monotonic, L, Common) ;
00413   }
00414     }
00415
00416     /* ---------------------------------------------------------------------- */
00417     /* free A1 and A2 if they exist */
00418     /* ---------------------------------------------------------------------- */
00419
00420     CHOLMOD(free_sparse) (&A1, Common) ;
00421     CHOLMOD(free_sparse) (&A2, Common) ;
00422     Common->status = MAX (Common->status, status) ;
00423     return (Common->status >= CHOLMOD_OK) ;
00424 }
00425 #endif
```