Amesos Package Browser (Single Doxygen Collection) Development
amesos_t_cholmod_rowfac.c
Go to the documentation of this file.
00001 /* ========================================================================== */
00002 /* === Cholesky/t_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 /* Template routine for cholmod_rowfac.  Supports any numeric xtype
00014  * (real, complex, or zomplex).
00015  *
00016  * workspace: Iwork (n), Flag (n), Xwork (n if real, 2*n if complex)
00017  */
00018 
00019 #include "amesos_cholmod_template.h"
00020 
00021 #ifdef MASK
00022 static int TEMPLATE (cholmod_rowfac_mask)
00023 #else
00024 static int TEMPLATE (cholmod_rowfac)
00025 #endif
00026 (
00027     /* ---- input ---- */
00028     cholmod_sparse *A,  /* matrix to factorize */
00029     cholmod_sparse *F,  /* used for A*A' case only. F=A' or A(:,f)' */
00030     double beta [2],  /* factorize beta*I+A or beta*I+AA' (beta [0] only) */
00031     size_t kstart,  /* first row to factorize */
00032     size_t kend,  /* last row to factorize is kend-1 */
00033 #ifdef MASK
00034     /* These inputs are used for cholmod_rowfac_mask only */
00035     Int *mask,    /* size A->nrow. if mask[i] then W(i) is set to zero */
00036     Int *RLinkUp, /* size A->nrow. link list of rows to compute */
00037 #endif
00038     /* ---- in/out --- */
00039     cholmod_factor *L,
00040     /* --------------- */
00041     cholmod_common *Common
00042 )
00043 {
00044     double yx [2], lx [2], fx [2], dk [1], di [1], fl = 0 ;
00045 #ifdef ZOMPLEX
00046     double yz [1], lz [1], fz [1] ;
00047 #endif
00048     double *Ax, *Az, *Lx, *Lz, *Wx, *Wz, *Fx, *Fz ;
00049     Int *Ap, *Anz, *Ai, *Lp, *Lnz, *Li, *Lnext, *Flag, *Stack, *Fp, *Fi, *Fnz,
00050   *Iwork ;
00051     Int i, p, k, t, pf, pfend, top, s, mark, pend, n, lnz, is_ll, multadds,
00052   use_dbound, packed, stype, Fpacked, sorted, nzmax, len, parent ;
00053 #ifndef REAL
00054     Int dk_imaginary ;
00055 #endif
00056 
00057     /* ---------------------------------------------------------------------- */
00058     /* get inputs */
00059     /* ---------------------------------------------------------------------- */
00060 
00061     PRINT1 (("\nin cholmod_rowfac, kstart %d kend %d stype %d\n",
00062     kstart, kend, A->stype)) ;
00063     DEBUG (CHOLMOD(dump_factor) (L, "Initial L", Common)) ;
00064 
00065     n = A->nrow ;
00066     stype = A->stype ;
00067 
00068     if (stype > 0)
00069     {
00070   /* symmetric upper case: F is not needed.  It may be NULL */
00071   Fp = NULL ;
00072   Fi = NULL ;
00073   Fx = NULL ;
00074   Fz = NULL ;
00075   Fnz = NULL ;
00076   Fpacked = TRUE ;
00077     }
00078     else
00079     {
00080   /* unsymmetric case: F is required. */
00081   Fp = F->p ;
00082   Fi = F->i ;
00083   Fx = F->x ;
00084   Fz = F->z ;
00085   Fnz = F->nz ;
00086   Fpacked = F->packed ;
00087     }
00088 
00089     Ap = A->p ;   /* size A->ncol+1, column pointers of A */
00090     Ai = A->i ;   /* size nz = Ap [A->ncol], row indices of A */
00091     Ax = A->x ;   /* size nz, numeric values of A */
00092     Az = A->z ;
00093     Anz = A->nz ;
00094     packed = A->packed ;
00095     sorted = A->sorted ;
00096 
00097     use_dbound = IS_GT_ZERO (Common->dbound) ;
00098 
00099     /* get the current factors L (and D for LDL'); allocate space if needed */
00100     is_ll = L->is_ll ;
00101     if (L->xtype == CHOLMOD_PATTERN)
00102     {
00103   /* ------------------------------------------------------------------ */
00104   /* L is symbolic only; allocate and initialize L (and D for LDL') */
00105   /* ------------------------------------------------------------------ */
00106 
00107   /* workspace: none */
00108   CHOLMOD(change_factor) (A->xtype, is_ll, FALSE, FALSE, TRUE, L, Common);
00109   if (Common->status < CHOLMOD_OK)
00110   {
00111       /* out of memory */
00112       return (FALSE) ;
00113   }
00114   ASSERT (L->minor == (size_t) n) ;
00115     }
00116     else if (kstart == 0 && kend == (size_t) n)
00117     {
00118   /* ------------------------------------------------------------------ */
00119   /* refactorization; reset L->nz and L->minor to restart factorization */
00120   /* ------------------------------------------------------------------ */
00121 
00122   L->minor = n ;
00123   Lnz = L->nz ;
00124   for (k = 0 ; k < n ; k++)
00125   {
00126       Lnz [k] = 1 ;
00127   }
00128     }
00129 
00130     ASSERT (is_ll == L->is_ll) ;
00131     ASSERT (L->xtype != CHOLMOD_PATTERN) ;
00132     DEBUG (CHOLMOD(dump_factor) (L, "L ready", Common)) ;
00133     DEBUG (CHOLMOD(dump_sparse) (A, "A ready", Common)) ;
00134     DEBUG (if (stype == 0) CHOLMOD(dump_sparse) (F, "F ready", Common)) ;
00135 
00136     /* inputs, can be modified on output: */
00137     Lp = L->p ;   /* size n+1 */
00138     ASSERT (Lp != NULL) ;
00139 
00140     /* outputs, contents defined on input for incremental case only: */
00141     Lnz = L->nz ; /* size n */
00142     Lnext = L->next ; /* size n+2 */
00143     Li = L->i ;   /* size L->nzmax, can change in size */
00144     Lx = L->x ;   /* size L->nzmax or 2*L->nzmax, can change in size */
00145     Lz = L->z ;   /* size L->nzmax for zomplex case, can change in size */
00146     nzmax = L->nzmax ;
00147     ASSERT (Lnz != NULL && Li != NULL && Lx != NULL) ;
00148 
00149     /* ---------------------------------------------------------------------- */
00150     /* get workspace */
00151     /* ---------------------------------------------------------------------- */
00152 
00153     Iwork = Common->Iwork ;
00154     Stack = Iwork ;   /* size n (i/i/l) */
00155     Flag = Common->Flag ; /* size n, Flag [i] < mark must hold */
00156     Wx = Common->Xwork ;  /* size n if real, 2*n if complex or 
00157          * zomplex.  Xwork [i] == 0 must hold. */
00158     Wz = Wx + n ;   /* size n for zomplex case only */
00159     mark = Common->mark ;
00160     ASSERT ((Int) Common->xworksize >= (L->xtype == CHOLMOD_REAL ? 1:2)*n) ;
00161 
00162     /* ---------------------------------------------------------------------- */
00163     /* compute LDL' or LL' factorization by rows */
00164     /* ---------------------------------------------------------------------- */
00165 
00166 #ifdef MASK
00167 #define NEXT(k) k = RLinkUp [k]
00168 #else
00169 #define NEXT(k) k++
00170 #endif
00171 
00172     for (k = kstart ; k < ((Int) kend) ; NEXT(k))
00173     {
00174   PRINT1 (("\n===============K "ID" Lnz [k] "ID"\n", k, Lnz [k])) ;
00175 
00176   /* ------------------------------------------------------------------ */
00177   /* compute pattern of kth row of L and scatter kth input column */
00178   /* ------------------------------------------------------------------ */
00179 
00180   /* column k of L is currently empty */
00181   ASSERT (Lnz [k] == 1) ;
00182   ASSERT (CHOLMOD(dump_work) (TRUE, TRUE, 2*n, Common)) ;
00183 
00184   top = n ;   /* Stack is empty */
00185   Flag [k] = mark ; /* do not include diagonal entry in Stack */
00186 
00187   /* use Li [Lp [i]+1] for etree */
00188 #define PARENT(i) (Lnz [i] > 1) ? (Li [Lp [i] + 1]) : EMPTY
00189 
00190   if (stype > 0)
00191   {
00192       /* scatter kth col of triu (beta*I+AA'), get pattern L(k,:) */
00193       p = Ap [k] ;
00194       pend = (packed) ? (Ap [k+1]) : (p + Anz [k]) ;
00195       /* W [i] = Ax [i] ; scatter column of A */
00196 #define SCATTER ASSIGN(Wx,Wz,i, Ax,Az,p)
00197       SUBTREE ;
00198 #undef SCATTER
00199   }
00200   else
00201   {
00202       /* scatter kth col of triu (beta*I+AA'), get pattern L(k,:) */
00203       pf = Fp [k] ;
00204       pfend = (Fpacked) ? (Fp [k+1]) : (pf + Fnz [k]) ;
00205       for ( ; pf < pfend ; pf++)
00206       {
00207     /* get nonzero entry F (t,k) */
00208     t = Fi [pf] ;
00209     /* fk = Fx [pf] */
00210     ASSIGN (fx, fz, 0, Fx, Fz, pf) ;
00211     p = Ap [t] ;
00212     pend = (packed) ? (Ap [t+1]) : (p + Anz [t]) ;
00213     multadds = 0 ;
00214     /* W [i] += Ax [p] * fx ; scatter column of A*A' */
00215 #define SCATTER MULTADD (Wx,Wz,i, Ax,Az,p, fx,fz,0) ; multadds++  ;
00216     SUBTREE ;
00217 #undef SCATTER
00218 #ifdef REAL
00219     fl += 2 * ((double) multadds) ;
00220 #else
00221     fl += 8 * ((double) multadds) ;
00222 #endif
00223       }
00224   }
00225 
00226 #undef PARENT
00227 
00228   /* ------------------------------------------------------------------ */
00229   /* if mask is present, set the corresponding entries in W to zero */
00230   /* ------------------------------------------------------------------ */
00231 
00232 #ifdef MASK
00233         /* remove the dead element of Wx */
00234         if (mask != NULL)
00235         {
00236 
00237 #if 0
00238       /* older version */
00239             for (p = n; p > top;)
00240             {
00241                 i = Stack [--p] ;
00242                 if ( mask [i] >= 0 )
00243     {
00244         CLEAR (Wx,Wz,i) ; /* set W(i) to zero */
00245     }
00246             }
00247 #endif
00248 
00249             for (s = top ; s < n ; s++)
00250             {
00251                 i = Stack [s] ;
00252                 if (mask [i] >= 0)
00253     {
00254         CLEAR (Wx,Wz,i) ; /* set W(i) to zero */
00255     }
00256             }
00257 
00258         }
00259 #endif
00260 
00261   /* nonzero pattern of kth row of L is now in Stack [top..n-1].
00262    * Flag [Stack [top..n-1]] is equal to mark, but no longer needed */
00263 
00264   mark = CHOLMOD(clear_flag) (Common) ;
00265 
00266   /* ------------------------------------------------------------------ */
00267   /* compute kth row of L and store in column form */
00268   /* ------------------------------------------------------------------ */
00269 
00270   /* Solve L (0:k-1, 0:k-1) * y (0:k-1) = b (0:k-1) where
00271    * b (0:k) = A (0:k,k) or A(0:k,:) * F(:,k) is in W and Stack.
00272    *
00273    * For LDL' factorization:
00274    * L (k, 0:k-1) = y (0:k-1) ./ D (0:k-1)
00275    * D (k) = b (k) - L (k, 0:k-1) * y (0:k-1)
00276    *
00277    * For LL' factorization:
00278    * L (k, 0:k-1) = y (0:k-1)
00279    * L (k,k) = sqrt (b (k) - L (k, 0:k-1) * L (0:k-1, k))
00280    */
00281 
00282   /* dk = W [k] + beta */
00283   ADD_REAL (dk,0, Wx,k, beta,0) ;
00284 
00285 #ifndef REAL
00286   /* In the unsymmetric case, the imaginary part of W[k] must be real,
00287    * since F is assumed to be the complex conjugate transpose of A.  In
00288    * the symmetric case, W[k] is the diagonal of A.  If the imaginary part
00289    * of W[k] is nonzero, then the Cholesky factorization cannot be
00290    * computed; A is not positive definite */
00291   dk_imaginary = (stype > 0) ? (IMAG_IS_NONZERO (Wx,Wz,k)) : FALSE ;
00292 #endif
00293 
00294   /* W [k] = 0.0 ; */
00295   CLEAR (Wx,Wz,k) ;
00296 
00297   for (s = top ; s < n ; s++)
00298   {
00299       /* get i for each nonzero entry L(k,i) */
00300       i = Stack [s] ;
00301 
00302       /* y = W [i] ; */
00303       ASSIGN (yx,yz,0, Wx,Wz,i) ;
00304 
00305       /* W [i] = 0.0 ; */
00306       CLEAR (Wx,Wz,i) ;
00307 
00308       lnz = Lnz [i] ;
00309       p = Lp [i] ;
00310       ASSERT (lnz > 0 && Li [p] == i) ;
00311       pend = p + lnz ;
00312 
00313       /* di = Lx [p] ; the diagonal entry L or D(i,i), which is real */
00314       ASSIGN_REAL (di,0, Lx,p) ;
00315 
00316       if (i >= (Int) L->minor || IS_ZERO (di [0]))
00317       {
00318     /* For the LL' factorization, L(i,i) is zero.  For the LDL',
00319      * D(i,i) is zero.  Skip column i of L, and set L(k,i) = 0. */
00320     CLEAR (lx,lz,0) ;
00321     p = pend ;
00322       }
00323       else if (is_ll)
00324       {
00325 #ifdef REAL
00326     fl += 2 * ((double) (pend - p - 1)) + 3 ;
00327 #else
00328     fl += 8 * ((double) (pend - p - 1)) + 6 ;
00329 #endif
00330     /* forward solve using L (i:(k-1),i) */
00331     /* divide by L(i,i), which must be real and nonzero */
00332     /* y /= di [0] */
00333     DIV_REAL (yx,yz,0, yx,yz,0, di,0) ;
00334     for (p++ ; p < pend ; p++)
00335     {
00336         /* W [Li [p]] -= Lx [p] * y ; */
00337         MULTSUB (Wx,Wz,Li[p], Lx,Lz,p, yx,yz,0) ;
00338     }
00339     /* do not scale L; compute dot product for L(k,k) */
00340     /* L(k,i) = conj(y) ; */
00341     ASSIGN_CONJ (lx,lz,0, yx,yz,0) ;
00342     /* d -= conj(y) * y ; */
00343     LLDOT (dk,0, yx,yz,0) ;
00344       }
00345       else
00346       {
00347 #ifdef REAL
00348     fl += 2 * ((double) (pend - p - 1)) + 3 ;
00349 #else
00350     fl += 8 * ((double) (pend - p - 1)) + 6 ;
00351 #endif
00352     /* forward solve using D (i,i) and L ((i+1):(k-1),i) */
00353     for (p++ ; p < pend ; p++)
00354     {
00355         /* W [Li [p]] -= Lx [p] * y ; */
00356         MULTSUB (Wx,Wz,Li[p], Lx,Lz,p, yx,yz,0) ;
00357     }
00358     /* Scale L (k,0:k-1) for LDL' factorization, compute D (k,k)*/
00359 #ifdef REAL
00360     /* L(k,i) = y/d */
00361     lx [0] = yx [0] / di [0] ;
00362     /* d -= L(k,i) * y */
00363     dk [0] -= lx [0] * yx [0] ;
00364 #else
00365     /* L(k,i) = conj(y) ; */
00366     ASSIGN_CONJ (lx,lz,0, yx,yz,0) ;
00367     /* L(k,i) /= di ; */
00368     DIV_REAL (lx,lz,0, lx,lz,0, di,0) ;
00369     /* d -= conj(y) * y / di */
00370     LDLDOT (dk,0, yx,yz,0, di,0) ;
00371 #endif
00372       }
00373 
00374       /* determine if column i of L can hold the new L(k,i) entry */
00375       if (p >= Lp [Lnext [i]])
00376       {
00377     /* column i needs to grow */
00378     PRINT1 (("Factor Colrealloc "ID", old Lnz "ID"\n", i, Lnz [i]));
00379     if (!CHOLMOD(reallocate_column) (i, lnz + 1, L, Common))
00380     {
00381         /* out of memory, L is now simplicial symbolic */
00382         for (i = 0 ; i < n ; i++)
00383         {
00384       /* W [i] = 0 ; */
00385       CLEAR (Wx,Wz,i) ;
00386         }
00387         ASSERT (CHOLMOD(dump_work) (TRUE, TRUE, n, Common)) ;
00388         return (FALSE) ;
00389     }
00390     Li = L->i ;   /* L->i, L->x, L->z may have moved */
00391     Lx = L->x ;
00392     Lz = L->z ;
00393     p = Lp [i] + lnz ;  /* contents of L->p changed */
00394     ASSERT (p < Lp [Lnext [i]]) ;
00395       }
00396 
00397       /* store L (k,i) in the column form matrix of L */
00398       Li [p] = k ;
00399       /* Lx [p] = L(k,i) ; */
00400       ASSIGN (Lx,Lz,p, lx,lz,0) ;
00401       Lnz [i]++ ;
00402   }
00403 
00404   /* ------------------------------------------------------------------ */
00405   /* ensure abs (d) >= dbound if dbound is given, and store it in L */
00406   /* ------------------------------------------------------------------ */
00407 
00408   p = Lp [k] ;
00409   Li [p] = k ;
00410 
00411   if (k >= (Int) L->minor)
00412   {
00413       /* the matrix is already not positive definite */
00414       dk [0] = 0 ;
00415   }
00416   else if (use_dbound)
00417   {
00418       /* modify the diagonal to force LL' or LDL' to exist */
00419       dk [0] = CHOLMOD(dbound) (is_ll ? fabs (dk [0]) : dk [0], Common) ;
00420   }
00421   else if ((is_ll ? (IS_LE_ZERO (dk [0])) : (IS_ZERO (dk [0])))
00422 #ifndef REAL
00423     || dk_imaginary
00424 #endif
00425     )
00426   {
00427       /* the matrix has just been found to be not positive definite */
00428       dk [0] = 0 ;
00429       L->minor = k ;
00430       ERROR (CHOLMOD_NOT_POSDEF, "not positive definite") ;
00431   }
00432 
00433   if (is_ll)
00434   {
00435       /* this is counted as one flop, below */
00436       dk [0] = sqrt (dk [0]) ;
00437   }
00438 
00439   /* Lx [p] = D(k,k) = d ; real part only */
00440   ASSIGN_REAL (Lx,p, dk,0) ;
00441   CLEAR_IMAG (Lx,Lz,p) ;
00442     }
00443 
00444 #undef NEXT
00445 
00446     if (is_ll) fl += MAX ((Int) kend - (Int) kstart, 0) ;   /* count sqrt's */
00447     Common->rowfacfl = fl ;
00448 
00449     DEBUG (CHOLMOD(dump_factor) (L, "final cholmod_rowfac", Common)) ;
00450     ASSERT (CHOLMOD(dump_work) (TRUE, TRUE, n, Common)) ;
00451     return (TRUE) ;
00452 }
00453 #undef PATTERN
00454 #undef REAL
00455 #undef COMPLEX
00456 #undef ZOMPLEX
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Defines