/* ========================================================================== */
/* === Core/cholmod_triplet ================================================= */
/* ========================================================================== */
/* -----------------------------------------------------------------------------
* 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
* -------------------------------------------------------------------------- */
/* Core utility routines for the cholmod_triplet object:
*
* A sparse matrix held in triplet form is the simplest one for a user to
* create. It consists of a list of nz entries in arbitrary order, held in
* three arrays: i, j, and x, each of length nk. The kth entry is in row i[k],
* column j[k], with value x[k]. There may be duplicate values; if A(i,j)
* appears more than once, its value is the sum of the entries with those row
* and column indices.
*
* Primary routines:
* -----------------
* cholmod_allocate_triplet allocate a triplet matrix
* cholmod_free_triplet free a triplet matrix
*
* Secondary routines:
* -------------------
* cholmod_reallocate_triplet reallocate a triplet matrix
* cholmod_sparse_to_triplet create a triplet matrix copy of a sparse matrix
* cholmod_triplet_to_sparse create a sparse matrix copy of a triplet matrix
* cholmod_copy_triplet create a copy of a triplet matrix
*
* The relationship between an m-by-n cholmod_sparse matrix A and a
* cholmod_triplet matrix (i, j, and x) is identical to how they are used in
* the MATLAB "sparse" and "find" functions:
*
* [i j x] = find (A)
* [m n] = size (A)
* A = sparse (i,j,x,m,n)
*
* with the exception that the cholmod_sparse matrix may be "unpacked", may
* have either sorted or unsorted columns (depending on the option selected),
* and may be symmetric with just the upper or lower triangular part stored.
* Likewise, the cholmod_triplet matrix may contain just the entries in the
* upper or lower triangular part of a symmetric matrix.
*
* MATLAB sparse matrices are always "packed", always have sorted columns,
* and always store both parts of a symmetric matrix. In some cases, MATLAB
* behaves like CHOLMOD by ignoring entries in the upper or lower triangular
* part of a matrix that is otherwise assumed to be symmetric (such as the
* input to chol). In CHOLMOD, that option is a characteristic of the object.
* In MATLAB, that option is based on how a matrix is used as the input to
* a function.
*
* The triplet matrix is provided to give the user a simple way of constructing
* a sparse matrix. There are very few operations supported for triplet
* matrices. The assumption is that they will be converted to cholmod_sparse
* matrix form first.
*
* Adding two triplet matrices simply involves concatenating the contents of
* the three arrays (i, j, and x). To permute a triplet matrix, just replace
* the row and column indices with their permuted values. For example, if
* P is a permutation vector, then P [k] = j means row/column j is the kth
* row/column in C=P*A*P'. In MATLAB notation, C=A(p,p). If Pinv is an array
* of size n and T is the triplet form of A, then:
*
* Ti = T->i ;
* Tj = T->j ;
* for (k = 0 ; k < n ; k++) Pinv [P [k]] = k ;
* for (k = 0 ; k < nz ; k++) Ti [k] = Pinv [Ti [k]] ;
* for (k = 0 ; k < nz ; k++) Tj [k] = Pinv [Tj [k]] ;
*
* overwrites T with the triplet form of C=P*A*P'. The conversion
*
* C = cholmod_triplet_to_sparse (T, 0, &Common) ;
*
* will then return the matrix C = P*A*P'.
*
* Note that T->stype > 0 means that entries in the lower triangular part of
* T are transposed into the upper triangular part when T is converted to
* sparse matrix (cholmod_sparse) form with cholmod_triplet_to_sparse. The
* opposite is true for T->stype < 0.
*
* Since the triplet matrix T is so simple to generate, it's quite easy
* to remove entries that you do not want, prior to converting T to the
* cholmod_sparse form. So if you include these entries in T, CHOLMOD
* assumes that there must be a reason (such as the one above). Thus,
* no entry in a triplet matrix is ever ignored.
*
* Other operations, such as extacting a submatrix, horizontal and vertical
* concatenation, multiply a triplet matrix times a dense matrix, are also
* simple. Multiplying two triplet matrices is not trivial; the simplest
* method is to convert them to cholmod_sparse matrices first.
*
* Supports all xtypes (pattern, real, complex, and zomplex).
*/
#include "cholmod_internal.h"
#include "cholmod_core.h"
/* ========================================================================== */
/* === TEMPLATE ============================================================= */
/* ========================================================================== */
#define PATTERN
#include "t_cholmod_triplet.c"
#define REAL
#include "t_cholmod_triplet.c"
#define COMPLEX
#include "t_cholmod_triplet.c"
#define ZOMPLEX
#include "t_cholmod_triplet.c"
/* ========================================================================== */
/* === cholmod_allocate_triplet ============================================= */
/* ========================================================================== */
/* allocate space for a triplet matrix
*
* workspace: none
*/
cholmod_triplet *CHOLMOD(allocate_triplet)
(
/* ---- input ---- */
size_t nrow, /* # of rows of T */
size_t ncol, /* # of columns of T */
size_t nzmax, /* max # of nonzeros of T */
int stype, /* stype of T */
int xtype, /* CHOLMOD_PATTERN, _REAL, _COMPLEX, or _ZOMPLEX */
/* --------------- */
cholmod_common *Common
)
{
cholmod_triplet *T ;
size_t nzmax0 ;
int ok = TRUE ;
/* ---------------------------------------------------------------------- */
/* get inputs */
/* ---------------------------------------------------------------------- */
RETURN_IF_NULL_COMMON (NULL) ;
if (xtype < CHOLMOD_PATTERN || xtype > CHOLMOD_ZOMPLEX)
{
ERROR (CHOLMOD_INVALID, "xtype invalid") ;
return (NULL) ;
}
/* ensure the dimensions do not cause integer overflow */
(void) CHOLMOD(add_size_t) (ncol, 2, &ok) ;
if (!ok || nrow > Int_max || ncol > Int_max || nzmax > Int_max)
{
ERROR (CHOLMOD_TOO_LARGE, "problem too large") ;
return (NULL) ;
}
Common->status = CHOLMOD_OK ;
/* ---------------------------------------------------------------------- */
/* allocate header */
/* ---------------------------------------------------------------------- */
T = CHOLMOD(malloc) (sizeof (cholmod_triplet), 1, Common) ;
if (Common->status < CHOLMOD_OK)
{
return (NULL) ; /* out of memory */
}
PRINT1 (("cholmod_allocate_triplet %d-by-%d nzmax %d xtype %d\n",
nrow, ncol, nzmax, xtype)) ;
nzmax = MAX (1, nzmax) ;
T->nrow = nrow ;
T->ncol = ncol ;
T->nzmax = nzmax ;
T->nnz = 0 ;
T->stype = stype ;
T->itype = ITYPE ;
T->xtype = xtype ;
T->dtype = DTYPE ;
T->j = NULL ;
T->i = NULL ;
T->x = NULL ;
T->z = NULL ;
/* ---------------------------------------------------------------------- */
/* allocate the matrix itself */
/* ---------------------------------------------------------------------- */
nzmax0 = 0 ;
CHOLMOD(realloc_multiple) (nzmax, 2, xtype, &(T->i), &(T->j),
&(T->x), &(T->z), &nzmax0, Common) ;
if (Common->status < CHOLMOD_OK)
{
CHOLMOD(free_triplet) (&T, Common) ;
return (NULL) ; /* out of memory */
}
return (T) ;
}
/* ========================================================================== */
/* === cholmod_free_triplet ================================================= */
/* ========================================================================== */
/* free a triplet matrix
*
* workspace: none
*/
int CHOLMOD(free_triplet)
(
/* ---- in/out --- */
cholmod_triplet **THandle, /* matrix to deallocate, NULL on output */
/* --------------- */
cholmod_common *Common
)
{
Int nz ;
cholmod_triplet *T ;
RETURN_IF_NULL_COMMON (FALSE) ;
if (THandle == NULL)
{
/* nothing to do */
return (TRUE) ;
}
T = *THandle ;
if (T == NULL)
{
/* nothing to do */
return (TRUE) ;
}
nz = T->nzmax ;
T->j = CHOLMOD(free) (nz, sizeof (Int), T->j, Common) ;
T->i = CHOLMOD(free) (nz, sizeof (Int), T->i, Common) ;
if (T->xtype == CHOLMOD_REAL)
{
T->x = CHOLMOD(free) (nz, sizeof (double), T->x, Common) ;
}
else if (T->xtype == CHOLMOD_COMPLEX)
{
T->x = CHOLMOD(free) (nz, 2*sizeof (double), T->x, Common) ;
}
else if (T->xtype == CHOLMOD_ZOMPLEX)
{
T->x = CHOLMOD(free) (nz, sizeof (double), T->x, Common) ;
T->z = CHOLMOD(free) (nz, sizeof (double), T->z, Common) ;
}
*THandle = CHOLMOD(free) (1, sizeof (cholmod_triplet), (*THandle), Common) ;
return (TRUE) ;
}
/* ========================================================================== */
/* === cholmod_reallocate_triplet =========================================== */
/* ========================================================================== */
/* Change the size of T->i, T->j, and T->x, or allocate them if their current
* size is zero. T->x is not modified if T->xtype is CHOLMOD_PATTERN.
*
* workspace: none
*/
int CHOLMOD(reallocate_triplet)
(
/* ---- input ---- */
size_t nznew, /* new # of entries in T */
/* ---- in/out --- */
cholmod_triplet *T, /* triplet matrix to modify */
/* --------------- */
cholmod_common *Common
)
{
/* ---------------------------------------------------------------------- */
/* get inputs */
/* ---------------------------------------------------------------------- */
RETURN_IF_NULL_COMMON (FALSE) ;
RETURN_IF_NULL (T, FALSE) ;
RETURN_IF_XTYPE_INVALID (T, CHOLMOD_PATTERN, CHOLMOD_ZOMPLEX, FALSE) ;
Common->status = CHOLMOD_OK ;
PRINT1 (("realloc triplet %d to %d, xtype: %d\n",
T->nzmax, nznew, T->xtype)) ;
/* ---------------------------------------------------------------------- */
/* resize the matrix */
/* ---------------------------------------------------------------------- */
CHOLMOD(realloc_multiple) (MAX (1,nznew), 2, T->xtype, &(T->i), &(T->j),
&(T->x), &(T->z), &(T->nzmax), Common) ;
return (Common->status == CHOLMOD_OK) ;
}
/* ========================================================================== */
/* === cholmod_triplet_to_sparse ============================================ */
/* ========================================================================== */
/* Convert a set of triplets into a cholmod_sparse matrix. In MATLAB notation,
* for unsymmetric matrices:
*
* A = sparse (Ti, Tj, Tx, nrow, ncol, nzmax) ;
*
* For the symmetric upper case:
*
* A = sparse (min(Ti,Tj), max(Ti,Tj), Tx, nrow, ncol, nzmax) ;
*
* For the symmetric lower case:
*
* A = sparse (max(Ti,Tj), min(Ti,Tj), Tx, nrow, ncol, nzmax) ;
*
* If Tx is NULL, then A->x is not allocated, and only the pattern of A is
* computed. A is returned in packed form, and can be of any stype
* (upper/lower/unsymmetric). It has enough space to hold the values in T,
* or nzmax, whichever is larger.
*
* workspace: Iwork (max (nrow,ncol))
* allocates a temporary copy of its output matrix.
*
* The resulting sparse matrix has the same xtype as the input triplet matrix.
*/
cholmod_sparse *CHOLMOD(triplet_to_sparse)
(
/* ---- input ---- */
cholmod_triplet *T, /* matrix to copy */
size_t nzmax, /* allocate at least this much space in output matrix */
/* --------------- */
cholmod_common *Common
)
{
cholmod_sparse *R, *A = NULL ;
Int *Wj, *Rp, *Ri, *Rnz, *Ti, *Tj ;
Int i, j, p, k, stype, nrow, ncol, nz, ok ;
size_t anz = 0 ;
/* ---------------------------------------------------------------------- */
/* check inputs */
/* ---------------------------------------------------------------------- */
RETURN_IF_NULL_COMMON (NULL) ;
RETURN_IF_NULL (T, NULL) ;
Ti = T->i ;
Tj = T->j ;
RETURN_IF_NULL (Ti, NULL) ;
RETURN_IF_NULL (Tj, NULL) ;
RETURN_IF_XTYPE_INVALID (T, CHOLMOD_PATTERN, CHOLMOD_ZOMPLEX, NULL) ;
stype = SIGN (T->stype) ;
if (stype && T->nrow != T->ncol)
{
/* inputs invalid */
ERROR (CHOLMOD_INVALID, "matrix invalid") ;
return (NULL) ;
}
Common->status = CHOLMOD_OK ;
DEBUG (CHOLMOD(dump_triplet) (T, "T", Common)) ;
/* ---------------------------------------------------------------------- */
/* get inputs */
/* ---------------------------------------------------------------------- */
nrow = T->nrow ;
ncol = T->ncol ;
nz = T->nnz ;
/* ---------------------------------------------------------------------- */
/* allocate workspace */
/* ---------------------------------------------------------------------- */
CHOLMOD(allocate_work) (0, MAX (nrow, ncol), 0, Common) ;
if (Common->status < CHOLMOD_OK)
{
return (NULL) ; /* out of memory */
}
/* ---------------------------------------------------------------------- */
/* allocate temporary matrix R */
/* ---------------------------------------------------------------------- */
R = CHOLMOD(allocate_sparse) (ncol, nrow, nz, FALSE, FALSE, -stype,
T->xtype, Common) ;
if (Common->status < CHOLMOD_OK)
{
return (NULL) ; /* out of memory */
}
Rp = R->p ;
Ri = R->i ;
Rnz = R->nz ;
/* ---------------------------------------------------------------------- */
/* count the entries in each row of A (also counting duplicates) */
/* ---------------------------------------------------------------------- */
for (i = 0 ; i < nrow ; i++)
{
Rnz [i] = 0 ;
}
if (stype > 0)
{
for (k = 0 ; k < nz ; k++)
{
i = Ti [k] ;
j = Tj [k] ;
if (i < 0 || i >= nrow || j < 0 || j >= ncol)
{
ERROR (CHOLMOD_INVALID, "index out of range") ;
break ;
}
/* A will be symmetric with just the upper triangular part stored.
* Create a matrix R that is lower triangular. Entries in the
* upper part of R are transposed to the lower part. */
Rnz [MIN (i,j)]++ ;
}
}
else if (stype < 0)
{
for (k = 0 ; k < nz ; k++)
{
i = Ti [k] ;
j = Tj [k] ;
if (i < 0 || i >= nrow || j < 0 || j >= ncol)
{
ERROR (CHOLMOD_INVALID, "index out of range") ;
break ;
}
/* A will be symmetric with just the lower triangular part stored.
* Create a matrix R that is upper triangular. Entries in the
* lower part of R are transposed to the upper part. */
Rnz [MAX (i,j)]++ ;
}
}
else
{
for (k = 0 ; k < nz ; k++)
{
i = Ti [k] ;
j = Tj [k] ;
if (i < 0 || i >= nrow || j < 0 || j >= ncol)
{
ERROR (CHOLMOD_INVALID, "index out of range") ;
break ;
}
/* constructing an unsymmetric matrix */
Rnz [i]++ ;
}
}
if (Common->status < CHOLMOD_OK)
{
/* triplet matrix is invalid */
CHOLMOD(free_sparse) (&R, Common) ;
return (NULL) ;
}
/* ---------------------------------------------------------------------- */
/* construct the row pointers */
/* ---------------------------------------------------------------------- */
p = 0 ;
for (i = 0 ; i < nrow ; i++)
{
Rp [i] = p ;
p += Rnz [i] ;
}
Rp [nrow] = p ;
/* use Wj (i/l/l) as temporary row pointers */
Wj = Common->Iwork ; /* size MAX (nrow,ncol) FUTURE WORK: (i/l/l) */
for (i = 0 ; i < nrow ; i++)
{
Wj [i] = Rp [i] ;
}
/* ---------------------------------------------------------------------- */
/* construct triplet matrix, using template routine */
/* ---------------------------------------------------------------------- */
switch (T->xtype)
{
case CHOLMOD_PATTERN:
anz = p_cholmod_triplet_to_sparse (T, R, Common) ;
break ;
case CHOLMOD_REAL:
anz = r_cholmod_triplet_to_sparse (T, R, Common) ;
break ;
case CHOLMOD_COMPLEX:
anz = c_cholmod_triplet_to_sparse (T, R, Common) ;
break ;
case CHOLMOD_ZOMPLEX:
anz = z_cholmod_triplet_to_sparse (T, R, Common) ;
break ;
}
/* ---------------------------------------------------------------------- */
/* A = R' (array transpose, not complex conjugate transpose) */
/* ---------------------------------------------------------------------- */
/* workspace: Iwork (R->nrow), which is A->ncol */
ASSERT (CHOLMOD(dump_sparse) (R, "R", Common) >= 0) ;
A = CHOLMOD(allocate_sparse) (nrow, ncol, MAX (anz, nzmax), TRUE, TRUE,
stype, T->xtype, Common) ;
if (stype)
{
ok = CHOLMOD(transpose_sym) (R, 1, NULL, A, Common) ;
}
else
{
ok = CHOLMOD(transpose_unsym) (R, 1, NULL, NULL, 0, A, Common) ;
}
CHOLMOD(free_sparse) (&R, Common) ;
if (Common->status < CHOLMOD_OK)
{
CHOLMOD(free_sparse) (&A, Common) ;
}
/* ---------------------------------------------------------------------- */
/* return result */
/* ---------------------------------------------------------------------- */
ASSERT (CHOLMOD(dump_sparse) (A, "A = triplet(T) result", Common) >= 0) ;
return (A) ;
}
/* ========================================================================== */
/* === cholmod_sparse_to_triplet ============================================ */
/* ========================================================================== */
/* Converts a sparse column-oriented matrix to triplet form.
* The resulting triplet matrix has the same xtype as the sparse matrix.
*
* workspace: none
*/
cholmod_triplet *CHOLMOD(sparse_to_triplet)
(
/* ---- input ---- */
cholmod_sparse *A, /* matrix to copy */
/* --------------- */
cholmod_common *Common
)
{
double *Ax, *Az, *Tx, *Tz ;
Int *Ap, *Ai, *Ti, *Tj, *Anz ;
cholmod_triplet *T ;
Int i, xtype, p, pend, k, j, nrow, ncol, nz, stype, packed, up, lo,
both ;
/* ---------------------------------------------------------------------- */
/* check inputs */
/* ---------------------------------------------------------------------- */
RETURN_IF_NULL_COMMON (NULL) ;
RETURN_IF_NULL (A, NULL) ;
RETURN_IF_XTYPE_INVALID (A, CHOLMOD_PATTERN, CHOLMOD_ZOMPLEX, NULL) ;
stype = SIGN (A->stype) ;
nrow = A->nrow ;
ncol = A->ncol ;
if (stype && nrow != ncol)
{
/* inputs invalid */
ERROR (CHOLMOD_INVALID, "matrix invalid") ;
return (NULL) ;
}
Ax = A->x ;
Az = A->z ;
xtype = A->xtype ;
Common->status = CHOLMOD_OK ;
ASSERT (CHOLMOD(dump_sparse) (A, "A", Common) >= 0) ;
/* ---------------------------------------------------------------------- */
/* allocate triplet matrix */
/* ---------------------------------------------------------------------- */
nz = CHOLMOD(nnz) (A, Common) ;
T = CHOLMOD(allocate_triplet) (nrow, ncol, nz, A->stype, A->xtype, Common) ;
if (Common->status < CHOLMOD_OK)
{
return (NULL) ; /* out of memory */
}
/* ---------------------------------------------------------------------- */
/* convert to a sparse matrix */
/* ---------------------------------------------------------------------- */
Ap = A->p ;
Ai = A->i ;
Anz = A->nz ;
packed = A->packed ;
Ti = T->i ;
Tj = T->j ;
Tx = T->x ;
Tz = T->z ;
T->stype = A->stype ;
both = (A->stype == 0) ;
up = (A->stype > 0) ;
lo = (A->stype < 0) ;
k = 0 ;
for (j = 0 ; j < ncol ; j++)
{
p = Ap [j] ;
pend = (packed) ? (Ap [j+1]) : (p + Anz [j]) ;
for ( ; p < pend ; p++)
{
i = Ai [p] ;
if (both || (up && i <= j) || (lo && i >= j))
{
Ti [k] = Ai [p] ;
Tj [k] = j ;
if (xtype == CHOLMOD_REAL)
{
Tx [k] = Ax [p] ;
}
else if (xtype == CHOLMOD_COMPLEX)
{
Tx [2*k ] = Ax [2*p ] ;
Tx [2*k+1] = Ax [2*p+1] ;
}
else if (xtype == CHOLMOD_ZOMPLEX)
{
Tx [k] = Ax [p] ;
Tz [k] = Az [p] ;
}
k++ ;
ASSERT (k <= nz) ;
}
}
}
T->nnz = k ;
/* ---------------------------------------------------------------------- */
/* return result */
/* ---------------------------------------------------------------------- */
ASSERT (CHOLMOD(dump_triplet) (T, "T", Common)) ;
return (T) ;
}
/* ========================================================================== */
/* === cholmod_copy_triplet ================================================= */
/* ========================================================================== */
/* Create an exact copy of a triplet matrix, except that entries in unused
* space are not copied (they might not be initialized, and copying them would
* cause program checkers such as purify and valgrind to complain).
* The output triplet matrix has the same xtype as the input triplet matrix.
*/
cholmod_triplet *CHOLMOD(copy_triplet)
(
/* ---- input ---- */
cholmod_triplet *T, /* matrix to copy */
/* --------------- */
cholmod_common *Common
)
{
double *Tx, *Tz, *Cx, *Cz ;
Int *Ci, *Cj, *Ti, *Tj ;
cholmod_triplet *C ;
Int xtype, k, nz ;
/* ---------------------------------------------------------------------- */
/* check inputs */
/* ---------------------------------------------------------------------- */
RETURN_IF_NULL_COMMON (NULL) ;
RETURN_IF_NULL (T, NULL) ;
RETURN_IF_XTYPE_INVALID (T, CHOLMOD_PATTERN, CHOLMOD_ZOMPLEX, NULL) ;
nz = T->nnz ;
Ti = T->i ;
Tj = T->j ;
Tx = T->x ;
Tz = T->z ;
xtype = T->xtype ;
RETURN_IF_NULL (Ti, NULL) ;
RETURN_IF_NULL (Tj, NULL) ;
Common->status = CHOLMOD_OK ;
DEBUG (CHOLMOD(dump_triplet) (T, "T input", Common)) ;
/* ---------------------------------------------------------------------- */
/* allocate copy */
/* ---------------------------------------------------------------------- */
C = CHOLMOD(allocate_triplet) (T->nrow, T->ncol, T->nzmax, T->stype,
xtype, Common) ;
if (Common->status < CHOLMOD_OK)
{
return (NULL) ; /* out of memory */
}
/* ---------------------------------------------------------------------- */
/* copy the triplet matrix */
/* ---------------------------------------------------------------------- */
Ci = C->i ;
Cj = C->j ;
Cx = C->x ;
Cz = C->z ;
C->nnz = nz ;
for (k = 0 ; k < nz ; k++)
{
Ci [k] = Ti [k] ;
}
for (k = 0 ; k < nz ; k++)
{
Cj [k] = Tj [k] ;
}
if (xtype == CHOLMOD_REAL)
{
for (k = 0 ; k < nz ; k++)
{
Cx [k] = Tx [k] ;
}
}
else if (xtype == CHOLMOD_COMPLEX)
{
for (k = 0 ; k < nz ; k++)
{
Cx [2*k ] = Tx [2*k ] ;
Cx [2*k+1] = Tx [2*k+1] ;
}
}
else if (xtype == CHOLMOD_ZOMPLEX)
{
for (k = 0 ; k < nz ; k++)
{
Cx [k] = Tx [k] ;
Cz [k] = Tz [k] ;
}
}
/* ---------------------------------------------------------------------- */
/* return the result */
/* ---------------------------------------------------------------------- */
ASSERT (CHOLMOD(dump_triplet) (C, "C triplet copy", Common)) ;
return (C) ;
}
syntax highlighted by Code2HTML, v. 0.9.1