dreadtriple.c.bak

SuperLU 2.2版本。对大型、稀疏、非对称的线性系统的直接求解

#include <stdio.h>
#include "superlu_ddefs.h"

#undef EXPAND_SYM

void
dreadtriple(FILE *fp, int_t *m, int_t *n, int_t *nonz,
	    double **nzval, int_t **rowind, int_t **colptr)
{
/*
 * Output parameters
 * =================
 *   (nzval, rowind, colptr): (*rowind)[*] contains the row subscripts of
 *      nonzeros in columns of matrix A; (*nzval)[*] the numerical values;
 *	column i of A is given by (*nzval)[k], k = (*rowind)[i],...,
 *      (*rowind)[i+1]-1.
 *
 */
    int_t    i, j, k, jsize, lasta, nnz, nz, new_nonz;
    double *a, *val;
    int_t    *asub, *xa, *row, *col;
    
    /* 	File format:
     *    First line:  #rows    #non-zero
     *    Triplet in the rest of lines:
     *                 row    col    value
     */

    /*fscanf(fp, "%d%d%d", m, n, nonz);*/
#ifdef _LONGINT
    fscanf(fp, "%ld%ld", n, nonz);
#else
    fscanf(fp, "%d%d", n, nonz);
#endif

#ifdef EXPAND_SYM
    new_nonz = 2 * *nonz - *n;
#else
    new_nonz = *nonz;
#endif
    *m = *n;
    printf("m %ld, n %ld, nonz %ld\n", *m, *n, *nonz);
    dallocateA_dist(*n, new_nonz, nzval, rowind, colptr); /* Allocate storage */
    a    = *nzval;
    asub = *rowind;
    xa   = *colptr;

    if ( !(val = (double *) SUPERLU_MALLOC(new_nonz * sizeof(double))) )
        ABORT("Malloc fails for val[]");
    if ( !(row = (int_t *) SUPERLU_MALLOC(new_nonz * sizeof(int_t))) )
        ABORT("Malloc fails for row[]");
    if ( !(col = (int_t *) SUPERLU_MALLOC(new_nonz * sizeof(int_t))) )
        ABORT("Malloc fails for col[]");

    for (j = 0; j < *n; ++j) xa[j] = 0;

    /* Read into the triplet array from a file */
    for (nnz = 0, nz = 0; nnz < *nonz; ++nnz) {
#ifdef _LONGINT
	fscanf(fp, "%ld%ld%lf\n", &row[nz], &col[nz], &val[nz]);
#else
	fscanf(fp, "%d%d%lf\n", &row[nz], &col[nz], &val[nz]);
#endif
#if 1
	/* Change to 0-based indexing. */
	--row[nz];
	--col[nz];
#endif
	if (row[nz] < 0 || row[nz] >= *m || col[nz] < 0 || col[nz] >= *n
	    /*|| val[nz] == 0.*/) {
	    fprintf(stderr, "nz %d, (%d, %d) = %e out of bound, removed\n", 
		    nz, row[nz], col[nz], val[nz]);
	    exit(-1);
	} else {
	    ++xa[col[nz]];
#ifdef EXPAND_SYM
	    if ( row[nz] != col[nz] ) { /* Excluding diagonal */
	      ++nz;
	      row[nz] = col[nz-1];
	      col[nz] = row[nz-1];
	      val[nz] = val[nz-1];
	      ++xa[col[nz]];
	    }
#endif	
	    ++nz;
	}
    }

    *nonz = nz;
#ifdef EXPAND_SYM
    printf("new_nonz after symmetric expansion:\t%d\n", *nonz);
#endif
    

    /* Initialize the array of column pointers */
    k = 0;
    jsize = xa[0];
    xa[0] = 0;
    for (j = 1; j < *n; ++j) {
	k += jsize;
	jsize = xa[j];
	xa[j] = k;
    }
    
    /* Copy the triplets into the column oriented storage */
    for (nz = 0; nz < *nonz; ++nz) {
	j = col[nz];
	k = xa[j];
	asub[k] = row[nz];
	a[k] = val[nz];
	++xa[j];
    }

    /* Reset the column pointers to the beginning of each column */
    for (j = *n; j > 0; --j)
	xa[j] = xa[j-1];
    xa[0] = 0;

    SUPERLU_FREE(val);
    SUPERLU_FREE(row);
    SUPERLU_FREE(col);

#ifdef CHK_INPUT
    for (i = 0; i < *n; i++) {
	printf("Col %d, xa %d\n", i, xa[i]);
	for (k = xa[i]; k < xa[i+1]; k++)
	    printf("%d\t%16.10f\n", asub[k], a[k]);
    }
#endif

}


void dreadrhs(int m, double *b)
{
    FILE *fp, *fopen();
    int i, j;

    if ( !(fp = fopen("b.dat", "r")) ) {
        fprintf(stderr, "dreadrhs: file does not exist\n");
	exit(-1);
    }
    for (i = 0; i < m; ++i)
      fscanf(fp, "%lf\n", &b[i]);
      /*fscanf(fp, "%d%lf\n", &j, &b[i]);*/
    /*        readpair_(j, &b[i]);*/
    fclose(fp);
}