kusano 7d535a
kusano 7d535a
/*! @file dreadrb.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
kusano 7d535a
#include "slu_ddefs.h"
kusano 7d535a
kusano 7d535a
kusano 7d535a
void
kusano 7d535a
dreadtriple(int *m, int *n, int *nonz,
kusano 7d535a
	    double **nzval, int **rowind, int **colptr)
kusano 7d535a
{
kusano 7d535a
/*
kusano 7d535a
 * Output parameters
kusano 7d535a
 * =================
kusano 7d535a
 *   (a,asub,xa): asub[*] contains the row subscripts of nonzeros
kusano 7d535a
 *	in columns of matrix A; a[*] the numerical values;
kusano 7d535a
 *	row i of A is given by a[k],k=xa[i],...,xa[i+1]-1.
kusano 7d535a
 *
kusano 7d535a
 */
kusano 7d535a
    int    j, k, jsize, nnz, nz;
kusano 7d535a
    double *a, *val;
kusano 7d535a
    int    *asub, *xa, *row, *col;
kusano 7d535a
    int    zero_base = 0;
kusano 7d535a
kusano 7d535a
    /*  Matrix format:
kusano 7d535a
     *    First line:  #rows, #cols, #non-zero
kusano 7d535a
     *    Triplet in the rest of lines:
kusano 7d535a
     *                 row, col, value
kusano 7d535a
     */
kusano 7d535a
kusano 7d535a
    scanf("%d%d", n, nonz);
kusano 7d535a
    *m = *n;
kusano 7d535a
    printf("m %d, n %d, nonz %d\n", *m, *n, *nonz);
kusano 7d535a
    dallocateA(*n, *nonz, nzval, rowind, colptr); /* Allocate storage */
kusano 7d535a
    a    = *nzval;
kusano 7d535a
    asub = *rowind;
kusano 7d535a
    xa   = *colptr;
kusano 7d535a
kusano 7d535a
    val = (double *) SUPERLU_MALLOC(*nonz * sizeof(double));
kusano 7d535a
    row = (int *) SUPERLU_MALLOC(*nonz * sizeof(int));
kusano 7d535a
    col = (int *) SUPERLU_MALLOC(*nonz * sizeof(int));
kusano 7d535a
kusano 7d535a
    for (j = 0; j < *n; ++j) xa[j] = 0;
kusano 7d535a
kusano 7d535a
    /* Read into the triplet array from a file */
kusano 7d535a
    for (nnz = 0, nz = 0; nnz < *nonz; ++nnz) {
kusano 7d535a
	scanf("%d%d%lf\n", &row[nz], &col[nz], &val[nz]);
kusano 7d535a
kusano 7d535a
        if ( nnz == 0 ) { /* first nonzero */
kusano 7d535a
	    if ( row[0] == 0 || col[0] == 0 ) {
kusano 7d535a
		zero_base = 1;
kusano 7d535a
		printf("triplet file: row/col indices are zero-based.\n");
kusano 7d535a
	    } else
kusano 7d535a
		printf("triplet file: row/col indices are one-based.\n");
kusano 7d535a
        }
kusano 7d535a
kusano 7d535a
        if ( !zero_base ) { 
kusano 7d535a
 	  /* Change to 0-based indexing. */
kusano 7d535a
	  --row[nz];
kusano 7d535a
	  --col[nz];
kusano 7d535a
        }
kusano 7d535a
kusano 7d535a
	if (row[nz] < 0 || row[nz] >= *m || col[nz] < 0 || col[nz] >= *n
kusano 7d535a
	    /*|| val[nz] == 0.*/) {
kusano 7d535a
	    fprintf(stderr, "nz %d, (%d, %d) = %e out of bound, removed\n",
kusano 7d535a
		    nz, row[nz], col[nz], val[nz]);
kusano 7d535a
	    exit(-1);
kusano 7d535a
	} else {
kusano 7d535a
	    ++xa[col[nz]];
kusano 7d535a
	    ++nz;
kusano 7d535a
	}
kusano 7d535a
    }
kusano 7d535a
kusano 7d535a
    *nonz = nz;
kusano 7d535a
kusano 7d535a
    /* Initialize the array of column pointers */
kusano 7d535a
    k = 0;
kusano 7d535a
    jsize = xa[0];
kusano 7d535a
    xa[0] = 0;
kusano 7d535a
    for (j = 1; j < *n; ++j) {
kusano 7d535a
	k += jsize;
kusano 7d535a
	jsize = xa[j];
kusano 7d535a
	xa[j] = k;
kusano 7d535a
    }
kusano 7d535a
    
kusano 7d535a
    /* Copy the triplets into the column oriented storage */
kusano 7d535a
    for (nz = 0; nz < *nonz; ++nz) {
kusano 7d535a
	j = col[nz];
kusano 7d535a
	k = xa[j];
kusano 7d535a
	asub[k] = row[nz];
kusano 7d535a
	a[k] = val[nz];
kusano 7d535a
	++xa[j];
kusano 7d535a
    }
kusano 7d535a
kusano 7d535a
    /* Reset the column pointers to the beginning of each column */
kusano 7d535a
    for (j = *n; j > 0; --j)
kusano 7d535a
	xa[j] = xa[j-1];
kusano 7d535a
    xa[0] = 0;
kusano 7d535a
kusano 7d535a
    SUPERLU_FREE(val);
kusano 7d535a
    SUPERLU_FREE(row);
kusano 7d535a
    SUPERLU_FREE(col);
kusano 7d535a
kusano 7d535a
#ifdef CHK_INPUT
kusano 7d535a
    {
kusano 7d535a
	int i;
kusano 7d535a
	for (i = 0; i < *n; i++) {
kusano 7d535a
	    printf("Col %d, xa %d\n", i, xa[i]);
kusano 7d535a
	    for (k = xa[i]; k < xa[i+1]; k++)
kusano 7d535a
		printf("%d\t%16.10f\n", asub[k], a[k]);
kusano 7d535a
	}
kusano 7d535a
    }
kusano 7d535a
#endif
kusano 7d535a
kusano 7d535a
}
kusano 7d535a
kusano 7d535a
kusano 7d535a
void dreadrhs(int m, double *b)
kusano 7d535a
{
kusano 7d535a
    FILE *fp, *fopen();
kusano 7d535a
    int i;
kusano 7d535a
    /*int j;*/
kusano 7d535a
kusano 7d535a
    if ( !(fp = fopen("b.dat", "r")) ) {
kusano 7d535a
        fprintf(stderr, "dreadrhs: file does not exist\n");
kusano 7d535a
	exit(-1);
kusano 7d535a
    }
kusano 7d535a
    for (i = 0; i < m; ++i)
kusano 7d535a
      fscanf(fp, "%lf\n", &b[i]);
kusano 7d535a
kusano 7d535a
    /*        readpair_(j, &b[i]);*/
kusano 7d535a
    fclose(fp);
kusano 7d535a
}