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

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

Parent Directory Parent Directory | Revision Log Revision Log


Revision 2275 - (show annotations)
Tue Feb 17 04:25:08 2009 UTC (12 years, 9 months ago) by artak
File MIME type: text/plain
File size: 2415 byte(s)
Extra tests are added for checking the quality of the solution from Jacobi and AMG.
1 #include <stdio.h>
2 #include "paso/Common.h"
3 #include "paso/Solver.h"
4 #include "paso/SystemMatrix.h"
5 #include "Paso_tests.h"
6 #include <math.h>
7
8 #define PI (3.141592653589793)
9
10 double Lsup(double* x, int n) {
11 double max=0;
12 int i;
13
14 for (i=0;i<n;i++) {
15 max=MAX(ABS(x[i]),max);
16 }
17
18 return max;
19 }
20
21 int main (int argc, char *argv[]) {
22 Paso_SystemMatrix *A = NULL;
23 double *b,*x,*x_ref;
24 dim_t i,n,level=1;
25 double *error;
26 double Lsuperror;
27
28 Paso_Options options;
29 Paso_Options_setDefaults(&options);
30
31
32 if (argc<2) {
33 fprintf(stderr,"Please enter the filename\n");
34 return -1;
35 }
36
37 if (argc>2) {
38 level=atoi(argv[2]);
39 }
40
41 A=MEMALLOC(1,Paso_SystemMatrix);
42
43 A=Paso_SystemMatrix_loadMM_toCSR(argv[1]);
44
45 if (A==NULL) {
46 fprintf(stderr,"CSR Matrix not Loaded\n");
47 return 0;
48 }
49 n=Paso_SystemMatrix_getTotalNumRows(A);
50 b=MEMALLOC(n,double);
51 x=MEMALLOC(n,double);
52 x_ref=MEMALLOC(n,double);
53 error=MEMALLOC(n,double);
54
55 if(argc==4) {
56 Paso_RHS_loadMM_toCSR(argv[3],b,n);
57 Paso_test_run(A,b,level);
58 }
59 else {
60 for(i=0;i<n;i++) {
61 x_ref[i]=cos(i*PI/n);
62 }
63 /* raw scaled vector update operation: out = alpha * A * in + beta * out */
64 Paso_SystemMatrix_MatrixVector_CSR_OFFSET0(DBLE(1),A,x_ref,DBLE(0),b);
65
66 if (level==3) {
67 options.method=PASO_PCG;
68 options.verbose=FALSE;
69 options.preconditioner=PASO_JACOBI;
70
71
72
73 Paso_solve(A,x,b,&options);
74
75
76 for(i=0;i<n;i++) {
77 error[i]=x[i]-x_ref[i];
78 }
79 Lsuperror=Lsup(error,n)/Lsup(x_ref,n);
80 fprintf(stdout,"Lsup error Jacobi %e\n",Lsuperror);
81
82 A->solver=NULL;
83 options.method=PASO_PCG;
84 options.preconditioner=PASO_AMG;
85
86 Paso_solve(A,x,b,&options);
87
88 for(i=0;i<n;i++) {
89 error[i]=x[i]-x_ref[i];
90 }
91 Lsuperror=Lsup(error,n)/Lsup(x_ref,n);
92 fprintf(stdout,"Lsup error AMG %e\n",Lsuperror);
93 }
94 else {
95 Paso_test_run(A,b,level);
96 }
97
98 }
99
100
101
102 MEMFREE(b);
103 MEMFREE(x);
104 MEMFREE(x_ref);
105 MEMFREE(error);
106 Paso_SystemMatrix_free(A);
107 return 1;
108 }
109
110

  ViewVC Help
Powered by ViewVC 1.1.26