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


#include "symswitches.c"




#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. This is the adapted symmetric version
// of it
//#define PP_PERM


// mixed strategies that finally switch to PP 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

// alternative Symmetric Maximum Weight Matching interface provided by PARDISO
//#define SMC64_RCM_PP
//#define SMC64_MMD_PP
//#define SMC64_AMF_PP
//#define SMC64_AMD_PP
//#define SMC64_METIS_E_PP
//#define SMC64_METIS_N_PP

//#define SMWM_RCM_PP
//#define SMWM_MMD_PP
//#define SMWM_AMF_PP
//#define SMWM_AMD_PP
//#define SMWM_METIS_E_PP
#define SMWM_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  sbcg
                        3  bcg
                        4  sqmr
                        8  gmres
                        9  fgmres */
    int SOLVER=4; /* sqmr */

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

    FILE *fp, *fo; 
    char rhstyp[3], title[73], key[9], 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,mynrhs=2,tmp0,tmp,tmp2,tmp3,ierr,flag,
         *invperm, *buff, *ws, *symmmd,flags,elbow,max_it,ipar[20],
         nrhsix, *rhsptr, *rhsind, sumit;
    REALS eps, DROP_TOL, CONDEST,condest,droptols[2],restol, val,vb;
    FLOAT *rhs,*sol,*w, *scale, *rhsval, *dbuff;
    float  systime, time_start,   time_stop,   secnds, secndsprep,
           timeAx_start, timeAx_stop, secndsAx, sumtime;
    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 in a symmetric positive definite matrix in Harwell-Boeing format. 
       By definition Harwell-Boeing format stores a matrix in compressed sparse
       COLUMN format. In the symmetric case this is equivalent to compressed 
       sparse ROW format which ILUPACK uses. 
       Note that only the upper triangular part is stored by rows.
                                 / 3.5 -1.0  0 \
       A symm. pos. def. matrix  |-1.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 
                                    upper triangular 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 "spdreadmatrix.c"
    // make matrix indefinite
    for (i=0; i<n;i++)
        for (j=A.ia[i]-1;j<A.ia[i+1]-1; j++) 
	    if (A.ja[j]==i+1) {
#if defined _SINGLE_REAL_ || defined _DOUBLE_REAL_
	       A.a[j]=-0.1*A.a[j];
	       if (!(i%3))
		  A.a[j]=0.0;
#else
	       A.a[j].r=-0.1*A.a[j].r;
	       A.a[j].i=-A.a[j].i;
#endif    
	    }
    
    // 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), "mainspd:sol");
    dbuff=(FLOAT *) MALLOC(3*(size_t)n*sizeof(FLOAT),      "mainspd:dbuff");
    ndbuff=3*(size_t)n;


    // set parameters to the default settings
    MYSYMAMGINIT(A, &param);
    param.dbuff=dbuff;
    param.ndbuff=ndbuff;
    
    // select reordering functions
    perm0=SYMPERMNULL;
    perm =SYMPERMNULL;
    permf=SYMPERMNULL;
#ifdef MINIMUM_DEGREE
    perm0=SYMPERMMMD;
    perm =SYMPERMMMD;
    permf=SYMPERMMMD;
    fprintf(fo,"mmds/mmds|");
#elif defined REVERSE_CM
    perm0=SYMPERMRCM;
    perm =SYMPERMRCM;
    permf=SYMPERMRCM;
    fprintf(fo,"rcms/rcms|");
#elif defined NESTED_DISSECTION
    perm0=SYMPERMND;
    perm =SYMPERMND;
    permf=SYMPERMND;
    fprintf(fo,"nds /nds |");
#elif defined IND_SET
    perm0=SYMPERMINDSET;
    perm =SYMPERMINDSET;
    permf=SYMPERMINDSET;
    fprintf(fo,"inds/inds|");
#elif defined AMF
    perm0=SYMPERMAMF;
    perm =SYMPERMAMF;
    permf=SYMPERMAMF;
    fprintf(fo,"amfs/amfs|");
#elif defined AMD
    perm0=SYMPERMAMD;
    perm =SYMPERMAMD;
    permf=SYMPERMAMD;
    fprintf(fo,"amds/amds|");
#elif defined PP_PERM
    perm0=SPDPERMPP;
    perm =SPDPERMPP;
    permf=SPDPERMPP;
    fprintf(fo,"PPs /PPs |");
#elif defined MMD_PP
    perm0=SYMPERMMMD;
    perm =SYMPERMMMD;
    permf=SPDPERMPP;
    fprintf(fo,"mmds/PPs |");
#elif defined AMF_PP
    perm0=SYMPERMAMF;
    perm =SYMPERMAMF;
    permf=SPDPERMPP;
    fprintf(fo,"amfs/PPs |");
#elif defined AMD_PP
    perm0=SYMPERMAMD;
    perm =SYMPERMAMD;
    permf=SPDPERMPP;
    fprintf(fo,"amds/PPs |");
#elif defined RCM_PP
    perm0=SYMPERMRCM;
    perm =SYMPERMRCM;
    permf=SPDPERMPP;
    fprintf(fo,"rcms/PPs |");
#elif defined FC_PP
    perm0=SYMPERMFC;
    perm =SYMPERMFC;
    permf=SPDPERMPP;
    fprintf(fo,"FCs /PPs |");
#elif defined METIS_E_PP
    perm0=SYMPERMMETISE;
    perm =SYMPERMMETISE;
    permf=SPDPERMPP;
    fprintf(fo,"mes /PQs |");
#elif defined METIS_N_PP
    perm0=SYMPERMMETISN;
    perm =SYMPERMMETISN;
    permf=SPDPERMPP;
    fprintf(fo,"mns /PPs |");
#elif defined SMWM_MMD_PP
    perm0=PERMSMWMMMD;
    perm =PERMSMWMMMD;
    permf=SPDPERMPP;
    fprintf(fo,"mmds/PPs |");
#elif defined SMWM_AMF_PP
    perm0=PERMSMWMAMF;
    perm =PERMSMWMAMF;
    permf=SPDPERMPP;
    fprintf(fo,"amfs/PPs |");
#elif defined SMWM_AMD_PP
    perm0=PERMSMWMAMD;
    perm =PERMSMWMAMD;
    permf=SPDPERMPP;
    fprintf(fo,"amds/PPs |");
#elif defined SMWM_RCM_PP
    perm0=PERMSMWMRCM;
    perm =PERMSMWMRCM;
    permf=SPDPERMPP;
    fprintf(fo,"rcms/PPs |");
#elif defined SMWM_METIS_E_PP
    perm0=PERMSMWMMETISE;
    perm =PERMSMWMMETISE;
    permf=SPDPERMPP;
    fprintf(fo,"mes /PQs |");
#elif defined SMWM_METIS_N_PP
    perm0=PERMSMWMMETISN;
    perm =PERMSMWMMETISN;
    permf=SPDPERMPP;
    fprintf(fo,"mns /PPs |");
#elif defined SMC64_MMD_PP
    perm0=PERMSMC64MMD;
    perm =PERMSMC64MMD;
    permf=SPDPERMPP;
    fprintf(fo,"mmds/PPs |");
#elif defined SMC64_AMF_PP
    perm0=PERMSMC64AMF;
    perm =PERMSMC64AMF;
    permf=SPDPERMPP;
    fprintf(fo,"amfs/PPs |");
#elif defined SMC64_AMD_PP
    perm0=PERMSMC64AMD;
    perm =PERMSMC64AMD;
    permf=SPDPERMPP;
    fprintf(fo,"amds/PPs |");
#elif defined SMC64_RCM_PP
    perm0=PERMSMC64RCM;
    perm =PERMSMC64RCM;
    permf=SPDPERMPP;
    fprintf(fo,"rcms/PPs |");
#elif defined SMC64_METIS_E_PP