maingnlsym.c

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

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

#include <blas.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
                        3  bcg
                        8  gmres
                        9  fgmres */
    int SOLVER=8; /* gmres */

    CSRMAT fullmat, A, AT, AS;
    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,flag, 
           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;


    MYGNLSYM(A,&AS,(int *)dbuff);
    
    /*
    for (i=0; i<AS.nr; i++) {
      printf("%3d\n",i+1); fflush(stdout);
      for (j=AS.ia[i]; j<AS.ia[i+1]; j++)
	printf("%8d",AS.ja[j-1]);
      printf("\n"); fflush(stdout);
      for (j=AS.ia[i]; j<AS.ia[i+1]; j++)
	printf("%8.1le",AS.a[j-1]);
      printf("\n"); fflush(stdout);
    }
    */

    // set parameters to the default settings
    MYSYMAMGINIT(AS, &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,"mwmd/PPs |");
#elif defined SMWM_AMF_PP
    perm0=PERMSMWMAMF;
    perm =PERMSMWMAMF;
    permf=SPDPERMPP;
    fprintf(fo,"mwaf/PPs |");
#elif defined SMWM_AMD_PP
    perm0=PERMSMWMAMD;
    perm =PERMSMWMAMD;
    permf=SPDPERMPP;
    fprintf(fo,"mwad/PPs |");
#elif defined SMWM_RCM_PP
    perm0=PERMSMWMRCM;
    perm =PERMSMWMRCM;
    permf=SPDPERMPP;
    fprintf(fo,"mwrc/PPs |");
#elif defined SMWM_METIS_E_PP
    perm0=PERMSMWMMETISE;
    perm =PERMSMWMMETISE;
    permf=SPDPERMPP;
    fprintf(fo,"mwme/PQs |");
#elif defined SMWM_METIS_N_PP
    perm0=PERMSMWMMETISN;
    perm =PERMSMWMMETISN;
    permf=SPDPERMPP;
    fprintf(fo,"mwmn/PPs |");
#elif defined SMC64_MMD_PP
    perm0=PERMSMC64MMD;
    perm =PERMSMC64MMD;
    permf=SPDPERMPP;
    fprintf(fo,"mc64md/PP|");
#elif defined SMC64_AMF_PP
    perm0=PERMSMC64AMF;
    perm =PERMSMC64AMF;
    permf=SPDPERMPP;
    fprintf(fo,"mc64af/PP|");
#elif defined SMC64_AMD_PP
    perm0=PERMSMC64AMD;
    perm =PERMSMC64AMD;
    permf=SPDPERMPP;
    fprintf(fo,"mc64ad/PP|");
#elif defined SMC64_RCM_PP
    perm0=PERMSMC64RCM;
    perm =PERMSMC64RCM;
    permf=SPDPERMPP;
    fprintf(fo,"mc64rc/PP|");
#elif defined SMC64_METIS_E_PP
    perm0=PERMSMC64METISE;
    perm =PERMSMC64METISE;
    permf=SPDPERMPP;
    fprintf(fo,"mc64me/PP|");
#elif defined SMC64_METIS_N_PP
    perm0=PERMSMC64METISN;
    perm =PERMSMC64METISN;
    permf=SPDPERMPP;
    fprintf(fo,"mc64mn/PP|");
#else
    fprintf(fo,"    /    |");
#endif



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

    // ONLY for mixed reordering strategies it is useful to
    // set the 'FINAL_PIVOTING' flag