/* ========================================================================== */
/* === Core/cholmod_band ==================================================== */
/* ========================================================================== */
/* -----------------------------------------------------------------------------
* 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 = tril (triu (A,k1), k2)
*
* C is a matrix consisting of the diagonals of A from k1 to k2.
*
* k=0 is the main diagonal of A, k=1 is the superdiagonal, k=-1 is the
* subdiagonal, and so on. If A is m-by-n, then:
*
* k1=-m C = tril (A,k2)
* k2=n C = triu (A,k1)
* k1=0 and k2=0 C = diag(A), except C is a matrix, not a vector
*
* Values of k1 and k2 less than -m are treated as -m, and values greater
* than n are treated as n.
*
* A can be of any symmetry (upper, lower, or unsymmetric); C is returned in
* the same form, and packed. If A->stype > 0, entries in the lower
* triangular part of A are ignored, and the opposite is true if
* A->stype < 0. If A has sorted columns, then so does C.
* C has the same size as A.
*
* If inplace is TRUE, then the matrix A is modified in place.
* Only packed matrices can be converted in place.
*
* C can be returned as a numerical valued matrix (if A has numerical values
* and mode > 0), as a pattern-only (mode == 0), or as a pattern-only but with
* the diagonal entries removed (mode < 0).
*
* workspace: none
*
* A can have an xtype of pattern or real. Complex and zomplex cases supported
* only if mode <= 0 (in which case the numerical values are ignored).
*/
#include "cholmod_internal.h"
#include "cholmod_core.h"
static cholmod_sparse *band /* returns C, or NULL if failure */
(
/* ---- input or in/out if inplace is TRUE --- */
cholmod_sparse *A,
/* ---- input ---- */
UF_long k1, /* ignore entries below the k1-st diagonal */
UF_long k2, /* ignore entries above the k2-nd diagonal */
int mode, /* >0: numerical, 0: pattern, <0: pattern (no diagonal) */
int inplace, /* if TRUE, then convert A in place */
/* --------------- */
cholmod_common *Common
)
{
double *Ax, *Cx ;
Int packed, nz, j, p, pend, i, ncol, nrow, jlo, jhi, ilo, ihi, sorted,
values, diag ;
Int *Ap, *Anz, *Ai, *Cp, *Ci ;
cholmod_sparse *C ;
/* ---------------------------------------------------------------------- */
/* 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) ;
packed = A->packed ;
diag = (mode >= 0) ;
if (inplace && !packed)
{
/* cannot operate on an unpacked matrix in place */
ERROR (CHOLMOD_INVALID, "cannot operate on unpacked matrix in-place") ;
return (NULL) ;
}
Common->status = CHOLMOD_OK ;
/* ---------------------------------------------------------------------- */
/* get inputs */
/* ---------------------------------------------------------------------- */
PRINT1 (("k1 %ld k2 %ld\n", k1, k2)) ;
Ap = A->p ;
Anz = A->nz ;
Ai = A->i ;
Ax = A->x ;
sorted = A->sorted ;
if (A->stype > 0)
{
/* ignore any entries in strictly lower triangular part of A */
k1 = MAX (k1, 0) ;
}
if (A->stype < 0)
{
/* ignore any entries in strictly upper triangular part of A */
k2 = MIN (k2, 0) ;
}
ncol = A->ncol ;
nrow = A->nrow ;
/* ensure k1 and k2 are in the range -nrow to +ncol to
* avoid possible integer overflow if k1 and k2 are huge */
k1 = MAX (-nrow, k1) ;
k1 = MIN (k1, ncol) ;
k2 = MAX (-nrow, k2) ;
k2 = MIN (k2, ncol) ;
/* consider columns jlo to jhi. columns outside this range are empty */
jlo = MAX (k1, 0) ;
jhi = MIN (k2+nrow, ncol) ;
if (k1 > k2)
{
/* nothing to do */
jlo = ncol ;
jhi = ncol ;
}
/* ---------------------------------------------------------------------- */
/* allocate C, or operate on A in place */
/* ---------------------------------------------------------------------- */
if (inplace)
{
/* convert A in place */
C = A ;
}
else
{
/* count the number of entries in the result C */
nz = 0 ;
if (sorted)
{
for (j = jlo ; j < jhi ; j++)
{
ilo = j-k2 ;
ihi = j-k1 ;
p = Ap [j] ;
pend = (packed) ? (Ap [j+1]) : (p + Anz [j]) ;
for ( ; p < pend ; p++)
{
i = Ai [p] ;
if (i > ihi)
{
break ;
}
if (i >= ilo && (diag || i != j))
{
nz++ ;
}
}
}
}
else
{
for (j = jlo ; j < jhi ; j++)
{
ilo = j-k2 ;
ihi = j-k1 ;
p = Ap [j] ;
pend = (packed) ? (Ap [j+1]) : (p + Anz [j]) ;
for ( ; p < pend ; p++)
{
i = Ai [p] ;
if (i >= ilo && i <= ihi && (diag || i != j))
{
nz++ ;
}
}
}
}
/* allocate C; A will not be modified. C is sorted if A is sorted */
C = CHOLMOD(allocate_sparse) (A->nrow, ncol, nz, sorted, TRUE,
A->stype, values ? A->xtype : CHOLMOD_PATTERN, Common) ;
if (Common->status < CHOLMOD_OK)
{
return (NULL) ; /* out of memory */
}
}
Cp = C->p ;
Ci = C->i ;
Cx = C->x ;
/* ---------------------------------------------------------------------- */
/* construct C */
/* ---------------------------------------------------------------------- */
/* columns 0 to jlo-1 are empty */
for (j = 0 ; j < jlo ; j++)
{
Cp [j] = 0 ;
}
nz = 0 ;
if (sorted)
{
if (values)
{
/* pattern and values */
ASSERT (diag) ;
for (j = jlo ; j < jhi ; j++)
{
ilo = j-k2 ;
ihi = j-k1 ;
p = Ap [j] ;
pend = (packed) ? (Ap [j+1]) : (p + Anz [j]) ;
Cp [j] = nz ;
for ( ; p < pend ; p++)
{
i = Ai [p] ;
if (i > ihi)
{
break ;
}
if (i >= ilo)
{
Ci [nz] = i ;
Cx [nz] = Ax [p] ;
nz++ ;
}
}
}
}
else
{
/* pattern only, perhaps with no diagonal */
for (j = jlo ; j < jhi ; j++)
{
ilo = j-k2 ;
ihi = j-k1 ;
p = Ap [j] ;
pend = (packed) ? (Ap [j+1]) : (p + Anz [j]) ;
Cp [j] = nz ;
for ( ; p < pend ; p++)
{
i = Ai [p] ;
if (i > ihi)
{
break ;
}
if (i >= ilo && (diag || i != j))
{
Ci [nz++] = i ;
}
}
}
}
}
else
{
if (values)
{
/* pattern and values */
ASSERT (diag) ;
for (j = jlo ; j < jhi ; j++)
{
ilo = j-k2 ;
ihi = j-k1 ;
p = Ap [j] ;
pend = (packed) ? (Ap [j+1]) : (p + Anz [j]) ;
Cp [j] = nz ;
for ( ; p < pend ; p++)
{
i = Ai [p] ;
if (i >= ilo && i <= ihi)
{
Ci [nz] = i ;
Cx [nz] = Ax [p] ;
nz++ ;
}
}
}
}
else
{
/* pattern only, perhaps with no diagonal */
for (j = jlo ; j < jhi ; j++)
{
ilo = j-k2 ;
ihi = j-k1 ;
p = Ap [j] ;
pend = (packed) ? (Ap [j+1]) : (p + Anz [j]) ;
Cp [j] = nz ;
for ( ; p < pend ; p++)
{
i = Ai [p] ;
if (i >= ilo && i <= ihi && (diag || i != j))
{
Ci [nz++] = i ;
}
}
}
}
}
/* columns jhi to ncol-1 are empty */
for (j = jhi ; j <= ncol ; j++)
{
Cp [j] = nz ;
}
/* ---------------------------------------------------------------------- */
/* reduce A in size if done in place */
/* ---------------------------------------------------------------------- */
if (inplace)
{
/* free the unused parts of A, and reduce A->i and A->x in size */
ASSERT (MAX (1,nz) <= A->nzmax) ;
CHOLMOD(reallocate_sparse) (nz, A, Common) ;
ASSERT (Common->status >= CHOLMOD_OK) ;
}
/* ---------------------------------------------------------------------- */
/* return the result C */
/* ---------------------------------------------------------------------- */
DEBUG (i = CHOLMOD(dump_sparse) (C, "band", Common)) ;
ASSERT (IMPLIES (mode < 0, i == 0)) ;
return (C) ;
}
/* ========================================================================== */
/* === cholmod_band ========================================================= */
/* ========================================================================== */
cholmod_sparse *CHOLMOD(band)
(
/* ---- input ---- */
cholmod_sparse *A, /* matrix to extract band matrix from */
UF_long k1, /* ignore entries below the k1-st diagonal */
UF_long k2, /* ignore entries above the k2-nd diagonal */
int mode, /* >0: numerical, 0: pattern, <0: pattern (no diag) */
/* --------------- */
cholmod_common *Common
)
{
return (band (A, k1, k2, mode, FALSE, Common)) ;
}
/* ========================================================================== */
/* === cholmod_band_inplace ================================================= */
/* ========================================================================== */
int CHOLMOD(band_inplace)
(
/* ---- input ---- */
UF_long k1, /* ignore entries below the k1-st diagonal */
UF_long k2, /* ignore entries above the k2-nd diagonal */
int mode, /* >0: numerical, 0: pattern, <0: pattern (no diag) */
/* ---- in/out --- */
cholmod_sparse *A, /* matrix from which entries not in band are removed */
/* --------------- */
cholmod_common *Common
)
{
return (band (A, k1, k2, mode, TRUE, Common) != NULL) ;
}
syntax highlighted by Code2HTML, v. 0.9.1