mainiluc.c

数学算法的实现库。可以实现常见的线性计算。

#include <stdio.h>
#include <stdlib.h>
#include <string.h>
#include <math.h>


#include <blas.h>
#include <sparspak.h>
#include <ilupack.h>

#include <ilupackmacros.h>




#define MAX_LINE        255
#define STDERR          stdout
#define STDOUT          stdout
#define PRINT_INFO
#define MAX(A,B)        (((A)>(B))?(A):(B))
#define MIN(A,B)        (((A)<(B))?(A):(B))

// maximum number of iterations independent on n
#define MAX_IT          500

// measure for terminating iterative solver
#define RESTOL_FUNC(A)  sqrt(A)
//#define RESTOL_FUNC(A)  (A)





// reorder the system according to the symmetric minimum degree ordering
//#define MINIMUM_DEGREE

// reorder the system according to the nested dissection ordering
//#define NESTED_DISSECTION 

// reorder the system according to the reverse Cuthill-McKee ordering
//#define REVERSE_CM

// reorder system using approximate minimum fill by Patrick Amestoy, 
// Tim Davis and Iain Duff
//#define AMF

// reorder system using approximate minimum degree by Patrick Amestoy
// Tim Davis and Iain Duff
//#define AMD

// reorder system using METIS multilevel nested dissection by Karypis & Kumar
// edge approach
//#define METIS_E
// nodal approach
//#define METIS_N

// Minimum Weight Matching interface MC64
//#define MC64_RCM
//#define MC64_MMD
//#define MC64_AMF
//#define MC64_AMD
//#define MC64_METIS_E
//#define MC64_METIS_N

// alternative Minimum Weight Matching interface provided by PARDISO
//#define MWM_RCM
//#define MWM_MMD
//#define MWM_AMF
//#define MWM_AMD
//#define MWM_METIS_E
#define MWM_METIS_N


// use an iterative solver from SPARSKIT defined by variable SOLVER
#define USE_SPARSKIT





int main(int argc, char **argv)
{
    /* SOLVER choice:   1  cg
                        2  cgnr
                        3  bcg
                        4  dbcg
                        5  bcgstab
                        6  tfqmr
                        7  fom
                        8  gmres
                        9  fgmres
                       10  dqgmres */
    int SOLVER=8; /* gmres */

    CSRMAT fullmat, ilutmat, A;

    FILE *fp, *fo; 
    char rhstyp[3], title[72], key[8], type[3], fname[100], foname[100];
    char line[MAX_LINE], *tmpstring, timeout[7], residout[7];
    int  i,j,k,fnamelen,n,nc,nz,nrhs,tmp0,tmp,tmp2,tmp3,ierr,ipar[20],flag,
         *p, *invq, *invperm, *ws, *symmmd,m,flags,
         nrhsix, *rhsptr, *rhsind;
    FLOAT *rhs,*sol,*w, *rowscale, *colscale, *rhsval, *dbuff;
    float  systime, time_start,   time_stop,   secnds, secndsprep;
    REALS eps,fpar[20], DROP_TOL, val, vb,*rbuff,nrm;
    int ELBOW, nB;
    ILUPACKPARAM param;
    size_t nibuff, ndbuff;

    /* the last argument passed serves as file name */
    if (argc!=4) {
      printf("usage '%s <drop tol.> <elbow space> <matrix>'\n",argv[0]);
       exit(0);
    }
    i=0;
    while (argv[argc-1][i]!='\0')
    {
          fname[i]=argv[argc-1][i];
          i++;
    } /* end while */
    fname[i]='\0';
    fnamelen=i;
    i=fnamelen-1;
    while (i>=0 && fname[i]!='/')
          i--;
    i++;
    j=0;
    while (i<fnamelen && fname[i]!='.')
          foname[j++]=fname[i++];
    while (j<16)
          foname[j++]=' ';
    foname[j]='\0';

    ELBOW=atoi(argv[argc-2]);
    DROP_TOL=atof(argv[argc-3]);


/*----------------------------------------------------------------------
|     Read a Harwell-Boeing matrix.
|     Use readmtc first time to determine sizes of arrays.
|     Read in values on the second call.
|---------------------------------------------------------------------*/

    nrhs = 0;
    tmp0 = 0;

    if ((fp=fopen(fname,"r"))==NULL) {
        fprintf(STDERR," file %s not found\n",fname);
        exit(0);
    }
    fclose(fp);

    READMTC(&tmp0,&tmp0,&tmp0,fname,fullmat.a,fullmat.ja,fullmat.ia,
	    rhs,&nrhs,rhstyp,&n,&nc,&nz,title,key,type,
	    &nrhsix,rhsptr,rhsind,rhsval,&ierr,fnamelen,2,72,8,3);
    if (ierr) {
        fprintf(STDERR," ierr = %d\n",ierr);
        fprintf(STDERR," error in reading the matrix, stop.\n");
	switch(ierr) {
	case 1:
	  fprintf(STDERR,"too many columns\n");
	  break;  
	case 2:
	  fprintf(STDERR,"too many nonzeros\n");
	  break;  
	case 3:
	  fprintf(STDERR,"too many columns and nonzeros\n");
	  break;  
	case 4:
	  fprintf(STDERR,"right hand side has incompatible type\n");
	  break;  
	case 5:
	  fprintf(STDERR,"too many right hand side entries\n");
	  break;  
	case 6:
	  fprintf(STDERR,"wrong type (real/complex)\n");
	  break;  
	}
        exit(ierr);
    }
    printf("Matrix: %s: size (%d,%d), nnz=%d(%4.1lfav.)\n", fname, n,nc,
	   nz,((double)nz)/n);

    if (fname[fnamelen-1]=='5')
      fo = fopen("out_mc64","aw");
    else
      fo = fopen("out_normal","aw");

    fprintf(fo,"%s|%7.1e|       |",foname,DROP_TOL);

    m=1;
    if (nrhs>0) {
       printf ("Number of right hand sides supplied: %d \n", nrhs) ;
       if (rhstyp[1]=='G' || rhstyp[1]=='g') {
	 m++;
	 printf ("Initial solution(s) offered\n") ;
       }
       else
	 printf ("\n") ;

       if (rhstyp[2]=='X' || rhstyp[2]=='x') {
	 m++;
	 printf ("Exact solution(s) provided\n") ;
       }
       else
	 printf ("\n") ;
    }
    else 
      printf("\n\n\n");
    printf("\n");

    rhsptr=NULL;
    rhsind=NULL;
    rhsval=NULL;
    if (rhstyp[0]=='M' || rhstyp[0]=='m') {
       rhsptr=(int *)  MALLOC((size_t)(nrhs+1)*sizeof(int),"main:rhsptr");
       rhsind=(int *)  MALLOC((size_t)nrhsix*sizeof(int),  "main:rhsind");
       rhsval=(FLOAT *)MALLOC((size_t)nrhsix*sizeof(FLOAT),"main:rhsval");
       // no dense right hand side
       m--;
       m*=n*MAX(1,nrhs);
       // in any case we need at least one vector for the r.h.s.
       m+=n;
    }
    else
      m*=n*MAX(1,nrhs);
    fullmat.ia=(int *)  MALLOC((size_t)(n+1)*sizeof(int),"main:fullmat.ia");
    fullmat.ja=(int *)  MALLOC((size_t)nz*sizeof(int),   "main:fullmat.ja");
    fullmat.a =(FLOAT *)MALLOC((size_t)nz*sizeof(FLOAT), "main:fullmat.a");
    fullmat.nr=n;
    fullmat.nc=n;
    rhs       =(FLOAT *) MALLOC((size_t)m*sizeof(FLOAT),  "main:rhs");
    // advance pointer to reserve space when uncompressing the right hand side
    if (rhstyp[0]=='M' || rhstyp[0]=='m')
       rhs+=n;

    sol  =(FLOAT *) MALLOC((size_t)n*sizeof(FLOAT),  "main:sol");
    dbuff=(FLOAT *) MALLOC(3*(size_t)n*sizeof(FLOAT),"main:dbuff");
    ndbuff=3*(size_t)n;

    tmp = 3;
    tmp2 = n;
    tmp3 = nz;

    if (rhstyp[0]=='M' || rhstyp[0]=='m')
       m-=n;
    READMTC(&tmp2,&tmp3,&tmp,fname,fullmat.a,fullmat.ja,fullmat.ia,
	    rhs,&m,rhstyp,&n,&nc,&nz,title,key,type,
	    &nrhsix,rhsptr,rhsind,rhsval,&ierr,fnamelen,2,72,8,3);
    if (rhstyp[0]=='M' || rhstyp[0]=='m')
       m+=n;



    if (ierr) {
        fprintf(STDERR," ierr = %d\n",ierr);
        fprintf(STDERR," error in reading the matrix, stop.\n");
	fprintf(fo,"out of memory\n");fclose(fo); 
        exit(ierr);
    }
    /*
    for (i=0; i<n;i++) {
      printf("%4d:\n",i+1);
      for (j=fullmat.ia[i];j<fullmat.ia[i+1]; j++)
	printf("%16d",fullmat.ja[j-1]);
      printf("\n");
      fflush(stdout);
      for (j=fullmat.ia[i];j<fullmat.ia[i+1]; j++) 
#if defined _DOUBLE_REAL_ || defined _SINGLE_REAL_
	printf("%16.1e",fullmat.a[j-1]);
#else
	printf("%8.1e%8.1e",fullmat.a[j-1].r,fullmat.a[j-1].i);
#endif
      printf("\n");
      fflush(stdout);
    }// end for i
    */
/*----------------------------------------------------------------------
|     Convert the matrix from csc to csr format.  First, allocate
|     space in (fullmat).  After conversion, free space occupied by
|     initial format.
|---------------------------------------------------------------------*/

    A.nr=A.nc=n;
    A.ia=(int *)  MALLOC((size_t)(n+1)*sizeof(int),"main:A.ia");
    A.ja=(int *)  MALLOC((size_t)nz*sizeof(int),   "main:A.ja");
    A.a =(FLOAT *)MALLOC((size_t)nz*sizeof(FLOAT), "main:A.a");

    tmp = 1;
    CSRCSC(&n,&tmp,&tmp,fullmat.a,fullmat.ja,fullmat.ia,A.a,A.ja,A.ia);
    /*
    for (i=0; i<n;i++) {
      printf("%4d:\n",i+1);
      for (j=A.ia[i];j<A.ia[i+1]; j++)
	printf("%16d",A.ja[j-1]);
      printf("\n");
      fflush(stdout);
      for (j=A.ia[i];j<A.ia[i+1]; j++) 
	printf("%8.1e%8.1e",A.a[j-1].r,A.a[j-1].i);
      printf("\n");
      fflush(stdout);
    }// end for i
    */
    free(fullmat.a);
    free(fullmat.ja);
    free(fullmat.ia);

    // release part of rhs that may store the uncompressed rhs
    if (rhstyp[0]=='M' || rhstyp[0]=='m')
       rhs-=n;

    // if no right hand side is provided, then set sol=1 and rhs=A*sol
    if (nrhs==0) { 
       for (i=0; i<n; i++) {
#if defined _DOUBLE_REAL_ || defined _SINGLE_REAL_
	   sol[i]=1;
#else
	   sol[i].r=1;
	   sol[i].i=0;
#endif
       } // end for i
       CSRMATVEC(A,sol,rhs);
    }
    else {
       if (rhstyp[0]=='M' || rhstyp[0]=='m') {
	  for (i=0; i<n; i++) {
#if defined _DOUBLE_REAL_ || defined _SINGLE_REAL_ 
	      rhs[i]=0;
#else
	      rhs[i].r=rhs[i].i=0;
#endif
	  }// end for i

	  // uncompress rhs
	  for (i=0; i<rhsptr[1]-rhsptr[0]; i++) {
	      j=rhsind[i]-1;
	      rhs[j]=rhsval[i];
	  } // end for i
       } // end if
    } // end if-else
    
    p   =(int *)MALLOC((size_t)n*sizeof(int),"mainilutp:p");
    invq=(int *)MALLOC((size_t)n*sizeof(int),"mainilutp:invq");
    rowscale=(FLOAT *)MALLOC((size_t)n*sizeof(FLOAT),"mainilutp:rowscale");
    colscale=(FLOAT *)MALLOC((size_t)n*sizeof(FLOAT),"mainilutp:colscale");
    for (i=0; i<n; i++) {
        p[i]=invq[i]=i+1;
#if defined _DOUBLE_REAL_ || defined _SINGLE_REAL_
        rowscale[i]=colscale[i]=1.0;
#else
        rowscale[i].r=colscale[i].r=1.0;
        rowscale[i].i=colscale[i].i=0.0;
#endif
    } // end for i
    ierr=0;
    nB=n;


    evaluate_time(&time_start,&systime);
    // set parameters to the default settings
    AMGINIT(A, &param);
    param.dbuff=dbuff;
    param.ndbuff=ndbuff;
    // turn off column scaling
    //param.ipar[7]&=~(2+16+128);
#if defined MINIMUM_DEGREE 
    ierr=PERMMMD(A,rowscale,colscale,p,invq,&nB,&param);
    fprintf(fo,"mmds/mmds|");
    printf("prescribe minimum degree ordering\n");
#elif defined NESTED_DISSECTION 
    ierr=PERMND(A,rowscale,colscale,p,invq,&nB,&param);
    fprintf(fo,"nds /nds |");
    printf("prescribe nested dissection ordering\n");
#elif defined REVERSE_CM 
    ierr=PERMRCM(A,rowscale,colscale,p,invq,&nB,&param);
    fprintf(fo,"rcms/rcms|");
    printf("prescribe reverse Cuthill-McKee ordering\n");
#elif defined AMF
    ierr=PERMAMF(A,rowscale,colscale,p,invq,&nB,&param);
    fprintf(fo,"amfs/amfs|");
    printf("prescribe approximate minimum fill ordering\n");
#elif defined AMD
    ierr=PERMAMD(A,rowscale,colscale,p,invq,&nB,&param);
    fprintf(fo,"amds/amds|");
    printf("prescribe approximate minimum degree ordering\n");
#elif defined METIS_E
    ierr=PERMMETISE(A,rowscale,colscale,p,invq,&nB,&param);
    fprintf(fo,"mes /mes |");
    printf("prescribe multilevel (edge) nested dissection ordering\n");
#elif defined METIS_N
    ierr=PERMMETISN(A,rowscale,colscale,p,invq,&nB,&param);
    fprintf(fo,"mns /mns |");
    printf("prescribe multilevel (node) nested dissection ordering\n");
#elif defined MWM_RCM
    ierr=PERMMWMRCM(A,rowscale,colscale,p,invq,&nB,&param);
    fprintf(fo,"mwrc/mwrc|");
    printf("prescribe MWM + RCM ordering\n");
#elif defined MWM_AMF
    ierr=PERMMWMAMF(A,rowscale,colscale,p,invq,&nB,&param);
    fprintf(fo,"mwaf/mwaf|");
    printf("prescribe MWM + AMF ordering\n");
#elif defined MWM_AMD
    ierr=PERMMWMAMD(A,rowscale,colscale,p,invq,&nB,&param);
    fprintf(fo,"mwad/mwad|");
    printf("prescribe MWM + AMD ordering\n");
#elif defined MWM_MMD
    ierr=PERMMWMMMD(A,rowscale,colscale,p,invq,&nB,&param);
    fprintf(fo,"mwmd/mwmd|");
    printf("prescribe MWM + MMD ordering\n");
#elif defined MWM_METIS_E
    ierr=PERMMWMMETISE(A,rowscale,colscale,p,invq,&nB,&param);
    fprintf(fo,"mwme/mwme|");
    printf("prescribe MWM + METIS multilvel (edge) ND ordering\n");
#elif defined MWM_METIS_N
    ierr=PERMMWMMETISN(A,rowscale,colscale,p,invq,&nB,&param);
    fprintf(fo,"mwmn/mwmn|");
    printf("prescribe MWM + METIS multilvel (node) ND ordering\n");
#elif defined MC64_RCM
    ierr=PERMMC64RCM(A,rowscale,colscale,p,invq,&nB,&param);
    fprintf(fo,"mcrc/mcrc|");
    printf("prescribe MC64 + RCM ordering\n");
#elif defined MC64_AMF
    ierr=PERMMC64AMF(A,rowscale,colscale,p,invq,&nB,&param);
    fprintf(fo,"mcaf/mcaf|");
    printf("prescribe MC64 + AMF ordering\n");
#elif defined MC64_AMD
    ierr=PERMMC64AMD(A,rowscale,colscale,p,invq,&nB,&param);
    fprintf(fo,"mcad/mcad|");
    printf("prescribe MC64 + AMD ordering\n");
#elif defined MC64_MMD
    ierr=PERMMC64MMD(A,rowscale,colscale,p,invq,&nB,&param);
    fprintf(fo,"mcmd/mcmd|");
    printf("prescribe MC64 + MMD ordering\n");
#elif defined MC64_METIS_E
    ierr=PERMMC64METISE(A,rowscale,colscale,p,invq,&nB,&param);
    fprintf(fo,"mcme/mcme|");
    printf("prescribe MC64 + METIS multilvel (edge) ND ordering\n");
#elif defined MC64_METIS_N
    ierr=PERMMC64METISN(A,rowscale,colscale,p,invq,&nB,&param);
    fprintf(fo,"mcmn/mcmn|");
    printf("prescribe MC64 + METIS multilvel (node) ND ordering\n");
#else
    fprintf(fo,"    /    |");
#endif


    // rescale right hand side
    for (i=0; i<n; i++) {
#if defined _DOUBLE_REAL_ || defined _SINGLE_REAL_
        rhs[i]*=rowscale[i];
#else
	val=rhs[i].r;
	rhs[i].r=val*rowscale[i].r-rhs[i].i*rowscale[i].i;
	rhs[i].i=val*rowscale[i].i+rhs[i].i*rowscale[i].r;
#endif
    } // end for i
    // reorder right hand side in the same way as the rows of A
    for (i=0; i<n; i++)
        param.dbuff[i]=rhs[p[i]-1];
    i=1;
    COPY(&n,param.dbuff,&i,rhs,&i);


    /*
    for (i=0; i<n; i++) 
      printf("%8d",p[i]);
    printf("\n");fflush(stdout);
    for (i=0; i<n; i++) {
        for (j=A.ia[i]; j<A.ia[i+1]; j++) {
	  printf("%8d",A.ja[j-1]);
	} // end for j
	printf("\n");fflush(stdout);
        for (j=A.ia[i]; j<A.ia[i+1]; j++) {
	  printf("%8.1e",A.a[j-1]);
	} // end for j
	printf("\n");fflush(stdout);
    } // end for i
    */

    // permutation for A
    CPERM(&A,invq);
    RPERM(&A,p);

    /*
    for (i=0; i<n; i++) {
        for (j=A.ia[i]; j<A.ia[i+1]; j++) {
	  printf("%8d",A.ja[j-1]);
	} // end for j
	printf("\n");fflush(stdout);
        for (j=A.ia[i]; j<A.ia[i+1]; j++) {
	  printf("%8.1e",A.a[j-1]);
	} // end for j
	printf("\n");fflush(stdout);
    } // end for i
    */

    evaluate_time(&time_stop,&systime);
    secndsprep=time_stop-time_start;
    printf("time: %8.1e [sec]\n",(double)secndsprep);



    /* ------------------------------------------------------------------- */
    /* -----------------------   Start MAIN part   ----------------------- */