kusano 7d535a
kusano 7d535a
/*! @file creadrb.c
kusano 7d535a
 * \brief Read a matrix stored in Rutherford-Boeing format
kusano 7d535a
 *
kusano 7d535a
 * 
kusano 7d535a
 * -- SuperLU routine (version 4.0) --
kusano 7d535a
 * Lawrence Berkeley National Laboratory.
kusano 7d535a
 * June 30, 2009
kusano 7d535a
 * 
kusano 7d535a
 *
kusano 7d535a
 * Purpose
kusano 7d535a
 * =======
kusano 7d535a
 *
kusano 7d535a
 * Read a COMPLEX PRECISION matrix stored in Rutherford-Boeing format 
kusano 7d535a
 * as described below.
kusano 7d535a
 *
kusano 7d535a
 * Line 1 (A72, A8)
kusano 7d535a
 *      Col. 1 - 72   Title (TITLE)
kusano 7d535a
 *      Col. 73 - 80  Matrix name / identifier (MTRXID)
kusano 7d535a
 *
kusano 7d535a
 * Line 2 (I14, 3(1X, I13))
kusano 7d535a
 *      Col. 1 - 14   Total number of lines excluding header (TOTCRD)
kusano 7d535a
 *      Col. 16 - 28  Number of lines for pointers (PTRCRD)
kusano 7d535a
 *      Col. 30 - 42  Number of lines for row (or variable) indices (INDCRD)
kusano 7d535a
 *      Col. 44 - 56  Number of lines for numerical values (VALCRD)
kusano 7d535a
 *
kusano 7d535a
 * Line 3 (A3, 11X, 4(1X, I13))
kusano 7d535a
 *      Col. 1 - 3    Matrix type (see below) (MXTYPE)
kusano 7d535a
 *      Col. 15 - 28  Compressed Column: Number of rows (NROW)
kusano 7d535a
 *                    Elemental: Largest integer used to index variable (MVAR)
kusano 7d535a
 *      Col. 30 - 42  Compressed Column: Number of columns (NCOL)
kusano 7d535a
 *                    Elemental: Number of element matrices (NELT)
kusano 7d535a
 *      Col. 44 - 56  Compressed Column: Number of entries (NNZERO)
kusano 7d535a
 *                    Elemental: Number of variable indeces (NVARIX)
kusano 7d535a
 *      Col. 58 - 70  Compressed Column: Unused, explicitly zero
kusano 7d535a
 *                    Elemental: Number of elemental matrix entries (NELTVL)
kusano 7d535a
 *
kusano 7d535a
 * Line 4 (2A16, A20)
kusano 7d535a
 *      Col. 1 - 16   Fortran format for pointers (PTRFMT)
kusano 7d535a
 *      Col. 17 - 32  Fortran format for row (or variable) indices (INDFMT)
kusano 7d535a
 *      Col. 33 - 52  Fortran format for numerical values of coefficient matrix
kusano 7d535a
 *                    (VALFMT)
kusano 7d535a
 *                    (blank in the case of matrix patterns)
kusano 7d535a
 *
kusano 7d535a
 * The three character type field on line 3 describes the matrix type.
kusano 7d535a
 * The following table lists the permitted values for each of the three
kusano 7d535a
 * characters. As an example of the type field, RSA denotes that the matrix
kusano 7d535a
 * is real, symmetric, and assembled.
kusano 7d535a
 *
kusano 7d535a
 * First Character:
kusano 7d535a
 *      R Real matrix
kusano 7d535a
 *      C Complex matrix
kusano 7d535a
 *      I integer matrix
kusano 7d535a
 *      P Pattern only (no numerical values supplied)
kusano 7d535a
 *      Q Pattern only (numerical values supplied in associated auxiliary value
kusano 7d535a
 *        file)
kusano 7d535a
 *
kusano 7d535a
 * Second Character:
kusano 7d535a
 *      S Symmetric
kusano 7d535a
 *      U Unsymmetric
kusano 7d535a
 *      H Hermitian
kusano 7d535a
 *      Z Skew symmetric
kusano 7d535a
 *      R Rectangular
kusano 7d535a
 *
kusano 7d535a
 * Third Character:
kusano 7d535a
 *      A Compressed column form
kusano 7d535a
 *      E Elemental form
kusano 7d535a
 *
kusano 7d535a
 * 
kusano 7d535a
 */
kusano 7d535a
kusano 7d535a
#include "slu_cdefs.h"
kusano 7d535a
kusano 7d535a
kusano 7d535a
/*! \brief Eat up the rest of the current line */
kusano 7d535a
static int cDumpLine(FILE *fp)
kusano 7d535a
{
kusano 7d535a
    register int c;
kusano 7d535a
    while ((c = fgetc(fp)) != '\n') ;
kusano 7d535a
    return 0;
kusano 7d535a
}
kusano 7d535a
kusano 7d535a
static int cParseIntFormat(char *buf, int *num, int *size)
kusano 7d535a
{
kusano 7d535a
    char *tmp;
kusano 7d535a
kusano 7d535a
    tmp = buf;
kusano 7d535a
    while (*tmp++ != '(') ;
kusano 7d535a
    sscanf(tmp, "%d", num);
kusano 7d535a
    while (*tmp != 'I' && *tmp != 'i') ++tmp;
kusano 7d535a
    ++tmp;
kusano 7d535a
    sscanf(tmp, "%d", size);
kusano 7d535a
    return 0;
kusano 7d535a
}
kusano 7d535a
kusano 7d535a
static int cParseFloatFormat(char *buf, int *num, int *size)
kusano 7d535a
{
kusano 7d535a
    char *tmp, *period;
kusano 7d535a
kusano 7d535a
    tmp = buf;
kusano 7d535a
    while (*tmp++ != '(') ;
kusano 7d535a
    *num = atoi(tmp); /*sscanf(tmp, "%d", num);*/
kusano 7d535a
    while (*tmp != 'E' && *tmp != 'e' && *tmp != 'D' && *tmp != 'd'
kusano 7d535a
           && *tmp != 'F' && *tmp != 'f') {
kusano 7d535a
        /* May find kP before nE/nD/nF, like (1P6F13.6). In this case the
kusano 7d535a
           num picked up refers to P, which should be skipped. */
kusano 7d535a
        if (*tmp=='p' || *tmp=='P') {
kusano 7d535a
           ++tmp;
kusano 7d535a
           *num = atoi(tmp); /*sscanf(tmp, "%d", num);*/
kusano 7d535a
        } else {
kusano 7d535a
           ++tmp;
kusano 7d535a
        }
kusano 7d535a
    }
kusano 7d535a
    ++tmp;
kusano 7d535a
    period = tmp;
kusano 7d535a
    while (*period != '.' && *period != ')') ++period ;
kusano 7d535a
    *period = '\0';
kusano 7d535a
    *size = atoi(tmp); /*sscanf(tmp, "%2d", size);*/
kusano 7d535a
kusano 7d535a
    return 0;
kusano 7d535a
}
kusano 7d535a
kusano 7d535a
static int ReadVector(FILE *fp, int n, int *where, int perline, int persize)
kusano 7d535a
{
kusano 7d535a
    register int i, j, item;
kusano 7d535a
    char tmp, buf[100];
kusano 7d535a
kusano 7d535a
    i = 0;
kusano 7d535a
    while (i < n) {
kusano 7d535a
        fgets(buf, 100, fp);    /* read a line at a time */
kusano 7d535a
        for (j=0; j
kusano 7d535a
            tmp = buf[(j+1)*persize];     /* save the char at that place */
kusano 7d535a
            buf[(j+1)*persize] = 0;       /* null terminate */
kusano 7d535a
            item = atoi(&buf[j*persize]); 
kusano 7d535a
            buf[(j+1)*persize] = tmp;     /* recover the char at that place */
kusano 7d535a
            where[i++] = item - 1;
kusano 7d535a
        }
kusano 7d535a
    }
kusano 7d535a
kusano 7d535a
    return 0;
kusano 7d535a
}
kusano 7d535a
kusano 7d535a
/*! \brief Read complex numbers as pairs of (real, imaginary) */
kusano 7d535a
static int cReadValues(FILE *fp, int n, complex *destination, int perline, int persize)
kusano 7d535a
{
kusano 7d535a
    register int i, j, k, s, pair;
kusano 7d535a
    register float realpart;
kusano 7d535a
    char tmp, buf[100];
kusano 7d535a
    
kusano 7d535a
    i = pair = 0;
kusano 7d535a
    while (i < n) {
kusano 7d535a
	fgets(buf, 100, fp);    /* read a line at a time */
kusano 7d535a
	for (j=0; j
kusano 7d535a
	    tmp = buf[(j+1)*persize];     /* save the char at that place */
kusano 7d535a
	    buf[(j+1)*persize] = 0;       /* null terminate */
kusano 7d535a
	    s = j*persize;
kusano 7d535a
	    for (k = 0; k < persize; ++k) /* No D_ format in C */
kusano 7d535a
		if ( buf[s+k] == 'D' || buf[s+k] == 'd' ) buf[s+k] = 'E';
kusano 7d535a
	    if ( pair == 0 ) {
kusano 7d535a
	  	/* The value is real part */
kusano 7d535a
		realpart = atof(&buf[s]);
kusano 7d535a
		pair = 1;
kusano 7d535a
	    } else {
kusano 7d535a
		/* The value is imaginary part */
kusano 7d535a
	        destination[i].r = realpart;
kusano 7d535a
		destination[i++].i = atof(&buf[s]);
kusano 7d535a
		pair = 0;
kusano 7d535a
	    }
kusano 7d535a
	    buf[(j+1)*persize] = tmp;     /* recover the char at that place */
kusano 7d535a
	}
kusano 7d535a
    }
kusano 7d535a
kusano 7d535a
    return 0;
kusano 7d535a
}
kusano 7d535a
kusano 7d535a
kusano 7d535a
void
kusano 7d535a
creadrb(int *nrow, int *ncol, int *nonz,
kusano 7d535a
        complex **nzval, int **rowind, int **colptr)
kusano 7d535a
{
kusano 7d535a
kusano 7d535a
    register int i, numer_lines = 0;
kusano 7d535a
    int tmp, colnum, colsize, rownum, rowsize, valnum, valsize;
kusano 7d535a
    char buf[100], type[4];
kusano 7d535a
    FILE *fp;
kusano 7d535a
kusano 7d535a
    fp = stdin;
kusano 7d535a
kusano 7d535a
    /* Line 1 */
kusano 7d535a
    fgets(buf, 100, fp);
kusano 7d535a
    fputs(buf, stdout);
kusano 7d535a
kusano 7d535a
    /* Line 2 */
kusano 7d535a
    for (i=0; i<4; i++) {
kusano 7d535a
        fscanf(fp, "%14c", buf); buf[14] = 0;
kusano 7d535a
        sscanf(buf, "%d", &tmp);
kusano 7d535a
        if (i == 3) numer_lines = tmp;
kusano 7d535a
    }
kusano 7d535a
    cDumpLine(fp);
kusano 7d535a
kusano 7d535a
    /* Line 3 */
kusano 7d535a
    fscanf(fp, "%3c", type);
kusano 7d535a
    fscanf(fp, "%11c", buf); /* pad */
kusano 7d535a
    type[3] = 0;
kusano 7d535a
#ifdef DEBUG
kusano 7d535a
    printf("Matrix type %s\n", type);
kusano 7d535a
#endif
kusano 7d535a
kusano 7d535a
    fscanf(fp, "%14c", buf); sscanf(buf, "%d", nrow);
kusano 7d535a
    fscanf(fp, "%14c", buf); sscanf(buf, "%d", ncol);
kusano 7d535a
    fscanf(fp, "%14c", buf); sscanf(buf, "%d", nonz);
kusano 7d535a
    fscanf(fp, "%14c", buf); sscanf(buf, "%d", &tmp);
kusano 7d535a
kusano 7d535a
    if (tmp != 0)
kusano 7d535a
        printf("This is not an assembled matrix!\n");
kusano 7d535a
    if (*nrow != *ncol)
kusano 7d535a
        printf("Matrix is not square.\n");
kusano 7d535a
    cDumpLine(fp);
kusano 7d535a
kusano 7d535a
    /* Allocate storage for the three arrays ( nzval, rowind, colptr ) */
kusano 7d535a
    callocateA(*ncol, *nonz, nzval, rowind, colptr);
kusano 7d535a
kusano 7d535a
    /* Line 4: format statement */
kusano 7d535a
    fscanf(fp, "%16c", buf);
kusano 7d535a
    cParseIntFormat(buf, &colnum, &colsize);
kusano 7d535a
    fscanf(fp, "%16c", buf);
kusano 7d535a
    cParseIntFormat(buf, &rownum, &rowsize);
kusano 7d535a
    fscanf(fp, "%20c", buf);
kusano 7d535a
    cParseFloatFormat(buf, &valnum, &valsize);
kusano 7d535a
    cDumpLine(fp);
kusano 7d535a
kusano 7d535a
#ifdef DEBUG
kusano 7d535a
    printf("%d rows, %d nonzeros\n", *nrow, *nonz);
kusano 7d535a
    printf("colnum %d, colsize %d\n", colnum, colsize);
kusano 7d535a
    printf("rownum %d, rowsize %d\n", rownum, rowsize);
kusano 7d535a
    printf("valnum %d, valsize %d\n", valnum, valsize);
kusano 7d535a
#endif
kusano 7d535a
kusano 7d535a
    ReadVector(fp, *ncol+1, *colptr, colnum, colsize);
kusano 7d535a
    ReadVector(fp, *nonz, *rowind, rownum, rowsize);
kusano 7d535a
    if ( numer_lines ) {
kusano 7d535a
        cReadValues(fp, *nonz, *nzval, valnum, valsize);
kusano 7d535a
    }
kusano 7d535a
kusano 7d535a
    fclose(fp);
kusano 7d535a
}