/[escript]/trunk/paso/profiling/Test.c
ViewVC logotype

Diff of /trunk/paso/profiling/Test.c

Parent Directory Parent Directory | Revision Log Revision Log | View Patch Patch

revision 2381 by artak, Tue Feb 17 04:25:08 2009 UTC revision 2382 by artak, Tue Apr 14 03:48:22 2009 UTC
# Line 1  Line 1 
1  #include <stdio.h>  #include <stdio.h>
2    #include <unistd.h>
3  #include "paso/Common.h"  #include "paso/Common.h"
4  #include "paso/Solver.h"  #include "paso/Solver.h"
5  #include "paso/SystemMatrix.h"  #include "paso/SystemMatrix.h"
# Line 7  Line 8 
8    
9  #define PI (3.141592653589793)  #define PI (3.141592653589793)
10    
11    /*
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  double Lsup(double* x, int n) {  double Lsup(double* x, int n) {
26      double max=0;      double max=0;
27      int i;      int i;
# Line 21  double Lsup(double* x, int n) { Line 36  double Lsup(double* x, int n) {
36  int main (int argc, char *argv[]) {  int main (int argc, char *argv[]) {
37      Paso_SystemMatrix *A = NULL;      Paso_SystemMatrix *A = NULL;
38      double *b,*x,*x_ref;      double *b,*x,*x_ref;
39      dim_t i,n,level=1;      dim_t i,n,level=0;
40      double *error;      double *error;
41      double Lsuperror;      double Lsuperror;
42        
43        int c;
44        char *filename,*solver,*prec,*rhs;
45        extern char *optarg;
46        extern int optopt;
47    
48      Paso_Options options;      Paso_Options options;
49      Paso_Options_setDefaults(&options);      Paso_Options_setDefaults(&options);
50    
51          
52      if (argc<2) {      
53          fprintf(stderr,"Please enter the filename\n");      options.verbose=TRUE;
         return -1;      
     }  
54    
55      if (argc>2) {      while ((c = getopt(argc, argv, "s:p:f:r:l:c:")) != -1) {
56        level=atoi(argv[2]);        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      }      }
133            
     A=MEMALLOC(1,Paso_SystemMatrix);  
     
     A=Paso_SystemMatrix_loadMM_toCSR(argv[1]);  
   
134      if (A==NULL) {      if (A==NULL) {
135        fprintf(stderr,"CSR Matrix not Loaded\n");        fprintf(stderr,"CSR Matrix not Loaded\n");
136        return 0;        return 0;
137      }      }
138      n=Paso_SystemMatrix_getTotalNumRows(A);    
139      b=MEMALLOC(n,double);      
140      x=MEMALLOC(n,double);      if (level==0) {
141      x_ref=MEMALLOC(n,double);          Paso_solve(A,x,b,&options);
     error=MEMALLOC(n,double);  
       
     if(argc==4) {  
         Paso_RHS_loadMM_toCSR(argv[3],b,n);  
         Paso_test_run(A,b,level);  
142      }      }
143      else {      else if (level==3) {
144          for(i=0;i<n;i++) {          options.method=PASO_PCG;
145            x_ref[i]=cos(i*PI/n);          options.verbose=TRUE;
146          }          options.preconditioner=PASO_JACOBI;
         /*  raw scaled vector update operation: out = alpha * A * in + beta * out */  
         Paso_SystemMatrix_MatrixVector_CSR_OFFSET0(DBLE(1),A,x_ref,DBLE(0),b);  
147    
148          if (level==3) {          Paso_solve(A,x,b,&options);
             options.method=PASO_PCG;  
             options.verbose=FALSE;  
             options.preconditioner=PASO_JACOBI;  
               
   
       
             Paso_solve(A,x,b,&options);  
149    
150                        for(i=0;i<n;i++) {
151              for(i=0;i<n;i++) {            error[i]=x[i]-x_ref[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;  
       
             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);  
152          }          }
153          else {          Lsuperror=Lsup(error,n)/Lsup(x_ref,n);
154              Paso_test_run(A,b,level);                      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          }          }
167            Lsuperror=Lsup(error,n)/Lsup(x_ref,n);
168            fprintf(stdout,"Lsup error AMG %e\n",Lsuperror);
169        }
170        else {
171            Paso_test_run(A,b,level);            
172      }      }
173    
174       if (A!=NULL) {
       
175      MEMFREE(b);      MEMFREE(b);
176      MEMFREE(x);      MEMFREE(x);
177      MEMFREE(x_ref);      MEMFREE(x_ref);
178      MEMFREE(error);      MEMFREE(error);
179      Paso_SystemMatrix_free(A);      Paso_SystemMatrix_free(A);
180       }
181        
182  return 1;  return 1;
183  }  }
184    
Legend:
Removed from v.2381  
changed lines
  Added in v.2382
  ViewVC Help
Powered by ViewVC 1.1.26