amesos_t_cholmod_rowfac.c

Go to the documentation of this file.

/* ========================================================================== */
/* === 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 "amesos_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