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

    CSRMAT fullmat, A, AT;
    FILE   *fp, *fo; 
    char   rhstyp[3], title[72], key[8], type[3], fname[100], foname[100];
    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,
           (*perm0)(),(*perm)(),(*permf)(), transpose=-1, sumit;
    REALS  DROP_TOL, condest, droptols[2], amgcancel, val,vb,
           CONDEST, restol;
    FLOAT  *rhs,*sol, *leftscale, *rightscale, *rhsval, *dbuff;
    float  systime, time_start, time_stop, secnds, sumtime;
    AMGLEVELMAT PRE, *next;
    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 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"
    AT=A;
    // 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(AT, &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;