/* ========================================================================== */
/* === MatrixOps/cholmod_norm =============================================== */
/* ========================================================================== */
/* -----------------------------------------------------------------------------
* 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
* -------------------------------------------------------------------------- */
/* r = norm (A), compute the infinity-norm, 1-norm, or 2-norm of a sparse or
* dense matrix. Can compute the 2-norm only for a dense column vector.
* Returns -1 if an error occurs.
*
* Pattern, real, complex, and zomplex sparse matrices are supported.
*/
#ifndef NMATRIXOPS
#include "cholmod_internal.h"
#include "cholmod_matrixops.h"
/* ========================================================================== */
/* === abs_value ============================================================ */
/* ========================================================================== */
/* Compute the absolute value of a real, complex, or zomplex value */
static double abs_value
(
int xtype,
double *Ax,
double *Az,
Int p,
cholmod_common *Common
)
{
double s = 0 ;
switch (xtype)
{
case CHOLMOD_PATTERN:
s = 1 ;
break ;
case CHOLMOD_REAL:
s = fabs (Ax [p]) ;
break ;
case CHOLMOD_COMPLEX:
s = Common->hypotenuse (Ax [2*p], Ax [2*p+1]) ;
break ;
case CHOLMOD_ZOMPLEX:
s = Common->hypotenuse (Ax [p], Az [p]) ;
break ;
}
return (s) ;
}
/* ========================================================================== */
/* === cholmod_norm_dense =================================================== */
/* ========================================================================== */
double CHOLMOD(norm_dense)
(
/* ---- input ---- */
cholmod_dense *X, /* matrix to compute the norm of */
int norm, /* type of norm: 0: inf. norm, 1: 1-norm, 2: 2-norm */
/* --------------- */
cholmod_common *Common
)
{
double xnorm, s, x, z ;
double *Xx, *Xz, *W ;
Int nrow, ncol, d, i, j, use_workspace, xtype ;
/* ---------------------------------------------------------------------- */
/* check inputs */
/* ---------------------------------------------------------------------- */
RETURN_IF_NULL_COMMON (EMPTY) ;
RETURN_IF_NULL (X, EMPTY) ;
RETURN_IF_XTYPE_INVALID (X, CHOLMOD_REAL, CHOLMOD_ZOMPLEX, EMPTY) ;
Common->status = CHOLMOD_OK ;
ncol = X->ncol ;
if (norm < 0 || norm > 2 || (norm == 2 && ncol > 1))
{
ERROR (CHOLMOD_INVALID, "invalid norm") ;
return (EMPTY) ;
}
/* ---------------------------------------------------------------------- */
/* get inputs */
/* ---------------------------------------------------------------------- */
nrow = X->nrow ;
d = X->d ;
Xx = X->x ;
Xz = X->z ;
xtype = X->xtype ;
/* ---------------------------------------------------------------------- */
/* allocate workspace, if needed */
/* ---------------------------------------------------------------------- */
W = NULL ;
use_workspace = (norm == 0 && ncol > 4) ;
if (use_workspace)
{
CHOLMOD(allocate_work) (0, 0, nrow, Common) ;
W = Common->Xwork ;
if (Common->status < CHOLMOD_OK)
{
/* oops, no workspace */
use_workspace = FALSE ;
}
}
/* ---------------------------------------------------------------------- */
/* compute the norm */
/* ---------------------------------------------------------------------- */
xnorm = 0 ;
if (use_workspace)
{
/* ------------------------------------------------------------------ */
/* infinity-norm = max row sum, using stride-1 access of X */
/* ------------------------------------------------------------------ */
DEBUG (for (i = 0 ; i < nrow ; i++) ASSERT (W [i] == 0)) ;
/* this is faster than stride-d, but requires O(nrow) workspace */
for (j = 0 ; j < ncol ; j++)
{
for (i = 0 ; i < nrow ; i++)
{
W [i] += abs_value (xtype, Xx, Xz, i+j*d, Common) ;
}
}
for (i = 0 ; i < nrow ; i++)
{
s = W [i] ;
if ((IS_NAN (s) || s > xnorm) && !IS_NAN (xnorm))
{
xnorm = s ;
}
W [i] = 0 ;
}
}
else if (norm == 0)
{
/* ------------------------------------------------------------------ */
/* infinity-norm = max row sum, using stride-d access of X */
/* ------------------------------------------------------------------ */
for (i = 0 ; i < nrow ; i++)
{
s = 0 ;
for (j = 0 ; j < ncol ; j++)
{
s += abs_value (xtype, Xx, Xz, i+j*d, Common) ;
}
if ((IS_NAN (s) || s > xnorm) && !IS_NAN (xnorm))
{
xnorm = s ;
}
}
}
else if (norm == 1)
{
/* ------------------------------------------------------------------ */
/* 1-norm = max column sum */
/* ------------------------------------------------------------------ */
for (j = 0 ; j < ncol ; j++)
{
s = 0 ;
for (i = 0 ; i < nrow ; i++)
{
s += abs_value (xtype, Xx, Xz, i+j*d, Common) ;
}
if ((IS_NAN (s) || s > xnorm) && !IS_NAN (xnorm))
{
xnorm = s ;
}
}
}
else
{
/* ------------------------------------------------------------------ */
/* 2-norm = sqrt (sum (X.^2)) */
/* ------------------------------------------------------------------ */
switch (xtype)
{
case CHOLMOD_REAL:
for (i = 0 ; i < nrow ; i++)
{
x = Xx [i] ;
xnorm += x*x ;
}
break ;
case CHOLMOD_COMPLEX:
for (i = 0 ; i < nrow ; i++)
{
x = Xx [2*i ] ;
z = Xx [2*i+1] ;
xnorm += x*x + z*z ;
}
break ;
case CHOLMOD_ZOMPLEX:
for (i = 0 ; i < nrow ; i++)
{
x = Xx [i] ;
z = Xz [i] ;
xnorm += x*x + z*z ;
}
break ;
}
xnorm = sqrt (xnorm) ;
}
/* ---------------------------------------------------------------------- */
/* return result */
/* ---------------------------------------------------------------------- */
return (xnorm) ;
}
/* ========================================================================== */
/* === cholmod_norm_sparse ================================================== */
/* ========================================================================== */
double CHOLMOD(norm_sparse)
(
/* ---- input ---- */
cholmod_sparse *A, /* matrix to compute the norm of */
int norm, /* type of norm: 0: inf. norm, 1: 1-norm */
/* --------------- */
cholmod_common *Common
)
{
double anorm, s ;
double *Ax, *Az, *W ;
Int *Ap, *Ai, *Anz ;
Int i, j, p, pend, nrow, ncol, packed, xtype ;
/* ---------------------------------------------------------------------- */
/* 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 ;
ncol = A->ncol ;
nrow = A->nrow ;
if (norm < 0 || norm > 1)
{
ERROR (CHOLMOD_INVALID, "invalid norm") ;
return (EMPTY) ;
}
if (A->stype && nrow != ncol)
{
ERROR (CHOLMOD_INVALID, "matrix invalid") ;
return (EMPTY) ;
}
/* ---------------------------------------------------------------------- */
/* get inputs */
/* ---------------------------------------------------------------------- */
Ap = A->p ;
Ai = A->i ;
Ax = A->x ;
Az = A->z ;
Anz = A->nz ;
packed = A->packed ;
xtype = A->xtype ;
/* ---------------------------------------------------------------------- */
/* allocate workspace, if needed */
/* ---------------------------------------------------------------------- */
W = NULL ;
if (A->stype || norm == 0)
{
CHOLMOD(allocate_work) (0, 0, nrow, Common) ;
W = Common->Xwork ;
if (Common->status < CHOLMOD_OK)
{
/* out of memory */
return (EMPTY) ;
}
DEBUG (for (i = 0 ; i < nrow ; i++) ASSERT (W [i] == 0)) ;
}
/* ---------------------------------------------------------------------- */
/* compute the norm */
/* ---------------------------------------------------------------------- */
anorm = 0 ;
if (A->stype > 0)
{
/* ------------------------------------------------------------------ */
/* A is symmetric with upper triangular part stored */
/* ------------------------------------------------------------------ */
/* infinity-norm = 1-norm = max row/col sum */
for (j = 0 ; j < ncol ; j++)
{
p = Ap [j] ;
pend = (packed) ? (Ap [j+1]) : (p + Anz [j]) ;
for ( ; p < pend ; p++)
{
i = Ai [p] ;
s = abs_value (xtype, Ax, Az, p, Common) ;
if (i == j)
{
W [i] += s ;
}
else if (i < j)
{
W [i] += s ;
W [j] += s ;
}
}
}
}
else if (A->stype < 0)
{
/* ------------------------------------------------------------------ */
/* A is symmetric with lower triangular part stored */
/* ------------------------------------------------------------------ */
/* infinity-norm = 1-norm = max row/col sum */
for (j = 0 ; j < ncol ; j++)
{
p = Ap [j] ;
pend = (packed) ? (Ap [j+1]) : (p + Anz [j]) ;
for ( ; p < pend ; p++)
{
i = Ai [p] ;
s = abs_value (xtype, Ax, Az, p, Common) ;
if (i == j)
{
W [i] += s ;
}
else if (i > j)
{
W [i] += s ;
W [j] += s ;
}
}
}
}
else if (norm == 0)
{
/* ------------------------------------------------------------------ */
/* A is unsymmetric, compute the infinity-norm */
/* ------------------------------------------------------------------ */
/* infinity-norm = max row sum */
for (j = 0 ; j < ncol ; j++)
{
p = Ap [j] ;
pend = (packed) ? (Ap [j+1]) : (p + Anz [j]) ;
for ( ; p < pend ; p++)
{
W [Ai [p]] += abs_value (xtype, Ax, Az, p, Common) ;
}
}
}
else
{
/* ------------------------------------------------------------------ */
/* A is unsymmetric, compute the 1-norm */
/* ------------------------------------------------------------------ */
/* 1-norm = max column sum */
for (j = 0 ; j < ncol ; j++)
{
p = Ap [j] ;
pend = (packed) ? (Ap [j+1]) : (p + Anz [j]) ;
if (xtype == CHOLMOD_PATTERN)
{
s = pend - p ;
}
else
{
s = 0 ;
for ( ; p < pend ; p++)
{
s += abs_value (xtype, Ax, Az, p, Common) ;
}
}
if ((IS_NAN (s) || s > anorm) && !IS_NAN (anorm))
{
anorm = s ;
}
}
}
/* ---------------------------------------------------------------------- */
/* compute the max row sum */
/* ---------------------------------------------------------------------- */
if (A->stype || norm == 0)
{
for (i = 0 ; i < nrow ; i++)
{
s = W [i] ;
if ((IS_NAN (s) || s > anorm) && !IS_NAN (anorm))
{
anorm = s ;
}
W [i] = 0 ;
}
}
/* ---------------------------------------------------------------------- */
/* return result */
/* ---------------------------------------------------------------------- */
return (anorm) ;
}
#endif
syntax highlighted by Code2HTML, v. 0.9.1