paso/profiling/Test.c

#include <stdio.h>
#include <unistd.h>
#include "paso/Common.h"
#include "paso/Solver.h"
#include "paso/SystemMatrix.h"
#include "Paso_tests.h"
#include <math.h>

#define PI (3.141592653589793)

/*
 Usage: PasoTests -f filename [-s solver] [-p preconditioner] [-l level] [-r rhs matrix] [-c coupling parameter for AMG]
        filename - matrix to be loaded in CSR Matrix-Market format 
        solver   - PCG, GMRES, PRES20, TFQMR and MINRES
        preconditioner - ILU0, RILU, JACOBI, GS and AMG
        level    - options are 1,2 and 3
                   0 - default option just solves with default of specified parameters 
                   1 - test all solvers with default preconditioner
                   2 - test all preconditioners with default solver
                   3 - compare solution obtained by using AMG and Jacobi precondioners
        rhs matris - right hand side vector in CSR Matrix Market format.
        coupling parameter for AMG - this is the threshold value used in AMG in courenening process. Default is 0.05. 
*/

double Lsup(double* x, int n) {
    double max=0;
    int i;
    
    for (i=0;i<n;i++) {
        max=MAX(ABS(x[i]),max);
    }
    
    return max;
}

int main (int argc, char *argv[]) {
    Paso_SystemMatrix *A = NULL;
    double *b,*x,*x_ref;
    dim_t i,n,level=0;
    double *error;
    double Lsuperror;
   
    int c;
    char *filename,*solver,*prec,*rhs;
    extern char *optarg;
    extern int optopt;

    Paso_Options options;
    Paso_Options_setDefaults(&options);

    
    
    options.verbose=TRUE;

    while ((c = getopt(argc, argv, "s:p:f:r:l:c:")) != -1) {
      switch(c) {
        case 's':
            solver=optarg;
            if (strcmp(solver,"PCG")==0)
                options.method=PASO_PCG;
            else if (strcmp(solver,"GMRES")==0)
                options.method=PASO_GMRES;
            else if (strcmp(solver,"PRES20")==0)    
                options.method=PASO_PRES20;
            else if (strcmp(solver,"BICGSTAB")==0)
                options.method=PASO_BICGSTAB;
            else if (strcmp(solver,"TFQMR")==0)
                options.method=PASO_TFQMR;
            else if (strcmp(solver,"MINRES")==0)
                options.method=PASO_MINRES;
        break;
        case 'p':
            prec=optarg;
            if (strcmp(prec,"JACOBI")==0) 
                options.preconditioner=PASO_JACOBI;
            else if (strcmp(prec,"RILU")==0) 
                options.preconditioner=PASO_RILU;
            else if (strcmp(prec,"ILU0")==0)
                options.preconditioner=PASO_ILU0;
            else if (strcmp(prec,"GS")==0)
                options.preconditioner=PASO_GS;
            else if (strcmp(prec,"AMG")==0) {
                options.preconditioner=PASO_AMG;
            }
        break;
        case 'f':
            filename = optarg;
            A=MEMALLOC(1,Paso_SystemMatrix);
            A=Paso_SystemMatrix_loadMM_toCSR(filename);
            n=Paso_SystemMatrix_getTotalNumRows(A);
            b=MEMALLOC(n,double);
            x=MEMALLOC(n,double);
            x_ref=MEMALLOC(n,double);
            error=MEMALLOC(n,double);
            for(i=0;i<n;i++) {
             x_ref[i]=cos(i*PI/n);
            }
            Paso_SystemMatrix_MatrixVector_CSR_OFFSET0(DBLE(1),A,x_ref,DBLE(0),b);
            break;
        case 'r':
            rhs=optarg;
            if (A==NULL) {
             fprintf(stderr,"Left hand side not loaded yet.\n");
             break;
            }
            n=Paso_SystemMatrix_getTotalNumRows(A);
            Paso_RHS_loadMM_toCSR(rhs,b,n);
            break;
        case 'l':
            level=atoi(optarg);
            break;
        case 'c':
            options.couplingParam=atof(optarg);
        break;
        case '?':
            printf("unknown arg %c\n", optopt);
            break;
        case ':':
            printf("Usage: PasoTests -f filename [-s solver] [-p preconditioner] [-l level] [-r rhs matrix] [-c coupling parameter for AMG]\n");
            printf("\t filename - matrix to be loaded in CSR Matrix-Market format\n");
            printf("\t solver   - PCG, GMRES, PRES20, TFQMR and MINRES\n");
            printf("\t preconditioner - ILU0, RILU, JACOBI, GS and AMG\n");
            printf("\t level    - options are 1,2 and 3\n");
            printf("\t\t 0 - default option just solves with default of specified parameters\n");
            printf("\t\t 1 - test all solvers with default preconditioner\n");            
            printf("\t\t 2 - test all preconditioners with default solver\n");
            printf("\t\t 3 - compare solution obtained by using AMG and Jacobi precondioners\n");            
            printf("rhs matris - right hand side vector in CSR Matrix Market format.\n");
            printf("coupling parameter for AMG - this is the threshold value used in AMG in courenening process. Default is 0.05.\n");            
            break;
        }
    }
    
    if (A==NULL) {
      fprintf(stderr,"CSR Matrix not Loaded\n");
      return 0;
    }
   
    
    if (level==0) {
        Paso_solve(A,x,b,&options);
    }
    else if (level==3) {
        options.method=PASO_PCG;
        options.verbose=TRUE;
        options.preconditioner=PASO_JACOBI;

        Paso_solve(A,x,b,&options);

        for(i=0;i<n;i++) {
          error[i]=x[i]-x_ref[i];
        }
        Lsuperror=Lsup(error,n)/Lsup(x_ref,n);
        fprintf(stdout,"Lsup error Jacobi %e\n",Lsuperror);
        
        A->solver=NULL;
        options.method=PASO_PCG;
        options.preconditioner=PASO_AMG;
        options.couplingParam=0.25;
        options.AMGlevels=3;
        
        Paso_solve(A,x,b,&options);
        
        for(i=0;i<n;i++) {
          error[i]=x[i]-x_ref[i];
        }
        Lsuperror=Lsup(error,n)/Lsup(x_ref,n);
        fprintf(stdout,"Lsup error AMG %e\n",Lsuperror);
    }
    else {
        Paso_test_run(A,b,level);            
    }

   if (A!=NULL) {
    MEMFREE(b);
    MEMFREE(x);
    MEMFREE(x_ref);
    MEMFREE(error);
    Paso_SystemMatrix_free(A);
   }
    
return 1;
}


1	artak	2160	#include <stdio.h>
2	artak	2382	#include <unistd.h>
3	artak	2160	#include "paso/Common.h"
4			#include "paso/Solver.h"
5			#include "paso/SystemMatrix.h"
6	artak	2163	#include "Paso_tests.h"
7	artak	2275	#include <math.h>
8	artak	2160
9	artak	2275	#define PI (3.141592653589793)
10
11	artak	2382	/*
12			Usage: PasoTests -f filename [-s solver] [-p preconditioner] [-l level] [-r rhs matrix] [-c coupling parameter for AMG]
13			filename - matrix to be loaded in CSR Matrix-Market format
14			solver - PCG, GMRES, PRES20, TFQMR and MINRES
15			preconditioner - ILU0, RILU, JACOBI, GS and AMG
16			level - options are 1,2 and 3
17			0 - default option just solves with default of specified parameters
18			1 - test all solvers with default preconditioner
19			2 - test all preconditioners with default solver
20			3 - compare solution obtained by using AMG and Jacobi precondioners
21			rhs matris - right hand side vector in CSR Matrix Market format.
22			coupling parameter for AMG - this is the threshold value used in AMG in courenening process. Default is 0.05.
23			*/
24
25	artak	2275	double Lsup(double* x, int n) {
26			double max=0;
27			int i;
28
29			for (i=0;i<n;i++) {
30			max=MAX(ABS(x[i]),max);
31			}
32
33			return max;
34			}
35
36	artak	2160	int main (int argc, char *argv[]) {
37			Paso_SystemMatrix *A = NULL;
38	artak	2275	double b,x,*x_ref;
39	artak	2382	dim_t i,n,level=0;
40	artak	2275	double *error;
41			double Lsuperror;
42	artak	2160
43	artak	2382	int c;
44			char filename,solver,prec,rhs;
45			extern char *optarg;
46			extern int optopt;
47
48	artak	2275	Paso_Options options;
49			Paso_Options_setDefaults(&options);
50
51	artak	2382
52
53			options.verbose=TRUE;
54	artak	2165
55	artak	2382	while ((c = getopt(argc, argv, "s:p:f:r:l:c:")) != -1) {
56			switch(c) {
57			case 's':
58			solver=optarg;
59			if (strcmp(solver,"PCG")==0)
60			options.method=PASO_PCG;
61			else if (strcmp(solver,"GMRES")==0)
62			options.method=PASO_GMRES;
63			else if (strcmp(solver,"PRES20")==0)
64			options.method=PASO_PRES20;
65			else if (strcmp(solver,"BICGSTAB")==0)
66			options.method=PASO_BICGSTAB;
67			else if (strcmp(solver,"TFQMR")==0)
68			options.method=PASO_TFQMR;
69			else if (strcmp(solver,"MINRES")==0)
70			options.method=PASO_MINRES;
71			break;
72			case 'p':
73			prec=optarg;
74			if (strcmp(prec,"JACOBI")==0)
75			options.preconditioner=PASO_JACOBI;
76			else if (strcmp(prec,"RILU")==0)
77			options.preconditioner=PASO_RILU;
78			else if (strcmp(prec,"ILU0")==0)
79			options.preconditioner=PASO_ILU0;
80			else if (strcmp(prec,"GS")==0)
81			options.preconditioner=PASO_GS;
82			else if (strcmp(prec,"AMG")==0) {
83			options.preconditioner=PASO_AMG;
84			}
85			break;
86			case 'f':
87			filename = optarg;
88			A=MEMALLOC(1,Paso_SystemMatrix);
89			A=Paso_SystemMatrix_loadMM_toCSR(filename);
90			n=Paso_SystemMatrix_getTotalNumRows(A);
91			b=MEMALLOC(n,double);
92			x=MEMALLOC(n,double);
93			x_ref=MEMALLOC(n,double);
94			error=MEMALLOC(n,double);
95			for(i=0;i<n;i++) {
96			x_ref[i]=cos(i*PI/n);
97			}
98			Paso_SystemMatrix_MatrixVector_CSR_OFFSET0(DBLE(1),A,x_ref,DBLE(0),b);
99			break;
100			case 'r':
101			rhs=optarg;
102			if (A==NULL) {
103			fprintf(stderr,"Left hand side not loaded yet.\n");
104			break;
105			}
106			n=Paso_SystemMatrix_getTotalNumRows(A);
107			Paso_RHS_loadMM_toCSR(rhs,b,n);
108			break;
109			case 'l':
110			level=atoi(optarg);
111			break;
112			case 'c':
113			options.couplingParam=atof(optarg);
114			break;
115			case '?':
116			printf("unknown arg %c\n", optopt);
117			break;
118			case ':':
119			printf("Usage: PasoTests -f filename [-s solver] [-p preconditioner] [-l level] [-r rhs matrix] [-c coupling parameter for AMG]\n");
120			printf("\t filename - matrix to be loaded in CSR Matrix-Market format\n");
121			printf("\t solver - PCG, GMRES, PRES20, TFQMR and MINRES\n");
122			printf("\t preconditioner - ILU0, RILU, JACOBI, GS and AMG\n");
123			printf("\t level - options are 1,2 and 3\n");
124			printf("\t\t 0 - default option just solves with default of specified parameters\n");
125			printf("\t\t 1 - test all solvers with default preconditioner\n");
126			printf("\t\t 2 - test all preconditioners with default solver\n");
127			printf("\t\t 3 - compare solution obtained by using AMG and Jacobi precondioners\n");
128			printf("rhs matris - right hand side vector in CSR Matrix Market format.\n");
129			printf("coupling parameter for AMG - this is the threshold value used in AMG in courenening process. Default is 0.05.\n");
130			break;
131			}
132	artak	2165	}
133	artak	2160
134			if (A==NULL) {
135	artak	2165	fprintf(stderr,"CSR Matrix not Loaded\n");
136	artak	2160	return 0;
137			}
138	artak	2382
139	artak	2168
140	artak	2382	if (level==0) {
141			Paso_solve(A,x,b,&options);
142	artak	2160	}
143	artak	2382	else if (level==3) {
144			options.method=PASO_PCG;
145			options.verbose=TRUE;
146			options.preconditioner=PASO_JACOBI;
147	artak	2275
148	artak	2382	Paso_solve(A,x,b,&options);
149	artak	2275
150	artak	2382	for(i=0;i<n;i++) {
151			error[i]=x[i]-x_ref[i];
152	artak	2275	}
153	artak	2382	Lsuperror=Lsup(error,n)/Lsup(x_ref,n);
154			fprintf(stdout,"Lsup error Jacobi %e\n",Lsuperror);
155
156			A->solver=NULL;
157			options.method=PASO_PCG;
158			options.preconditioner=PASO_AMG;
159			options.couplingParam=0.25;
160			options.AMGlevels=3;
161
162			Paso_solve(A,x,b,&options);
163
164			for(i=0;i<n;i++) {
165			error[i]=x[i]-x_ref[i];
166	artak	2275	}
167	artak	2382	Lsuperror=Lsup(error,n)/Lsup(x_ref,n);
168			fprintf(stdout,"Lsup error AMG %e\n",Lsuperror);
169	artak	2168	}
170	artak	2382	else {
171			Paso_test_run(A,b,level);
172			}
173	artak	2275
174	artak	2382	if (A!=NULL) {
175	artak	2160	MEMFREE(b);
176			MEMFREE(x);
177	artak	2275	MEMFREE(x_ref);
178			MEMFREE(error);
179	artak	2160	Paso_SystemMatrix_free(A);
180	artak	2382	}
181
182	artak	2160	return 1;
183	artak	2275	}
184
185