/[escript]/trunk/paso/src/Solver.c
ViewVC logotype

Contents of /trunk/paso/src/Solver.c

Parent Directory Parent Directory | Revision Log Revision Log


Revision 584 - (show annotations)
Thu Mar 9 23:03:38 2006 UTC (13 years, 1 month ago) by gross
Original Path: trunk/paso/src/Solvers/Solver.c
File MIME type: text/plain
File size: 10961 byte(s)
eigenvalues: compiles and passes tests on altix now
1 /* $Id$ */
2
3 /**************************************************************/
4
5 /* Paso: SystemMatrix: controls iterative linear system solvers */
6
7 /**************************************************************/
8
9 /* Copyrights by ACcESS Australia 2003/04 */
10 /* Author: gross@access.edu.au */
11
12 /**************************************************************/
13
14 #include "Paso.h"
15 #include "SystemMatrix.h"
16 #include "Solver.h"
17
18 /***********************************************************************************/
19
20 /* free space */
21
22 void Paso_Solver_free(Paso_SystemMatrix* A) {
23 Paso_Preconditioner_free(A->solver);
24 A->solver=NULL;
25 }
26 /* call the iterative solver: */
27
28 void Paso_Solver(Paso_SystemMatrix* A,double* x,double* b,Paso_Options* options,Paso_Performance* pp) {
29 double norm2_of_b,tol,tolerance,time_iter,*r=NULL,norm2_of_residual,last_norm2_of_residual,norm_max_of_b,
30 norm2_of_b_local,norm_max_of_b_local,norm2_of_residual_local,norm_max_of_residual_local,norm_max_of_residual,
31 last_norm_max_of_residual,*scaling;
32 dim_t i,totIter,cntIter,method;
33 bool_t finalizeIteration;
34 err_t errorCode;
35 dim_t n_col = A->num_cols * A-> col_block_size;
36 dim_t n_row = A->num_rows * A-> row_block_size;
37
38 tolerance=MAX(options->tolerance,EPSILON);
39 Paso_resetError();
40 method=Paso_Options_getSolver(options->method,PASO_PASO,options->symmetric);
41 /* check matrix type */
42 if ((A->type & MATRIX_FORMAT_CSC) || (A->type & MATRIX_FORMAT_OFFSET1) || (A->type & MATRIX_FORMAT_SYM) ) {
43 Paso_setError(TYPE_ERROR,"Iterative solver requires CSR format with unsymmetric storage scheme and index offset 0.");
44 }
45 if (A->col_block_size != A->row_block_size) {
46 Paso_setError(TYPE_ERROR,"Iterative solver requires row and column block sizes to be equal.");
47 }
48 if (A->num_cols != A->num_rows) {
49 Paso_setError(TYPE_ERROR,"Iterative solver requires requires a square matrix.");
50 return;
51 }
52 Performance_startMonitor(pp,PERFORMANCE_ALL);
53 if (Paso_noError()) {
54 /* get normalization */
55 scaling=Paso_SystemMatrix_borrowNormalization(A);
56 if (Paso_noError()) {
57 /* get the norm of the right hand side */
58 norm2_of_b=0.;
59 norm_max_of_b=0.;
60 #pragma omp parallel private(norm2_of_b_local,norm_max_of_b_local)
61 {
62 norm2_of_b_local=0.;
63 norm_max_of_b_local=0.;
64 #pragma omp for private(i) schedule(static)
65 for (i = 0; i < n_row ; ++i) {
66 norm2_of_b_local += b[i] * b[i];
67 norm_max_of_b_local = MAX(ABS(scaling[i]*b[i]),norm_max_of_b_local);
68 }
69 #pragma omp critical
70 {
71 norm2_of_b += norm2_of_b_local;
72 norm_max_of_b = MAX(norm_max_of_b_local,norm_max_of_b);
73 }
74 }
75 norm2_of_b=sqrt(norm2_of_b);
76
77 /* if norm2_of_b==0 we are ready: x=0 */
78 if (norm2_of_b <=0.) {
79 #pragma omp parallel for private(i) schedule(static)
80 for (i = 0; i < n_col; i++) x[i]=0.;
81 if (options->verbose) printf("right hand side is identical zero.\n");
82 } else {
83 if (options->verbose) {
84 printf("l2/lmax-norm of right hand side is %e/%e.\n",norm2_of_b,norm_max_of_b);
85 printf("l2/lmax-stopping criterion is %e/%e.\n",norm2_of_b*tolerance,norm_max_of_b*tolerance);
86 switch (method) {
87 case PASO_BICGSTAB:
88 printf("Iterative method is BiCGStab.\n");
89 break;
90 case PASO_PCG:
91 printf("Iterative method is PCG.\n");
92 break;
93 case PASO_PRES20:
94 printf("Iterative method is PRES20.\n");
95 break;
96 case PASO_GMRES:
97 if (options->restart>0) {
98 printf("Iterative method is GMRES(%d,%d)\n",options->truncation,options->restart);
99 } else {
100 printf("Iterative method is GMRES(%d)\n",options->truncation);
101 }
102 break;
103 }
104 }
105 /* construct the preconditioner */
106
107 Performance_startMonitor(pp,PERFORMANCE_PRECONDITIONER_INIT);
108 Paso_Solver_setPreconditioner(A,options);
109 Performance_stopMonitor(pp,PERFORMANCE_PRECONDITIONER_INIT);
110 if (! Paso_noError()) return;
111
112 /* get an initial guess by evaluating the preconditioner */
113 time_iter=Paso_timer();
114 #pragma omp parallel
115 Paso_Solver_solvePreconditioner(A,x,b);
116 if (! Paso_noError()) return;
117 /* start the iteration process :*/
118 r=TMPMEMALLOC(n_row,double);
119 if (Paso_checkPtr(r)) return;
120 totIter = 0;
121 finalizeIteration = FALSE;
122 last_norm2_of_residual=norm2_of_b;
123 last_norm_max_of_residual=norm_max_of_b;
124 while (! finalizeIteration) {
125 cntIter = options->iter_max - totIter;
126 finalizeIteration = TRUE;
127 /* Set initial residual. */
128 norm2_of_residual = 0;
129 norm_max_of_residual = 0;
130 #pragma omp parallel private(norm2_of_residual_local,norm_max_of_residual_local)
131 {
132 #pragma omp for private(i) schedule(static)
133 for (i = 0; i < n_row; i++) r[i]=b[i];
134
135 Paso_SystemMatrix_MatrixVector_CSR_OFFSET0(DBLE(-1), A, x, DBLE(1), r);
136
137 norm2_of_residual_local = 0;
138 norm_max_of_residual_local = 0;
139 #pragma omp for private(i) schedule(static)
140 for (i = 0; i < n_row; i++) {
141 norm2_of_residual_local+= r[i] * r[i];
142 norm_max_of_residual_local=MAX(ABS(scaling[i]*r[i]),norm_max_of_residual_local);
143 }
144 #pragma omp critical
145 {
146 norm2_of_residual += norm2_of_residual_local;
147 norm_max_of_residual = MAX(norm_max_of_residual_local,norm_max_of_residual);
148 }
149 }
150 norm2_of_residual =sqrt(norm2_of_residual);
151
152 if (options->verbose) printf("Step %5d: l2/lmax-norm of residual is %e/%e",totIter,norm2_of_residual,norm_max_of_residual);
153 if (totIter>0 && norm2_of_residual>=last_norm2_of_residual && norm_max_of_residual>=last_norm_max_of_residual) {
154 if (options->verbose) printf(" divergence!\n");
155 Paso_setError(WARNING, "No improvement during iteration. Iterative solver gives up.");
156 } else {
157 /* */
158 if (norm2_of_residual>tolerance*norm2_of_b || norm_max_of_residual>tolerance*norm_max_of_b ) {
159 tol=tolerance*MIN(norm2_of_b,0.1*norm2_of_residual/norm_max_of_residual*norm_max_of_b);
160 if (options->verbose) printf(" (new tolerance = %e).\n",tol);
161 last_norm2_of_residual=norm2_of_residual;
162 last_norm_max_of_residual=norm_max_of_residual;
163 /* call the solver */
164 switch (method) {
165 case PASO_BICGSTAB:
166 errorCode = Paso_Solver_BiCGStab(A, r, x, &cntIter, &tol,pp);
167 break;
168 case PASO_PCG:
169 errorCode = Paso_Solver_PCG(A, r, x, &cntIter, &tol,pp);
170 break;
171 case PASO_PRES20:
172 errorCode = Paso_Solver_GMRES(A, r, x, &cntIter, &tol,5,20,pp);
173 break;
174 case PASO_GMRES:
175 errorCode = Paso_Solver_GMRES(A, r, x, &cntIter, &tol,options->truncation,options->restart,pp);
176 break;
177 }
178 totIter += cntIter;
179
180 /* error handling */
181 if (errorCode==NO_ERROR) {
182 finalizeIteration = FALSE;
183 } else if (errorCode==SOLVER_MAXITER_REACHED) {
184 Paso_setError(WARNING,"maximum number of iteration step reached.\nReturned solution does not fulfil stopping criterion.");
185 if (options->verbose) printf("Maximum number of iterations reached.!\n");
186 } else if (errorCode == SOLVER_INPUT_ERROR ) {
187 Paso_setError(SYSTEM_ERROR,"illegal dimension in iterative solver.");
188 if (options->verbose) printf("Internal error!\n");
189 } else if ( errorCode == SOLVER_BREAKDOWN ) {
190 if (cntIter <= 1) {
191 Paso_setError(ZERO_DIVISION_ERROR, "fatal break down in iterative solver.");
192 if (options->verbose) printf("Uncurable break down!\n");
193 } else {
194 if (options->verbose) printf("Breakdown at iter %d (residual = %e). Restarting ...\n", cntIter+totIter, tol);
195 finalizeIteration = FALSE;
196 }
197 } else if (errorCode == SOLVER_MEMORY_ERROR) {
198 Paso_setError(MEMORY_ERROR,"memory allocation failed.");
199 if (options->verbose) printf("Memory allocation failed!\n");
200 } else if (errorCode !=SOLVER_NO_ERROR ) {
201 Paso_setError(SYSTEM_ERROR,"unidentified error in iterative solver.");
202 if (options->verbose) printf("Unidentified error!\n");
203 }
204 } else {
205 if (options->verbose) printf(". convergence! \n");
206 }
207 }
208 } /* while */
209 MEMFREE(r);
210 time_iter=Paso_timer()-time_iter;
211 if (options->verbose) {
212 printf("timing: solver: %.4e sec\n",time_iter);
213 if (totIter>0) printf("timing: per iteration step: %.4e sec\n",time_iter/totIter);
214 }
215 }
216 }
217 }
218 Performance_stopMonitor(pp,PERFORMANCE_ALL);
219 }

Properties

Name Value
svn:eol-style native
svn:keywords Author Date Id Revision

  ViewVC Help
Powered by ViewVC 1.1.26