mainexpert.c

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

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

#include <blas.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 the columns and rows of a system differently using a new unsymmetric
// reordering strategy by Yousef Saad
//#define PQ_PERM


// mixed strategies that finally switch to PQ if necessary
//#define MMD_PQ
//#define AMF_PQ
//#define AMD_PQ
//#define RCM_PQ
//#define FC_PQ
//#define METIS_E_PQ
//#define METIS_N_PQ

// Minimum Weight Matching interface MC64
//#define MC64_RCM_PQ
//#define MC64_MMD_PQ
//#define MC64_AMF_PQ
//#define MC64_AMD_PQ
//#define MC64_METIS_E_PQ
//#define MC64_METIS_N_PQ

// alternative Minimum Weight Matching interface provided by PARDISO
//#define MWM_RCM_PQ
//#define MWM_MMD_PQ
//#define MWM_AMD_PQ
//#define MWM_AMF_PQ
//#define MWM_METIS_E_PQ
#define MWM_METIS_N_PQ


// 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=8; /* gmres */


    CSRMAT fullmat, A;
    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,l,fnamelen,n,m,nc,nz,nrhs,tmp0,tmp,tmp2,tmp3,ierr,flags,
         *invperm, *buff, *ws, *symmmd, *p, *invp, *q, *invq,nLU,
         nrhsix, *rhsptr, *rhsind;
    REALS eps, DROP_TOL, condest, CONDEST,droptols[2],amgcancel,val,
           vb,*rbuff,restol;
    FLOAT *rhs,*sol,*w, *leftscale, *rightscale, *rhsval, *dbuff;
    float  systime, time_start,   time_stop,   secnds,
           timeAx_start, timeAx_stop, secndsAx;
    int ELBOW, elbow, nnzL,nnzU, nAx=0, nrestart, max_it;
    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;
    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]);
    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,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|%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);
    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);



    // set parameters to the default settings
    AMGINIT(A, &param);
    param.dbuff=dbuff;
    param.ndbuff=ndbuff;
    
    perm0=PERMNULL;
    perm =PERMNULL;
    permf=PERMNULL;
#ifdef MINIMUM_DEGREE
    perm0=PERMMMD;
    perm =PERMMMD;
    permf=PERMMMD;
    //printf("prescribe minimum degree ordering\n");
    fprintf(fo,"mmds/mmds|");
#elif defined REVERSE_CM
    perm0=PERMRCM;
    perm =PERMRCM;
    permf=PERMRCM;
    fprintf(fo,"rcms/rcms|");
    //printf("prescribe reverse Cuthill-McKee ordering\n");
#elif defined NESTED_DISSECTION
    perm0=PERMND;
    perm =PERMND;
    permf=PERMND;
    fprintf(fo,"nds /nds |");
    //printf("prescribe nested dissection ordering\n");
#elif defined IND_SET
    perm0=PERMINDSET;
    perm =PERMINDSET;
    permf=PERMINDSET;
    fprintf(fo,"inds/inds|");
    //printf("prescribe independent set ordering\n");
#elif defined AMF
    perm0=PERMAMF;
    perm =PERMAMF;
    permf=PERMAMF;
    fprintf(fo,"amfs/amfs|");
#elif defined AMD
    perm0=PERMAMD;
    perm =PERMAMD;
    permf=PERMAMD;
    fprintf(fo,"amds/amds|");
#elif defined PQ_PERM
    perm0=PERMPQ;
    perm =PERMPQ;
    permf=PERMPQ;
    fprintf(fo,"PQs /PQs |");
    //printf("prescribe PQ ordering\n");
#elif defined MMD_PQ
    perm0=PERMMMD;
    perm =PERMMMD;
    permf=PERMPQ;
    fprintf(fo,"mmds/PQs |");
    //printf("prescribe MMD/PQ ordering\n");
#elif defined AMF_PQ
    perm0=PERMAMF;
    perm =PERMAMF;
    permf=PERMPQ;
    fprintf(fo,"amfs/PQs |");
#elif defined AMD_PQ
    perm0=PERMAMD;
    perm =PERMAMD;
    permf=PERMPQ;
    fprintf(fo,"amds/PQs |");
#elif defined RCM_PQ
    perm0=PERMRCM;
    perm =PERMRCM;
    permf=PERMPQ;
    fprintf(fo,"rcms/PQs |");
    //printf("prescribe MMD/PQ ordering\n");
#elif defined FC_PQ
    perm0=PERMFC;
    perm =PERMFC;
    permf=PERMPQ;
    fprintf(fo,"FCs /PQs |");
#elif defined METIS_E_PQ
    perm0=PERMMETISE;
    perm =PERMMETISE;
    permf=PERMPQ;
    fprintf(fo,"mes /PQs |");
#elif defined METIS_N_PQ
    perm0=PERMMETISN;
    perm =PERMMETISN;
    permf=PERMPQ;
    fprintf(fo,"mns /PQs |");
#elif defined MC64_RCM_PQ
    perm0=PERMMC64RCM;
    perm =PERMRCM;
    permf=PERMPQ;
    fprintf(fo,"mc64rc/PQ|");
    //printf("prescribe MC64 ordering\n");
#elif defined MC64_MMD_PQ
    perm0=PERMMC64MMD;
    perm =PERMMMD;
    permf=PERMPQ;
    fprintf(fo,"mc64md/PQ|");
    //printf("prescribe MC64 ordering\n");
#elif defined MC64_AMF_PQ
    perm0=PERMMC64AMF;
    perm =PERMAMF;
    permf=PERMPQ;
    fprintf(fo,"mc64af/PQ|");
    //printf("prescribe MC64 ordering\n");
#elif defined MC64_AMD_PQ
    perm0=PERMMC64AMD;
    perm =PERMAMD;
    permf=PERMPQ;
    fprintf(fo,"mc64ad/PQ|");
    //printf("prescribe MC64 ordering\n");
#elif defined MC64_METIS_E_PQ
    perm0=PERMMC64METISE;
    perm =PERMMETISE;
    permf=PERMPQ;
    fprintf(fo,"mc64me/PQ|");
    //printf("prescribe MC64 ordering\n");
#elif defined MC64_METIS_N_PQ
    perm0=PERMMC64METISN;
    perm =PERMMETISN;
    permf=PERMPQ;
    fprintf(fo,"mc64mn/PQ|");