/* ========================================================================== */
/* === Check/cholmod_read =================================================== */
/* ========================================================================== */
/* -----------------------------------------------------------------------------
* CHOLMOD/Check Module. Copyright (C) 2005-2006, Timothy A. Davis.
* The CHOLMOD/Check 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
* -------------------------------------------------------------------------- */
/* Read a sparse matrix in triplet or dense form. A triplet matrix can be
* returned as compressed-column sparse matrix. The file format is compatible
* with all variations of the Matrix Market "coordinate" and "array" format
* (http://www.nist.gov/MatrixMarket). The format supported by these routines
* also allow other formats, where the Matrix Market header is optional.
*
* Although the Matrix Market header is optional, I recommend that users stick
* with the strict Matrix Market format. The optional format appears here to
* support the reading of symmetric matrices stored with just their upper
* triangular parts present, for testing and development of the A->stype > 0
* format in CHOLMOD. That format is not included in the Matrix Market format.
*
* If the first line of the file starts with %%MatrixMarket, then it is
* interpretted as a file in Matrix Market format. This line must have
* the following format:
*
* %%MatrixMarket matrix <fmt> <type> <storage>
*
* <fmt> is one of: coordinate or array. The former is a sparse matrix in
* triplet form. The latter is a dense matrix in column-major form.
*
* <type> is one of: real, complex, pattern, or integer.
* The functions here convert the "integer" and "pattern" types to real.
*
* <storage> is one of: general, hermitian, symmetric, or skew-symmetric
*
* The strings are case-insensitive. Only the first character is
* significant (or the first two for skew-symmetric).
*
* <type> is ignored for all matrices; the actual type (real, complex,
* or pattern) is inferred from the number of tokens in each line of the
* file. For a "coordinate" matrix: 2: pattern, 3: real, 4: complex; for
* a dense "array" matrix: 1: real, 2: complex. This is compatible with
* the Matrix Market format, since pattern matrices must have two tokens
* per line, real matrices must have 3, and complex matrices must have 4.
* A storage of "general" implies an stype of zero (see below).
* "symmetric" and "hermitian" imply an stype of -1. Skew-symmetric and
* complex symmetric matrices are always returned with both upper and lower
* triangular parts present, with an stype of zero, since CHOLMOD does not
* have a method for representing skew-symmetric and complex symmetric
* matrices. Real symmetric and complex Hermitian matrices may optionally
* be returned with both parts present.
*
* Any other lines starting with "%" are treated as comments, and are ignored.
* Blank lines are ignored. The Matrix Market header is optional in this
* routine (it is not optional in the Matrix Market format).
*
* Note that complex matrices are always returned in CHOLMOD_COMPLEX format,
* not CHOLMOD_ZOMPLEX.
*
* -----------------------------------------------------------------------------
* Triplet matrices:
* -----------------------------------------------------------------------------
*
* The first data line of a triplet matrix contains 3 or 4 integers:
*
* nrow ncol nnz stype
*
* where stype is optional (stype does not appear in the Matrix Market format).
* The matrix is nrow-by-ncol. The following nnz lines (excluding comments
* and blank lines) each contain a single entry. Duplicates are permitted,
* and are summed in the output matrix.
*
* The stype is first derived from the Matrix Market header. If the stype
* appears as the fourth integer in the first data line, it is determined from
* that line.
*
* If stype is present, it denotes the storage format for the matrix.
* stype = 0 denotes an unsymmetric matrix (same as Matrix Market "general").
* stype = -1 denotes a real symmetric or complex Hermitian matrix whose lower
* triangular entries are stored. Entries may be present in the upper
* triangular part, but these are ignored (same as Matrix Market
* "real symmetric" and "complex Hermitian").
* stype = 1 denotes a real symmetric or complex Hermitian matrix whose upper
* triangular entries are stored. Entries may be present in the lower
* triangular part, but these are ignored. This option is not present
* in the Matrix Market format.
*
* If stype is not present (no Matrix Market header and not in the first data
* line) it is inferred from the rest of the data. If the matrix is
* rectangular, or has entries in both the upper and lower triangular parts,
* then it is assumed to be unsymmetric (stype=0). If only entries in the
* lower triangular part are present, the matrix is assumed to have stype = -1.
* If only entries in the upper triangular part are present, the matrix is
* assumed to have stype = 1.
*
* After the first data line (with nrow, ncol, nnz, and optionally stype),
* each nonzero consists of one line with 2, 3, or 4 entries. All lines must
* have the same number of entries. The first two entries are the row and
* column indices of the nonzero. If 3 entries are present, the 3rd entry is
* the numerical value, and the matrix is real. If 4 entries are present,
* the 3rd and 4th entries in the line are the real and imaginary parts of
* a complex value.
*
* The matrix can be either 0-based or 1-based. It is first assumed to be
* one-based (all matrices in the Matrix Market are one-based), with row indices
* in the range 1 to ncol and column indices in the range 1 to nrow. If a row
* or column index of zero is found, the matrix is assumed to be zero-based
* (with row indices in the range 0 to ncol-1 and column indices in the range 0
* to nrow-1).
*
* If Common->prefer_binary is set to its default value of FALSE, then
* for symmetric pattern-only matrices, the kth diagonal (if present) is set to
* one plus the degree of the row/column k, and the off-diagonal entries are set
* to -1. A symmetric pattern-only matrix with a zero-free diagonal is thus
* converted into a symmetric positive definite matrix. All entries are set to
* one for an unsymmetric pattern-only matrix. This differs from the
* Matrix Market format (A = mmread ('file') returns a binary pattern for A for
* symmetric pattern-only matrices). If Common->prefer_binary is TRUE, then
* this function returns a binary matrix (just like mmread('file')).
*
* -----------------------------------------------------------------------------
* Dense matrices:
* -----------------------------------------------------------------------------
*
* A dense matrix is specified by the Matrix Market "array" format. The
* Matrix Market header is optional; if not present, the matrix is assumed to
* be in the Matrix Market "general" format. The first data line contains just
* two integers:
*
* nrow ncol
*
* The <type> can be real, integer, or complex (not pattern). These functions
* convert an integer type to real. The entries in the matrix are stored in
* column-major format, with one line per entry. Two entries are present in
* each line for complex matrices, one for real and integer matrices. In
* rectangular and unsymmetric matrices, all entries are present. For real
* symmetric or complex Hermitian matrices, only entries in the lower triangular
* part appear. For skew-symmetric matrices, only entries in the strictly
* lower triangular part appear.
*
* Since CHOLMOD does not have a data structure for presenting dense symmetric/
* Hermitian matrices, these functions always return a dense matrix in its
* general form, with both upper and lower parts present.
*/
#ifndef NCHECK
#include "cholmod_internal.h"
#include "cholmod_check.h"
#include <string.h>
#include <ctype.h>
/* The MatrixMarket format specificies a maximum line length of 1024 */
#define MAXLINE 1030
/* ========================================================================== */
/* === get_line ============================================================= */
/* ========================================================================== */
/* Read one line of the file, return TRUE if successful, FALSE if EOF. */
static int get_line (FILE *f, char *buf)
{
buf [0] = '\0' ;
buf [1] = '\0' ;
buf [MAXLINE] = '\0' ;
return (fgets (buf, MAXLINE, f) != NULL) ;
}
/* ========================================================================== */
/* === fix_inf ============================================================== */
/* ========================================================================== */
/* Replace huge values with +/- Inf's, since scanf and printf don't deal
* with Inf's properly.
*/
static double fix_inf (double x)
{
if ((x >= HUGE_DOUBLE) || (x <= -HUGE_DOUBLE))
{
/* treat this as +/- Inf (assume 2*x leads to overflow) */
x = 2*x ;
}
return (x) ;
}
/* ========================================================================== */
/* === is_blank_line ======================================================== */
/* ========================================================================== */
/* TRUE if s is a blank line or comment, FALSE otherwise */
static int is_blank_line
(
char *s
)
{
int c, k ;
if (s [0] == '%')
{
/* a comment line */
return (TRUE) ;
}
for (k = 0 ; k <= MAXLINE ; k++)
{
c = s [k] ;
if (c == '\0')
{
/* end of line */
break ;
}
if (!isspace (c))
{
/* non-space character */
return (FALSE) ;
}
}
return (TRUE) ;
}
/* ========================================================================== */
/* === read_header ========================================================== */
/* ========================================================================== */
/* Read the header. This consists of zero or more comment lines (blank, or
* starting with a "%" in the first column), followed by a single data line
* containing up to four numerical values.
*
* The first line may optionally be a Matrix Market header line, of the form
*
* %%MatrixMarket matrix <fmt> <type> <storage>
*
* The first data line of a sparse matrix in triplet form consists of 3 or 4
* numerical values:
*
* nrow ncol nnz stype
*
* where stype is optional (it does not appear in the Matrix Market file
* format). The first line of a dense matrix in column-major form consists of
* two numerical values:
*
* nrow ncol
*
* The stype of the matrix is determine either from the Matrix Market header,
* or (optionally) from the first data line. stypes of 0 to -3 directly
* correlate with the Matrix Market format; stype = 1 is an extension to that
* format.
*
* 999: unknown (will be inferred from the data)
* 1: real symmetric or complex Hermitian with upper part stored
* (not in the Matrix Market format)
* 0: unsymmetric (same as Matrix Market "general")
* -1: real symmetric or complex Hermitian, with lower part stored
* (Matrix Market "real symmetric" or "complex hermitian")
* -2: real or complex skew symmetric (lower part stored, can only be
* specified by Matrix Market header)
* -3: complex symmetric (lower part stored)
* specified by Matrix Market header)
*
* The Matrix Market header is optional. If stype appears in the first data
* line, it is determine by that data line. Otherwise, if the Matrix Market
* header appears, stype is determined from that header. If stype does not
* appear, it is set to "unknown" (999).
*/
#define STYPE_UNKNOWN 999
#define STYPE_SYMMETRIC_UPPER 1
#define STYPE_UNSYMMETRIC 0
#define STYPE_SYMMETRIC_LOWER -1
#define STYPE_SKEW_SYMMETRIC -2
#define STYPE_COMPLEX_SYMMETRIC_LOWER -3
static int read_header /* returns TRUE if successful, FALSE on error */
(
/* ---- input ---- */
FILE *f, /* file to read from */
/* ---- output --- */
char *buf, /* a character array of size MAXLINE+1 */
int *mtype, /* CHOLMOD_TRIPLET or CHOLMOD_DENSE */
size_t *nrow, /* number of rows in the matrix */
size_t *ncol, /* number of columns in the matrix */
size_t *nnz, /* number of entries in a triplet matrix (0 for dense)*/
int *stype /* stype (see above) */
)
{
char *p ;
int first = TRUE, got_mm_header = FALSE, c, c2, is_complex, nitems ;
double l1, l2, l3, l4 ;
*mtype = CHOLMOD_TRIPLET ;
*nrow = 0 ;
*ncol = 0 ;
*nnz = 0 ;
*stype = STYPE_UNKNOWN ;
for ( ; ; )
{
/* ------------------------------------------------------------------ */
/* get the next line */
/* ------------------------------------------------------------------ */
if (!get_line (f, buf))
{
/* premature end of file */
return (FALSE) ;
}
if (first && (strncmp (buf, "%%MatrixMarket", 14) == 0))
{
/* -------------------------------------------------------------- */
/* read a Matrix Market header */
/* -------------------------------------------------------------- */
got_mm_header = TRUE ;
p = buf ;
/* -------------------------------------------------------------- */
/* get "matrix" token */
/* -------------------------------------------------------------- */
while (*p && !isspace (*p)) p++ ;
while (*p && isspace (*p)) p++ ;
c = tolower (*p) ;
if (c != 'm')
{
/* bad format */
return (FALSE) ;
}
/* -------------------------------------------------------------- */
/* get the fmt token ("coord" or "array") */
/* -------------------------------------------------------------- */
while (*p && !isspace (*p)) p++ ;
while (*p && isspace (*p)) p++ ;
c = tolower (*p) ;
if (c == 'c')
{
*mtype = CHOLMOD_TRIPLET ;
}
else if (c == 'a')
{
*mtype = CHOLMOD_DENSE ;
}
else
{
/* bad format, neither "coordinate" nor "array" */
return (FALSE) ;
}
/* -------------------------------------------------------------- */
/* get type token (real, pattern, complex, integer) */
/* -------------------------------------------------------------- */
while (*p && !isspace (*p)) p++ ;
while (*p && isspace (*p)) p++ ;
c = tolower (*p) ;
if (!(c == 'r' || c == 'p' || c == 'c' || c == 'i'))
{
/* bad format */
return (FALSE) ;
}
is_complex = (c == 'c') ;
/* -------------------------------------------------------------- */
/* get storage token (general, hermitian, symmetric, skew) */
/* -------------------------------------------------------------- */
while (*p && !isspace (*p)) p++ ;
while (*p && isspace (*p)) p++ ;
c = tolower (*p) ;
c2 = tolower (*(p+1)) ;
if (c == 'g')
{
/* "general" storage (unsymmetric matrix), both parts present */
*stype = STYPE_UNSYMMETRIC ;
}
else if (c == 's' && c2 == 'y')
{
/* "symmetric" */
if (is_complex)
{
/* complex symmetric, lower triangular part present */
*stype = STYPE_COMPLEX_SYMMETRIC_LOWER ;
}
else
{
/* real symmetric, lower triangular part present */
*stype = STYPE_SYMMETRIC_LOWER ;
}
}
else if (c == 'h')
{
/* "hermitian" matrix, lower triangular part present */
*stype = STYPE_SYMMETRIC_LOWER ;
}
else if (c == 's' && c2 == 'k')
{
/* "skew-symmetric" (real or complex), lower part present */
*stype = STYPE_SKEW_SYMMETRIC ;
}
else
{
/* bad format */
return (FALSE) ;
}
}
else if (is_blank_line (buf))
{
/* -------------------------------------------------------------- */
/* blank line or comment line */
/* -------------------------------------------------------------- */
continue ;
}
else
{
/* -------------------------------------------------------------- */
/* read the first data line and return */
/* -------------------------------------------------------------- */
/* format: nrow ncol nnz stype */
l1 = EMPTY ;
l2 = EMPTY ;
l3 = 0 ;
l4 = 0 ;
nitems = sscanf (buf, "%lg %lg %lg %lg\n", &l1, &l2, &l3, &l4) ;
if (nitems < 2 || nitems > 4 || l1 > Int_max || l2 > Int_max)
{
/* invalid matrix */
return (FALSE) ;
}
*nrow = l1 ;
*ncol = l2 ;
if (nitems == 2)
{
/* a dense matrix */
if (!got_mm_header)
{
*mtype = CHOLMOD_DENSE ;
*stype = STYPE_UNSYMMETRIC ;
}
}
if (nitems == 3 || nitems == 4)
{
/* a sparse triplet matrix */
*nnz = l3 ;
if (!got_mm_header)
{
*mtype = CHOLMOD_TRIPLET ;
}
}
if (nitems == 4)
{
/* an stype specified here can only be 1, 0, or -1 */
if (l4 < 0)
{
*stype = STYPE_SYMMETRIC_LOWER ;
}
else if (l4 > 0)
{
*stype = STYPE_SYMMETRIC_UPPER ;
}
else
{
*stype = STYPE_UNSYMMETRIC ;
}
}
if (*nrow != *ncol)
{
/* a rectangular matrix must be unsymmetric */
*stype = STYPE_UNSYMMETRIC ;
}
return (TRUE) ;
}
first = FALSE ;
}
}
/* ========================================================================== */
/* === read_triplet ========================================================= */
/* ========================================================================== */
/* Header has already been read in, including first line (nrow ncol nnz stype).
* Read the triplets. */
static cholmod_triplet *read_triplet
(
/* ---- input ---- */
FILE *f, /* file to read from, must already be open */
size_t nrow, /* number of rows */
size_t ncol, /* number of columns */
size_t nnz, /* number of triplets in file to read */
int stype, /* stype from header, or "unknown" */
int prefer_unsym, /* if TRUE, always return T->stype of zero */
/* ---- workspace */
char *buf, /* of size MAXLINE+1 */
/* --------------- */
cholmod_common *Common
)
{
double x, z ;
double *Tx ;
Int *Ti, *Tj, *Rdeg, *Cdeg ;
cholmod_triplet *T ;
double l1, l2 ;
Int nitems, xtype, unknown, k, nshould, is_lower, is_upper, one_based, i, j,
imax, jmax, ignore, skew_symmetric, p, complex_symmetric ;
size_t s, nnz2, extra ;
int ok = TRUE ;
/* ---------------------------------------------------------------------- */
/* quick return for empty matrix */
/* ---------------------------------------------------------------------- */
if (nrow == 0 || ncol == 0 || nnz == 0)
{
/* return an empty matrix */
return (CHOLMOD(allocate_triplet) (nrow, ncol, 0, 0, CHOLMOD_REAL,
Common)) ;
}
/* ---------------------------------------------------------------------- */
/* special stype cases: unknown, skew symmetric, and complex symmetric */
/* ---------------------------------------------------------------------- */
unknown = (stype == STYPE_UNKNOWN) ;
skew_symmetric = (stype == STYPE_SKEW_SYMMETRIC) ;
complex_symmetric = (stype == STYPE_COMPLEX_SYMMETRIC_LOWER) ;
extra = 0 ;
if (stype < STYPE_SYMMETRIC_LOWER
|| (prefer_unsym && stype != STYPE_UNSYMMETRIC))
{
/* 999: unknown might be converted to unsymmetric */
/* 1: symmetric upper converted to unsym. if prefer_unsym is TRUE */
/* -1: symmetric lower converted to unsym. if prefer_unsym is TRUE */
/* -2: real or complex skew symmetric converted to unsymmetric */
/* -3: complex symmetric converted to unsymmetric */
stype = STYPE_UNSYMMETRIC ;
extra = nnz ;
}
nnz2 = CHOLMOD(add_size_t) (nnz, extra, &ok) ;
/* ---------------------------------------------------------------------- */
/* allocate workspace */
/* ---------------------------------------------------------------------- */
/* s = nrow + ncol */
s = CHOLMOD(add_size_t) (nrow, ncol, &ok) ;
if (!ok || nrow > Int_max || ncol > Int_max || nnz > Int_max)
{
ERROR (CHOLMOD_TOO_LARGE, "problem too large") ;
return (NULL) ;
}
CHOLMOD(allocate_work) (0, s, 0, Common) ;
Rdeg = Common->Iwork ; /* size nrow */
Cdeg = Rdeg + nrow ; /* size ncol */
/* ---------------------------------------------------------------------- */
/* read the triplets */
/* ---------------------------------------------------------------------- */
is_lower = TRUE ;
is_upper = TRUE ;
one_based = TRUE ;
imax = 0 ;
jmax = 0 ;
Tx = NULL ;
Ti = NULL ;
Tj = NULL ;
xtype = 999 ;
nshould = 0 ;
for (k = 0 ; k < (Int) nnz ; k++)
{
/* ------------------------------------------------------------------ */
/* get the next triplet, skipping blank lines and comment lines */
/* ------------------------------------------------------------------ */
l1 = EMPTY ;
l2 = EMPTY ;
x = 0 ;
z = 0 ;
for ( ; ; )
{
if (!get_line (f, buf))
{
/* premature end of file - not enough triplets read in */
ERROR (CHOLMOD_INVALID, "premature EOF") ;
return (NULL) ;
}
if (is_blank_line (buf))
{
/* blank line or comment */
continue ;
}
nitems = sscanf (buf, "%lg %lg %lg %lg\n", &l1, &l2, &x, &z) ;
x = fix_inf (x) ;
z = fix_inf (z) ;
break ;
}
nitems = (nitems == EOF) ? 0 : nitems ;
i = l1 ;
j = l2 ;
/* ------------------------------------------------------------------ */
/* for first triplet: determine type and allocate triplet matrix */
/* ------------------------------------------------------------------ */
if (k == 0)
{
if (nitems < 2 || nitems > 4)
{
/* invalid matrix */
ERROR (CHOLMOD_INVALID, "invalid format") ;
return (NULL) ;
}
else if (nitems == 2)
{
/* this will be converted into a real matrix later */
xtype = CHOLMOD_PATTERN ;
}
else if (nitems == 3)
{
xtype = CHOLMOD_REAL ;
}
else if (nitems == 4)
{
xtype = CHOLMOD_COMPLEX ;
}
/* the rest of the lines should have the same number of entries */
nshould = nitems ;
/* allocate triplet matrix */
T = CHOLMOD(allocate_triplet) (nrow, ncol, nnz2, stype,
(xtype == CHOLMOD_PATTERN ? CHOLMOD_REAL : xtype), Common) ;
if (Common->status < CHOLMOD_OK)
{
/* out of memory */
return (NULL) ;
}
Ti = T->i ;
Tj = T->j ;
Tx = T->x ;
T->nnz = nnz ;
}
/* ------------------------------------------------------------------ */
/* save the entry in the triplet matrix */
/* ------------------------------------------------------------------ */
if (nitems != nshould || i < 0 || j < 0)
{
/* wrong format, premature end-of-file, or negative indices */
CHOLMOD(free_triplet) (&T, Common) ;
ERROR (CHOLMOD_INVALID, "invalid matrix file") ;
return (NULL) ;
}
Ti [k] = i ;
Tj [k] = j ;
if (i < j)
{
/* this entry is in the upper triangular part */
is_lower = FALSE ;
}
if (i > j)
{
/* this entry is in the lower triangular part */
is_upper = FALSE ;
}
if (xtype == CHOLMOD_REAL)
{
Tx [k] = x ;
}
else if (xtype == CHOLMOD_COMPLEX)
{
Tx [2*k ] = x ; /* real part */
Tx [2*k+1] = z ; /* imaginary part */
}
if (i == 0 || j == 0)
{
one_based = FALSE ;
}
imax = MAX (i, imax) ;
jmax = MAX (j, jmax) ;
}
/* ---------------------------------------------------------------------- */
/* convert to zero-based */
/* ---------------------------------------------------------------------- */
if (one_based)
{
/* input matrix is one-based; convert matrix to zero-based */
for (k = 0 ; k < (Int) nnz ; k++)
{
Ti [k]-- ;
Tj [k]-- ;
}
}
if (one_based ?
(imax > (Int) nrow || jmax > (Int) ncol) :
(imax >= (Int) nrow || jmax >= (Int) ncol))
{
/* indices out of range */
CHOLMOD(free_triplet) (&T, Common) ;
ERROR (CHOLMOD_INVALID, "indices out of range") ;
return (NULL) ;
}
/* ---------------------------------------------------------------------- */
/* determine the stype, if not yet known */
/* ---------------------------------------------------------------------- */
if (unknown)
{
if (is_lower && is_upper)
{
/* diagonal matrix, symmetric with upper part present */
stype = STYPE_SYMMETRIC_UPPER ;
}
else if (is_lower && !is_upper)
{
/* symmetric, lower triangular part present */
stype = STYPE_SYMMETRIC_LOWER ;
}
else if (!is_lower && is_upper)
{
/* symmetric, upper triangular part present */
stype = STYPE_SYMMETRIC_UPPER ;
}
else
{
/* unsymmetric */
stype = STYPE_UNSYMMETRIC ;
extra = 0 ;
}
}
/* ---------------------------------------------------------------------- */
/* add the remainder of symmetric, skew-symmetric or Hermitian matrices */
/* ---------------------------------------------------------------------- */
/* note that this step is not done for real symmetric or complex Hermitian
* matrices, unless prefer_unsym is TRUE */
if (extra > 0)
{
p = nnz ;
for (k = 0 ; k < (Int) nnz ; k++)
{
i = Ti [k] ;
j = Tj [k] ;
if (i != j)
{
Ti [p] = j ;
Tj [p] = i ;
if (xtype == CHOLMOD_REAL)
{
if (skew_symmetric)
{
Tx [p] = -Tx [k] ;
}
else
{
Tx [p] = Tx [k] ;
}
}
else if (xtype == CHOLMOD_COMPLEX)
{
if (skew_symmetric)
{
Tx [2*p ] = -Tx [2*k ] ;
Tx [2*p+1] = -Tx [2*k+1] ;
}
else if (complex_symmetric)
{
Tx [2*p ] = Tx [2*k ] ;
Tx [2*p+1] = Tx [2*k+1] ;
}
else /* Hermitian */
{
Tx [2*p ] = Tx [2*k ] ;
Tx [2*p+1] = -Tx [2*k+1] ;
}
}
p++ ;
}
}
T->nnz = p ;
nnz = p ;
}
T->stype = stype ;
/* ---------------------------------------------------------------------- */
/* create values for a pattern-only matrix */
/* ---------------------------------------------------------------------- */
if (xtype == CHOLMOD_PATTERN)
{
if (stype == STYPE_UNSYMMETRIC || Common->prefer_binary)
{
/* unsymmetric case, or binary case */
for (k = 0 ; k < (Int) nnz ; k++)
{
Tx [k] = 1 ;
}
}
else
{
/* compute the row and columm degrees (excluding the diagonal) */
for (i = 0 ; i < (Int) nrow ; i++)
{
Rdeg [i] = 0 ;
}
for (j = 0 ; j < (Int) ncol ; j++)
{
Cdeg [j] = 0 ;
}
for (k = 0 ; k < (Int) nnz ; k++)
{
i = Ti [k] ;
j = Tj [k] ;
if ((stype < 0 && i > j) || (stype > 0 && i < j))
{
/* both a(i,j) and a(j,i) appear in the matrix */
Rdeg [i]++ ;
Cdeg [j]++ ;
Rdeg [j]++ ;
Cdeg [i]++ ;
}
}
/* assign the numerical values */
for (k = 0 ; k < (Int) nnz ; k++)
{
i = Ti [k] ;
j = Tj [k] ;
Tx [k] = (i == j) ? (1 + MAX (Rdeg [i], Cdeg [j])) : (-1) ;
}
}
}
/* ---------------------------------------------------------------------- */
/* return the new triplet matrix */
/* ---------------------------------------------------------------------- */
return (T) ;
}
/* ========================================================================== */
/* === read_dense =========================================================== */
/* ========================================================================== */
/* Header has already been read in, including first line (nrow ncol).
* Read a dense matrix. */
static cholmod_dense *read_dense
(
/* ---- input ---- */
FILE *f, /* file to read from, must already be open */
size_t nrow, /* number of rows */
size_t ncol, /* number of columns */
int stype, /* stype from header */
/* ---- workspace */
char *buf, /* of size MAXLINE+1 */
/* --------------- */
cholmod_common *Common
)
{
double x, z ;
double *Xx ;
cholmod_dense *X ;
Int nitems, xtype, nshould, i, j, k, kup, first ;
/* ---------------------------------------------------------------------- */
/* quick return for empty matrix */
/* ---------------------------------------------------------------------- */
if (nrow == 0 || ncol == 0)
{
/* return an empty dense matrix */
return (CHOLMOD(zeros) (nrow, ncol, CHOLMOD_REAL, Common)) ;
}
/* ---------------------------------------------------------------------- */
/* read the entries */
/* ---------------------------------------------------------------------- */
first = TRUE ;
for (j = 0 ; j < (Int) ncol ; j++)
{
/* ------------------------------------------------------------------ */
/* get the row index of the first entry in the file for column j */
/* ------------------------------------------------------------------ */
if (stype == STYPE_UNSYMMETRIC)
{
i = 0 ;
}
else if (stype == STYPE_SKEW_SYMMETRIC)
{
i = j+1 ;
}
else /* real symmetric or complex Hermitian lower */
{
i = j ;
}
/* ------------------------------------------------------------------ */
/* get column j */
/* ------------------------------------------------------------------ */
for ( ; i < (Int) nrow ; i++)
{
/* -------------------------------------------------------------- */
/* get the next entry, skipping blank lines and comment lines */
/* -------------------------------------------------------------- */
x = 0 ;
z = 0 ;
for ( ; ; )
{
if (!get_line (f, buf))
{
/* premature end of file - not enough entries read in */
ERROR (CHOLMOD_INVALID, "premature EOF") ;
return (NULL) ;
}
if (is_blank_line (buf))
{
/* blank line or comment */
continue ;
}
nitems = sscanf (buf, "%lg %lg\n", &x, &z) ;
x = fix_inf (x) ;
z = fix_inf (z) ;
break ;
}
nitems = (nitems == EOF) ? 0 : nitems ;
/* -------------------------------------------------------------- */
/* for first entry: determine type and allocate dense matrix */
/* -------------------------------------------------------------- */
if (first)
{
first = FALSE ;
if (nitems < 1 || nitems > 2)
{
/* invalid matrix */
ERROR (CHOLMOD_INVALID, "invalid format") ;
return (NULL) ;
}
else if (nitems == 1)
{
/* a real matrix */
xtype = CHOLMOD_REAL ;
}
else if (nitems == 2)
{
/* a complex matrix */
xtype = CHOLMOD_COMPLEX ;
}
/* the rest of the lines should have same number of entries */
nshould = nitems ;
/* allocate the result */
X = CHOLMOD(zeros) (nrow, ncol, xtype, Common) ;
if (Common->status < CHOLMOD_OK)
{
/* out of memory */
return (NULL) ;
}
Xx = X->x ;
}
/* -------------------------------------------------------------- */
/* save the entry in the dense matrix */
/* -------------------------------------------------------------- */
if (nitems != nshould)
{
/* wrong format or premature end-of-file */
CHOLMOD(free_dense) (&X, Common) ;
ERROR (CHOLMOD_INVALID, "invalid matrix file") ;
return (NULL) ;
}
k = i + j*nrow ;
kup = j + i*nrow ;
if (xtype == CHOLMOD_REAL)
{
/* real matrix */
Xx [k] = x ;
if (k != kup)
{
if (stype == STYPE_SYMMETRIC_LOWER)
{
/* real symmetric matrix */
Xx [kup] = x ;
}
else if (stype == STYPE_SKEW_SYMMETRIC)
{
/* real skew symmetric matrix */
Xx [kup] = -x ;
}
}
}
else if (xtype == CHOLMOD_COMPLEX)
{
Xx [2*k ] = x ; /* real part */
Xx [2*k+1] = z ; /* imaginary part */
if (k != kup)
{
if (stype == STYPE_SYMMETRIC_LOWER)
{
/* complex Hermitian */
Xx [2*kup ] = x ; /* real part */
Xx [2*kup+1] = -z ; /* imaginary part */
}
else if (stype == STYPE_SKEW_SYMMETRIC)
{
/* complex skew symmetric */
Xx [2*kup ] = -x ; /* real part */
Xx [2*kup+1] = -z ; /* imaginary part */
}
if (stype == STYPE_COMPLEX_SYMMETRIC_LOWER)
{
/* complex symmetric */
Xx [2*kup ] = x ; /* real part */
Xx [2*kup+1] = z ; /* imaginary part */
}
}
}
}
}
/* ---------------------------------------------------------------------- */
/* return the new dense matrix */
/* ---------------------------------------------------------------------- */
return (X) ;
}
/* ========================================================================== */
/* === cholmod_read_triplet ================================================= */
/* ========================================================================== */
/* Read in a triplet matrix from a file. */
cholmod_triplet *CHOLMOD(read_triplet)
(
/* ---- input ---- */
FILE *f, /* file to read from, must already be open */
/* --------------- */
cholmod_common *Common
)
{
char buf [MAXLINE+1] ;
size_t nrow, ncol, nnz ;
int stype, mtype ;
/* ---------------------------------------------------------------------- */
/* check inputs */
/* ---------------------------------------------------------------------- */
RETURN_IF_NULL_COMMON (NULL) ;
RETURN_IF_NULL (f, NULL) ;
Common->status = CHOLMOD_OK ;
/* ---------------------------------------------------------------------- */
/* read the header and first data line */
/* ---------------------------------------------------------------------- */
if (!read_header (f, buf, &mtype, &nrow, &ncol, &nnz, &stype) ||
mtype != CHOLMOD_TRIPLET)
{
/* invalid matrix - this function can only read in a triplet matrix */
ERROR (CHOLMOD_INVALID, "invalid format") ;
return (NULL) ;
}
/* ---------------------------------------------------------------------- */
/* read the triplet matrix */
/* ---------------------------------------------------------------------- */
return (read_triplet (f, nrow, ncol, nnz, stype, FALSE, buf, Common)) ;
}
/* ========================================================================== */
/* === cholmod_read_sparse ================================================== */
/* ========================================================================== */
/* Read a sparse matrix from a file. See cholmod_read_triplet for a discussion
* of the file format.
*
* If Common->prefer_upper is TRUE (the default case), a symmetric matrix is
* returned stored in upper-triangular form (A->stype == 1).
*/
cholmod_sparse *CHOLMOD(read_sparse)
(
/* ---- input ---- */
FILE *f, /* file to read from, must already be open */
/* --------------- */
cholmod_common *Common
)
{
cholmod_sparse *A, *A2 ;
cholmod_triplet *T ;
/* ---------------------------------------------------------------------- */
/* check inputs */
/* ---------------------------------------------------------------------- */
RETURN_IF_NULL_COMMON (NULL) ;
RETURN_IF_NULL (f, NULL) ;
Common->status = CHOLMOD_OK ;
/* ---------------------------------------------------------------------- */
/* convert to a sparse matrix in compressed-column form */
/* ---------------------------------------------------------------------- */
T = CHOLMOD(read_triplet) (f, Common) ;
A = CHOLMOD(triplet_to_sparse) (T, 0, Common) ;
CHOLMOD(free_triplet) (&T, Common) ;
if (Common->prefer_upper && A != NULL && A->stype == -1)
{
/* A=A' */
A2 = CHOLMOD(transpose) (A, 2, Common) ;
CHOLMOD(free_sparse) (&A, Common) ;
A = A2 ;
}
return (A) ;
}
/* ========================================================================== */
/* === cholmod_read_dense =================================================== */
/* ========================================================================== */
/* Read a dense matrix from a file. */
cholmod_dense *CHOLMOD(read_dense)
(
/* ---- input ---- */
FILE *f, /* file to read from, must already be open */
/* --------------- */
cholmod_common *Common
)
{
char buf [MAXLINE+1] ;
size_t nrow, ncol, nnz ;
int stype, mtype ;
/* ---------------------------------------------------------------------- */
/* check inputs */
/* ---------------------------------------------------------------------- */
RETURN_IF_NULL_COMMON (NULL) ;
RETURN_IF_NULL (f, NULL) ;
Common->status = CHOLMOD_OK ;
/* ---------------------------------------------------------------------- */
/* read the header and first data line */
/* ---------------------------------------------------------------------- */
if (!read_header (f, buf, &mtype, &nrow, &ncol, &nnz, &stype) ||
mtype != CHOLMOD_DENSE)
{
/* invalid matrix - this function can only read in a dense matrix */
ERROR (CHOLMOD_INVALID, "invalid format") ;
return (NULL) ;
}
/* ---------------------------------------------------------------------- */
/* read the dense matrix */
/* ---------------------------------------------------------------------- */
return (read_dense (f, nrow, ncol, stype, buf, Common)) ;
}
/* ========================================================================== */
/* === cholmod_read_matrix ================================================== */
/* ========================================================================== */
/* Read a triplet matrix, sparse matrix or a dense matrix from a file. Returns
* a void pointer to either a cholmod_triplet, cholmod_sparse, or cholmod_dense
* object. The type of object is passed back to the caller as the mtype
* argument. */
void *CHOLMOD(read_matrix)
(
/* ---- input ---- */
FILE *f, /* file to read from, must already be open */
int prefer, /* If 0, a sparse matrix is always return as a
* cholmod_triplet form. It can have any stype
* (symmetric-lower, unsymmetric, or
* symmetric-upper).
* If 1, a sparse matrix is returned as an unsymmetric
* cholmod_sparse form (A->stype == 0), with both
* upper and lower triangular parts present.
* This is what the MATLAB mread mexFunction does,
* since MATLAB does not have an stype.
* If 2, a sparse matrix is returned with an stype of 0
* or 1 (unsymmetric, or symmetric with upper part
* stored).
* This argument has no effect for dense matrices.
*/
/* ---- output---- */
int *mtype, /* CHOLMOD_TRIPLET, CHOLMOD_SPARSE or CHOLMOD_DENSE */
/* --------------- */
cholmod_common *Common
)
{
void *G = NULL ;
cholmod_sparse *A, *A2 ;
cholmod_triplet *T ;
char buf [MAXLINE+1] ;
size_t nrow, ncol, nnz ;
int stype ;
/* ---------------------------------------------------------------------- */
/* check inputs */
/* ---------------------------------------------------------------------- */
RETURN_IF_NULL_COMMON (NULL) ;
RETURN_IF_NULL (f, NULL) ;
RETURN_IF_NULL (mtype, NULL) ;
Common->status = CHOLMOD_OK ;
/* ---------------------------------------------------------------------- */
/* read the header to determine the mtype */
/* ---------------------------------------------------------------------- */
if (!read_header (f, buf, mtype, &nrow, &ncol, &nnz, &stype))
{
/* invalid matrix */
ERROR (CHOLMOD_INVALID, "invalid format") ;
return (NULL) ;
}
/* ---------------------------------------------------------------------- */
/* read a matrix */
/* ---------------------------------------------------------------------- */
if (*mtype == CHOLMOD_TRIPLET)
{
/* read in the triplet matrix */
T = read_triplet (f, nrow, ncol, nnz, stype, prefer == 1, buf, Common) ;
if (prefer == 0)
{
/* return matrix in its original triplet form */
G = T ;
}
else
{
/* return matrix in a compressed-column form */
A = CHOLMOD(triplet_to_sparse) (T, 0, Common) ;
CHOLMOD(free_triplet) (&T, Common) ;
if (A != NULL && prefer == 2 && A->stype == -1)
{
/* convert A from symmetric-lower to symmetric-upper */
A2 = CHOLMOD(transpose) (A, 2, Common) ;
CHOLMOD(free_sparse) (&A, Common) ;
A = A2 ;
}
*mtype = CHOLMOD_SPARSE ;
G = A ;
}
}
else if (*mtype == CHOLMOD_DENSE)
{
/* return a dense matrix */
G = read_dense (f, nrow, ncol, stype, buf, Common) ;
}
return (G) ;
}
#endif
syntax highlighted by Code2HTML, v. 0.9.1