/* ========================================================================== */
/* === MatrixOps/cholmod_symmetry =========================================== */
/* ========================================================================== */
/* -----------------------------------------------------------------------------
* CHOLMOD/MatrixOps Module. Copyright (C) 2005-2006, Timothy A. Davis
* The CHOLMOD/MatrixOps Module is licensed under Version 2.0 of the GNU
* General Public License. See gpl.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
* -------------------------------------------------------------------------- */
/* Determines if a sparse matrix is rectangular, unsymmetric, symmetric,
* skew-symmetric, or Hermitian. It does so by looking at its numerical values
* of both upper and lower triangular parts of a CHOLMOD "unsymmetric"
* matrix, where A->stype == 0. The transpose of A is NOT constructed.
*
* If not unsymmetric, it also determines if the matrix has a diagonal whose
* entries are all real and positive (and thus a candidate for sparse Cholesky
* if A->stype is changed to a nonzero value).
*
* Note that a Matrix Market "general" matrix is either rectangular or
* unsymmetric.
*
* The row indices in the column of each matrix MUST be sorted for this function
* to work properly (A->sorted must be TRUE). This routine returns EMPTY if
* A->stype is not zero, or if A->sorted is FALSE. The exception to this rule
* is if A is rectangular.
*
* If option == 0, then this routine returns immediately when it finds a
* non-positive diagonal entry (or one with nonzero imaginary part). If the
* matrix is not a candidate for sparse Cholesky, it returns the value
* CHOLMOD_MM_UNSYMMETRIC, even if the matrix might in fact be symmetric or
* Hermitian.
*
* This routine is useful inside the MATLAB backslash, which must look at an
* arbitrary matrix (A->stype == 0) and determine if it is a candidate for
* sparse Cholesky. In that case, option should be 0.
*
* This routine is also useful when writing a MATLAB matrix to a file in
* Rutherford/Boeing or Matrix Market format. Those formats require a
* determination as to the symmetry of the matrix, and thus this routine should
* not return upon encountering the first non-positive diagonal. In this case,
* option should be 1.
*
* If option is 2, this function can be used to compute the numerical and
* pattern symmetry, where 0 is a completely unsymmetric matrix, and 1 is a
* perfectly symmetric matrix. This option is used when computing the following
* statistics for the matrices in the UF Sparse Matrix Collection.
*
* numerical symmetry: number of matched offdiagonal nonzeros over
* the total number of offdiagonal entries. A real entry A(i,j), i ~= j,
* is matched if A (j,i) == A (i,j), but this is only counted if both
* A(j,i) and A(i,j) are nonzero. This does not depend on Z.
* (If A is complex, then the above test is modified; A (i,j) is matched
* if conj (A (j,i)) == A (i,j)).
*
* Then numeric symmetry = xmatched / nzoffdiag, or 1 if nzoffdiag = 0.
*
* pattern symmetry: number of matched offdiagonal entries over the
* total number of offdiagonal entries. An entry A(i,j), i ~= j, is
* matched if A (j,i) is also an entry.
*
* Then pattern symmetry = pmatched / nzoffdiag, or 1 if nzoffdiag = 0.
*
* The symmetry of a matrix with no offdiagonal entries is equal to 1.
*
* A workspace of size ncol integers is allocated; EMPTY is returned if this
* allocation fails.
*
* Summary of return values:
*
* EMPTY (-1) out of memory, stype not zero, A not sorted
* CHOLMOD_MM_RECTANGULAR 1 A is rectangular
* CHOLMOD_MM_UNSYMMETRIC 2 A is unsymmetric
* CHOLMOD_MM_SYMMETRIC 3 A is symmetric, but with non-pos. diagonal
* CHOLMOD_MM_HERMITIAN 4 A is Hermitian, but with non-pos. diagonal
* CHOLMOD_MM_SKEW_SYMMETRIC 5 A is skew symmetric
* CHOLMOD_MM_SYMMETRIC_POSDIAG 6 A is symmetric with positive diagonal
* CHOLMOD_MM_HERMITIAN_POSDIAG 7 A is Hermitian with positive diagonal
*
* See also the spsym mexFunction, which is a MATLAB interface for this code.
*
* If the matrix is a candidate for sparse Cholesky, it will return a result
* CHOLMOD_MM_SYMMETRIC_POSDIAG if real, or CHOLMOD_MM_HERMITIAN_POSDIAG if
* complex. Otherwise, it will return a value less than this. This is true
* regardless of the value of the option parameter.
*/
#ifndef NMATRIXOPS
#include "cholmod_internal.h"
#include "cholmod_matrixops.h"
/* ========================================================================== */
/* === get_value ============================================================ */
/* ========================================================================== */
/* Get the pth value in the matrix. */
static void get_value
(
double *Ax, /* real values, or real/imag. for CHOLMOD_COMPLEX type */
double *Az, /* imaginary values for CHOLMOD_ZOMPLEX type */
Int p, /* get the pth entry */
Int xtype, /* A->xtype: pattern, real, complex, or zomplex */
double *x, /* the real part */
double *z /* the imaginary part */
)
{
switch (xtype)
{
case CHOLMOD_PATTERN:
*x = 1 ;
*z = 0 ;
break ;
case CHOLMOD_REAL:
*x = Ax [p] ;
*z = 0 ;
break ;
case CHOLMOD_COMPLEX:
*x = Ax [2*p] ;
*z = Ax [2*p+1] ;
break ;
case CHOLMOD_ZOMPLEX:
*x = Ax [p] ;
*z = Az [p] ;
break ;
}
}
/* ========================================================================== */
/* === cholmod_symmetry ===================================================== */
/* ========================================================================== */
/* Determine the symmetry of a matrix, and check its diagonal.
*
* option 0: Do not count # of matched pairs. Quick return if the
* the matrix has a zero, negative, or imaginary diagonal entry.
*
* option 1: Do not count # of matched pairs. Do not return quickly if
* the matrix has a zero, negative, or imaginary diagonal entry.
* The result 1 to 7 is accurately computed:
*
* EMPTY (-1) out of memory, stype not zero, A not sorted
* CHOLMOD_MM_RECTANGULAR 1 A is rectangular
* CHOLMOD_MM_UNSYMMETRIC 2 A is unsymmetric
* CHOLMOD_MM_SYMMETRIC 3 A is symmetric, with non-pos. diagonal
* CHOLMOD_MM_HERMITIAN 4 A is Hermitian, with non-pos. diagonal
* CHOLMOD_MM_SKEW_SYMMETRIC 5 A is skew symmetric
* CHOLMOD_MM_SYMMETRIC_POSDIAG 6 is symmetric with positive diagonal
* CHOLMOD_MM_HERMITIAN_POSDIAG 7 A is Hermitian with positive diagonal
*
* The routine returns as soon as the above is determined (that is, it
* can return as soon as it determines the matrix is unsymmetric).
*
* option 2: All of the above, but also compute the number of matched off-
* diagonal entries (of two types). xmatched is the number of
* nonzero entries for which A(i,j) = conj(A(j,i)). pmatched is
* the number of entries (i,j) for which A(i,j) and A(j,i) are both in
* the pattern of A (the value doesn't matter). nzoffdiag is the total
* number of off-diagonal entries in the pattern. nzdiag is the number of
* diagonal entries in the pattern.
*
* With option 0 or 1, or if the matrix is rectangular, xmatched, pmatched,
* nzoffdiag, and nzdiag are not computed.
*
* Note that a matched pair, A(i,j) and A(j,i) for i != j, is counted twice
* (once per entry).
*/
int CHOLMOD(symmetry)
(
/* ---- input ---- */
cholmod_sparse *A,
int option, /* option 0, 1, or 2 (see above) */
/* ---- output --- */ /* outputs ignored if any are NULL */
Int *p_xmatched, /* # of matched numerical entries */
Int *p_pmatched, /* # of matched entries in pattern */
Int *p_nzoffdiag, /* # of off diagonal entries */
Int *p_nzdiag, /* # of diagonal entries */
/* --------------- */
cholmod_common *Common
)
{
double aij_real, aij_imag, aji_real, aji_imag ;
double *Ax, *Az ;
Int *Ap, *Ai, *Anz, *munch ;
Int packed, nrow, ncol, xtype, is_symmetric, is_skew, is_hermitian, posdiag,
j, p, pend, i, piend, result, xmatched, pmatched, nzdiag, i2, found ;
/* ---------------------------------------------------------------------- */
/* check inputs */
/* ---------------------------------------------------------------------- */
RETURN_IF_NULL_COMMON (EMPTY) ;
RETURN_IF_NULL (A, EMPTY) ;
RETURN_IF_XTYPE_INVALID (A, CHOLMOD_PATTERN, CHOLMOD_ZOMPLEX, EMPTY) ;
Common->status = CHOLMOD_OK ;
ASSERT (CHOLMOD(dump_sparse) (A, "cholmod_symmetry", Common) >= 0) ;
if (p_xmatched == NULL || p_pmatched == NULL
|| p_nzoffdiag == NULL || p_nzdiag == NULL)
{
/* option 2 is not performed if any output parameter is NULL */
option = MAX (option, 1) ;
}
/* ---------------------------------------------------------------------- */
/* get inputs */
/* ---------------------------------------------------------------------- */
Ap = A->p ;
Ai = A->i ;
Ax = A->x ;
Az = A->z ;
Anz = A->nz ;
packed = A->packed ;
ncol = A->ncol ;
nrow = A->nrow ;
xtype = A->xtype ;
/* ---------------------------------------------------------------------- */
/* check if rectangular, unsorted, or stype is not zero */
/* ---------------------------------------------------------------------- */
if (nrow != ncol)
{
/* matrix is rectangular */
return (CHOLMOD_MM_RECTANGULAR) ;
}
if (!(A->sorted) || A->stype != 0)
{
/* this function cannot determine the type or symmetry */
return (EMPTY) ;
}
/* ---------------------------------------------------------------------- */
/* allocate workspace */
/* ---------------------------------------------------------------------- */
/* this function requires uninitialized Int workspace of size ncol */
CHOLMOD(allocate_work) (0, ncol, 0, Common) ;
if (Common->status < CHOLMOD_OK)
{
/* out of memory */
return (EMPTY) ;
}
munch = Common->Iwork ; /* the munch array is size ncol */
/* ---------------------------------------------------------------------- */
/* determine symmetry of a square matrix */
/* ---------------------------------------------------------------------- */
/* a complex or zomplex matrix is Hermitian until proven otherwise */
is_hermitian = (xtype >= CHOLMOD_COMPLEX) ;
/* any matrix is symmetric until proven otherwise */
is_symmetric = TRUE ;
/* a non-pattern matrix is skew-symmetric until proven otherwise */
is_skew = (xtype != CHOLMOD_PATTERN) ;
/* a matrix has positive diagonal entries until proven otherwise */
posdiag = TRUE ;
/* munch pointers start at the top of each column */
for (j = 0 ; j < ncol ; j++)
{
munch [j] = Ap [j] ;
}
xmatched = 0 ;
pmatched = 0 ;
nzdiag = 0 ;
for (j = 0 ; j < ncol ; j++) /* examine each column of A */
{
/* ------------------------------------------------------------------ */
/* look at the entire munch column j */
/* ------------------------------------------------------------------ */
/* start at the munch point of column j, and go to end of the column */
p = munch [j] ;
pend = (packed) ? (Ap [j+1]) : (Ap [j] + Anz [j]) ;
for ( ; p < pend ; p++)
{
/* get the row index of A(i,j) */
i = Ai [p] ;
if (i < j)
{
/* ---------------------------------------------------------- */
/* A(i,j) in triu(A), but matching A(j,i) not in tril(A) */
/* ---------------------------------------------------------- */
/* entry A(i,j) is unmatched; it appears in the upper triangular
* part, but not the lower triangular part. The matrix is
* unsymmetric. */
is_hermitian = FALSE ;
is_symmetric = FALSE ;
is_skew = FALSE ;
}
else if (i == j)
{
/* ---------------------------------------------------------- */
/* the diagonal A(j,j) is present; check its value */
/* ---------------------------------------------------------- */
get_value (Ax, Az, p, xtype, &aij_real, &aij_imag) ;
if (aij_real != 0. || aij_imag != 0.)
{
/* diagonal is nonzero; matrix is not skew-symmetric */
nzdiag++ ;
is_skew = FALSE ;
}
if (aij_real <= 0. || aij_imag != 0.)
{
/* diagonal negative or imaginary; not chol candidate */
posdiag = FALSE ;
}
}
else /* i > j */
{
/* ---------------------------------------------------------- */
/* consider column i, up to and including row j */
/* ---------------------------------------------------------- */
/* munch the entry at top of column i up to and incl row j */
piend = (packed) ? (Ap [i+1]) : (Ap [i] + Anz [i]) ;
found = FALSE ;
for ( ; munch [i] < piend ; munch [i]++)
{
i2 = Ai [munch [i]] ;
if (i2 < j)
{
/* -------------------------------------------------- */
/* A(i2,i) in triu(A) but A(i,i2) not in tril(A) */
/* -------------------------------------------------- */
/* The matrix is unsymmetric. */
is_hermitian = FALSE ;
is_symmetric = FALSE ;
is_skew = FALSE ;
}
else if (i2 == j)
{
/* -------------------------------------------------- */
/* both A(i,j) and A(j,i) exist in the matrix */
/* -------------------------------------------------- */
/* this is one more matching entry in the pattern */
pmatched += 2 ;
found = TRUE ;
/* get the value of A(i,j) */
get_value (Ax, Az, p, xtype, &aij_real, &aij_imag) ;
/* get the value of A(j,i) */
get_value (Ax, Az, munch [i],
xtype, &aji_real, &aji_imag) ;
/* compare A(i,j) with A(j,i) */
if (aij_real != aji_real || aij_imag != aji_imag)
{
/* the matrix cannot be symmetric */
is_symmetric = FALSE ;
}
if (aij_real != -aji_real || aij_imag != aji_imag)
{
/* the matrix cannot be skew-symmetric */
is_skew = FALSE ;
}
if (aij_real != aji_real || aij_imag != -aji_imag)
{
/* the matrix cannot be Hermitian */
is_hermitian = FALSE ;
}
else
{
/* A(i,j) and A(j,i) are numerically matched */
xmatched += 2 ;
}
}
else /* i2 > j */
{
/* -------------------------------------------------- */
/* entry A(i2,i) is not munched; consider it later */
/* -------------------------------------------------- */
break ;
}
}
if (!found)
{
/* A(i,j) in tril(A) but A(j,i) not in triu(A).
* The matrix is unsymmetric. */
is_hermitian = FALSE ;
is_symmetric = FALSE ;
is_skew = FALSE ;
}
}
if (option < 2 && !(is_symmetric || is_skew || is_hermitian))
{
/* matrix is unsymmetric; terminate the test */
return (CHOLMOD_MM_UNSYMMETRIC) ;
}
}
/* ------------------------------------------------------------------ */
/* quick return if not Cholesky candidate */
/* ------------------------------------------------------------------ */
if (option < 1 && (!posdiag || nzdiag < ncol))
{
/* Diagonal entry not present, or present but negative or with
* nonzero imaginary part. Quick return for option 0. */
return (CHOLMOD_MM_UNSYMMETRIC) ;
}
}
/* ---------------------------------------------------------------------- */
/* return the results */
/* ---------------------------------------------------------------------- */
if (option >= 2)
{
*p_xmatched = xmatched ;
*p_pmatched = pmatched ;
*p_nzoffdiag = CHOLMOD(nnz) (A, Common) - nzdiag ;
*p_nzdiag = nzdiag ;
}
result = CHOLMOD_MM_UNSYMMETRIC ;
if (is_hermitian)
{
/* complex Hermitian matrix, with either pos. or non-pos. diagonal */
result = posdiag ? CHOLMOD_MM_HERMITIAN_POSDIAG : CHOLMOD_MM_HERMITIAN ;
}
else if (is_symmetric)
{
/* real or complex symmetric matrix, with pos. or non-pos. diagonal */
result = posdiag ? CHOLMOD_MM_SYMMETRIC_POSDIAG : CHOLMOD_MM_SYMMETRIC ;
}
else if (is_skew)
{
/* real or complex skew-symmetric matrix */
result = CHOLMOD_MM_SKEW_SYMMETRIC ;
}
return (result) ;
}
#endif
syntax highlighted by Code2HTML, v. 0.9.1