mainspdexpert.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 the system according to some independent set ordering
//#define IND_SET

// 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 edges
//#define METIS_E

// reorder system using MeTiS multilevel nested dissection by nodes
//#define METIS_N

// reorder the columns and rows of a system differently using a new unsymmetric
// reordering strategy by Yousef Saad
//#define PP_PERM

// mixed strategies that finally switch to PQ if necessary
//#define MMD_PP
//#define AMF_PP
//#define AMD_PP
//#define RCM_PP
//#define FC_PP
//#define METIS_E_PP
#define METIS_N_PP


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


// variant of PILUC that uses a repeated multiple factorization approach
//#define USE_MPILUC






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

    CSRMAT A, ilutmat;
    AMGLEVELMAT PRE, *next;
    int nlev=0, nprev, nB;
    int (*perm0)(),(*perm)(),(*permf)();

    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,m,fnamelen,n,nc,nz,nrhs,tmp0,tmp,tmp2,tmp3,ierr,flag,
      *invperm, *buff, *ws, *symmmd,flags,elbow,max_it,
         nrhsix, *rhsptr, *rhsind;
    REALS  eps, DROP_TOL, CONDEST,condest,droptols[2],restol,
          val,vb, *rbuff;
    FLOAT *rhs,*sol,*w,  *scale, *rhsval, *dbuff;
    float  systime, time_start,   time_stop,   secnds, secndsprep,
           timeAx_start, timeAx_stop, secndsAx;
    int ELBOW, nnzU, l, nAx;
    ILUPACKPARAM param;
    size_t   nibuff, ndbuff;

    /* the last argument passed serves as file name */
    if (argc!=5) {
      printf("usage '%s <drop tol.> <bound for L^{-1},U^{-1}> <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;
    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]);
    CONDEST=atof(argv[argc-3]);
    DROP_TOL=atof(argv[argc-4]);


/*----------------------------------------------------------------------
|     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,A.a,A.ja,A.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|%7.1e|",foname,DROP_TOL,CONDEST);

    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);
    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");
    A.nr=n;
    A.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;
    

    // set parameters to the default settings
    SPDAMGINIT(A, &param);

    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,A.a,A.ja,A.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=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++) 
#if defined _DOUBLE_REAL_ || defined _SINGLE_REAL_
	printf("%16.1e",A.a[j-1]);
#else
	printf("%8.1e%8.1e",A.a[j-1].r,A.a[j-1].i);
#endif
      printf("\n");
      fflush(stdout);
    }// end for i
    */


    SPDAMGINIT(A, &param);
    param.dbuff=dbuff;
    param.ndbuff=ndbuff;
    
    perm0=SPDPERMNULL;
    perm =SPDPERMNULL;
    permf=SPDPERMNULL;
#ifdef MINIMUM_DEGREE
    perm0=SPDPERMMMD;
    perm =SPDPERMMMD;
    permf=SPDPERMMMD;
    //printf("prescribe minimum degree ordering\n");
    fprintf(fo,"mmds/mmds|");
#elif defined REVERSE_CM
    perm0=SPDPERMRCM;
    perm =SPDPERMRCM;
    permf=SPDPERMRCM;
    fprintf(fo,"rcms/rcms|");
    //printf("prescribe reverse Cuthill-McKee ordering\n");
#elif defined NESTED_DISSECTION
    perm0=SPDPERMND;
    perm =SPDPERMND;
    permf=SPDPERMND;
    fprintf(fo,"nds /nds |");
    //printf("prescribe nested dissection ordering\n");
#elif defined IND_SET
    perm0=SPDPERMINDSET;
    perm =SPDPERMINDSET;
    permf=SPDPERMINDSET;
    fprintf(fo,"inds/inds|");
    //printf("prescribe independent set ordering\n");
#elif defined AMF
    perm0=SPDPERMAMF;
    perm =SPDPERMAMF;
    permf=SPDPERMAMF;
    fprintf(fo,"amfs/amfs|");
    //printf("prescribe independent set ordering\n");
#elif defined PP_PERM
    perm0=SPDPERMPP;
    perm =SPDPERMPP;
    permf=SPDPERMPP;
    fprintf(fo,"PPs /PPs |");
    //printf("prescribe PP ordering\n");
#elif defined MMD_PP
    perm0=SPDPERMMMD;
    perm =SPDPERMMMD;
    permf=SPDPERMPP;
    fprintf(fo,"mmds/PPs |");
    //printf("prescribe MMD/PP ordering\n");
#elif defined AMF_PP
    perm0=SPDPERMAMF;
    perm =SPDPERMAMF;
    permf=SPDPERMPP;
    fprintf(fo,"amfs/PPs |");
    //printf("prescribe MMD/PP ordering\n");
#elif defined RCM_PP
    perm0=SPDPERMRCM;
    perm =SPDPERMRCM;
    permf=SPDPERMPP;
    fprintf(fo,"rcms/PPs |");
    //printf("prescribe MMD/PP ordering\n");
#else
    fprintf(fo,"    /    |");
#endif


    // modify the default settings
    SPDAMGGETPARAMS(param, &flags, &elbow, droptols, &condest,
		    &restol, &max_it);

    // ONLY for mixed reordering strategies it is useful to
    // set the 'FINAL_PIVOTING' flag
#if !defined RCM_PP && !defined AMF_PP && !defined AMD_PP && !defined MMD_PP && !defined FC_PP && !defined METIS_E_PP && !defined METIS_N_PP 
    flags&=~FINAL_PIVOTING;
#endif
    // change flags if mpiluc is desired
#ifdef USE_MPILUC
    flags|=MULTI_PILUC;
#endif

    // flags
    //flags|=TISMENETSKY_SC;
    //flags&=~ENSURE_SPD;

    // overwrite the default drop tolerances
    droptols[0]=1.0;
    droptols[1]=DROP_TOL; 

    // overwrite the default value for elbow space
    elbow=ELBOW;

    // overwrite default values for condest
    condest=CONDEST;

    // overwrite the default value for the residual norm
    restol=1e-8;

    // overwrite the default value for the maximum number of iterations
    max_it=MIN(1.1*n+10,MAX_IT);

    // rewrite the updated parameters