/* ========================================================================== */
/* === Cholesky/t_cholmod_rowfac ============================================ */
/* ========================================================================== */
/* -----------------------------------------------------------------------------
* CHOLMOD/Cholesky Module. Copyright (C) 2005-2006, Timothy A. Davis
* The CHOLMOD/Cholesky Module is licensed under Version 2.1 of the GNU
* Lesser General Public License. See lesser.txt for a text of the license.
* CHOLMOD is also available under other licenses; contact authors for details.
* http://www.cise.ufl.edu/research/sparse
* -------------------------------------------------------------------------- */
/* Template routine for cholmod_rowfac. Supports any numeric xtype
* (real, complex, or zomplex).
*
* workspace: Iwork (n), Flag (n), Xwork (n if real, 2*n if complex)
*/
#include "cholmod_template.h"
#ifdef MASK
static int TEMPLATE (cholmod_rowfac_mask)
#else
static int TEMPLATE (cholmod_rowfac)
#endif
(
/* ---- input ---- */
cholmod_sparse *A, /* matrix to factorize */
cholmod_sparse *F, /* used for A*A' case only. F=A' or A(:,f)' */
double beta [2], /* factorize beta*I+A or beta*I+AA' (beta [0] only) */
size_t kstart, /* first row to factorize */
size_t kend, /* last row to factorize is kend-1 */
#ifdef MASK
/* These inputs are used for cholmod_rowfac_mask only */
Int *mask, /* size A->nrow. if mask[i] then W(i) is set to zero */
Int *RLinkUp, /* size A->nrow. link list of rows to compute */
#endif
/* ---- in/out --- */
cholmod_factor *L,
/* --------------- */
cholmod_common *Common
)
{
double yx [2], lx [2], fx [2], dk [1], di [1], fl = 0 ;
#ifdef ZOMPLEX
double yz [1], lz [1], fz [1] ;
#endif
double *Ax, *Az, *Lx, *Lz, *Wx, *Wz, *Fx, *Fz ;
Int *Ap, *Anz, *Ai, *Lp, *Lnz, *Li, *Lnext, *Flag, *Stack, *Fp, *Fi, *Fnz,
*Iwork ;
Int i, p, k, t, pf, pfend, top, s, mark, pend, n, lnz, is_ll, multadds,
use_dbound, packed, stype, Fpacked, sorted, nzmax, len, parent ;
#ifndef REAL
Int dk_imaginary ;
#endif
/* ---------------------------------------------------------------------- */
/* get inputs */
/* ---------------------------------------------------------------------- */
PRINT1 (("\nin cholmod_rowfac, kstart %d kend %d stype %d\n",
kstart, kend, A->stype)) ;
DEBUG (CHOLMOD(dump_factor) (L, "Initial L", Common)) ;
n = A->nrow ;
stype = A->stype ;
if (stype > 0)
{
/* symmetric upper case: F is not needed. It may be NULL */
Fp = NULL ;
Fi = NULL ;
Fx = NULL ;
Fz = NULL ;
Fnz = NULL ;
Fpacked = TRUE ;
}
else
{
/* unsymmetric case: F is required. */
Fp = F->p ;
Fi = F->i ;
Fx = F->x ;
Fz = F->z ;
Fnz = F->nz ;
Fpacked = F->packed ;
}
Ap = A->p ; /* size A->ncol+1, column pointers of A */
Ai = A->i ; /* size nz = Ap [A->ncol], row indices of A */
Ax = A->x ; /* size nz, numeric values of A */
Az = A->z ;
Anz = A->nz ;
packed = A->packed ;
sorted = A->sorted ;
use_dbound = IS_GT_ZERO (Common->dbound) ;
/* get the current factors L (and D for LDL'); allocate space if needed */
is_ll = L->is_ll ;
if (L->xtype == CHOLMOD_PATTERN)
{
/* ------------------------------------------------------------------ */
/* L is symbolic only; allocate and initialize L (and D for LDL') */
/* ------------------------------------------------------------------ */
/* workspace: none */
CHOLMOD(change_factor) (A->xtype, is_ll, FALSE, FALSE, TRUE, L, Common);
if (Common->status < CHOLMOD_OK)
{
/* out of memory */
return (FALSE) ;
}
ASSERT (L->minor == (size_t) n) ;
}
else if (kstart == 0 && kend == (size_t) n)
{
/* ------------------------------------------------------------------ */
/* refactorization; reset L->nz and L->minor to restart factorization */
/* ------------------------------------------------------------------ */
L->minor = n ;
Lnz = L->nz ;
for (k = 0 ; k < n ; k++)
{
Lnz [k] = 1 ;
}
}
ASSERT (is_ll == L->is_ll) ;
ASSERT (L->xtype != CHOLMOD_PATTERN) ;
DEBUG (CHOLMOD(dump_factor) (L, "L ready", Common)) ;
DEBUG (CHOLMOD(dump_sparse) (A, "A ready", Common)) ;
DEBUG (if (stype == 0) CHOLMOD(dump_sparse) (F, "F ready", Common)) ;
/* inputs, can be modified on output: */
Lp = L->p ; /* size n+1 */
ASSERT (Lp != NULL) ;
/* outputs, contents defined on input for incremental case only: */
Lnz = L->nz ; /* size n */
Lnext = L->next ; /* size n+2 */
Li = L->i ; /* size L->nzmax, can change in size */
Lx = L->x ; /* size L->nzmax or 2*L->nzmax, can change in size */
Lz = L->z ; /* size L->nzmax for zomplex case, can change in size */
nzmax = L->nzmax ;
ASSERT (Lnz != NULL && Li != NULL && Lx != NULL) ;
/* ---------------------------------------------------------------------- */
/* get workspace */
/* ---------------------------------------------------------------------- */
Iwork = Common->Iwork ;
Stack = Iwork ; /* size n (i/i/l) */
Flag = Common->Flag ; /* size n, Flag [i] < mark must hold */
Wx = Common->Xwork ; /* size n if real, 2*n if complex or
* zomplex. Xwork [i] == 0 must hold. */
Wz = Wx + n ; /* size n for zomplex case only */
mark = Common->mark ;
ASSERT ((Int) Common->xworksize >= (L->xtype == CHOLMOD_REAL ? 1:2)*n) ;
/* ---------------------------------------------------------------------- */
/* compute LDL' or LL' factorization by rows */
/* ---------------------------------------------------------------------- */
#ifdef MASK
#define NEXT(k) k = RLinkUp [k]
#else
#define NEXT(k) k++
#endif
for (k = kstart ; k < ((Int) kend) ; NEXT(k))
{
PRINT1 (("\n===============K "ID" Lnz [k] "ID"\n", k, Lnz [k])) ;
/* ------------------------------------------------------------------ */
/* compute pattern of kth row of L and scatter kth input column */
/* ------------------------------------------------------------------ */
/* column k of L is currently empty */
ASSERT (Lnz [k] == 1) ;
ASSERT (CHOLMOD(dump_work) (TRUE, TRUE, 2*n, Common)) ;
top = n ; /* Stack is empty */
Flag [k] = mark ; /* do not include diagonal entry in Stack */
/* use Li [Lp [i]+1] for etree */
#define PARENT(i) (Lnz [i] > 1) ? (Li [Lp [i] + 1]) : EMPTY
if (stype > 0)
{
/* scatter kth col of triu (beta*I+AA'), get pattern L(k,:) */
p = Ap [k] ;
pend = (packed) ? (Ap [k+1]) : (p + Anz [k]) ;
/* W [i] = Ax [i] ; scatter column of A */
#define SCATTER ASSIGN(Wx,Wz,i, Ax,Az,p)
SUBTREE ;
#undef SCATTER
}
else
{
/* scatter kth col of triu (beta*I+AA'), get pattern L(k,:) */
pf = Fp [k] ;
pfend = (Fpacked) ? (Fp [k+1]) : (pf + Fnz [k]) ;
for ( ; pf < pfend ; pf++)
{
/* get nonzero entry F (t,k) */
t = Fi [pf] ;
/* fk = Fx [pf] */
ASSIGN (fx, fz, 0, Fx, Fz, pf) ;
p = Ap [t] ;
pend = (packed) ? (Ap [t+1]) : (p + Anz [t]) ;
multadds = 0 ;
/* W [i] += Ax [p] * fx ; scatter column of A*A' */
#define SCATTER MULTADD (Wx,Wz,i, Ax,Az,p, fx,fz,0) ; multadds++ ;
SUBTREE ;
#undef SCATTER
#ifdef REAL
fl += 2 * ((double) multadds) ;
#else
fl += 8 * ((double) multadds) ;
#endif
}
}
#undef PARENT
/* ------------------------------------------------------------------ */
/* if mask is present, set the corresponding entries in W to zero */
/* ------------------------------------------------------------------ */
#ifdef MASK
/* remove the dead element of Wx */
if (mask != NULL)
{
#if 0
/* older version */
for (p = n; p > top;)
{
i = Stack [--p] ;
if ( mask [i] >= 0 )
{
CLEAR (Wx,Wz,i) ; /* set W(i) to zero */
}
}
#endif
for (s = top ; s < n ; s++)
{
i = Stack [s] ;
if (mask [i] >= 0)
{
CLEAR (Wx,Wz,i) ; /* set W(i) to zero */
}
}
}
#endif
/* nonzero pattern of kth row of L is now in Stack [top..n-1].
* Flag [Stack [top..n-1]] is equal to mark, but no longer needed */
mark = CHOLMOD(clear_flag) (Common) ;
/* ------------------------------------------------------------------ */
/* compute kth row of L and store in column form */
/* ------------------------------------------------------------------ */
/* Solve L (0:k-1, 0:k-1) * y (0:k-1) = b (0:k-1) where
* b (0:k) = A (0:k,k) or A(0:k,:) * F(:,k) is in W and Stack.
*
* For LDL' factorization:
* L (k, 0:k-1) = y (0:k-1) ./ D (0:k-1)
* D (k) = b (k) - L (k, 0:k-1) * y (0:k-1)
*
* For LL' factorization:
* L (k, 0:k-1) = y (0:k-1)
* L (k,k) = sqrt (b (k) - L (k, 0:k-1) * L (0:k-1, k))
*/
/* dk = W [k] + beta */
ADD_REAL (dk,0, Wx,k, beta,0) ;
#ifndef REAL
/* In the unsymmetric case, the imaginary part of W[k] must be real,
* since F is assumed to be the complex conjugate transpose of A. In
* the symmetric case, W[k] is the diagonal of A. If the imaginary part
* of W[k] is nonzero, then the Cholesky factorization cannot be
* computed; A is not positive definite */
dk_imaginary = (stype > 0) ? (IMAG_IS_NONZERO (Wx,Wz,k)) : FALSE ;
#endif
/* W [k] = 0.0 ; */
CLEAR (Wx,Wz,k) ;
for (s = top ; s < n ; s++)
{
/* get i for each nonzero entry L(k,i) */
i = Stack [s] ;
/* y = W [i] ; */
ASSIGN (yx,yz,0, Wx,Wz,i) ;
/* W [i] = 0.0 ; */
CLEAR (Wx,Wz,i) ;
lnz = Lnz [i] ;
p = Lp [i] ;
ASSERT (lnz > 0 && Li [p] == i) ;
pend = p + lnz ;
/* di = Lx [p] ; the diagonal entry L or D(i,i), which is real */
ASSIGN_REAL (di,0, Lx,p) ;
if (i >= (Int) L->minor || IS_ZERO (di [0]))
{
/* For the LL' factorization, L(i,i) is zero. For the LDL',
* D(i,i) is zero. Skip column i of L, and set L(k,i) = 0. */
CLEAR (lx,lz,0) ;
p = pend ;
}
else if (is_ll)
{
#ifdef REAL
fl += 2 * ((double) (pend - p - 1)) + 3 ;
#else
fl += 8 * ((double) (pend - p - 1)) + 6 ;
#endif
/* forward solve using L (i:(k-1),i) */
/* divide by L(i,i), which must be real and nonzero */
/* y /= di [0] */
DIV_REAL (yx,yz,0, yx,yz,0, di,0) ;
for (p++ ; p < pend ; p++)
{
/* W [Li [p]] -= Lx [p] * y ; */
MULTSUB (Wx,Wz,Li[p], Lx,Lz,p, yx,yz,0) ;
}
/* do not scale L; compute dot product for L(k,k) */
/* L(k,i) = conj(y) ; */
ASSIGN_CONJ (lx,lz,0, yx,yz,0) ;
/* d -= conj(y) * y ; */
LLDOT (dk,0, yx,yz,0) ;
}
else
{
#ifdef REAL
fl += 2 * ((double) (pend - p - 1)) + 3 ;
#else
fl += 8 * ((double) (pend - p - 1)) + 6 ;
#endif
/* forward solve using D (i,i) and L ((i+1):(k-1),i) */
for (p++ ; p < pend ; p++)
{
/* W [Li [p]] -= Lx [p] * y ; */
MULTSUB (Wx,Wz,Li[p], Lx,Lz,p, yx,yz,0) ;
}
/* Scale L (k,0:k-1) for LDL' factorization, compute D (k,k)*/
#ifdef REAL
/* L(k,i) = y/d */
lx [0] = yx [0] / di [0] ;
/* d -= L(k,i) * y */
dk [0] -= lx [0] * yx [0] ;
#else
/* L(k,i) = conj(y) ; */
ASSIGN_CONJ (lx,lz,0, yx,yz,0) ;
/* L(k,i) /= di ; */
DIV_REAL (lx,lz,0, lx,lz,0, di,0) ;
/* d -= conj(y) * y / di */
LDLDOT (dk,0, yx,yz,0, di,0) ;
#endif
}
/* determine if column i of L can hold the new L(k,i) entry */
if (p >= Lp [Lnext [i]])
{
/* column i needs to grow */
PRINT1 (("Factor Colrealloc "ID", old Lnz "ID"\n", i, Lnz [i]));
if (!CHOLMOD(reallocate_column) (i, lnz + 1, L, Common))
{
/* out of memory, L is now simplicial symbolic */
for (i = 0 ; i < n ; i++)
{
/* W [i] = 0 ; */
CLEAR (Wx,Wz,i) ;
}
ASSERT (CHOLMOD(dump_work) (TRUE, TRUE, n, Common)) ;
return (FALSE) ;
}
Li = L->i ; /* L->i, L->x, L->z may have moved */
Lx = L->x ;
Lz = L->z ;
p = Lp [i] + lnz ; /* contents of L->p changed */
ASSERT (p < Lp [Lnext [i]]) ;
}
/* store L (k,i) in the column form matrix of L */
Li [p] = k ;
/* Lx [p] = L(k,i) ; */
ASSIGN (Lx,Lz,p, lx,lz,0) ;
Lnz [i]++ ;
}
/* ------------------------------------------------------------------ */
/* ensure abs (d) >= dbound if dbound is given, and store it in L */
/* ------------------------------------------------------------------ */
p = Lp [k] ;
Li [p] = k ;
if (k >= (Int) L->minor)
{
/* the matrix is already not positive definite */
dk [0] = 0 ;
}
else if (use_dbound)
{
/* modify the diagonal to force LL' or LDL' to exist */
dk [0] = CHOLMOD(dbound) (is_ll ? fabs (dk [0]) : dk [0], Common) ;
}
else if ((is_ll ? (IS_LE_ZERO (dk [0])) : (IS_ZERO (dk [0])))
#ifndef REAL
|| dk_imaginary
#endif
)
{
/* the matrix has just been found to be not positive definite */
dk [0] = 0 ;
L->minor = k ;
ERROR (CHOLMOD_NOT_POSDEF, "not positive definite") ;
}
if (is_ll)
{
/* this is counted as one flop, below */
dk [0] = sqrt (dk [0]) ;
}
/* Lx [p] = D(k,k) = d ; real part only */
ASSIGN_REAL (Lx,p, dk,0) ;
CLEAR_IMAG (Lx,Lz,p) ;
}
#undef NEXT
if (is_ll) fl += MAX ((Int) kend - (Int) kstart, 0) ; /* count sqrt's */
Common->rowfacfl = fl ;
DEBUG (CHOLMOD(dump_factor) (L, "final cholmod_rowfac", Common)) ;
ASSERT (CHOLMOD(dump_work) (TRUE, TRUE, n, Common)) ;
return (TRUE) ;
}
#undef PATTERN
#undef REAL
#undef COMPLEX
#undef ZOMPLEX
syntax highlighted by Code2HTML, v. 0.9.1