Does anyone can help me on spooles package's iterative method ?

Does anyone can help me on spooles package's iterative method ?
Dear everyone :
I need your help !
I am using spooles , and I have successfully
use it in direct methods which is in directory
/misc/drivers , and I have run into problems when
I use the code in /Iter/drivers , for the demos in it
, I have run successfully , but I want to use the same
matrix of A and b as in direct method when I do in direct,
in fact , I am use a fortran code for main program , I
only to take use of spooles as a subroutine , to solve Ax=b for
me, the code iter.c in Iter/drivers , is a main c program , and I want to
change it as a subroutine for my fortran main program .
and I don't know how to change it to a suitable code . because in the
/misc/drivers ,for the code allInOne.c , the input of matrix A and
right-hand side matrix b , I have successuflly write the matrix for them
.
unfortunatelly , in /Iter/drivers . I don't know how to change them to
as a .inpmtxf and .densemtxf .
could you please see it for me ?
Do you know how to contact the author of spooles,
as for me , I have written email to him long ago, but I
can't successfully, the email is rejected by the internet.

the following is the iterative code changed from iter.c :

/*
-----------------------------------------------------
test the factor method for a grid matrix
(0) read in matrix from source file
(1) conver data matrix to InpMtx object if necessary
(2) create Graph and ETree object if necessary
(3) read in/create an ETree object
(4) create a solution matrix object
(5) multiply the solution with the matrix
to get a right hand side matrix object
(6) factor the matrix
(7) solve the system

created -- 98dec30, jwu
-----------------------------------------------------
*/

/* iter.c */
/* this program is from iter.c , I make it as a fortran's subroutine to iterative
solve Ax=b */
#include "../Iter.h"
#include "../../Graph.h"
#include
#define METHODS 10
#include "../../misc.h"
#include "../../FrontMtx.h"
#include "../../SymbFac.h"
#include "../../DenseMtx.h" //add by ljp
#include "../../timings.h" //add by ljp
#include "../../InpMtx.h"
// add by ljp

int
InpMtx_readFromAIJ2file(InpMtx *mtxA, char *srcFileName);
int
InpMtx_createGraph(InpMtx *mtxA, Graph *graph);

/*--------------------------------------------------------------------*/
void ljp_iter_9acous(char *parfilename,int iw ) ;
void ljp_iter_9acous_( char *PARFILENAME, int *IW )
{ char *parfilename ;
int iw;
iw = *IW ;
ljp_iter_9acous(parfilename ,iw) ;
}
// char parfilenamee[256]="in2";
// call ljp_iter_9acous_();
//void ljp_iter_9acous( char *parfilename )
void ljp_iter_9acous(char *parfilename ,int iw)
{
// add by fb
char parfilename[256]="in2";
char subname[256]="ljp_iter_9acous";
// add by fb

// ***** add by ljp ***
FILE *solFileName; // add by ljp
FILE *rhsFile_B, *matrixFile_A, *outfile3_X ;
char matrixFileName_A[80] , rhsFileName_origin[80] , solFile_X[80] ;
char matrixFileName_original[80] , solFile[80] ;

char msgFileName_A[80];
char msgFileName_B[80];
FILE *matrixFile;
FILE *rhsFile;
FILE *msgFile_A , *msgFile_B ; // for A and B matrix
FILE *inputFile ;
int nrow,ncol,nent,ient,irow,jcol,jrow,jrhs;
// ***** add by ljp ***

char etreeFileName[80], mtxFileName[80], *cpt, rhsFileName[80],
srcFileName[80], ctemp[81], msgFileName[80], slnFileName[80] ;
Chv *chv, *rootchv ;
ChvManager *chvmanager ;
DenseMtx *mtxB, *mtxQ, *mtxX, *mtxZ ;
double one[2] = { 1.0, 0.0 } ;
FrontMtx *frontmtx ;
InpMtx *mtxA ;
SubMtxManager *mtxmanager ;
double cputotal, droptol, conv_tol, factorops ;
double cpus[9] ;
Drand drand ;
double nops, tau, t1, t2 ;
ETree *frontETree ;
Graph *graph ;
FILE *msgFile, *inFile ;
int error, loc, msglvl, neqns, nzf, iformat,
pivotingflag, rc, seed, sparsityflag, symmetryflag,
method[METHODS], type, nrhs, etreeflag ;
int stats[6] ;
int nnzA, Ik, itermax, zversion, iterout ;
IV *newToOldIV, *oldToNewIV ;
IVL *symbfacIVL ;
int i, j, k, m, n, imethod, maxdomainsize, maxzeros, maxsize;
int nouter,ninner ;

// **************** add by ljp ************
FILE *outfile, *outfile1, *outfile2, *outfile3;
// for entry in B matrix ,same to direct LU method
char filename_prefix[80]="B_rhs_9_acoust_iter.input.iw";
char filename_suffix[80];
sprintf(filename_suffix,"%d",iw); // iw is 1,2,3...
strcpy(rhsFileName_origin,filename_prefix);
strcat(rhsFileName_origin,filename_suffix);
// 4-4-night ljp rhsFile =fopen(rhsFileName_origin,"r");
// for entry in B matrix ,same to direct LU method

// for entry in A matrix ,same to direct LU method
char filename_prefix1[80]="matrix_A_9_acoust_iter.input.iw";
char filename_suffix1[80];
sprintf(filename_suffix1,"%d",iw);
strcpy(matrixFileName_original,filename_prefix1);
strcat(matrixFileName_original,filename_suffix1);
// 4-4-night ljp matrixFile =fopen(matrixFileName_original,"r");
// for entry in A matrix ,same to direct LU method

/ /进入解决x矩阵,直接LU相同method
char filename_prefix3[80]="ljp_solve_x_9_acoust_iter.output.iw";
char filename_suffix3[80];
sprintf(filename_suffix3,"%d",iw);
strcpy(solFile,filename_prefix3);
strcat(solFile,filename_suffix3);
outfile3 =fopen(solFile,"w");
/ /进入解决x矩阵,直接LU相同method

char filename_prefix2[80]="ljp_message_9_acoust_iter.output.iw";
char filename_suffix2[80];
sprintf(filename_suffix2,"%d",iw);
strcpy(msgFileName,filename_prefix2);
strcat(msgFileName,filename_suffix2);
// 4-4-night ljp outfile2 =fopen(msgFileName,"w");

/ /通过ljpfor B.densemtxf
char rhsFileName2[80];
char filename_suffix_B[80]="iw.B_rhs_9_acoust_iter.densemtxf";
sprintf(rhsFileName,"%d",iw);
strcpy(rhsFileName2,filename_suffix_B);
strcat(rhsFileName, rhsFileName2);
rhsFile_B =fopen(rhsFileName, "w"); //4-3-2009 ljp
// rhsFile_B is : 1iw.B_rhs_9_acoust_iter.densemtxf , etc
/ /通过ljpfor B.densemtxf

/*
/ /通过ljpfor B.densemtxf
char rhsFileName1[80];
char filename_suffix_B[80]="iw.B_rhs_9_acoust_iter.densemtxf";
sprintf(rhsFileName1,"%d",iw);
strcpy(rhsFileName,filename_suffix_B);
strcat(rhsFileName1, rhsFileName);
rhsFile_B =fopen(rhsFileName1,"w"); //4-3-2009 ljp
// rhsFile_B is : 1iw.B_rhs_9_acoust_iter.densemtxf , etc
/ /通过ljpfor B.densemtxf
*/

/ /通过ljpfor A.inputmtxf
char srcFileName2[80];
char filename_suffix1A[80]="iw.matrix_A_9_acoust_iter.inpmtxf";
sprintf(srcFileName,"%d",iw);
strcpy(srcFileName2,filename_suffix1A);
strcat(srcFileName, srcFileName2);
matrixFile_A =fopen(srcFileName ,"w"); //4-3-2009 ljp
// matrixFile_A is : 1iw.matrix_A_9_acoust_iter.inpmtxf , etc
/ /通过ljpfor A.inputmtxf

/ /通过ljpfor X.densemtxf
char filename_prefix3_X[80];
char filename_suffix3_X[80]="iw.ljp_solve_x_9_acoust_iter.densemtxf";
sprintf(filename_prefix3_X,"%d",iw);
strcpy(solFile_X,filename_suffix3_X);
strcat(filename_prefix3_X, solFile_X);
//4-4-2009 night slnFileName solFileName = fopen(filename_prefix3_X, "w") ;
// slnFileName = fopen(filename_prefix3_X, "w") ; //4-3-2009 ljp
solFileName = fopen(filename_prefix3_X, "w") ; //4-3-2009 ljp
// solFileName is : 1iw.ljp_solve_x_9_acoust_iter.densemtxf
/ /通过ljpfor X.densemtxf

// printf("rhsFile_B =%s\n",rhsFile_B);
// printf("matrixFile_A =%s\n",matrixFile_A);
// printf("solFileName =%s\n",solFileName);

// **************** add by ljp ************
// **************** add by ljp ************
/*
--------------------------------------------
STEP 1: read the entries from the input file
and create the InpMtx object
--------------------------------------------
*/

inputFile = fopen(matrixFileName_original, "r") ; // the original entry matrix A
msgFile = fopen(msgFileName, "w") ;

fscanf(inputFile, "%d %d %d", &nrow, &ncol, &nent) ;
neqns = nrow ;

printf("pass 6"); // by ljp
printf("nrow=%d\n", nrow); // by ljp
printf("ncol=%d\n", ncol); // by ljp
printf("nent=%d\n", nent); // by ljp
printf("neqns=%d\n",neqns); // by ljp
mtxA = InpMtx_new() ;
printf("pass 6-1"); // by ljp
type=2; // add by ljp
InpMtx_init(mtxA, INPMTX_BY_ROWS, type, nent, neqns) ;
if ( type == SPOOLES_REAL )
{
double value ;
for ( ient = 0 ; ient < nent ; ient++ )
{
fscanf(inputFile, "%d %d %le", &irow, &jcol, &value) ;
InpMtx_inputRealEntry(mtxA, irow, jcol, value) ;
}
}
else
{
printf("pass 6-2"); // by ljp
double imag, real ;
for ( ient = 0 ; ient < nent ; ient++ )
{
fscanf(inputFile, "%d %d %le %le", &irow, &jcol, &real, &imag) ;
InpMtx_inputComplexEntry(mtxA, irow, jcol, real, imag) ;
}
}

fclose(inputFile) ;
InpMtx_changeStorageMode(mtxA, INPMTX_BY_VECTORS) ; // 4-5-2009
printf("pass 6-2-1"); // by ljp
if ( msglvl > 0 ) {
fprintf(msgFile, "\n\n input matrix") ;
InpMtx_writeForHumanEye(mtxA, msgFile) ;
fflush(msgFile) ;
}
printf("pass 6-3"); // by ljp
InpMtx_writeForMatlab(mtxA, "A", msgFile) ;
InpMtx_writeToFormattedFile(mtxA, matrixFile_A) ; // the A.inpmtxf

/*
-----------------------------------------
STEP 2: read the right hand side matrix Y
-----------------------------------------
*/

inputFile = fopen(rhsFileName_origin, "r") ; // the original entry matrix B
printf("pass 6-3-1"); // by ljp
fscanf(inputFile, "%d %d", &nrow, &nrhs) ;
mtxB = DenseMtx_new() ;
// DenseMtx_init(mtxB, type, 0, 0, neqns, 1, 1, neqns) ;
printf("pass 6-3-2"); // by ljp
DenseMtx_init(mtxB, type, 0, 0, neqns, nrhs, 1, neqns) ;
DenseMtx_zero(mtxB) ;
if ( type == SPOOLES_REAL )
{
printf("pass 6-3-3-1 is type == SPOOLES_REAL"); // by ljp
double value ;
for ( irow = 0 ; irow < nrow ; irow++ )
{
fscanf(inputFile, "%d", &jrow) ;
for ( jrhs = 0 ; jrhs < nrhs ; jrhs++ )
{
fscanf(inputFile, "%le", &value) ;
DenseMtx_setRealEntry(mtxB, jrow, jrhs, value) ;
}
}
}
else
{
printf("pass 6-3-3 is type ==SPOOLES_Complex"); // by ljp
double imag, real ;
for ( irow = 0 ; irow < nrow ; irow++ )
{
fscanf(inputFile, "%d", &jrow) ;
for ( jrhs = 0 ; jrhs < nrhs ; jrhs++ )
{
fscanf(inputFile, "%le %le", &real, &imag) ;
DenseMtx_setComplexEntry(mtxB, jrow, jrhs, real, imag) ;
}
}
}
printf("pass 6-3-4"); // by ljp
fclose(inputFile) ;
printf("pass 6-8"); // by ljp
// if ( msglvl > 0 ) {
fprintf(msgFile, "\n\n rhs matrix in original ordering") ;
DenseMtx_writeForHumanEye(mtxB, msgFile) ;
fflush(msgFile) ;
// }
DenseMtx_writeForMatlab(mtxB, "b", msgFile) ;

DenseMtx_writeToFormattedFile(mtxB, rhsFile_B) ; // the B.densemtxf matrix
// DenseMtx_writeToFormattedFile(mtxB, msgFile_B) ;

// ***** add by ljp ***
// ***** add by ljp ***

// now begin read the inpmtxf and densemtxf matrix
// now begin read the inpmtxf and densemtxf matrix
/* read input file */
inFile = fopen(parfilename, "r");
if (inFile == (FILE *)NULL) {
fprintf(stderr, "\n fatal error in %s: unable to open file %s\n",
subname, parfilename) ; // ????ljp
//by ljp return(-1) ;
}
for (i=0; i 2) {
fprintf(stderr, "\n fatal error in %s: "
"invalid source matrix format\n",subname) ;
//by ljp return(-1) ;
}
}
else if (k==1)
sscanf(ctemp, "%s", srcFileName);
else if (k==2)
sscanf(ctemp, "%s", mtxFileName);
else if (k==3) {
sscanf(ctemp, "%d", &etreeflag);
if (etreeflag < 0 || etreeflag > 4) {
fprintf(stderr, "\n fatal error in %s: "
"invalid etree file status\n",subname) ;
//by ljp return(-1) ;
}
}
else if (k==4)
sscanf(ctemp, "%s", etreeFileName);
else if (k==5)
sscanf(ctemp, "%s", rhsFileName);
else if (k==6)
sscanf(ctemp, "%s", slnFileName);
else if (k==7){
sscanf(ctemp, "%s", msgFileName);
if ( strcmp(msgFileName, "stdout") == 0 ) {
msgFile = stdout ;
}
如果(msgFile = fopen (msgFileName, "a")) == NULL ) {
fprintf(stderr, "\n fatal error in %s"
"\n unable to open file %s\n", subname, ctemp) ;
//by ljp return(-1) ;
}
}
else if (k==8)
sscanf(ctemp, "%d %d %d %d %d %d",
&msglvl, &seed, &nrhs, &Ik, &itermax, &iterout);
else if (k==9)
sscanf(ctemp, "%d %d %d", &symmetryflag, &sparsityflag, &pivotingflag);
else if (k==10)
sscanf(ctemp, "%lf %lf %lf", &tau, &droptol, &conv_tol);
else if (k==11) {
/*
for (j=0; j 1) {
fprintf(msgFile,"*** Multiple right-hand-side vectors is not allowed yet.\n");
fprintf(msgFile,"*** nrhs is reset to 1.\n");
nrhs =1;
}

// the parfilename file's parameters name is not correct now ????
printf("pass 7"); // by ljp
printf("type=%d\n",type); // by ljp
printf (" srcFileName = % s \ n ", srcFileName);/ /通过ljp
printf("mtxFileName=%s\n", mtxFileName); // by ljp
printf("etreeFileName=%s\n",etreeFileName); // by ljp
printf("rhsFileName=%s\n", rhsFileName); // by ljp
printf("msglvl=%d\n",msglvl); // by ljp
printf("seed=%d\n", seed); // by ljp
printf("symmetryflag=%d\n",symmetryflag); // by ljp
printf("sparsityflag=%d\n", sparsityflag); // by ljp
printf("pivotingflag=%d\n",pivotingflag); // by ljp
printf("tau=%e\n",tau ); // by ljp
printf("droptol=%e\n",droptol ); // by ljp
printf("conv_tol=%e\n",conv_tol); // by ljp
printf("pass 7-2"); // by ljp

// now begin read the inpmtxf and densemtxf matrix
// now begin read the inpmtxf and densemtxf matrix
fprintf(msgFile,
"\n %s "
"\n srcFileName -- %s"
"\n mtxFileName -- %s"
"\n etreeFileName -- %s"
"\n rhsFileName -- %s"
"\n msglvl -- %d"
"\n seed -- %d"
"\n symmetryflag -- %d"
"\n sparsityflag -- %d"
"\n pivotingflag -- %d"
"\n tau -- %e"
"\n droptol -- %e"
"\n conv_tol -- %e"
"\n method -- ",
subname, srcFileName, mtxFileName, etreeFileName, rhsFileName,
msglvl, seed, symmetryflag, sparsityflag, pivotingflag,
tau, droptol, conv_tol) ;

for (k=0; k0 && strcmp(mtxFileName, "none") != 0 ) {
rc = InpMtx_writeToFile(mtxA, mtxFileName) ;
if ( rc != 1 )
fprintf(msgFile, "\n return value %d from InpMtx_writeToFile(%p,%s)",
rc, mtxA, mtxFileName) ;
}

printf("pass 7-8"); // by ljp
fprintf(msgFile, "\n CPU %8.3f : read in (+ convert to) mtxA from file %s",
t2 - t1, mtxFileName) ;
if (rc != 1) {
goto end_read;
}
type = mtxA->inputMode ;
neqns = 1 + IVmax(mtxA->nent, InpMtx_ivec1(mtxA), &loc) ;
printf("pass 7-9"); // by ljp
if ( INPMTX_IS_BY_ROWS(mtxA) ) {
fprintf(msgFile, "\n matrix coordinate type is rows") ;
} else if ( INPMTX_IS_BY_COLUMNS(mtxA) ) {
fprintf(msgFile, "\n matrix coordinate type is columns") ;
} else if ( INPMTX_IS_BY_CHEVRONS(mtxA) ) {
fprintf(msgFile, "\n matrix coordinate type is chevrons") ;
} else {
fprintf(msgFile, "\n\n, error, bad coordinate type") ;
rc=-1;
goto end_read;
}
if ( INPMTX_IS_RAW_DATA(mtxA) ) {
fprintf(msgFile, "\n matrix storage mode is raw data\n") ;
} else if ( INPMTX_IS_SORTED(mtxA) ) {
fprintf(msgFile, "\n matrix storage mode is sorted\n") ;
} else if ( INPMTX_IS_BY_VECTORS(mtxA) ) {
fprintf(msgFile, "\n matrix storage mode is by vectors\n") ;
} else {
fprintf(msgFile, "\n\n, error, bad storage mode") ;
rc=-1;
goto end_read;
}

printf("pass 8-1"); // by ljp
if ( msglvl > 1 ) {
fprintf(msgFile, "\n\n after reading InpMtx object from file %s",
mtxFileName) ;
if ( msglvl == 2 ) {
InpMtx_writeStats(mtxA, msgFile) ;
} else {
InpMtx_writeForHumanEye(mtxA, msgFile) ;
}
fflush(msgFile) ;
}
printf("pass 8-2"); // by ljp
/*
Get the nonzeros in matrix A and print it
*/
nnzA = InpMtx_nent( mtxA );
fprintf(msgFile, "\n\n Input matrix size %d NNZ %d",
neqns, nnzA) ;

/*
--------------------------------------------------------
generate the linear system
1.生成解决方案矩阵和充满随机numbers
2. generate rhs matrix and fill with zeros
3. compute matrix-matrix multiply
--------------------------------------------------------
*/
MARKTIME(t1) ;
mtxX = DenseMtx_new() ;
DenseMtx_init(mtxX, type, 0, -1, neqns, nrhs, 1, neqns) ;
mtxB = DenseMtx_new() ;

printf("pass 8-2"); // by ljp
if (strcmp(rhsFileName, "none")) {
// printf("pass 8-2-1"); // by ljp
rc = DenseMtx_readFromFile(mtxB, rhsFileName) ;
if ( rc != 1 )
fprintf(msgFile, "\n return value %d from DenseMtx_readFromFile(%p,%s)",
rc, mtxB, rhsFileName) ;
DenseMtx_zero(mtxX) ;
}
// printf("pass 8-2-2"); // by ljp
else {
DenseMtx_init(mtxB, type, 1, -1, neqns, nrhs, 1, neqns) ;
DenseMtx_fillRandomEntries(mtxX, &drand) ;
DenseMtx_zero(mtxB) ;
switch ( symmetryflag ) {
//printf("pass 8-2-3"); // by ljp
case SPOOLES_SYMMETRIC :
InpMtx_sym_mmm(mtxA, mtxB, one, mtxX) ;
break ;
case SPOOLES_HERMITIAN :
// printf("pass 8-2-4"); // by ljp
InpMtx_herm_mmm(mtxA, mtxB, one, mtxX) ;
break ;
case SPOOLES_NONSYMMETRIC :
// printf("pass 8-2-5"); // by ljp
InpMtx_nonsym_mmm(mtxA, mtxB, one, mtxX) ;
break ;
default :
break ;
}
}

printf("pass 8-3"); // by ljp
MARKTIME(t2) ;
流(msgFile,“\ n CPU % 8.3 f:设置solution and rhs ",
t2 - t1) ;
if ( msglvl > 2 ) {
fprintf(msgFile, "\n\n original mtxX") ;
DenseMtx_writeForHumanEye(mtxX, msgFile) ;
fprintf(msgFile, "\n\n original mtxB") ;
DenseMtx_writeForHumanEye(mtxB, msgFile) ;
fflush(msgFile) ;
}
if (rc != 1) {
InpMtx_free(mtxA);
DenseMtx_free(mtxX);
DenseMtx_free(mtxB);
goto end_init;
}

/*
------------------------
read in/create the ETree object
------------------------
*/

printf("pass 8-4"); // by ljp
MARKTIME(t1) ;
if (etreeflag == 0) { /* read in ETree from file */
if ( strcmp(etreeFileName, "none") == 0 )
fprintf(msgFile, "\n no file to read from") ;
frontETree = ETree_new() ;
rc = ETree_readFromFile(frontETree, etreeFileName) ;
if (rc!=1)
fprintf(msgFile, "\n return value %d from ETree_readFromFile(%p,%s)",
rc, frontETree, etreeFileName) ;
}
else {
graph = Graph_new() ;
rc = InpMtx_createGraph(mtxA, graph);
if (rc!=1) {
fprintf(msgFile, "\n return value %d from InpMtx_createGraph(%p,%p)",
rc, mtxA, graph) ;
Graph_free(graph);
goto end_tree;
}
if (etreeflag == 1) { /* Via BestOfNDandMS */
maxdomainsize = 500; maxzeros = 1000; maxsize = 64 ;
frontETree = orderViaBestOfNDandMS(graph, maxdomainsize, maxzeros,
maxsize, seed, msglvl, msgFile) ;
}
else if (etreeflag == 2) { /* Via MMD */
frontETree = orderViaMMD(graph, seed, msglvl, msgFile) ;
}
else if (etreeflag == 3) { /* Via MS */
maxdomainsize = 500;
frontETree = orderViaMS(graph, maxdomainsize, seed, msglvl, msgFile) ;
}
else if (etreeflag == 4) { /* Via ND */
maxdomainsize = 500;
frontETree = orderViaND(graph, maxdomainsize, seed, msglvl, msgFile) ;
}
Graph_free(graph);

/* optionally write out the ETree object */
if ( strcmp(etreeFileName, "none") != 0 ) {
fprintf(msgFile, "\n\n writing out ETree to file %s",
etreeFileName) ;
ETree_writeToFile(frontETree, etreeFileName) ;
}
}
MARKTIME(t2) ;
fprintf(msgFile, "\n CPU %8.3f : read in/create frontETree from file %s",
t2 - t1, etreeFileName) ;
if ( rc != 1 ) {
ETree_free(frontETree);
goto end_tree;
}

ETree_leftJustify(frontETree) ;
if ( msglvl > 1 ) {
fprintf(msgFile, "\n\n after reading ETree object from file %s",
etreeFileName) ;
if ( msglvl == 2 ) {
ETree_writeStats(frontETree, msgFile) ;
} else {
ETree_writeForHumanEye(frontETree, msgFile) ;
}
}
fflush(msgFile) ;
/*
--------------------------------------------------
get the permutations, permute the matrix and the
front tree, and compute the symbolic factorization
--------------------------------------------------
*/
MARKTIME(t1) ;
oldToNewIV = ETree_oldToNewVtxPerm(frontETree) ;
newToOldIV = ETree_newToOldVtxPerm(frontETree) ;
MARKTIME(t2) ;
fprintf(msgFile, "\n CPU %8.3f : get permutations", t2 - t1) ;
MARKTIME(t1) ;
ETree_permuteVertices(frontETree, oldToNewIV) ;
MARKTIME(t2) ;
fprintf(msgFile, "\n CPU %8.3f : permute front tree", t2 - t1) ;
MARKTIME(t1) ;
InpMtx_permute(mtxA, IV_entries(oldToNewIV), IV_entries(oldToNewIV)) ;
MARKTIME(t2) ;
fprintf(msgFile, "\n CPU %8.3f : permute mtxA", t2 - t1) ;
if ( symmetryflag == SPOOLES_SYMMETRIC
|| symmetryflag == SPOOLES_HERMITIAN ) {
MARKTIME(t1) ;
InpMtx_mapToUpperTriangle(mtxA) ;
MARKTIME(t2) ;
fprintf(msgFile, "\n CPU %8.3f : map to upper triangle", t2 - t1) ;
}
if ( ! INPMTX_IS_BY_CHEVRONS(mtxA) ) {
MARKTIME(t1) ;
InpMtx_changeCoordType(mtxA, INPMTX_BY_CHEVRONS) ;
MARKTIME(t2) ;
fprintf(msgFile, "\n CPU %8.3f : change coordinate type", t2 - t1) ;
}
if ( INPMTX_IS_RAW_DATA(mtxA) ) {
MARKTIME(t1) ;
InpMtx_changeStorageMode(mtxA, INPMTX_SORTED) ;
MARKTIME(t2) ;
fprintf(msgFile, "\n CPU %8.3f : sort entries ", t2 - t1) ;
}
if ( INPMTX_IS_SORTED(mtxA) ) {
MARKTIME(t1) ;
InpMtx_changeStorageMode(mtxA, INPMTX_BY_VECTORS) ;
MARKTIME(t2) ;
fprintf(msgFile, "\n CPU %8.3f : convert to vectors ", t2 - t1) ;
}
MARKTIME(t1) ;
symbfacIVL = SymbFac_initFromInpMtx(frontETree, mtxA) ;
MARKTIME(t2) ;
fprintf(msgFile, "\n CPU %8.3f : symbolic factorization", t2 - t1) ;
MARKTIME(t1) ;
DenseMtx_permuteRows(mtxB, oldToNewIV) ;
MARKTIME(t2) ;
fprintf(msgFile, "\n CPU %8.3f : permute rhs", t2 - t1) ;

/*
------------------------------
initialize the FrontMtx object
------------------------------
*/
MARKTIME(t1) ;
frontmtx = FrontMtx_new() ;
mtxmanager = SubMtxManager_new() ;
SubMtxManager_init(mtxmanager, NO_LOCK, 0) ;
FrontMtx_init(frontmtx, frontETree, symbfacIVL,
type, symmetryflag, sparsityflag, pivotingflag,
NO_LOCK, 0, NULL, mtxmanager, msglvl, msgFile) ;
MARKTIME(t2) ;
fprintf(msgFile, "\n\n CPU %8.3f : initialize the front matrix",
t2 - t1) ;
if ( msglvl > 1 ) {
fprintf(msgFile,
"\n nendD = %d, nentL = %d, nentU = %d",
frontmtx->nentD, frontmtx->nentL, frontmtx->nentU) ;
SubMtxManager_writeForHumanEye(mtxmanager, msgFile) ;
}
if ( msglvl > 2 ) {
fprintf(msgFile, "\n front matrix initialized") ;
FrontMtx_writeForHumanEye(frontmtx, msgFile) ;
fflush(msgFile) ;
}
/*
-----------------
factor the matrix
-----------------
*/
nzf = ETree_nFactorEntries(frontETree, symmetryflag) ;
factorops = ETree_nFactorOps (frontETree、类型、信谊metryflag) ;
fprintf(msgFile,
"\n %d factor entries, %.0f factor ops, %8.3f ratio",
nzf, factorops, factorops/nzf) ;
IVzero(6, stats) ;
DVzero(9, cpus) ;
chvmanager = ChvManager_new() ;
ChvManager_init(chvmanager, NO_LOCK, 1) ;
MARKTIME(t1) ;
rootchv = FrontMtx_factorInpMtx(frontmtx, mtxA, tau, droptol,
chvmanager, &error, cpus,
stats, msglvl, msgFile) ;
MARKTIME(t2) ;
fprintf(msgFile, "\n\n CPU %8.3f : factor matrix, %8.3f mflops",
t2 - t1, 1.e-6*factorops/(t2-t1)) ;
if ( rootchv != NULL ) {
fprintf(msgFile, "\n\n factorization did not complete") ;
for ( chv = rootchv ; chv != NULL ; chv = chv->next ) {
fprintf(stdout, "\n chv %d, nD = %d, nL = %d, nU = %d",
chv->id, chv->nD, chv->nL, chv->nU) ;
}
}
if ( error >= 0 ) {
fprintf(msgFile, "\n\n error encountered at front %d\n", error) ;
rc=error ;
goto end_front;
}
fprintf(msgFile,
"\n %8d pivots, %8d pivot tests, %8d delayed rows and columns",
stats[0], stats[1], stats[2]) ;
if ( frontmtx->rowadjIVL != NULL ) {
fprintf(msgFile,
"\n %d entries in rowadjIVL", frontmtx->rowadjIVL->tsize) ;
}
if ( frontmtx->coladjIVL != NULL ) {
fprintf(msgFile,
", %d entries in coladjIVL", frontmtx->coladjIVL->tsize) ;
}
if ( frontmtx->upperblockIVL != NULL ) {
fprintf(msgFile,
"\n %d fronts, %d entries in upperblockIVL",
frontmtx->nfront, frontmtx->upperblockIVL->tsize) ;
}
if ( frontmtx->lowerblockIVL != NULL ) {
fprintf(msgFile,
", %d entries in lowerblockIVL",
frontmtx->lowerblockIVL->tsize) ;
}
fprintf(msgFile,
"\n %d entries in D, %d entries in L, %d entries in U",
stats[3], stats[4], stats[5]) ;
fprintf(msgFile, "\n %d locks", frontmtx->nlocks) ;
if ( FRONTMTX_IS_SYMMETRIC(frontmtx)
|| FRONTMTX_IS_HERMITIAN(frontmtx) ) {
int nneg, npos, nzero ;

FrontMtx_inertia(frontmtx, &nneg, &nzero, &npos) ;
fprintf(msgFile,
"\n %d negative, %d zero and %d positive eigenvalues",
nneg, nzero, npos) ;
fflush(msgFile) ;
}
cputotal = cpus[8] ;
if ( cputotal > 0.0 ) {
fprintf(msgFile,
"\n initialize fronts %8.3f %6.2f"
"\n load original entries %8.3f %6.2f"
"\n update fronts %8.3f %6.2f"
"\n assemble postponed data %8.3f %6.2f"
"\n factor fronts %8.3f %6.2f"
"\n extract postponed data %8.3f %6.2f"
"\n store factor entries %8.3f %6.2f"
"\n miscellaneous %8.3f %6.2f"
"\n total time %8.3f",
cpus[0], 100.*cpus[0]/cputotal,
cpus[1], 100.*cpus[1]/cputotal,
cpus[2], 100.*cpus[2]/cputotal,
cpus[3], 100.*cpus[3]/cputotal,
cpus[4], 100.*cpus[4]/cputotal,
cpus[5], 100.*cpus[5]/cputotal,
cpu [6], 100.*cpus[6]/cputotal,
cpus[7], 100.*cpus[7]/cputotal, cputotal) ;
}
if ( msglvl > 1 ) {
SubMtxManager_writeForHumanEye(mtxmanager, msgFile) ;
ChvManager_writeForHumanEye(chvmanager, msgFile) ;
}
if ( msglvl > 2 ) {
fprintf(msgFile, "\n\n front factor matrix") ;
FrontMtx_writeForHumanEye(frontmtx, msgFile) ;
}

/*
------------------------------
post-process the factor matrix
------------------------------
*/
MARKTIME(t1) ;
FrontMtx_postProcess(frontmtx, msglvl, msgFile) ;
MARKTIME(t2) ;
fprintf(msgFile, "\n\n CPU %8.3f : post-process the matrix", t2 - t1) ;
if ( msglvl > 2 ) {
fprintf(msgFile, "\n\n front factor matrix after post-processing") ;
FrontMtx_writeForHumanEye(frontmtx, msgFile) ;
}
fprintf(msgFile, "\n\n after post-processing") ;
if ( msglvl > 1 ) SubMtxManager_writeForHumanEye(frontmtx->manager, msgFile) ;
/*
----------------
solve the system
----------------
*/
neqns = mtxB->nrow ;
mtxZ = DenseMtx_new() ;
DenseMtx_init(mtxZ, type, 0, 0, neqns, nrhs, 1, neqns) ;
zversion=INPMTX_IS_COMPLEX_ENTRIES(mtxA);

for (k=0; k 2 ) {
fprintf(msgFile, "\n\n rhs") ;
DenseMtx_writeForHumanEye(mtxB, msgFile) ;
fflush(stdout) ;
}
fprintf(msgFile, "\n\n itemax %d", itermax) ;
DVzero(6, cpus) ;
MARKTIME(t1) ;
switch ( method[k] ) {
case BiCGStabR :
if (zversion)
rc=zbicgstabr(neqns, type, symmetryflag, mtxA, frontmtx, mtxZ, mtxB,
itermax, conv_tol, msglvl, msgFile);
else
rc=bicgstabr(neqns, type, symmetryflag, mtxA, frontmtx, mtxZ, mtxB,
itermax, conv_tol, msglvl, msgFile);

break;
case BiCGStabL :
if (zversion)
rc=zbicgstabl(neqns, type, symmetryflag, mtxA, frontmtx, mtxZ, mtxB,
itermax, conv_tol, msglvl, msgFile);
else
rc=bicgstabl(neqns, type, symmetryflag, mtxA, frontmtx, mtxZ, mtxB,
itermax, conv_tol, msglvl, msgFile);
break;
case TFQMRR :
if (zversion)
rc=ztfqmrr(neqns, type, symmetryflag, mtxA, frontmtx, mtxZ, mtxB,
itermax, conv_tol, msglvl, msgFile);
else
rc=tfqmrr(neqns, type, symmetryflag, mtxA, frontmtx, mtxZ, mtxB,
itermax, conv_tol, msglvl, msgFile);
break;
case TFQMRL :
if (zversion)
rc=ztfqmrl(neqns, type, symmetryflag, mtxA, frontmtx, mtxZ, mtxB,
itermax, conv_tol, msglvl, msgFile);
else
rc=tfqmrl(neqns, type, symmetryflag, mtxA, frontmtx, mtxZ, mtxB,
itermax, conv_tol, msglvl, msgFile);
break;
case PCGR :
if (zversion)
rc=zpcgr(neqns, type, symmetryflag, mtxA, frontmtx, mtxZ, mtxB,
itermax, conv_tol, msglvl, msgFile);
else
rc=pcgr(neqns, type, symmetryflag, mtxA, frontmtx, mtxZ, mtxB,
itermax, conv_tol, msglvl, msgFile);
break;
case PCGL :
if (zversion)
rc=zpcgl(neqns, type, symmetryflag, mtxA, frontmtx, mtxZ, mtxB,
itermax, conv_tol, msglvl, msgFile);
else
rc=pcgl(neqns, type, symmetryflag, mtxA, frontmtx, mtxZ, mtxB,
itermax, conv_tol, msglvl, msgFile);
break;
case MLBiCGStabR :
mtxQ = DenseMtx_new() ;
DenseMtx_init(mtxQ, type, 0, -1, neqns, Ik, 1, neqns) ;
Drand_setUniform(&drand, 0.0, 1.0) ;
DenseMtx_fillRandomEntries(mtxQ, &drand) ;
if (zversion)
rc=zmlbicgstabr(neqns, type, symmetryflag, mtxA, frontmtx, mtxQ, mtxZ,
mtxB, itermax, conv_tol, msglvl, msgFile);
else
rc=mlbicgstabr(neqns, type, symmetryflag, mtxA, frontmtx, mtxQ, mtxZ,
mtxB, itermax, conv_tol, msglvl, msgFile);
DenseMtx_free(mtxQ) ;
break;
案例MLBiCGStabL:
mtxQ = DenseMtx_new() ;
DenseMtx_init(mtxQ, type, 0, -1, neqns, Ik, 1, neqns) ;
Drand_setUniform(&drand, 0.0, 1.0) ;
DenseMtx_fillRandomEntries(mtxQ, &drand) ;
if (zversion)
rc=zmlbicgstabl(neqns, type, symmetryflag, mtxA, frontmtx, mtxQ, mtxZ,
mtxB, itermax, conv_tol, msglvl, msgFile);
else
rc=mlbicgstabl(neqns, type, symmetryflag, mtxA, frontmtx, mtxQ, mtxZ,
mtxB, itermax, conv_tol, msglvl, msgFile);
DenseMtx_free(mtxQ) ;
break;
case BGMRESR:
if (zversion)
fprintf(msgFile, "\n\n *** BGMRESR complex version is not available "
"at this moment. ") ;
else
rc=bgmresr(neqns, type, symmetryflag, mtxA, frontmtx, mtxZ,
mtxB, iterout, itermax, &nouter, &ninner, conv_tol,
msglvl, msgFile);
break;
case BGMRESL:
if (zversion)
fprintf(msgFile, "\n\n *** BGMRESR complex version is not available "
"at this moment. ") ;
else
rc=bgmresl(neqns, type, symmetryflag, mtxA, frontmtx, mtxZ,
mtxB, iterout, itermax, &nouter, &ninner, conv_tol,
msglvl, msgFile);
break;
default:
fprintf(msgFile, "\n\n *** Invalid method number ") ;
}

MARKTIME(t2) ;
fprintf(msgFile, "\n\n CPU %8.3f : solve the system", t2 - t1) ;
if ( msglvl > 2 ) {
fprintf(msgFile, "\n\n computed solution") ;
DenseMtx_writeForHumanEye(mtxZ, msgFile) ;
fflush(stdout) ;
}

/*
-------------------------------------------------------------
permute the computed solution back into the original ordering
-------------------------------------------------------------
*/
MARKTIME(t1) ;
DenseMtx_permuteRows(mtxZ, newToOldIV) ;
MARKTIME(t2) ;
fprintf(msgFile, "\n CPU %8.3f : permute solution", t2 - t1) ;
if ( msglvl > 2 ) {
fprintf(msgFile, "\n\n permuted solution") ;
DenseMtx_writeForHumanEye(mtxZ, msgFile) ;
fflush(stdout) ;
}
/*
-------------
save solution
x should write_for_matlab
-------------
*/
if ( strcmp(slnFileName, "none") != 0 )
{
DenseMtx_writeToFile(mtxZ, slnFileName) ;
}
// ljp *************** add 3-24-2009
DenseMtx_writeForMatlab(mtxZ, "z", solFileName) ;
// ljp *************** add 3-24-2009

/*
-----------------
compute the error
-----------------
*/
if (!strcmp(rhsFileName, "none")) {
DenseMtx_sub(mtxZ, mtxX) ;
if (method[k] <8) {
mtxQ = DenseMtx_new() ;
DenseMtx_init(mtxQ, type, 0, -1, neqns, 1, 1, neqns) ;
rc=DenseMtx_initAsSubmatrix (mtxQ, mtxZ, 0, neqns-1, 0, 0);
fprintf(msgFile, "\n\n maxabs error = %12.4e", DenseMtx_maxabs(mtxQ)) ;
DenseMtx_free(mtxQ) ;
}
else
fprintf(msgFile, "\n\n maxabs error = %12.4e", DenseMtx_maxabs(mtxZ)) ;

if ( msglvl > 1 ) {
fprintf(msgFile, "\n\n error") ;
DenseMtx_writeForHumanEye(mtxZ, msgFile) ;
fflush(stdout) ;
}
if ( msglvl > 1 )
SubMtxManager_writeForHumanEye(frontmtx->manager, msgFile) ;
}
fprintf(msgFile, "\n--------- End of Method %d -------\n",method[k]) ;

}
/*
------------------------
free the working storage
------------------------
*/
DenseMtx_free(mtxZ) ;

end_front:
ChvManager_free(chvmanager) ;
SubMtxManager_free(mtxmanager) ;
FrontMtx_free(frontmtx) ;
IVL_free(symbfacIVL) ;
IV_free(oldToNewIV) ;
IV_free(newToOldIV) ;

end_tree:
ETree_free(frontETree) ;

end_init:
DenseMtx_free(mtxB) ;
DenseMtx_free(mtxX) ;

end_read:
InpMtx_free(mtxA) ;

fprintf(msgFile, "\n") ;
fclose(msgFile) ;

// ljp return(rc) ;
}

/*--------------------------------------------------------------------*/
/*
-------------------------------------------------------------------
InpMtx_readFromAIJ2file reads in an AIJ2 matrix from file and
converts to InpMtx object.

created -- 98dec01, jwu
-------------------------------------------------------------------
*/
int
InpMtx_readFromAIJ2file(InpMtx *mtxA, char *srcFileName)
{

FILE *srcFile;
int rc, nrow, ncol, irow, icol, jrow, jcol, ient, nentInRow;
double ent;
/*
----------------------------
open the input file and read
#rows #columns
----------------------------
*/
if ( (srcFile = fopen(srcFileName, "r")) == NULL ) {
fprintf(stderr, "\n fatal error in InpMtx_readFromAIJ2file"
"\n unable to open file %s\n",
srcFileName) ;
return(-1) ;
}
rc = fscanf(srcFile, "%d %d", &nrow, &ncol) ;
if ( rc != 2 ) {
fprintf(stderr, "\n fatal error in InpMtx_readFromAIJ2file"
"\n %d of 2 fields read on first line of file %s",
rc, srcFileName) ;
return(-1) ;
}

/*
--------------------------------------------------
initialize the object
set coordType = INPMTX_BY_ROWS --> row coordinates
--------------------------------------------------
*/
InpMtx_init(mtxA, INPMTX_BY_ROWS, SPOOLES_REAL, 10*nrow, 0) ;
/*
-------------------------------------------------
read in the entries and load them into the object
-------------------------------------------------
*/
ient=0;
for ( irow = 0 ; irow < nrow ; irow++ ) {
rc = fscanf(srcFile, "%d %d", &jrow, &nentInRow) ;
if ( rc != 2 ) {
fprintf(stderr, "\n fatal error in InpMtx_readFromAIJ2file"
"\n %d of 2 fields read on entry %d of file %s",
rc, ient, srcFileName) ;
return(-1) ;
}
for ( icol = 0 ; icol < nentInRow ; icol++ ) {
rc = fscanf(srcFile, "%d %lf", &jcol, &ent) ;
if ( rc != 2 ) {
fprintf(stderr, "\n fatal error in InpMtx_readFromAIJ2file"
"\n %d of 2 fields read on entry %d of file %s",
rc, ient, srcFileName) ;
return(-1) ;
}
InpMtx_inputRealEntry(mtxA, jrow-1, jcol-1, ent) ;
ient++;
}
}
/*
-----------------------------
sort and compress the entries
-----------------------------
*/
InpMtx_changeStorageMode(mtxA, 3) ;

fclose(srcFile);
return(1) ; }

/*--------------------------------------------------------------------*/

/* createGraph.c */
/*--------------------------------------------------------------------*/
int
InpMtx_createGraph (InpMtx *mtxA, Graph *graph)
/*
----------------------------------------------------
read in a InpMtx object and create the Graph object

created -- 97feb14, cca
----------------------------------------------------
*/
{
int count, msglvl, nvtx, rc ;
IVL *adjIVL ;

InpMtx_changeStorageMode(mtxA, 3) ;
nvtx = 1 + IV_max(&mtxA->ivec1IV) ;
count = 1 + IV_max(&mtxA->ivec2IV) ;
if ( nvtx < count ) {
nvtx = count ;
}
/*
------------------------------------
create the full adjacency IVL object
------------------------------------
*/
adjIVL = InpMtx_fullAdjacency(mtxA) ;
/*
---------------------
fill the Graph object
---------------------
*/
Graph_init2(graph, 0, nvtx, 0, adjIVL->tsize, nvtx, adjIVL->tsize,
adjIVL, NULL, NULL) ;

return(1) ; }

/*--------------------------------------------------------------------*/

thank you very much !
thanks a lot
Best regards
yours
Jianping liao