/* ========================================================================== */
/* === MATLAB/cholmod_matlab ================================================ */
/* ========================================================================== */
/* Utility routines for the CHOLMOD MATLAB mexFunctions.
*
* If CHOLMOD runs out of memory, MATLAB will terminate the mexFunction
* immediately since it uses mxMalloc (see sputil_config, below). Likewise,
* if mxCreate* or mxMalloc (as called in this file) fails, MATLAB will also
* terminate the mexFunction. When this occurs, MATLAB frees all allocated
* memory, so we don't have to worry about memory leaks. If this were not the
* case, the routines in this file would suffer from memory leaks whenever an
* error occurred.
*/
#include "cholmod_matlab.h"
#ifndef INT64_T
#define INT64_T long long
#endif
/* This file pointer is used for the mread and mwrite mexFunctions. It must
* be a global variable, because the file pointer is not passed to the
* sputil_error_handler function when an error occurs. */
FILE *sputil_file = NULL ;
/* ========================================================================== */
/* === sputil_config ======================================================== */
/* ========================================================================== */
/* Define function pointers and other parameters for a mexFunction */
void sputil_config (int spumoni, cholmod_common *cm)
{
/* cholmod_solve must return a real or zomplex X for MATLAB */
cm->prefer_zomplex = TRUE ;
/* use mxMalloc and related memory management routines */
cm->malloc_memory = mxMalloc ;
cm->free_memory = mxFree ;
cm->realloc_memory = mxRealloc ;
cm->calloc_memory = mxCalloc ;
/* printing and error handling */
if (spumoni == 0)
{
/* do not print anything from within CHOLMOD */
cm->print = -1 ;
cm->print_function = NULL ;
}
else
{
/* spumoni = 1: print warning and error messages. cholmod_print_*
* routines will print a one-line summary of each object printed.
* spumoni = 2: also print a short summary of each object.
*/
cm->print = spumoni + 2 ;
cm->print_function = mexPrintf ;
}
/* error handler */
cm->error_handler = sputil_error_handler ;
/* complex arithmetic */
cm->complex_divide = cholmod_divcomplex ;
cm->hypotenuse = cholmod_hypot ;
#ifndef NPARTITION
#if defined(METIS_VERSION)
#if (METIS_VERSION >= METIS_VER(4,0,2))
/* METIS 4.0.2 uses function pointers for malloc and free */
METIS_malloc = cm->malloc_memory ;
METIS_free = cm->free_memory ;
#endif
#endif
#endif
/* Turn off METIS memory guard. It is not needed, because mxMalloc will
* safely terminate the mexFunction and free any workspace without killing
* all of MATLAB. This assumes cholmod_make was used to compile CHOLMOD
* for MATLAB. */
cm->metis_memory = 0.0 ;
}
/* ========================================================================== */
/* === sputil_error_handler ================================================= */
/* ========================================================================== */
void sputil_error_handler (int status, char *file, int line, char *message)
{
if (status < CHOLMOD_OK)
{
/*
mexPrintf ("ERROR: file %s line %d, status %d\n", file, line, status) ;
*/
if (sputil_file != NULL)
{
fclose (sputil_file) ;
sputil_file = NULL ;
}
mexErrMsgTxt (message) ;
}
/*
else
{
mexPrintf ("Warning: file %s line %d, status %d\n", file, line, status);
}
*/
}
/* ========================================================================== */
/* === sputil_get_sparse ==================================================== */
/* ========================================================================== */
/* Create a shallow CHOLMOD copy of a MATLAB sparse matrix. No memory is
* allocated. The resulting matrix A must not be modified.
*/
cholmod_sparse *sputil_get_sparse
(
const mxArray *Amatlab, /* MATLAB version of the matrix */
cholmod_sparse *A, /* CHOLMOD version of the matrix */
double *dummy, /* a pointer to a valid scalar double */
int stype /* -1: lower, 0: unsymmetric, 1: upper */
)
{
int *Ap ;
A->nrow = mxGetM (Amatlab) ;
A->ncol = mxGetN (Amatlab) ;
A->p = mxGetJc (Amatlab) ;
A->i = mxGetIr (Amatlab) ;
Ap = A->p ;
A->nzmax = Ap [A->ncol] ;
A->packed = TRUE ;
A->sorted = TRUE ;
A->nz = NULL ;
A->itype = CHOLMOD_INT ;
A->dtype = CHOLMOD_DOUBLE ;
A->stype = stype ;
if (mxIsLogical (Amatlab))
{
A->x = NULL ;
A->z = NULL ;
A->xtype = CHOLMOD_PATTERN ;
}
else if (mxIsEmpty (Amatlab))
{
/* this is not dereferenced, but the existence (non-NULL) of these
* pointers is checked in CHOLMOD */
A->x = dummy ;
A->z = dummy ;
A->xtype = mxIsComplex (Amatlab) ? CHOLMOD_ZOMPLEX : CHOLMOD_REAL ;
}
else if (mxIsDouble (Amatlab))
{
A->x = mxGetPr (Amatlab) ;
A->z = mxGetPi (Amatlab) ;
A->xtype = mxIsComplex (Amatlab) ? CHOLMOD_ZOMPLEX : CHOLMOD_REAL ;
}
else
{
/* only logical and complex/real double matrices supported */
sputil_error (ERROR_INVALID_TYPE, 0) ;
}
return (A) ;
}
/* ========================================================================== */
/* === sputil_get_dense ===================================================== */
/* ========================================================================== */
/* Create a shallow CHOLMOD copy of a MATLAB dense matrix. No memory is
* allocated. Only double (real and zomplex) matrices are supported. The
* resulting matrix B must not be modified.
*/
cholmod_dense *sputil_get_dense
(
const mxArray *Amatlab, /* MATLAB version of the matrix */
cholmod_dense *A, /* CHOLMOD version of the matrix */
double *dummy /* a pointer to a valid scalar double */
)
{
A->nrow = mxGetM (Amatlab) ;
A->ncol = mxGetN (Amatlab) ;
A->d = A->nrow ;
A->nzmax = A->nrow * A->ncol ;
A->dtype = CHOLMOD_DOUBLE ;
if (mxIsEmpty (Amatlab))
{
A->x = dummy ;
A->z = dummy ;
}
else if (mxIsDouble (Amatlab))
{
A->x = mxGetPr (Amatlab) ;
A->z = mxGetPi (Amatlab) ;
}
else
{
/* only full double matrices supported by sputil_get_dense */
sputil_error (ERROR_INVALID_TYPE, 0) ;
}
A->xtype = mxIsComplex (Amatlab) ? CHOLMOD_ZOMPLEX : CHOLMOD_REAL ;
return (A) ;
}
/* ========================================================================== */
/* === sputil_get_sparse_pattern ============================================ */
/* ========================================================================== */
/* Create a CHOLMOD_PATTERN sparse matrix for a MATLAB matrix, depending on the
* type:
*
* (1) MATLAB full real double: duplicate CHOLMOD_REAL sparse matrix.
* (2) MATLAB full complex double: duplicate CHOLMOD_ZOMPLEX sparse matrix.
* (3) MATLAB full logical: duplicate CHOLMOD_PATTERN sparse matrix.
* (4) MATLAB sparse real double: shallow CHOLMOD_REAL copy.
* (5) MATLAB sparse complex double: shallow CHOLMOD_ZOMPLEX copy.
* (6) MATLAB sparse logical: shallow CHOLMOD_PATTERN copy.
*
* A shallow copy or duplicate is returned; the shallow copy must not be freed.
* For a shallow copy, the return value A is the same as Ashallow. For a
* complete duplicate, A and Ashallow will differ.
*/
cholmod_sparse *sputil_get_sparse_pattern
(
const mxArray *Amatlab, /* MATLAB version of the matrix */
cholmod_sparse *Ashallow, /* shallow CHOLMOD version of the matrix */
double *dummy, /* a pointer to a valid scalar double */
cholmod_common *cm
)
{
cholmod_sparse *A ;
if (!mxIsSparse (Amatlab))
{
/* ------------------------------------------------------------------ */
/* A = sparse (X) where X is full */
/* ------------------------------------------------------------------ */
if (mxIsDouble (Amatlab))
{
/* -------------------------------------------------------------- */
/* convert full double X into sparse matrix A (pattern only) */
/* -------------------------------------------------------------- */
cholmod_dense Xmatrix, *X ;
X = sputil_get_dense (Amatlab, &Xmatrix, dummy) ;
A = cholmod_dense_to_sparse (X, FALSE, cm) ;
}
else if (mxIsLogical (Amatlab))
{
/* -------------------------------------------------------------- */
/* convert full logical MATLAB matrix into CHOLMOD_PATTERN */
/* -------------------------------------------------------------- */
/* (this is copied and modified from t_cholmod_dense.c) */
char *x ;
int *Ap, *Ai ;
int nrow, ncol, i, j, nz, nzmax, p ;
/* -------------------------------------------------------------- */
/* count the number of nonzeros in the result */
/* -------------------------------------------------------------- */
nrow = mxGetM (Amatlab) ;
ncol = mxGetN (Amatlab) ;
x = (char *) mxGetData (Amatlab) ;
nzmax = nrow * ncol ;
for (nz = 0, j = 0 ; j < nzmax ; j++)
{
if (x [j])
{
nz++ ;
}
}
/* -------------------------------------------------------------- */
/* allocate the result A */
/* -------------------------------------------------------------- */
A = cholmod_allocate_sparse (nrow, ncol, nz, TRUE, TRUE, 0,
CHOLMOD_PATTERN, cm) ;
if (cm->status < CHOLMOD_OK)
{
return (NULL) ; /* out of memory */
}
Ap = A->p ;
Ai = A->i ;
/* -------------------------------------------------------------- */
/* copy the full logical matrix into the sparse matrix A */
/* -------------------------------------------------------------- */
p = 0 ;
for (j = 0 ; j < ncol ; j++)
{
Ap [j] = p ;
for (i = 0 ; i < nrow ; i++)
{
if (x [i+j*nrow])
{
Ai [p++] = i ;
}
}
}
/* ASSERT (p == nz) ; */
Ap [ncol] = nz ;
}
else
{
/* only double and logical matrices supported */
sputil_error (ERROR_INVALID_TYPE, 0) ;
}
}
else
{
/* ------------------------------------------------------------------ */
/* create a shallow copy of sparse matrix A (default stype is zero) */
/* ------------------------------------------------------------------ */
A = sputil_get_sparse (Amatlab, Ashallow, dummy, 0) ;
A->x = NULL ;
A->z = NULL ;
A->xtype = CHOLMOD_PATTERN ;
}
return (A) ;
}
/* ========================================================================== */
/* === sputil_put_sparse ==================================================== */
/* ========================================================================== */
/* Creates a true MATLAB version of a CHOLMOD sparse matrix. The CHOLMOD sparse
* matrix is destroyed. Both real and zomplex matrices are supported.
*/
mxArray *sputil_put_sparse
(
cholmod_sparse **Ahandle, /* CHOLMOD version of the matrix */
cholmod_common *cm
)
{
mxArray *Amatlab ;
cholmod_sparse *A ;
A = *Ahandle ;
Amatlab = mxCreateSparse (0, 0, 0,
(A->xtype != CHOLMOD_REAL) ? mxCOMPLEX: mxREAL) ;
mxSetM (Amatlab, A->nrow) ;
mxSetN (Amatlab, A->ncol) ;
mxSetNzmax (Amatlab, A->nzmax) ;
mxFree (mxGetJc (Amatlab)) ;
mxFree (mxGetIr (Amatlab)) ;
mxFree (mxGetPr (Amatlab)) ;
mxSetJc (Amatlab, A->p) ;
mxSetIr (Amatlab, A->i) ;
mxSetPr (Amatlab, A->x) ;
mexMakeMemoryPersistent (A->p) ;
mexMakeMemoryPersistent (A->i) ;
mexMakeMemoryPersistent (A->x) ;
if (A->xtype != CHOLMOD_REAL)
{
mxFree (mxGetPi (Amatlab)) ;
mxSetPi (Amatlab, A->z) ;
mexMakeMemoryPersistent (A->z) ;
}
A->p = NULL ;
A->i = NULL ;
A->x = NULL ;
A->z = NULL ;
cholmod_free_sparse (Ahandle, cm) ;
return (Amatlab) ;
}
/* ========================================================================== */
/* === sputil_put_dense ===================================================== */
/* ========================================================================== */
/* Creates a true MATLAB version of a CHOLMOD dense matrix. The CHOLMOD dense
* matrix is destroyed. Both real and zomplex matrices are supported.
*/
mxArray *sputil_put_dense
(
cholmod_dense **Ahandle, /* CHOLMOD version of the matrix */
cholmod_common *cm
)
{
mxArray *Amatlab ;
cholmod_dense *A ;
A = *Ahandle ;
Amatlab = mxCreateDoubleMatrix (0, 0,
(A->xtype != CHOLMOD_REAL) ? mxCOMPLEX: mxREAL) ;
mxSetM (Amatlab, A->nrow) ;
mxSetN (Amatlab, A->ncol) ;
mxFree (mxGetPr (Amatlab)) ;
mxSetPr (Amatlab, A->x) ;
mexMakeMemoryPersistent (A->x) ;
if (A->xtype != CHOLMOD_REAL)
{
mxFree (mxGetPi (Amatlab)) ;
mxSetPi (Amatlab, A->z) ;
mexMakeMemoryPersistent (A->z) ;
}
A->x = NULL ;
A->z = NULL ;
cholmod_free_dense (Ahandle, cm) ;
return (Amatlab) ;
}
/* ========================================================================== */
/* === sputil_put_int ======================================================= */
/* ========================================================================== */
/* Convert an int vector into a double mxArray */
mxArray *sputil_put_int
(
int *P, /* vector to convert */
int n, /* length of P */
int one_based /* 1 if convert from 0-based to 1-based, 0 otherwise */
)
{
double *p ;
mxArray *Q ;
int i ;
Q = mxCreateDoubleMatrix (1, n, mxREAL) ;
p = mxGetPr (Q) ;
for (i = 0 ; i < n ; i++)
{
p [i] = (double) (P [i] + one_based) ;
}
return (Q) ;
}
/* ========================================================================== */
/* === sputil_error ========================================================= */
/* ========================================================================== */
/* An integer is out of range, or other error has occurred. */
void sputil_error
(
int error, /* kind of error */
int is_index /* TRUE if a matrix index, FALSE if a matrix dimension */
)
{
if (error == ERROR_TOO_SMALL)
{
mexErrMsgTxt (is_index ?
"sparse: index into matrix must be positive" :
"sparse: sparse matrix sizes must be non-negative integers") ;
}
else if (error == ERROR_HUGE)
{
mexErrMsgTxt (is_index ?
"sparse: index into matrix is too large" :
"sparse: sparse matrix size is too large") ;
}
else if (error == ERROR_NOT_INTEGER)
{
mexErrMsgTxt (is_index ?
"sparse: index into matrix must be an integer" :
"sparse: sparse matrix size must be an integer") ;
}
else if (error == ERROR_TOO_LARGE)
{
mexErrMsgTxt ("sparse: index exceeds matrix dimensions") ;
}
else if (error == ERROR_USAGE)
{
mexErrMsgTxt (
"Usage:\n"
"A = sparse (S)\n"
"A = sparse (i,j,s,m,n,nzmax)\n"
"A = sparse (i,j,s,m,n)\n"
"A = sparse (i,j,s)\n"
"A = sparse (m,n)\n") ;
}
else if (error == ERROR_LENGTH)
{
mexErrMsgTxt ("sparse: vectors must be the same lengths") ;
}
else if (error == ERROR_INVALID_TYPE)
{
mexErrMsgTxt ("matrix class not supported") ;
}
}
/* ========================================================================== */
/* === sputil_double_to_int ================================================= */
/* ========================================================================== */
/* convert a double into an integer */
int sputil_double_to_int /* returns integer value of x */
(
double x, /* double value to convert */
int is_index, /* TRUE if a matrix index, FALSE if a matrix dimension */
int n /* if a matrix index, x cannot exceed this dimension,
* except that -1 is treated as infinity */
)
{
int i ;
if (x > INT_MAX)
{
/* x is way too big for an integer */
sputil_error (ERROR_HUGE, is_index) ;
}
else if (x < 0)
{
/* x must be non-negative */
sputil_error (ERROR_TOO_SMALL, is_index) ;
}
i = (int) x ;
if (x != (double) i)
{
/* x must be an integer */
sputil_error (ERROR_NOT_INTEGER, is_index) ;
}
if (is_index)
{
if (i < 1)
{
sputil_error (ERROR_TOO_SMALL, is_index) ;
}
else if (i > n && n != EMPTY)
{
sputil_error (ERROR_TOO_LARGE, is_index) ;
}
}
return (i) ;
}
/* ========================================================================== */
/* === sputil_nelements ===================================================== */
/* ========================================================================== */
/* return the number of elements in an mxArray. Trigger an error on integer
* overflow (in case the argument is sparse) */
int sputil_nelements (const mxArray *arg)
{
double size ;
const int *dims ;
int k, ndims ;
ndims = mxGetNumberOfDimensions (arg) ;
dims = mxGetDimensions (arg) ;
size = 1 ;
for (k = 0 ; k < ndims ; k++)
{
size *= dims [k] ;
}
return (sputil_double_to_int (size, FALSE, 0)) ;
}
/* ========================================================================== */
/* === sputil_get_double ==================================================== */
/* ========================================================================== */
double sputil_get_double (const mxArray *arg)
{
if (sputil_nelements (arg) < 1)
{
/* [] is not a scalar, but its value is zero so that
* sparse ([],[],[]) is a 0-by-0 matrix */
return (0) ;
}
return (mxGetScalar (arg)) ;
}
/* ========================================================================== */
/* === sputil_get_integer =================================================== */
/* ========================================================================== */
/* return an argument as a non-negative integer scalar, or -1 if error */
int sputil_get_integer
(
const mxArray *arg, /* MATLAB argument to convert */
int is_index, /* TRUE if an index, FALSE if a matrix dimension */
int n /* maximum value, if an index */
)
{
double x = sputil_get_double (arg) ;
if (mxIsInf (x) || mxIsNaN (x))
{
/* arg is Inf or NaN, return -1 */
return (EMPTY) ;
}
return (sputil_double_to_int (x, is_index, n)) ;
}
/* ========================================================================== */
/* === sputil_trim ========================================================== */
/* ========================================================================== */
/* Remove columns k to n-1 from a sparse matrix S, leaving columns 0 to k-1.
* S must be packed (there can be no S->nz array). This condition is not
* checked, since only packed matrices are passed to this routine. */
void sputil_trim
(
cholmod_sparse *S,
int k,
cholmod_common *cm
)
{
int *Sp ;
int ncol ;
size_t n1, nznew ;
if (S == NULL)
{
return ;
}
ncol = S->ncol ;
if (k < 0 || k >= ncol)
{
/* do not modify S */
return ;
}
/* reduce S->p in size. This cannot fail. */
n1 = ncol + 1 ;
S->p = cholmod_realloc (k+1, sizeof (int), S->p, &n1, cm) ;
/* get the new number of entries in S */
Sp = S->p ;
nznew = Sp [k] ;
/* reduce S->i, S->x, and S->z (if present) to size nznew */
cholmod_reallocate_sparse (nznew, S, cm) ;
/* S now has only k columns */
S->ncol = k ;
}
/* ========================================================================== */
/* === sputil_extract_zeros ================================================= */
/* ========================================================================== */
/* Create a sparse binary (real double) matrix Z that contains the pattern
* of explicit zeros in the sparse real/zomplex double matrix A. */
cholmod_sparse *sputil_extract_zeros
(
cholmod_sparse *A,
cholmod_common *cm
)
{
int *Ap, *Ai, *Zp, *Zi ;
double *Ax, *Az, *Zx ;
int j, p, nzeros = 0, is_complex, pz, nrow, ncol ;
cholmod_sparse *Z ;
if (A == NULL || A->xtype == CHOLMOD_PATTERN || A->xtype == CHOLMOD_COMPLEX)
{
/* only sparse real/zomplex double matrices supported */
sputil_error (ERROR_INVALID_TYPE, 0) ;
}
Ap = A->p ;
Ai = A->i ;
Ax = A->x ;
Az = A->z ;
ncol = A->ncol ;
nrow = A->nrow ;
is_complex = (A->xtype == CHOLMOD_ZOMPLEX) ;
/* count the number of zeros in a sparse matrix A */
for (j = 0 ; j < ncol ; j++)
{
for (p = Ap [j] ; p < Ap [j+1] ; p++)
{
if (CHOLMOD_IS_ZERO (Ax [p]) &&
((is_complex) ? CHOLMOD_IS_ZERO (Az [p]) : TRUE))
{
nzeros++ ;
}
}
}
/* allocate the Z matrix with space for all the zero entries */
Z = cholmod_spzeros (nrow, ncol, nzeros, CHOLMOD_REAL, cm) ;
/* extract the zeros from A and store them in Z as binary values */
if (nzeros > 0)
{
Zp = Z->p ;
Zi = Z->i ;
Zx = Z->x ;
pz = 0 ;
for (j = 0 ; j < ncol ; j++)
{
Zp [j] = pz ;
for (p = Ap [j] ; p < Ap [j+1] ; p++)
{
if (CHOLMOD_IS_ZERO (Ax [p]) &&
((is_complex) ? CHOLMOD_IS_ZERO (Az [p]) : TRUE))
{
Zi [pz] = Ai [p] ;
Zx [pz] = 1 ;
pz++ ;
}
}
}
Zp [ncol] = pz ;
}
return (Z) ;
}
/* ========================================================================== */
/* === sputil_drop_zeros ==================================================== */
/* ========================================================================== */
/* Drop zeros from a CHOLMOD sparse matrix (zomplex or real). This is very
* similar to CHOLMOD/MatrixOps/cholmod_drop, except that this routine has
* no tolerance parameter and it can handle zomplex matrices. NaN's are left
* in the matrix. If this is used on the sparse matrix version of the factor
* L, then the update/downdate methods cannot be applied to L (ldlupdate).
*/
void sputil_drop_zeros
(
cholmod_sparse *S
)
{
double sik, zik ;
int *Sp, *Si ;
double *Sx, *Sz ;
int pdest, k, ncol, p, pend ;
if (S == NULL)
{
return ;
}
Sp = S->p ;
Si = S->i ;
Sx = S->x ;
Sz = S->z ;
pdest = 0 ;
ncol = S->ncol ;
if (S->xtype == CHOLMOD_ZOMPLEX)
{
for (k = 0 ; k < ncol ; k++)
{
p = Sp [k] ;
pend = Sp [k+1] ;
Sp [k] = pdest ;
for ( ; p < pend ; p++)
{
sik = Sx [p] ;
zik = Sz [p] ;
if (CHOLMOD_IS_NONZERO (sik) || CHOLMOD_IS_NONZERO (zik))
{
if (p != pdest)
{
Si [pdest] = Si [p] ;
Sx [pdest] = sik ;
Sz [pdest] = zik ;
}
pdest++ ;
}
}
}
}
else
{
for (k = 0 ; k < ncol ; k++)
{
p = Sp [k] ;
pend = Sp [k+1] ;
Sp [k] = pdest ;
for ( ; p < pend ; p++)
{
sik = Sx [p] ;
if (CHOLMOD_IS_NONZERO (sik))
{
if (p != pdest)
{
Si [pdest] = Si [p] ;
Sx [pdest] = sik ;
}
pdest++ ;
}
}
}
}
Sp [ncol] = pdest ;
}
/* ========================================================================== */
/* === sputil_copy_ij ======================================================= */
/* ========================================================================== */
/* copy i or j arguments into an int vector. For small integer types, i and
* and j can be returned with negative entries; this error condition is caught
* later, in cholmod_triplet_to_sparse.
*/
int sputil_copy_ij /* returns the dimension, n */
(
int is_scalar, /* TRUE if argument is a scalar, FALSE otherwise */
int scalar, /* scalar value of the argument */
void *vector, /* vector value of the argument */
mxClassID category, /* type of vector */
int nz, /* length of output vector I */
int n, /* maximum dimension, EMPTY if not yet known */
int *I /* vector of length nz to copy into */
)
{
int i, k, ok, ok2, ok3, n2 ;
if (is_scalar)
{
n2 = scalar ;
if (n == EMPTY)
{
n = scalar ;
}
i = scalar - 1 ;
for (k = 0 ; k < nz ; k++)
{
I [k] = i ;
}
}
else
{
/* copy double input into int vector (convert to 0-based) */
ok = TRUE ;
ok2 = TRUE ;
ok3 = TRUE ;
n2 = 0 ;
switch (category)
{
case mxLOGICAL_CLASS:
for (k = 0 ; k < nz ; k++)
{
mxLogical y = ((mxLogical *) vector) [k] ;
i = (int) y ;
I [k] = i - 1 ;
n2 = MAX (n2, i) ;
}
break ;
case mxCHAR_CLASS:
for (k = 0 ; k < nz ; k++)
{
mxChar y = ((mxChar *) vector) [k] ;
i = (int) y ;
I [k] = i - 1 ;
n2 = MAX (n2, i) ;
}
break ;
case mxINT8_CLASS:
for (k = 0 ; k < nz ; k++)
{
char y = ((char *) vector) [k] ;
i = (int) y ;
I [k] = i - 1 ;
n2 = MAX (n2, i) ;
}
break ;
case mxUINT8_CLASS:
for (k = 0 ; k < nz ; k++)
{
unsigned char y = ((unsigned char *) vector) [k] ;
i = (int) y ;
I [k] = i - 1 ;
n2 = MAX (n2, i) ;
}
break ;
case mxINT16_CLASS:
for (k = 0 ; k < nz ; k++)
{
short y = ((short *) vector) [k] ;
i = (int) y ;
I [k] = i - 1 ;
n2 = MAX (n2, i) ;
}
break ;
case mxUINT16_CLASS:
for (k = 0 ; k < nz ; k++)
{
unsigned short y = ((unsigned short *) vector) [k] ;
i = (int) y ;
I [k] = i - 1 ;
n2 = MAX (n2, i) ;
}
break ;
case mxINT32_CLASS:
for (k = 0 ; k < nz ; k++)
{
int y = ((int *) vector) [k] ;
I [k] = y - 1 ;
n2 = MAX (n2, y) ;
}
break ;
case mxUINT32_CLASS:
for (k = 0 ; ok3 && k < nz ; k++)
{
unsigned int y = ((unsigned int *) vector) [k] ;
i = (int) y ;
ok3 = (y < INT_MAX) ;
I [k] = i - 1 ;
n2 = MAX (n2, i) ;
}
break ;
case mxINT64_CLASS:
for (k = 0 ; ok2 && ok3 && k < nz ; k++)
{
INT64_T y = ((INT64_T *) vector) [k] ;
i = (int) y ;
ok2 = (y > 0) ;
ok3 = (y < INT_MAX) ;
I [k] = i - 1 ;
n2 = MAX (n2, i) ;
}
break ;
case mxUINT64_CLASS:
for (k = 0 ; ok2 && ok3 && k < nz ; k++)
{
unsigned INT64_T y = ((unsigned INT64_T *) vector) [k] ;
i = (int) y ;
ok2 = (y > 0) ;
ok3 = (y < INT_MAX) ;
I [k] = i - 1 ;
n2 = MAX (n2, i) ;
}
break ;
case mxSINGLE_CLASS:
for (k = 0 ; ok && ok2 && ok3 && k < nz ; k++)
{
float y = ((float *) vector) [k] ;
i = (int) y ;
ok = (y == (float) i) ;
ok2 = (y > 0) ;
ok3 = (y < INT_MAX) ;
I [k] = i - 1 ;
n2 = MAX (n2, i) ;
}
break ;
case mxDOUBLE_CLASS:
for (k = 0 ; ok && ok2 && ok3 && k < nz ; k++)
{
double y = ((double *) vector) [k] ;
i = (int) y ;
ok = (y == (double) i) ;
ok2 = (y > 0) ;
ok3 = (y < INT_MAX) ;
I [k] = i - 1 ;
n2 = MAX (n2, i) ;
}
break ;
default:
sputil_error (ERROR_INVALID_TYPE, FALSE) ;
break ;
}
if (!ok)
{
sputil_error (ERROR_NOT_INTEGER, TRUE) ;
}
if (!ok2)
{
sputil_error (ERROR_TOO_SMALL, TRUE) ;
}
if (!ok3)
{
sputil_error (ERROR_HUGE, TRUE) ;
}
}
return ((n == EMPTY) ? n2 : n) ;
}
/* ========================================================================== */
/* === sputil_dense_to_sparse =============================================== */
/* ========================================================================== */
/* Convert a dense matrix of any numeric type into a
* sparse double or sparse logical matrix.
*/
#define COUNT_NZ \
{ \
for (j = 0 ; j < ncol ; j++) \
{ \
for (i = 0 ; i < nrow ; i++) \
{ \
xij = X [i + j*nrow] ; \
if (CHOLMOD_IS_NONZERO (xij)) \
{ \
nz++ ; \
} \
} \
} \
}
#define COPY_DENSE_TO_SPARSE(stype) \
{ \
stype *Sx ; \
Sp = mxGetJc (S) ; \
Si = mxGetIr (S) ; \
Sx = (stype *) mxGetData (S) ; \
nz = 0 ; \
for (j = 0 ; j < ncol ; j++) \
{ \
Sp [j] = nz ; \
for (i = 0 ; i < nrow ; i++) \
{ \
xij = X [i + j*nrow] ; \
if (CHOLMOD_IS_NONZERO (xij)) \
{ \
Si [nz] = i ; \
Sx [nz] = (stype) xij ; \
nz++ ; \
} \
} \
} \
Sp [ncol] = nz ; \
}
#define DENSE_TO_SPARSE(type) \
{ \
type *X, xij ; \
X = (type *) mxGetData (arg) ; \
COUNT_NZ ; \
S = mxCreateSparse (nrow, ncol, nz, mxREAL) ; \
COPY_DENSE_TO_SPARSE (double) ; \
}
mxArray *sputil_dense_to_sparse (const mxArray *arg)
{
mxArray *S = NULL ;
int *Sp, *Si ;
int nrow, ncol, nz, i, j ;
nrow = mxGetM (arg) ;
ncol = mxGetN (arg) ;
nz = 0 ;
if (mxIsComplex (arg))
{
/* ------------------------------------------------------------------ */
/* convert a complex dense matrix into a complex sparse matrix */
/* ------------------------------------------------------------------ */
double xij, zij ;
double *X, *Z, *Sx, *Sz ;
if (mxGetClassID (arg) != mxDOUBLE_CLASS)
{
/* A complex matrix can have any class (int8, int16, single, etc),
* but this function only supports complex double. This condition
* is not checked in the caller. */
sputil_error (ERROR_INVALID_TYPE, FALSE) ;
}
X = mxGetPr (arg) ;
Z = mxGetPi (arg) ;
for (j = 0 ; j < ncol ; j++)
{
for (i = 0 ; i < nrow ; i++)
{
xij = X [i + j*nrow] ;
zij = Z [i + j*nrow] ;
if (CHOLMOD_IS_NONZERO (xij) || CHOLMOD_IS_NONZERO (zij))
{
nz++ ;
}
}
}
S = mxCreateSparse (nrow, ncol, nz, mxCOMPLEX) ;
Sp = mxGetJc (S) ;
Si = mxGetIr (S) ;
Sx = mxGetPr (S) ;
Sz = mxGetPi (S) ;
nz = 0 ;
for (j = 0 ; j < ncol ; j++)
{
Sp [j] = nz ;
for (i = 0 ; i < nrow ; i++)
{
xij = X [i + j*nrow] ;
zij = Z [i + j*nrow] ;
if (CHOLMOD_IS_NONZERO (xij) || CHOLMOD_IS_NONZERO (zij))
{
Si [nz] = i ;
Sx [nz] = xij ;
Sz [nz] = zij ;
nz++ ;
}
}
}
Sp [ncol] = nz ;
}
else
{
/* ------------------------------------------------------------------ */
/* convert real matrix (any class) to sparse double or logical */
/* ------------------------------------------------------------------ */
switch (mxGetClassID (arg))
{
case mxLOGICAL_CLASS:
{
mxLogical *X, xij ;
X = (mxLogical *) mxGetData (arg) ;
COUNT_NZ ;
S = mxCreateSparseLogicalMatrix (nrow, ncol, nz) ;
COPY_DENSE_TO_SPARSE (mxLogical) ;
}
break ;
case mxCHAR_CLASS:
DENSE_TO_SPARSE (mxChar) ;
break ;
case mxINT8_CLASS:
DENSE_TO_SPARSE (char) ;
break ;
case mxUINT8_CLASS:
DENSE_TO_SPARSE (unsigned char) ;
break ;
case mxINT16_CLASS:
DENSE_TO_SPARSE (short) ;
break ;
case mxUINT16_CLASS:
DENSE_TO_SPARSE (unsigned short) ;
break ;
case mxINT32_CLASS:
DENSE_TO_SPARSE (int) ;
break ;
case mxUINT32_CLASS:
DENSE_TO_SPARSE (unsigned int) ;
break ;
case mxINT64_CLASS:
DENSE_TO_SPARSE (INT64_T) ;
break ;
case mxUINT64_CLASS:
DENSE_TO_SPARSE (unsigned INT64_T) ;
break ;
case mxSINGLE_CLASS:
DENSE_TO_SPARSE (float) ;
break ;
case mxDOUBLE_CLASS:
DENSE_TO_SPARSE (double) ;
break ;
default:
sputil_error (ERROR_INVALID_TYPE, FALSE) ;
break ;
}
}
return (S) ;
}
/* ========================================================================== */
/* === sputil_triplet_to_sparse ============================================= */
/* ========================================================================== */
/* Convert a triplet form into a sparse matrix. If complex, s must be double.
* If real, s can be of any class. Optionally creates a
*/
cholmod_sparse *sputil_triplet_to_sparse
(
int nrow, int ncol, int nz, int nzmax,
int i_is_scalar, int i, void *i_vector, mxClassID i_class,
int j_is_scalar, int j, void *j_vector, mxClassID j_class,
int s_is_scalar, double x, double z, void *x_vector, double *z_vector,
mxClassID s_class, int s_complex,
cholmod_common *cm
)
{
double dummy = 0 ;
cholmod_triplet *T ;
cholmod_sparse *S ;
double *Tx, *Tz ;
int *Ti, *Tj ;
int k, x_patch ;
/* ---------------------------------------------------------------------- */
/* allocate the triplet form */
/* ---------------------------------------------------------------------- */
/* Note that nrow and ncol may be EMPTY; this is not an error condition.
* Allocate the numerical part of T only if s is a scalar. */
x_patch = (!s_is_scalar && (s_class == mxDOUBLE_CLASS || s_complex)) ;
T = cholmod_allocate_triplet (MAX (0,nrow), MAX (0,ncol), nz, 0,
x_patch ? CHOLMOD_PATTERN :
(s_complex ? CHOLMOD_ZOMPLEX : CHOLMOD_REAL), cm) ;
Ti = T->i ;
Tj = T->j ;
Tx = T->x ;
Tz = T->z ;
/* ---------------------------------------------------------------------- */
/* fill the triplet form */
/* ---------------------------------------------------------------------- */
if (s_is_scalar)
{
/* ------------------------------------------------------------------ */
/* fill T->x and T->z with a scalar value */
/* ------------------------------------------------------------------ */
for (k = 0 ; k < nz ; k++)
{
Tx [k] = x ;
}
if (s_complex)
{
for (k = 0 ; k < nz ; k++)
{
Tz [k] = z ;
}
}
}
else
{
/* ------------------------------------------------------------------ */
/* copy x/z_vector into T->x and T->z, and convert to double */
/* ------------------------------------------------------------------ */
if (s_complex)
{
/* Patch in s as the numerical values of the triplet matrix.
* Note that T->x and T->z must not be free'd when done. */
T->x = (x_vector == NULL) ? &dummy : x_vector ;
T->z = (z_vector == NULL) ? &dummy : z_vector ;
T->xtype = CHOLMOD_ZOMPLEX ;
}
else switch (s_class)
{
case mxLOGICAL_CLASS:
for (k = 0 ; k < nz ; k++)
{
mxLogical y = ((mxLogical *) x_vector) [k] ;
Tx [k] = (double) y ;
}
break ;
case mxCHAR_CLASS:
for (k = 0 ; k < nz ; k++)
{
mxChar y = ((mxChar *) x_vector) [k] ;
Tx [k] = (double) y ;
}
break ;
case mxINT8_CLASS:
for (k = 0 ; k < nz ; k++)
{
char y = ((char *) x_vector) [k] ;
Tx [k] = (double) y ;
}
break ;
case mxUINT8_CLASS:
for (k = 0 ; k < nz ; k++)
{
unsigned char y = ((unsigned char *) x_vector) [k] ;
Tx [k] = (double) y ;
}
break ;
case mxINT16_CLASS:
for (k = 0 ; k < nz ; k++)
{
short y = ((short *) x_vector) [k] ;
Tx [k] = (double) y ;
}
break ;
case mxUINT16_CLASS:
for (k = 0 ; k < nz ; k++)
{
unsigned short y = ((unsigned short *) x_vector) [k] ;
Tx [k] = (double) y ;
}
break ;
case mxINT32_CLASS:
for (k = 0 ; k < nz ; k++)
{
int y = ((int *) x_vector) [k] ;
Tx [k] = (double) y ;
}
break ;
case mxUINT32_CLASS:
for (k = 0 ; k < nz ; k++)
{
unsigned int y = ((unsigned int *) x_vector) [k] ;
Tx [k] = (double) y ;
}
break ;
case mxINT64_CLASS:
for (k = 0 ; k < nz ; k++)
{
INT64_T y = ((INT64_T *) x_vector) [k] ;
Tx [k] = (double) y ;
}
break ;
case mxUINT64_CLASS:
for (k = 0 ; k < nz ; k++)
{
unsigned INT64_T y = ((unsigned INT64_T *) x_vector)[k];
Tx [k] = (double) y ;
}
break ;
case mxSINGLE_CLASS:
for (k = 0 ; k < nz ; k++)
{
float y = ((float *) x_vector) [k] ;
Tx [k] = (double) y ;
}
break ;
case mxDOUBLE_CLASS:
/* Patch in s as the numerical values of the triplet matrix.
* Note that T->x must not be free'd when done. */
T->x = (x_vector == NULL) ? &dummy : x_vector ;
T->xtype = CHOLMOD_REAL ;
break ;
default:
sputil_error (ERROR_INVALID_TYPE, FALSE) ;
break ;
}
}
/* copy i in to the integer vector T->i */
nrow = sputil_copy_ij (i_is_scalar, i, i_vector, i_class, nz, nrow, Ti) ;
/* copy j in to the integer vector T->j */
ncol = sputil_copy_ij (j_is_scalar, j, j_vector, j_class, nz, ncol, Tj) ;
/* nrow and ncol are known */
T->nrow = nrow ;
T->ncol = ncol ;
T->nnz = nz ;
/* ---------------------------------------------------------------------- */
/* convert triplet to sparse matrix */
/* ---------------------------------------------------------------------- */
/* If the triplet matrix T is invalid, or if CHOLMOD runs out of memory,
* then S is NULL. */
S = cholmod_triplet_to_sparse (T, nzmax, cm) ;
/* ---------------------------------------------------------------------- */
/* free workspace */
/* ---------------------------------------------------------------------- */
/* do not free T->x or T->z if it points to input x_vector */
if (x_patch)
{
T->x = NULL ;
T->z = NULL ;
T->xtype = CHOLMOD_PATTERN ;
}
cholmod_free_triplet (&T, cm) ;
return (S) ;
}
/* ========================================================================== */
/* === sputil_copy_sparse =================================================== */
/* ========================================================================== */
/* copy a sparse matrix, S = sparse(A), dropping any zero entries and ensuring
* the nzmax(S) == nnz(S). Explicit zero entries in A "cannot" occur, in
* the current version of MATLAB ... but a user mexFunction might generate a
* matrix with explicit zeros. This function ensures S=sparse(A) drops those
* explicit zeros. */
mxArray *sputil_copy_sparse (const mxArray *A)
{
double aij, zij ;
mxArray *S ;
double *Ax, *Az, *Sx, *Sz ;
mxLogical *Al, *Sl ;
int *Ap, *Ai, *Sp, *Si ;
int anz, snz, p, j, nrow, ncol, pend ;
if (! (mxGetClassID (A) == mxLOGICAL_CLASS
|| mxGetClassID (A) == mxDOUBLE_CLASS))
{
/* Only sparse logical and real/complex double matrices supported.
* This condition is not checked in the caller. */
sputil_error (ERROR_INVALID_TYPE, 0) ;
}
nrow = mxGetM (A) ;
ncol = mxGetN (A) ;
Ap = mxGetJc (A) ;
Ai = mxGetIr (A) ;
anz = Ap [ncol] ;
snz = 0 ;
if (mxIsLogical (A))
{
/* ------------------------------------------------------------------ */
/* copy a sparse logical matrix */
/* ------------------------------------------------------------------ */
/* count the number of nonzeros in A */
Al = mxGetLogicals (A) ;
for (p = 0 ; p < anz ; p++)
{
if (Al [p])
{
snz++ ;
}
}
/* allocate S */
S = mxCreateSparseLogicalMatrix (nrow, ncol, snz) ;
Sp = mxGetJc (S) ;
Si = mxGetIr (S) ;
Sl = mxGetLogicals (S) ;
/* copy A into S, dropping zero entries */
snz = 0 ;
for (j = 0 ; j < ncol ; j++)
{
Sp [j] = snz ;
pend = Ap [j+1] ;
for (p = Ap [j] ; p < pend ; p++)
{
if (Al [p])
{
Si [snz] = Ai [p] ;
Sl [snz] = 1 ;
snz++ ;
}
}
}
}
else if (mxIsComplex (A))
{
/* ------------------------------------------------------------------ */
/* copy a sparse complex double matrix */
/* ------------------------------------------------------------------ */
/* count the number of nonzeros in A */
Ax = mxGetPr (A) ;
Az = mxGetPi (A) ;
for (p = 0 ; p < anz ; p++)
{
aij = Ax [p] ;
zij = Az [p] ;
if (CHOLMOD_IS_NONZERO (aij) || CHOLMOD_IS_NONZERO (zij))
{
snz++ ;
}
}
/* allocate S */
S = mxCreateSparse (nrow, ncol, snz, mxCOMPLEX) ;
Sp = mxGetJc (S) ;
Si = mxGetIr (S) ;
Sx = mxGetPr (S) ;
Sz = mxGetPi (S) ;
/* copy A into S, dropping zero entries */
snz = 0 ;
for (j = 0 ; j < ncol ; j++)
{
Sp [j] = snz ;
pend = Ap [j+1] ;
for (p = Ap [j] ; p < pend ; p++)
{
aij = Ax [p] ;
zij = Az [p] ;
if (CHOLMOD_IS_NONZERO (aij) || CHOLMOD_IS_NONZERO (zij))
{
Si [snz] = Ai [p] ;
Sx [snz] = aij ;
Sz [snz] = zij ;
snz++ ;
}
}
}
}
else
{
/* ------------------------------------------------------------------ */
/* copy a sparse real double matrix */
/* ------------------------------------------------------------------ */
/* count the number of nonzeros in A */
Ax = mxGetPr (A) ;
for (p = 0 ; p < anz ; p++)
{
aij = Ax [p] ;
if (CHOLMOD_IS_NONZERO (aij))
{
snz++ ;
}
}
/* allocate S */
S = mxCreateSparse (nrow, ncol, snz, mxREAL) ;
Sp = mxGetJc (S) ;
Si = mxGetIr (S) ;
Sx = mxGetPr (S) ;
/* copy A into S, dropping zero entries */
snz = 0 ;
for (j = 0 ; j < ncol ; j++)
{
Sp [j] = snz ;
pend = Ap [j+1] ;
for (p = Ap [j] ; p < pend ; p++)
{
aij = Ax [p] ;
if (CHOLMOD_IS_NONZERO (aij))
{
Si [snz] = Ai [p] ;
Sx [snz] = aij ;
snz++ ;
}
}
}
}
Sp [ncol] = snz ;
return (S) ;
}
/* ========================================================================== */
/* === sputil_sparse_to_dense =============================================== */
/* ========================================================================== */
/* convert a sparse double or logical array to a dense double array */
mxArray *sputil_sparse_to_dense (const mxArray *S)
{
mxArray *X ;
mxLogical *Sl ;
double *Sx, *Sz, *Xx, *Xz ;
int *Sp, *Si ;
int nrow, ncol, i, j, p, pend, j2 ;
if (! (mxGetClassID (S) == mxLOGICAL_CLASS
&& mxGetClassID (S) == mxDOUBLE_CLASS))
{
/* only sparse logical and real/complex double matrices supported */
sputil_error (ERROR_INVALID_TYPE, 0) ;
}
nrow = mxGetM (S) ;
ncol = mxGetN (S) ;
Sp = mxGetJc (S) ;
Si = mxGetIr (S) ;
if (mxIsLogical (S))
{
/* logical */
Sl = (mxLogical *) mxGetData (S) ;
X = mxCreateDoubleMatrix (nrow, ncol, mxREAL) ;
Xx = mxGetPr (X) ;
for (j = 0 ; j < ncol ; j++)
{
pend = Sp [j+1] ;
j2 = j*nrow ;
for (p = Sp [j] ; p < pend ; p++)
{
Xx [Si [p] + j2] = (double) (Sl [p]) ;
}
}
}
else if (mxIsComplex (S))
{
/* complex */
Sx = mxGetPr (S) ;
Sz = mxGetPi (S) ;
X = mxCreateDoubleMatrix (nrow, ncol, mxCOMPLEX) ;
Xx = mxGetPr (X) ;
Xz = mxGetPi (X) ;
for (j = 0 ; j < ncol ; j++)
{
pend = Sp [j+1] ;
j2 = j*nrow ;
for (p = Sp [j] ; p < pend ; p++)
{
i = Si [p] ;
Xx [i + j2] = Sx [p] ;
Xz [i + j2] = Sz [p] ;
}
}
}
else
{
/* real */
Sx = mxGetPr (S) ;
X = mxCreateDoubleMatrix (nrow, ncol, mxREAL) ;
Xx = mxGetPr (X) ;
for (j = 0 ; j < ncol ; j++)
{
pend = Sp [j+1] ;
j2 = j*nrow ;
for (p = Sp [j] ; p < pend ; p++)
{
Xx [Si [p] + j2] = Sx [p] ;
}
}
}
return (X) ;
}
/* ========================================================================== */
/* === sputil_check_ijvector ================================================ */
/* ========================================================================== */
/* Check a sparse i or j input argument */
void sputil_check_ijvector (const mxArray *arg)
{
if (mxIsComplex (arg))
{
/* i and j cannot be complex */
sputil_error (ERROR_NOT_INTEGER, TRUE) ;
}
if (mxIsSparse (arg))
{
/* the i and j arguments for sparse(i,j,s,...) can be sparse, but if so
* they must have no zero entries. */
double mn, m, nz ;
int *p, n ;
m = (double) mxGetM (arg) ;
n = mxGetN (arg) ;
mn = m*n ;
p = mxGetJc (arg) ;
nz = p [n] ;
if (mn != nz)
{
/* i or j contains at least one zero, which is invalid */
sputil_error (ERROR_TOO_SMALL, TRUE) ;
}
}
}
/* ========================================================================== */
/* === sputil_sparse ======================================================== */
/* ========================================================================== */
/* Implements the sparse2 mexFunction */
void sputil_sparse
(
int nargout,
mxArray *pargout [ ],
int nargin,
const mxArray *pargin [ ]
)
{
double x, z ;
double *z_vector ;
void *i_vector, *j_vector, *x_vector ;
mxArray *s_array ;
cholmod_sparse *S, *Z ;
cholmod_common Common, *cm ;
int nrow, ncol, k, nz, i_is_scalar, j_is_scalar, s_is_sparse,
s_is_scalar, ilen, jlen, slen, nzmax, i, j, s_complex ;
mxClassID i_class, j_class, s_class ;
/* ---------------------------------------------------------------------- */
/* start CHOLMOD and set defaults */
/* ---------------------------------------------------------------------- */
cm = &Common ;
cholmod_start (cm) ;
sputil_config (SPUMONI, cm) ;
/* ---------------------------------------------------------------------- */
/* get inputs */
/* ---------------------------------------------------------------------- */
if (nargout > 2 || nargin > 6 || nargin == 4 || nargin == 0)
{
sputil_error (ERROR_USAGE, FALSE) ;
}
/* ---------------------------------------------------------------------- */
/* convert inputs into a sparse matrix S */
/* ---------------------------------------------------------------------- */
S = NULL ;
Z = NULL ;
if (nargin == 1)
{
/* ------------------------------------------------------------------ */
/* S = sparse (A) where A is sparse or full */
/* ------------------------------------------------------------------ */
if (mxIsSparse (pargin [0]))
{
/* -------------------------------------------------------------- */
/* S = sparse (A) where A is sparse (double, complex, or logical) */
/* -------------------------------------------------------------- */
pargout [0] = sputil_copy_sparse (pargin [0]) ;
}
else
{
/* -------------------------------------------------------------- */
/* S = sparse (A) where A is full (real or complex) */
/* -------------------------------------------------------------- */
/* A can be of any numeric type (mxLogical, int8, ..., double),
* except that if A is complex, it must also be double. */
pargout [0] = sputil_dense_to_sparse (pargin [0]) ;
}
}
else if (nargin == 2)
{
/* ------------------------------------------------------------------ */
/* S = sparse (m,n) */
/* ------------------------------------------------------------------ */
int *Sp ;
nrow = sputil_get_integer (pargin [0], FALSE, 0) ;
ncol = sputil_get_integer (pargin [1], FALSE, 0) ;
pargout [0] = mxCreateSparse (nrow, ncol, 1, mxREAL) ;
Sp = mxGetJc (pargout [0]) ;
Sp [0] = 0 ;
}
else
{
/* ------------------------------------------------------------------ */
/* S = sparse (i,j,s), sparse (i,j,s,m,n) or sparse (i,j,s,m,n,nzmax) */
/* ------------------------------------------------------------------ */
/* i, j, and s can be of any numeric type */
/* ensure i and j are valid (i and j cannot be complex) */
sputil_check_ijvector (pargin [0]) ;
sputil_check_ijvector (pargin [1]) ;
/* convert s from sparse to dense (double), if necessary */
s_is_sparse = mxIsSparse (pargin [2]) ;
if (s_is_sparse)
{
/* s must be double (real/complex) or logical */
s_array = sputil_sparse_to_dense (pargin [2]) ;
}
else
{
s_array = (mxArray *) pargin [2] ;
}
/* s is now full. It can be any class, except if complex it must also
* be double */
s_class = mxGetClassID (s_array) ;
s_complex = mxIsComplex (s_array) ;
if (s_complex && s_class != mxDOUBLE_CLASS)
{
/* for complex case, only double class is supported */
sputil_error (ERROR_INVALID_TYPE, 0) ;
}
/* get sizes of inputs */
ilen = sputil_nelements (pargin [0]) ;
jlen = sputil_nelements (pargin [1]) ;
slen = sputil_nelements (s_array) ;
/* if i, j, s are scalars, they "float" to sizes of non-scalar args */
i_is_scalar = (ilen == 1) ;
j_is_scalar = (jlen == 1) ;
s_is_scalar = (slen == 1) ;
/* find the length */
if (!i_is_scalar)
{
/* if i is not a scalar, let it determine the length */
nz = ilen ;
}
else if (!j_is_scalar)
{
/* otherwise, if j is not a scalar, let it determine the length */
nz = jlen ;
}
else
{
/* finally, i and j are both scalars, so let s determine length */
nz = slen ;
}
/* make sure the sizes are compatible */
if (!((i_is_scalar || ilen == nz) &&
(j_is_scalar || jlen == nz) &&
(s_is_scalar || slen == nz)))
{
sputil_error (ERROR_LENGTH, FALSE) ;
}
if (nargin > 4)
{
nrow = sputil_get_integer (pargin [3], FALSE, 0) ;
ncol = sputil_get_integer (pargin [4], FALSE, 0) ;
}
else
{
/* nrow and ncol will be discovered by scanning i and j */
nrow = EMPTY ;
ncol = EMPTY ;
}
if (nargin > 5)
{
nzmax = sputil_get_integer (pargin [5], FALSE, 0) ;
nzmax = MAX (nzmax, nz) ;
}
else
{
nzmax = nz ;
}
/* ------------------------------------------------------------------ */
/* convert triplet form to sparse form */
/* ------------------------------------------------------------------ */
i = i_is_scalar ? sputil_get_integer (pargin [0], TRUE, nrow) : 0 ;
i_vector = mxGetData (pargin [0]) ;
i_class = mxGetClassID (pargin [0]) ;
j = j_is_scalar ? sputil_get_integer (pargin [1], TRUE, ncol) : 0 ;
j_vector = mxGetData (pargin [1]) ;
j_class = mxGetClassID (pargin [1]) ;
x_vector = mxGetData (s_array) ;
z_vector = mxGetPi (s_array) ;
x = sputil_get_double (s_array) ;
z = (s_complex && z_vector != NULL) ? (z_vector [0]) : 0 ;
S = sputil_triplet_to_sparse (nrow, ncol, nz, nzmax,
i_is_scalar, i, i_vector, i_class,
j_is_scalar, j, j_vector, j_class,
s_is_scalar, x, z, x_vector, z_vector,
s_class, s_complex,
cm) ;
/* set nzmax(A) to nnz(S), unless nzmax is specified on input */
if (nargin <= 5 && S != NULL)
{
cholmod_reallocate_sparse (cholmod_nnz (S, cm), S, cm) ;
}
if (nargout > 1)
{
/* return a binary pattern of the explicit zero entries, for the
* S = sparse(i,j,x, ...) form. */
Z = sputil_extract_zeros (S, cm) ;
}
/* drop explicit zeros from S */
sputil_drop_zeros (S) ;
if (s_is_sparse)
{
mxDestroyArray (s_array) ;
}
}
/* ---------------------------------------------------------------------- */
/* convert S into a MATLAB sparse matrix */
/* ---------------------------------------------------------------------- */
k = 0 ;
if (S != NULL)
{
if (mxIsLogical (pargin [2]))
{
/* copy S into a MATLAB sparse logical matrix */
mxLogical *s_logical ;
pargout [0] = mxCreateSparseLogicalMatrix (0, 0, 0) ;
s_logical = cholmod_malloc (S->nzmax, sizeof (mxLogical), cm) ;
for (k = 0 ; k < (int) (S->nzmax) ; k++)
{
s_logical [k] = 1 ;
}
mxFree (mxGetData (pargout [0])) ;
mxSetData (pargout [0], s_logical) ;
mexMakeMemoryPersistent (s_logical) ;
k++ ;
}
else if (mxIsComplex (pargin [2]))
{
/* copy S into a MATLAB sparse complex double matrix */
pargout [0] = mxCreateSparse (0, 0, 0, mxCOMPLEX) ;
mxFree (mxGetPr (pargout [0])) ;
mxFree (mxGetPi (pargout [0])) ;
mxSetPr (pargout [0], S->x) ;
mxSetPi (pargout [0], S->z) ;
mexMakeMemoryPersistent (S->x) ;
mexMakeMemoryPersistent (S->z) ;
k += 2 ;
S->x = NULL ;
S->z = NULL ;
}
else
{
/* copy S into a MATLAB sparse real double matrix */
pargout [0] = mxCreateSparse (0, 0, 0, mxREAL) ;
mxSetPr (pargout [0], S->x) ;
mexMakeMemoryPersistent (S->x) ;
k++ ;
S->x = NULL ;
}
mxSetM (pargout [0], S->nrow) ;
mxSetN (pargout [0], S->ncol) ;
mxSetNzmax (pargout [0], S->nzmax) ;
mxFree (mxGetJc (pargout [0])) ;
mxFree (mxGetIr (pargout [0])) ;
mxSetJc (pargout [0], S->p) ;
mxSetIr (pargout [0], S->i) ;
mexMakeMemoryPersistent (S->p) ;
mexMakeMemoryPersistent (S->i) ;
k += 2 ;
/* free cholmod_sparse S, except for what has been given to MATLAB */
S->p = NULL ;
S->i = NULL ;
cholmod_free_sparse (&S, cm) ;
}
/* ---------------------------------------------------------------------- */
/* return Z to MATLAB, if requested */
/* ---------------------------------------------------------------------- */
if (nargout > 1)
{
if (Z == NULL)
{
/* Z not computed; return an empty matrix */
Z = cholmod_spzeros (nrow, ncol, 0, CHOLMOD_REAL, cm) ;
}
pargout [1] = sputil_put_sparse (&Z, cm) ;
}
cholmod_finish (cm) ;
cholmod_print_common (" ", cm) ;
/*
if (cm->malloc_count != k) mexErrMsgTxt ("!") ;
*/
}
syntax highlighted by Code2HTML, v. 0.9.1