amesos_cholmod_aat.c

Go to the documentation of this file.

/* ========================================================================== */
/* === Core/cholmod_aat ===================================================== */
/* ========================================================================== */

/* -----------------------------------------------------------------------------
 * CHOLMOD/Core Module.  Copyright (C) 2005-2006,
 * Univ. of Florida.  Author: Timothy A. Davis
 * The CHOLMOD/Core 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
 * -------------------------------------------------------------------------- */

/* C = A*A' or C = A(:,f)*A(:,f)'
 *
 * A can be packed or unpacked, sorted or unsorted, but must be stored with
 * both upper and lower parts (A->stype of zero).  C is returned as packed,
 * C->stype of zero (both upper and lower parts present), and unsorted.  See
 * cholmod_ssmult in the MatrixOps Module for a more general matrix-matrix
 * multiply.
 *
 * You can trivially convert C into a symmetric upper/lower matrix by
 * changing C->stype = 1 or -1 after calling this routine.
 *
 * workspace:
 *  Flag (A->nrow),
 *  Iwork (max (A->nrow, A->ncol)) if fset present,
 *  Iwork (A->nrow) if no fset,
 *  W (A->nrow) if mode > 0,
 *  allocates temporary copy for A'.
 *
 * A can be pattern or real.  Complex or zomplex cases are supported only
 * if the mode is <= 0 (in which case the numerical values are ignored).
 */

#include "amesos_cholmod_internal.h"
#include "amesos_cholmod_core.h"

cholmod_sparse *CHOLMOD(aat)
(
    /* ---- input ---- */
    cholmod_sparse *A,  /* input matrix; C=A*A' is constructed */
    Int *fset,    /* subset of 0:(A->ncol)-1 */
    size_t fsize, /* size of fset */
    int mode,   /* >0: numerical, 0: pattern, <0: pattern (no diag)
       * -2: pattern only, no diagonal, add 50% + n extra
       * space to C */
    /* --------------- */
    cholmod_common *Common
)
{
    double fjt ;
    double *Ax, *Fx, *Cx, *W ;
    Int *Ap, *Anz, *Ai, *Fp, *Fi, *Cp, *Ci, *Flag ;
    cholmod_sparse *C, *F ;
    Int packed, j, i, pa, paend, pf, pfend, n, mark, cnz, t, p, values, diag,
  extra ;

    /* ---------------------------------------------------------------------- */
    /* check inputs */
    /* ---------------------------------------------------------------------- */

    RETURN_IF_NULL_COMMON (NULL) ;
    RETURN_IF_NULL (A, NULL) ;
    values = (mode > 0) && (A->xtype != CHOLMOD_PATTERN) ;
    RETURN_IF_XTYPE_INVALID (A, CHOLMOD_PATTERN,
      values ? CHOLMOD_REAL : CHOLMOD_ZOMPLEX, NULL) ;
    if (A->stype)
    {
  ERROR (CHOLMOD_INVALID, "matrix cannot be symmetric") ;
  return (NULL) ;
    }
    Common->status = CHOLMOD_OK ;

    /* ---------------------------------------------------------------------- */
    /* allocate workspace */
    /* ---------------------------------------------------------------------- */

    diag = (mode >= 0) ;
    n = A->nrow ;
    CHOLMOD(allocate_work) (n, MAX (A->ncol, A->nrow), values ? n : 0, Common) ;
    if (Common->status < CHOLMOD_OK)
    {
  return (NULL) ;     /* out of memory */
    }
    ASSERT (CHOLMOD(dump_work) (TRUE, TRUE, values ? n : 0, Common)) ;

    /* ---------------------------------------------------------------------- */
    /* get inputs */
    /* ---------------------------------------------------------------------- */

    ASSERT (CHOLMOD(dump_sparse) (A, "A", Common) >= 0) ;

    /* get the A matrix */
    Ap  = A->p ;
    Anz = A->nz ;
    Ai  = A->i ;
    Ax  = A->x ;
    packed = A->packed ;

    /* get workspace */
    W = Common->Xwork ;   /* size n, unused if values is FALSE */
    Flag = Common->Flag ; /* size n, Flag [0..n-1] < mark on input*/

    /* ---------------------------------------------------------------------- */
    /* F = A' or A(:,f)' */
    /* ---------------------------------------------------------------------- */

    /* workspace: Iwork (nrow if no fset; MAX (nrow,ncol) if fset)*/
    F = CHOLMOD(ptranspose) (A, values, NULL, fset, fsize, Common) ;
    if (Common->status < CHOLMOD_OK)
    {
  return (NULL) ;     /* out of memory */
    }

    Fp = F->p ;
    Fi = F->i ;
    Fx = F->x ;

    /* ---------------------------------------------------------------------- */
    /* count the number of entries in the result C */
    /* ---------------------------------------------------------------------- */

    cnz = 0 ;
    for (j = 0 ; j < n ; j++)
    {
  /* clear the Flag array */
  mark = CHOLMOD(clear_flag) (Common) ;

  /* exclude the diagonal, if requested */
  if (!diag)
  {
      Flag [j] = mark ;
  }

  /* for each nonzero F(t,j) in column j, do: */
  pfend = Fp [j+1] ;
  for (pf = Fp [j] ; pf < pfend ; pf++)
  {
      /* F(t,j) is nonzero */
      t = Fi [pf] ;

      /* add the nonzero pattern of A(:,t) to the pattern of C(:,j) */
      pa = Ap [t] ;
      paend = (packed) ? (Ap [t+1]) : (pa + Anz [t]) ;
      for ( ; pa < paend ; pa++)
      {
    i = Ai [pa] ;
    if (Flag [i] != mark)
    {
        Flag [i] = mark ;
        cnz++ ;
    }
      }
  }
  if (cnz < 0)
  {
      break ;     /* integer overflow case */
  }
    }

    extra = (mode == -2) ? (cnz/2 + n) : 0 ;

    mark = CHOLMOD(clear_flag) (Common) ;

    /* ---------------------------------------------------------------------- */
    /* check for integer overflow */
    /* ---------------------------------------------------------------------- */

    if (cnz < 0 || (cnz + extra) < 0)
    {
  ERROR (CHOLMOD_TOO_LARGE, "problem too large") ;
  CHOLMOD(clear_flag) (Common) ;
  CHOLMOD(free_sparse) (&F, Common) ;
  return (NULL) ;     /* problem too large */
    }

    /* ---------------------------------------------------------------------- */
    /* allocate C */
    /* ---------------------------------------------------------------------- */

    C = CHOLMOD(allocate_sparse) (n, n, cnz + extra, FALSE, TRUE, 0,
      values ? A->xtype : CHOLMOD_PATTERN, Common) ;
    if (Common->status < CHOLMOD_OK)
    {
  CHOLMOD(free_sparse) (&F, Common) ;
  return (NULL) ;     /* out of memory */
    }

    Cp = C->p ;
    Ci = C->i ;
    Cx = C->x ;

    /* ---------------------------------------------------------------------- */
    /* C = A*A' */
    /* ---------------------------------------------------------------------- */

    cnz = 0 ;

    if (values)
    {

  /* pattern and values */
  for (j = 0 ; j < n ; j++)
  {
      /* clear the Flag array */
      mark = CHOLMOD(clear_flag) (Common) ;

      /* start column j of C */
      Cp [j] = cnz ;

      /* for each nonzero F(t,j) in column j, do: */
      pfend = Fp [j+1] ;
      for (pf = Fp [j] ; pf < pfend ; pf++)
      {
    /* F(t,j) is nonzero */
    t = Fi [pf] ;
    fjt = Fx [pf] ;

    /* add the nonzero pattern of A(:,t) to the pattern of C(:,j)
     * and scatter the values into W */
    pa = Ap [t] ;
    paend = (packed) ? (Ap [t+1]) : (pa + Anz [t]) ;
    for ( ; pa < paend ; pa++)
    {
        i = Ai [pa] ;
        if (Flag [i] != mark)
        {
      Flag [i] = mark ;
      Ci [cnz++] = i ;
        }
        W [i] += Ax [pa] * fjt ;
    }
      }

      /* gather the values into C(:,j) */
      for (p = Cp [j] ; p < cnz ; p++)
      {
    i = Ci [p] ;
    Cx [p] = W [i] ;
    W [i] = 0 ;
      }
  }

    }
    else
    {

  /* pattern only */
  for (j = 0 ; j < n ; j++)
  {
      /* clear the Flag array */
      mark = CHOLMOD(clear_flag) (Common) ;

      /* exclude the diagonal, if requested */
      if (!diag)
      {
    Flag [j] = mark ;
      }

      /* start column j of C */
      Cp [j] = cnz ;

      /* for each nonzero F(t,j) in column j, do: */
      pfend = Fp [j+1] ;
      for (pf = Fp [j] ; pf < pfend ; pf++)
      {
    /* F(t,j) is nonzero */
    t = Fi [pf] ;

    /* add the nonzero pattern of A(:,t) to the pattern of C(:,j) */
    pa = Ap [t] ;
    paend = (packed) ? (Ap [t+1]) : (pa + Anz [t]) ;
    for ( ; pa < paend ; pa++)
    {
        i = Ai [pa] ;
        if (Flag [i] != mark)
        {
      Flag [i] = mark ;
      Ci [cnz++] = i ;
        }
    }
      }
  }
    }

    Cp [n] = cnz ;
    ASSERT (IMPLIES (mode != -2, MAX (1,cnz) == C->nzmax)) ;

    /* ---------------------------------------------------------------------- */
    /* clear workspace and free temporary matrices and return result */
    /* ---------------------------------------------------------------------- */

    CHOLMOD(free_sparse) (&F, Common) ;
    CHOLMOD(clear_flag) (Common) ;
    ASSERT (CHOLMOD(dump_work) (TRUE, TRUE, values ? n : 0, Common)) ;
    DEBUG (i = CHOLMOD(dump_sparse) (C, "aat", Common)) ;
    ASSERT (IMPLIES (mode < 0, i == 0)) ;
    return (C) ;
}