Amesos Package Browser (Single Doxygen Collection) Development
amesos_cholmod_rowfac.c
Go to the documentation of this file.
00001 /* ========================================================================== */
00002 /* === Cholesky/cholmod_rowfac ============================================== */
00003 /* ========================================================================== */
00004 
00005 /* -----------------------------------------------------------------------------
00006  * CHOLMOD/Cholesky Module.  Copyright (C) 2005-2006, Timothy A. Davis
00007  * The CHOLMOD/Cholesky Module is licensed under Version 2.1 of the GNU
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 /* Full or incremental numerical LDL' or LL' factorization (simplicial, not
00014  * supernodal) cholmod_factorize is the "easy" wrapper for this code, but it
00015  * does not provide access to incremental factorization.
00016  *
00017  * cholmod_rowfac computes the full or incremental LDL' or LL' factorization of
00018  * A+beta*I (where A is symmetric) or A*F+beta*I (where A and F are unsymmetric
00019  * and only the upper triangular part of A*F+beta*I is used).  It computes
00020  * L (and D, for LDL') one row at a time.  beta is real.
00021  *
00022  * A is nrow-by-ncol or nrow-by-nrow.  In "packed" form it is a conventional
00023  * column-oriented sparse matrix.  Row indices of column j are in
00024  * Ai [Ap [j] ... Ap [j+1]-1] and values in the same locations of Ax.
00025  * will be faster if A has sorted columns.  In "unpacked" form the column
00026  * of A ends at Ap [j] + Anz [j] - 1 instead of Ap [j+1] - 1.
00027  *
00028  * Row indices in each column of A can be sorted or unsorted, but the routine
00029  * routine works fastest if A is sorted, or if only triu(A) is provided
00030  * for the symmetric case.
00031  *
00032  * The unit-diagonal nrow-by-nrow output matrix L is returned in "unpacked"
00033  * column form, with row indices of column j in Li [Lp [j] ...
00034  * Lp [j] + Lnz [j] - 1] and values in the same location in Lx.  The row
00035  * indices in each column of L are in sorted order.  The unit diagonal of L
00036  * is not stored.
00037  *
00038  * L can be a simplicial symbolic or numeric (L->is_super must be FALSE).
00039  * A symbolic factor is converted immediately into a numeric factor containing
00040  * the identity matrix.
00041  *
00042  * For a full factorization, kstart = 0 and kend = nrow.  The existing nonzero
00043  * entries (numerical values in L->x and L->z for the zomplex case, and indices
00044  * in L->i), if any, are overwritten.
00045  *
00046  * To compute an incremental factorization, select kstart and kend as the range
00047  * of rows of L you wish to compute.  A correct factorization will be computed
00048  * only if all descendants of all nodes k = kstart to kend-1 in the etree have
00049  * been factorized by a prior call to this routine, and if rows kstart to kend-1
00050  * have not been factorized.  This condition is NOT checked on input.
00051  *
00052  * ---------------
00053  * Symmetric case:
00054  * ---------------
00055  *
00056  *  The factorization (in MATLAB notation) is:
00057  *
00058  *  S = beta*I + A
00059  *  S = triu (S) + triu (S,1)'
00060  *  L*D*L' = S, or L*L' = S
00061  *
00062  *  A is a conventional sparse matrix in compressed column form.  Only the
00063  *  diagonal and upper triangular part of A is accessed; the lower
00064  *  triangular part is ignored and assumed to be equal to the upper
00065  *  triangular part.  For an incremental factorization, only columns kstart
00066  *  to kend-1 of A are accessed.  F is not used.
00067  *
00068  * ---------------
00069  * Unsymmetric case:
00070  * ---------------
00071  *
00072  *  The factorization (in MATLAB notation) is:
00073  *
00074  *  S = beta*I + A*F
00075  *  S = triu (S) + triu (S,1)'
00076  *  L*D*L' = S, or L*L' = S
00077  *
00078  *  The typical case is F=A'.  Alternatively, if F=A(:,f)', then this
00079  *  routine factorizes S = beta*I + A(:,f)*A(:,f)'.
00080  *
00081  *  All of A and F are accessed, but only the upper triangular part of A*F
00082  *  is used.  F must be of size A->ncol by A->nrow.  F is used for the
00083  *  unsymmetric case only.  F can be packed or unpacked and it need not be
00084  *  sorted.
00085  *
00086  *  For a complete factorization of beta*I + A*A',
00087  *  this routine performs a number of flops exactly equal to:
00088  *
00089  *  sum (for each column j of A) of (Anz (j)^2 + Anz (j)), to form S
00090  *  +
00091  *  sum (for each column j of L) of (Lnz (j)^2 + 3*Lnz (j)), to factorize S
00092  *
00093  *  where Anz (j) is the number of nonzeros in column j of A, and Lnz (j)
00094  *  is the number of nonzero in column j of L below the diagonal.
00095  *
00096  *
00097  * workspace: Flag (nrow), W (nrow if real, 2*nrow if complex/zomplex),
00098  * Iwork (nrow)
00099  *
00100  * Supports any xtype, except a pattern-only input matrix A cannot be
00101  * factorized.
00102  */
00103 
00104 #ifndef NCHOLESKY
00105 
00106 #include "amesos_cholmod_internal.h"
00107 #include "amesos_cholmod_cholesky.h"
00108 
00109 /* ========================================================================== */
00110 /* === subtree ============================================================== */
00111 /* ========================================================================== */
00112 
00113 /* Compute the nonzero pattern of the sparse triangular solve Lx=b, where L in
00114  * this case is L(0:k-1,0:k-1), and b is a column of A.  This is done by
00115  * traversing the kth row-subtree of the elimination tree of L, starting from
00116  * each nonzero entry in b.  The pattern is returned postordered, and is valid
00117  * for a subsequent numerical triangular solve of Lx=b.  The elimination tree
00118  * can be provided in a Parent array, or extracted from the pattern of L itself.
00119  *
00120  * The pattern of x = inv(L)*b is returned in Stack [top...].
00121  * Also scatters b, or a multiple of b, into the work vector W.
00122  *
00123  * The SCATTER macro is defines how the numerical values of A or A*A' are to be
00124  * scattered.
00125  *
00126  * PARENT(i) is a macro the defines how the etree is accessed.  It is either:
00127  *  #define PARENT(i) Parent [i]
00128  *  #define PARENT(i) (Lnz [i] > 1) ? (Li [Lp [i] + 1]) : EMPTY
00129  */
00130 
00131 #define SUBTREE \
00132     for ( ; p < pend ; p++) \
00133     { \
00134   i = Ai [p] ; \
00135   if (i <= k) \
00136   { \
00137       /* scatter the column of A, or A*A' into Wx and Wz */ \
00138       SCATTER ; \
00139       /* start at node i and traverse up the subtree, stop at node k */ \
00140       for (len = 0 ; i < k && i != EMPTY && Flag [i] < mark ; i = parent) \
00141       { \
00142     /* L(k,i) is nonzero, and seen for the first time */ \
00143     Stack [len++] = i ;     /* place i on the stack */ \
00144     Flag [i] = mark ;     /* mark i as visited */ \
00145     parent = PARENT (i) ;   /* traverse up the etree to the parent */ \
00146       } \
00147       /* move the path down to the bottom of the stack */ \
00148       while (len > 0) \
00149       { \
00150     Stack [--top] = Stack [--len] ; \
00151       } \
00152   } \
00153   else if (sorted) \
00154   { \
00155       break ; \
00156   } \
00157     }
00158 
00159 
00160 /* ========================================================================== */
00161 /* === TEMPLATE ============================================================= */
00162 /* ========================================================================== */
00163 
00164 #define REAL
00165 #include "amesos_t_cholmod_rowfac.c"
00166 #define COMPLEX
00167 #include "amesos_t_cholmod_rowfac.c"
00168 #define ZOMPLEX
00169 #include "amesos_t_cholmod_rowfac.c"
00170 
00171 #define MASK
00172 #define REAL
00173 #include "amesos_t_cholmod_rowfac.c"
00174 #define COMPLEX
00175 #include "amesos_t_cholmod_rowfac.c"
00176 #define ZOMPLEX
00177 #include "amesos_t_cholmod_rowfac.c"
00178 #undef MASK
00179 
00180 
00181 /* ========================================================================== */
00182 /* === cholmod_row_subtree ================================================== */
00183 /* ========================================================================== */
00184 
00185 /* Compute the nonzero pattern of the solution to the lower triangular system
00186  * L(0:k-1,0:k-1) * x = A (0:k-1,k) if A is symmetric, or
00187  * L(0:k-1,0:k-1) * x = A (0:k-1,:) * A (:,k)' if A is unsymmetric.
00188  * This gives the nonzero pattern of row k of L (excluding the diagonal).
00189  * The pattern is returned postordered.
00190  *
00191  * The symmetric case requires A to be in symmetric-upper form.
00192  *
00193  * The result is returned in R, a pre-allocated sparse matrix of size nrow-by-1,
00194  * with R->nzmax >= nrow.  R is assumed to be packed (Rnz [0] is not updated);
00195  * the number of entries in R is given by Rp [0].
00196  *
00197  * FUTURE WORK:  a very minor change to this routine could allow it to compute
00198  * the nonzero pattern of x for any system Lx=b.  The SUBTREE macro would need
00199  * to change, to eliminate its dependence on k.
00200  *
00201  * workspace: Flag (nrow)
00202  */
00203 
00204 int CHOLMOD(row_subtree)
00205 (
00206     /* ---- input ---- */
00207     cholmod_sparse *A,  /* matrix to analyze */
00208     cholmod_sparse *F,  /* used for A*A' case only. F=A' or A(:,f)' */
00209     size_t krow,  /* row k of L */
00210     Int *Parent,  /* elimination tree */
00211     /* ---- output --- */
00212     cholmod_sparse *R,  /* pattern of L(k,:), 1-by-n with R->nzmax >= n */
00213     /* --------------- */
00214     cholmod_common *Common
00215 )
00216 {
00217     Int *Rp, *Stack, *Flag, *Ap, *Ai, *Anz, *Fp, *Fi, *Fnz ;
00218     Int p, pend, parent, t, stype, nrow, k, pf, pfend, Fpacked, packed,
00219   sorted, top, len, i, mark ;
00220 
00221     /* ---------------------------------------------------------------------- */
00222     /* check inputs */
00223     /* ---------------------------------------------------------------------- */
00224 
00225     RETURN_IF_NULL_COMMON (FALSE) ;
00226     RETURN_IF_NULL (A, FALSE) ;
00227     RETURN_IF_NULL (R, FALSE) ;
00228     RETURN_IF_NULL (Parent, FALSE) ;
00229     RETURN_IF_XTYPE_INVALID (A, CHOLMOD_PATTERN, CHOLMOD_ZOMPLEX, FALSE) ;
00230     RETURN_IF_XTYPE_INVALID (R, CHOLMOD_PATTERN, CHOLMOD_ZOMPLEX, FALSE) ;
00231     stype = A->stype ;
00232     if (stype == 0)
00233     {
00234   RETURN_IF_NULL (F, FALSE) ;
00235   RETURN_IF_XTYPE_INVALID (F, CHOLMOD_PATTERN, CHOLMOD_ZOMPLEX, FALSE) ;
00236     }
00237     if (krow >= A->nrow)
00238     {
00239   ERROR (CHOLMOD_INVALID, "subtree: k invalid") ;
00240   return (FALSE) ;
00241     }
00242     if (R->ncol != 1 || A->nrow != R->nrow || A->nrow > R->nzmax)
00243     {
00244   ERROR (CHOLMOD_INVALID, "subtree: R invalid") ;
00245   return (FALSE) ;
00246     }
00247     Common->status = CHOLMOD_OK ;
00248 
00249     /* ---------------------------------------------------------------------- */
00250     /* allocate workspace */
00251     /* ---------------------------------------------------------------------- */
00252 
00253     nrow = A->nrow ;
00254     CHOLMOD(allocate_work) (nrow, 0, 0, Common) ;
00255     if (Common->status < CHOLMOD_OK)
00256     {
00257   return (FALSE) ;
00258     }
00259     ASSERT (CHOLMOD(dump_work) (TRUE, TRUE, 0, Common)) ;
00260 
00261     /* ---------------------------------------------------------------------- */
00262     /* get inputs */
00263     /* ---------------------------------------------------------------------- */
00264 
00265     if (stype > 0)
00266     {
00267   /* symmetric upper case: F is not needed.  It may be NULL */
00268   Fp = NULL ;
00269   Fi = NULL ;
00270   Fnz = NULL ;
00271   Fpacked = TRUE ;
00272     }
00273     else if (stype == 0)
00274     {
00275   /* unsymmetric case: F is required. */
00276   Fp = F->p ;
00277   Fi = F->i ;
00278   Fnz = F->nz ;
00279   Fpacked = F->packed ;
00280     }
00281     else
00282     {
00283   /* symmetric lower triangular form not supported */
00284   ERROR (CHOLMOD_INVALID, "symmetric lower not supported") ;
00285   return (FALSE) ;
00286     }
00287 
00288     Ap = A->p ;
00289     Ai = A->i ;
00290     Anz = A->nz ;
00291     packed = A->packed ;
00292     sorted = A->sorted ;
00293 
00294     k = krow ;
00295     Stack = R->i ;
00296 
00297     /* ---------------------------------------------------------------------- */
00298     /* get workspace */
00299     /* ---------------------------------------------------------------------- */
00300 
00301     Flag = Common->Flag ; /* size nrow, Flag [i] < mark must hold */
00302     mark = CHOLMOD(clear_flag) (Common) ;
00303 
00304     /* ---------------------------------------------------------------------- */
00305     /* compute the pattern of L(k,:) */
00306     /* ---------------------------------------------------------------------- */
00307 
00308     top = nrow ;    /* Stack is empty */
00309     Flag [k] = mark ;   /* do not include diagonal entry in Stack */
00310 
00311 #define SCATTER     /* do not scatter numerical values */
00312 #define PARENT(i) Parent [i]  /* use Parent for etree */
00313 
00314     if (stype != 0)
00315     {
00316   /* scatter kth col of triu (A), get pattern L(k,:) */
00317   p = Ap [k] ;
00318   pend = (packed) ? (Ap [k+1]) : (p + Anz [k]) ;
00319   SUBTREE ;
00320     }
00321     else
00322     {
00323   /* scatter kth col of triu (beta*I+AA'), get pattern L(k,:) */
00324   pf = Fp [k] ;
00325   pfend = (Fpacked) ? (Fp [k+1]) : (pf + Fnz [k]) ;
00326   for ( ; pf < pfend ; pf++)
00327   {
00328       /* get nonzero entry F (t,k) */
00329       t = Fi [pf] ;
00330       p = Ap [t] ;
00331       pend = (packed) ? (Ap [t+1]) : (p + Anz [t]) ;
00332       SUBTREE ;
00333   }
00334     }
00335 
00336 #undef SCATTER
00337 #undef PARENT
00338 
00339     /* shift the stack upwards, to the first part of R */
00340     len = nrow - top ;
00341     for (i = 0 ; i < len ; i++)
00342     {
00343   Stack [i] = Stack [top + i] ;
00344     }
00345 
00346     Rp = R->p ;
00347     Rp [0] = 0 ;
00348     Rp [1] = len ;
00349     R->sorted = FALSE ;
00350 
00351     CHOLMOD(clear_flag) (Common) ;
00352     ASSERT (CHOLMOD(dump_work) (TRUE, TRUE, 0, Common)) ;
00353     return (TRUE) ;
00354 }
00355 
00356 
00357 /* ========================================================================== */
00358 /* === cholmod_row_lsubtree ================================================= */
00359 /* ========================================================================== */
00360 
00361 /* Identical to cholmod_row_subtree, except that the elimination tree is
00362  * obtained from L itself, as the first off-diagonal entry in each column.
00363  * L must be simplicial, not supernodal */
00364 
00365 int CHOLMOD(row_lsubtree)
00366 (
00367     /* ---- input ---- */
00368     cholmod_sparse *A,  /* matrix to analyze */
00369     Int *Fi, size_t fnz,    /* nonzero pattern of kth row of A', not required
00370            * for the symmetric case.  Need not be sorted. */
00371     size_t krow,  /* row k of L */
00372     cholmod_factor *L,  /* the factor L from which parent(i) is derived */
00373     /* ---- output --- */
00374     cholmod_sparse *R,  /* pattern of L(k,:), 1-by-n with R->nzmax >= n */
00375     /* --------------- */
00376     cholmod_common *Common
00377 )
00378 {
00379     Int *Rp, *Stack, *Flag, *Ap, *Ai, *Anz, *Lp, *Li, *Lnz ;
00380     Int p, pend, parent, t, stype, nrow, k, pf, packed, sorted, top, len, i,
00381   mark ;
00382 
00383     /* ---------------------------------------------------------------------- */
00384     /* check inputs */
00385     /* ---------------------------------------------------------------------- */
00386 
00387     RETURN_IF_NULL_COMMON (FALSE) ;
00388     RETURN_IF_NULL (A, FALSE) ;
00389     RETURN_IF_NULL (R, FALSE) ;
00390     RETURN_IF_NULL (L, FALSE) ;
00391     RETURN_IF_XTYPE_INVALID (A, CHOLMOD_PATTERN, CHOLMOD_ZOMPLEX, FALSE) ;
00392     RETURN_IF_XTYPE_INVALID (R, CHOLMOD_PATTERN, CHOLMOD_ZOMPLEX, FALSE) ;
00393     RETURN_IF_XTYPE_INVALID (L, CHOLMOD_REAL, CHOLMOD_ZOMPLEX, FALSE) ;
00394     stype = A->stype ;
00395     if (stype == 0)
00396     {
00397   RETURN_IF_NULL (Fi, FALSE) ;
00398     }
00399     if (krow >= A->nrow)
00400     {
00401   ERROR (CHOLMOD_INVALID, "lsubtree: k invalid") ;
00402   return (FALSE) ;
00403     }
00404     if (R->ncol != 1 || A->nrow != R->nrow || A->nrow > R->nzmax)
00405     {
00406   ERROR (CHOLMOD_INVALID, "lsubtree: R invalid") ;
00407   return (FALSE) ;
00408     }
00409     if (L->is_super)
00410     {
00411   ERROR (CHOLMOD_INVALID, "lsubtree: L invalid (cannot be supernodal)") ;
00412   return (FALSE) ;
00413     }
00414     Common->status = CHOLMOD_OK ;
00415 
00416     /* ---------------------------------------------------------------------- */
00417     /* allocate workspace */
00418     /* ---------------------------------------------------------------------- */
00419 
00420     nrow = A->nrow ;
00421     CHOLMOD(allocate_work) (nrow, 0, 0, Common) ;
00422     if (Common->status < CHOLMOD_OK)
00423     {
00424   return (FALSE) ;
00425     }
00426     ASSERT (CHOLMOD(dump_work) (TRUE, TRUE, 0, Common)) ;
00427 
00428     /* ---------------------------------------------------------------------- */
00429     /* get inputs */
00430     /* ---------------------------------------------------------------------- */
00431 
00432     if (stype < 0)
00433     {
00434   /* symmetric lower triangular form not supported */
00435   ERROR (CHOLMOD_INVALID, "symmetric lower not supported") ;
00436   return (FALSE) ;
00437     }
00438 
00439     Ap = A->p ;
00440     Ai = A->i ;
00441     Anz = A->nz ;
00442     packed = A->packed ;
00443     sorted = A->sorted ;
00444 
00445     k = krow ;
00446     Stack = R->i ;
00447 
00448     Lp = L->p ;
00449     Li = L->i ;
00450     Lnz = L->nz ;
00451 
00452     /* ---------------------------------------------------------------------- */
00453     /* get workspace */
00454     /* ---------------------------------------------------------------------- */
00455 
00456     Flag = Common->Flag ; /* size nrow, Flag [i] < mark must hold */
00457     mark = CHOLMOD(clear_flag) (Common) ;
00458 
00459     /* ---------------------------------------------------------------------- */
00460     /* compute the pattern of L(k,:) */
00461     /* ---------------------------------------------------------------------- */
00462 
00463     top = nrow ;    /* Stack is empty */
00464     Flag [k] = mark ;   /* do not include diagonal entry in Stack */
00465 
00466 #define SCATTER     /* do not scatter numerical values */
00467 #define PARENT(i) (Lnz [i] > 1) ? (Li [Lp [i] + 1]) : EMPTY
00468 
00469     if (stype != 0)
00470     {
00471   /* scatter kth col of triu (A), get pattern L(k,:) */
00472   p = Ap [k] ;
00473   pend = (packed) ? (Ap [k+1]) : (p + Anz [k]) ;
00474   SUBTREE ;
00475     }
00476     else
00477     {
00478   /* scatter kth col of triu (beta*I+AA'), get pattern L(k,:) */
00479   for (pf = 0 ; pf < (Int) fnz ; pf++)
00480   {
00481       /* get nonzero entry F (t,k) */
00482       t = Fi [pf] ;
00483       p = Ap [t] ;
00484       pend = (packed) ? (Ap [t+1]) : (p + Anz [t]) ;
00485       SUBTREE ;
00486   }
00487     }
00488 
00489 #undef SCATTER
00490 #undef PARENT
00491 
00492     /* shift the stack upwards, to the first part of R */
00493     len = nrow - top ;
00494     for (i = 0 ; i < len ; i++)
00495     {
00496   Stack [i] = Stack [top + i] ;
00497     }
00498 
00499     Rp = R->p ;
00500     Rp [0] = 0 ;
00501     Rp [1] = len ;
00502     R->sorted = FALSE ;
00503 
00504     CHOLMOD(clear_flag) (Common) ;
00505     ASSERT (CHOLMOD(dump_work) (TRUE, TRUE, 0, Common)) ;
00506     return (TRUE) ;
00507 }
00508 
00509 
00510 /* ========================================================================== */
00511 /* === cholmod_rowfac ======================================================= */
00512 /* ========================================================================== */
00513 
00514 /* This is the incremental factorization for general purpose usage. */
00515 
00516 int CHOLMOD(rowfac)
00517 (
00518     /* ---- input ---- */
00519     cholmod_sparse *A,  /* matrix to factorize */
00520     cholmod_sparse *F,  /* used for A*A' case only. F=A' or A(:,f)' */
00521     double beta [2],  /* factorize beta*I+A or beta*I+AA' */
00522     size_t kstart,  /* first row to factorize */
00523     size_t kend,  /* last row to factorize is kend-1 */
00524     /* ---- in/out --- */
00525     cholmod_factor *L,
00526     /* --------------- */
00527     cholmod_common *Common
00528 )
00529 {
00530     return (CHOLMOD(rowfac_mask) (A, F, beta, kstart, kend, NULL, NULL, L,
00531   Common)) ;
00532 }
00533 
00534 
00535 /* ========================================================================== */
00536 /* === cholmod_rowfac_mask ================================================== */
00537 /* ========================================================================== */
00538 
00539 /* This is meant for use in LPDASA only. */
00540 
00541 int CHOLMOD(rowfac_mask)
00542 (
00543     /* ---- input ---- */
00544     cholmod_sparse *A,  /* matrix to factorize */
00545     cholmod_sparse *F,  /* used for A*A' case only. F=A' or A(:,f)' */
00546     double beta [2],  /* factorize beta*I+A or beta*I+AA' */
00547     size_t kstart,  /* first row to factorize */
00548     size_t kend,  /* last row to factorize is kend-1 */
00549     Int *mask,    /* size A->nrow. if mask[i] >= 0 row i is set to zero */
00550     Int *RLinkUp, /* size A->nrow. link list of rows to compute */
00551     /* ---- in/out --- */
00552     cholmod_factor *L,
00553     /* --------------- */
00554     cholmod_common *Common
00555 )
00556 {
00557     Int n ;
00558     size_t s ;
00559     int ok = TRUE ;
00560 
00561     /* ---------------------------------------------------------------------- */
00562     /* check inputs */
00563     /* ---------------------------------------------------------------------- */
00564 
00565     RETURN_IF_NULL_COMMON (FALSE) ;
00566     RETURN_IF_NULL (A, FALSE) ;
00567     RETURN_IF_NULL (L, FALSE) ;
00568     RETURN_IF_XTYPE_INVALID (A, CHOLMOD_REAL, CHOLMOD_ZOMPLEX, FALSE) ;
00569     RETURN_IF_XTYPE_INVALID (L, CHOLMOD_PATTERN, CHOLMOD_ZOMPLEX, FALSE) ;
00570     if (L->xtype != CHOLMOD_PATTERN && A->xtype != L->xtype)
00571     {
00572   ERROR (CHOLMOD_INVALID, "xtype of A and L do not match") ;
00573   return (FALSE) ;
00574     }
00575     if (L->is_super)
00576     {
00577   ERROR (CHOLMOD_INVALID, "can only do simplicial factorization");
00578   return (FALSE) ;
00579     }
00580     if (A->stype == 0)
00581     {
00582   RETURN_IF_NULL (F, FALSE) ;
00583   if (A->xtype != F->xtype)
00584   {
00585       ERROR (CHOLMOD_INVALID, "xtype of A and F do not match") ;
00586       return (FALSE) ;
00587   }
00588     }
00589     if (A->stype < 0)
00590     {
00591   /* symmetric lower triangular form not supported */
00592   ERROR (CHOLMOD_INVALID, "symmetric lower not supported") ;
00593   return (FALSE) ;
00594     }
00595     if (kend > L->n)
00596     {
00597   ERROR (CHOLMOD_INVALID, "kend invalid") ;
00598   return (FALSE) ;
00599     }
00600     if (A->nrow != L->n)
00601     {
00602   ERROR (CHOLMOD_INVALID, "dimensions of A and L do not match") ;
00603   return (FALSE) ;
00604     }
00605     Common->status = CHOLMOD_OK ;
00606     Common->rowfacfl = 0 ;
00607 
00608     /* ---------------------------------------------------------------------- */
00609     /* allocate workspace */
00610     /* ---------------------------------------------------------------------- */
00611 
00612     /* Xwork is of size n for the real case, 2*n for complex/zomplex */
00613     n = L->n  ;
00614 
00615     /* s = ((A->xtype != CHOLMOD_REAL) ? 2:1)*n */
00616     s = CHOLMOD(mult_size_t) (n, ((A->xtype != CHOLMOD_REAL) ? 2:1), &ok) ;
00617     if (!ok)
00618     {
00619   ERROR (CHOLMOD_TOO_LARGE, "problem too large") ;
00620   return (FALSE) ;
00621     }
00622 
00623     CHOLMOD(allocate_work) (n, n, s, Common) ;
00624     if (Common->status < CHOLMOD_OK)
00625     {
00626   return (FALSE) ;
00627     }
00628     ASSERT (CHOLMOD(dump_work) (TRUE, TRUE, A->nrow, Common)) ;
00629 
00630     /* ---------------------------------------------------------------------- */
00631     /* factorize the matrix, using template routine */
00632     /* ---------------------------------------------------------------------- */
00633 
00634     if (RLinkUp == NULL)
00635     {
00636 
00637   switch (A->xtype)
00638   {
00639       case CHOLMOD_REAL:
00640     ok = r_cholmod_rowfac (A, F, beta, kstart, kend, L, Common) ;
00641     break ;
00642 
00643       case CHOLMOD_COMPLEX:
00644     ok = c_cholmod_rowfac (A, F, beta, kstart, kend, L, Common) ;
00645     break ;
00646 
00647       case CHOLMOD_ZOMPLEX:
00648     ok = z_cholmod_rowfac (A, F, beta, kstart, kend, L, Common) ;
00649     break ;
00650   }
00651 
00652     }
00653     else
00654     {
00655 
00656   switch (A->xtype)
00657   {
00658       case CHOLMOD_REAL:
00659     ok = r_cholmod_rowfac_mask (A, F, beta, kstart, kend,
00660         mask, RLinkUp, L, Common) ;
00661     break ;
00662 
00663       case CHOLMOD_COMPLEX:
00664     ok = c_cholmod_rowfac_mask (A, F, beta, kstart, kend,
00665         mask, RLinkUp, L, Common) ;
00666     break ;
00667 
00668       case CHOLMOD_ZOMPLEX:
00669     ok = z_cholmod_rowfac_mask (A, F, beta, kstart, kend,
00670         mask, RLinkUp, L, Common) ;
00671     break ;
00672   }
00673     }
00674 
00675     return (ok) ;
00676 }
00677 #endif
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Defines