main3.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

// Maximum 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 Maximum Weight Matching interface provided by PARDISO
//#define MWM_RCM_PQ
//#define MWM_MMD_PQ
//#define MWM_AMF_PQ
//#define MWM_AMD_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
                        3  bcg
                        8  gmres
                        9  fgmres */
    int SOLVER=8; /* gmres */

    CSRMAT fullmat, A, AT;
    FILE   *fp, *fo; 
    char   rhstyp[3], title[73], key[9], type[3], fname[100], foname[100], extension[3];
    int    i,j,k,l,fnamelen,n,m,nc,nz,nrhs,tmp0,tmp,tmp2,tmp3,ierr, nLU,
           nrhsix, *rhsptr, *rhsind, ELBOW, flags, nnzL,nnzU, nAx=0, 
           elbow, nrestart, max_it, nlev=0, mynrhs=2, nju,njlu,nalu,
           (*perm0)(),(*perm)(),(*permf)(), transpose=0, sumit, *ibuff, *ju, *jlu;
    REALS  DROP_TOL, condest, droptols[2], amgcancel, val,vb,
           CONDEST, restol;
    FLOAT  *rhs,*sol, *leftscale, *rightscale, *rhsval, *dbuff, *alu;
    float  systime, time_start, time_stop, secnds, sumtime;
    AMGLEVELMAT PRE, *next;
    ILUPACKPARAM param;
    size_t ndbuff, nibuff;

    
    /* 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 in a matrix in Harwell-Boeing format. By definition Harwell-Boeing
       format stores a matrix in compressed sparse COLUMN format. However,
       ILUPACK uses compressed sparse ROW format. Therefore it is necessary
       to transpose the matrix manually. 
                 /3.5 -1.0  0 \
       A matrix  | 0   2.0  0 | is stored as follows
                 \ 0    0  1.5/ 

       A.ia:   1  3  4  5           pointer to the start of every compressed 
                                    row plus pointer to the first space 
				    behind the compressed rows

       A.ja:   1    2    2    3     nonzero column indices

       A.a:   3.5 -1.0  2.0  1.5    nonzero numerical values

       The read part finally yields the following data structures
        -  A:  matrix in compressed sparse row format
	    o  A.nr, A.nc: number of rows and columns of A
            o  A.ia:  pointer array
            o  A.ja:  nonzero column index array 
	    o  A.a:   nonzero numerical values
	-  rhs:  right hand side(s) and additional data like exact solution
	         or initial guess
	-  n:  same as A.nr,A.nc
	-  nz:  number of nonzero entries
     */
#include "readmatrix.c"
    // if right hand sides are provided, then run AMGSOLVER for any of these
    // right hand sides. Otherwise use own set of right hand sides

    // allocate memory for the solution vector and some buffer
    sol  =(FLOAT *)MALLOC(mynrhs*(size_t)n*sizeof(FLOAT),  "main:sol");
    dbuff=(FLOAT *)MALLOC(3*(size_t)n*sizeof(FLOAT),"main:dbuff");
    ndbuff=3*(size_t)n;


    // set parameters to the default settings
    AMGINIT(A, &param);
    param.dbuff=dbuff;
    param.ndbuff=ndbuff;


    // select reordering functions
    perm0=PERMNULL;
    perm =PERMNULL;
    permf=PERMNULL;
#ifdef MINIMUM_DEGREE
    perm0=PERMMMD;
    perm =PERMMMD;
    permf=PERMMMD;
    fprintf(fo,"mmds/mmds|");
#elif defined REVERSE_CM
    perm0=PERMRCM;
    perm =PERMRCM;
    permf=PERMRCM;
    fprintf(fo,"rcms/rcms|");
#elif defined NESTED_DISSECTION
    perm0=PERMND;
    perm =PERMND;
    permf=PERMND;
    fprintf(fo,"nds /nds |");
#elif defined IND_SET
    perm0=PERMINDSET;
    perm =PERMINDSET;
    permf=PERMINDSET;
    fprintf(fo,"inds/inds|");
#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 |");
#elif defined MMD_PQ
    perm0=PERMMMD; 
    perm =PERMMMD; 
    permf=PERMPQ;
    fprintf(fo,"mmds/PQs |");
#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 |");
#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|");
#elif defined MC64_MMD_PQ
    perm0=PERMMC64MMD;
    perm =PERMMMD;
    permf=PERMPQ;
    fprintf(fo,"mc64md/PQ|");
#elif defined MC64_AMF_PQ
    perm0=PERMMC64AMF;
    perm =PERMAMF;
    permf=PERMPQ;
    fprintf(fo,"mc64af/PQ|");
#elif defined MC64_AMD_PQ
    perm0=PERMMC64AMD;
    perm =PERMAMD;
    permf=PERMPQ;
    fprintf(fo,"mc64ad/PQ|");
#elif defined MC64_METIS_E_PQ
    perm0=PERMMC64METISE;
    perm =PERMMETISE;
    permf=PERMPQ;
    fprintf(fo,"mc64me/PQ|");
#elif defined MC64_METIS_N_PQ
    perm0=PERMMC64METISN;
    perm =PERMMETISN;
    permf=PERMPQ;
    fprintf(fo,"mc64mn/PQ|");
#elif defined MWM_RCM_PQ
    perm0=PERMMWMRCM;
    perm =PERMRCM;
    permf=PERMPQ;
    fprintf(fo,"mwrc/PQs |");
#elif defined MWM_MMD_PQ
    perm0=PERMMWMMMD;
    perm =PERMMMD;
    permf=PERMPQ;
    fprintf(fo,"mwmd/PQs |");
#elif defined MWM_AMF_PQ
    perm0=PERMMWMAMF;
    perm =PERMAMF;
    permf=PERMPQ;
    fprintf(fo,"mwaf/PQs |");
#elif defined MWM_AMD_PQ
    perm0=PERMMWMAMD;
    perm =PERMAMD;
    permf=PERMPQ;
    fprintf(fo,"mwad/PQs |");
#elif defined MWM_METIS_E_PQ
    perm0=PERMMWMMETISE;
    perm =PERMMETISE;
    permf=PERMPQ;
    fprintf(fo,"mwme/PQs |");
#elif defined MWM_METIS_N_PQ
    perm0=PERMMWMMETISN;
    perm =PERMMETISN;
    permf=PERMPQ;
    fprintf(fo,"mwmn/PQs |");
#else
    fprintf(fo,"    /    |");
#endif



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

    // ONLY for mixed reordering strategies it is useful to
    // set the 'FINAL_PIVOTING' flag
#if !defined RCM_PQ && !defined AMF_PQ && !defined AMD_PQ && !defined MMD_PQ && !defined FC_PQ && !defined METIS_E_PQ && !defined METIS_N_PQ && !defined MWM_RCM_PQ && !defined MWM_MMD_PQ && !defined MWM_AMF_PQ && !defined MWM_AMD_PQ && !defined MWM_METIS_E_PQ && !defined MWM_METIS_N_PQ && !defined MC64_RCM_PQ && !defined MC64_MMD_PQ && !defined MC64_AMF_PQ && !defined MC64_AMD_PQ && !defined MC64_METIS_E_PQ && !defined MC64_METIS_N_PQ
    flags&=~FINAL_PIVOTING;
#endif
    // change flags if mpiluc is desired
#ifdef USE_MPILUC
    flags|=MULTI_PILUC;
#endif

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

    // choose the iterative solver, default: 8 (GMRES)
    param.ipar[5]=SOLVER;

    // 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=1.0e-8;

    // turn of column scaling
    //param.ipar[7]&=~(2+16+128);
    
    // use diagonal compensation
    //flags|=DIAGONAL_COMPENSATION;
    //flags|=DIAGONAL_COMPENSATION|TISMENETSKY_SC;
    //flags|=TISMENETSKY_SC;

    // use Standard Schur-complement
    // flags&=~SIMPLE_SC;

    // indicate that we want to compute more than one ILU
    flags|=RE_FACTOR;

    // keep the (1,2) block and the (1,2) block in the multilevel ILU
    // factorization, this may lead to more fill-in
    // flags&=~COARSE_REDUCE;

    // limit maximum number of entries in each column of L, row of U
    //param.ipar[3]=ELBOW*(nz/n)+5;
    // limit maximum number of entries in each row of S
    //param.ipar[9]=ELBOW*(nz/n)+5;

    // rewrite the updated parameters
    AMGSETPARAMS(A, &param, flags, elbow, droptols, condest,
		 restol, max_it, nrestart);

    // manually change the drop tolerance for the approximate Schur complement
    //param.fpar[8]=DROP_TOL;


    // print some messages that give information about flags and reorderings
#include "messages.c"

    evaluate_time(&time_start,&systime);
    ierr=AMGFACTOR(&A, &PRE, &nlev, &param, perm0,perm, permf);
    // update buffer size
    nibuff=param.nibuff;
    ndbuff=param.ndbuff;
    evaluate_time(&time_stop,&systime);
    secnds=time_stop-time_start;


    switch (ierr)
    {
           case  0: /* perfect! */
	            break;
           case -1: /* Error. input matrix may be wrong.
                       (The elimination process has generated a
			row in L or U whose length is .gt.  n.) */
	            printf("Error. input matrix may be wrong at level %d\n",nlev);
		    fprintf(fo,"input matrix may be wrong\n");fclose(fo); 
		    break;
           case -2: /* The matrix L overflows the array alu */
	            printf("The matrix L overflows the array alu at level %d\n",nlev);
		    fprintf(fo,"out of memory\n");fclose(fo); 
		    break;
           case -3: /* The matrix U overflows the array alu */
	            printf("The matrix U overflows the array alu at level %d\n",nlev);
		    fprintf(fo,"out of memory\n");fclose(fo); 
		    break;
           case -4: /* Illegal value for lfil */
	            printf("Illegal value for lfil at level %d\n",nlev);
		    fprintf(fo,"Illegal value for lfil\n");fclose(fo); 
		    break;
           case -5: /* zero row encountered */
	            printf("zero row encountered at level %d\n",nlev);
		    fprintf(fo,"zero row encountered\n");fclose(fo); 
		    break;
           case -6: /* zero column encountered */
	            printf("zero column encountered at level %d\n",nlev);
		    fprintf(fo,"zero column encountered\n");fclose(fo); 
		    break;
           case -7: /* buffers too small */
	            printf("buffers are too small\n");
		    // This error message would not be necessary if AMGsetup is
		    // called with the correct values of nibuff, ndbuff
		    printf("increase buffer size to at least %ld (float), %ld (int)\n",
			   ndbuff, nibuff);
		    fflush(stdout);
		    fprintf(fo,"buffers are too small\n");fclose(fo); 
           default: /* zero pivot encountered at step number ierr */
	            printf("zero pivot encountered at step number %d of level %d\n",ierr,nlev);
		    fprintf(fo,"zero pivot encountered\n");fclose(fo); 
		    break;
    } /* end switch */
    if (ierr) {
       fprintf(fo,"Iterative solver(s) cannot be applied\n");
       fflush(fo);
       exit(ierr);
    }

    // print some statistics about the levels, their size and the 
    // computation time
#include "printperformance.c"