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

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

Parent Directory Parent Directory | Revision Log Revision Log


Revision 3259 - (show annotations)
Mon Oct 11 01:48:14 2010 UTC (8 years, 11 months ago) by jfenwick
File MIME type: text/plain
File size: 14956 byte(s)
Merging dudley and scons updates from branches

1 /*******************************************************
2 *
3 * Copyright (c) 2003-2010 by University of Queensland
4 * Earth Systems Science Computational Center (ESSCC)
5 * http://www.uq.edu.au/esscc
6 *
7 * Primary Business: Queensland, Australia
8 * Licensed under the Open Software License version 3.0
9 * http://www.opensource.org/licenses/osl-3.0.php
10 *
11 *******************************************************/
12
13
14 /**************************************************************/
15
16 /* Paso: SystemMatrix: controls iterative linear system solvers */
17
18 /**************************************************************/
19
20 /* Copyrights by ACcESS Australia 2003/04 */
21 /* Author: Lutz Gross, l.gross@uq.edu.au */
22
23 /**************************************************************/
24
25 #include "Paso.h"
26 #include "SystemMatrix.h"
27 #include "Solver.h"
28 #include "esysUtils/blocktimer.h"
29
30 /***********************************************************************************/
31
32 /* free space */
33
34 void Paso_Solver_free(Paso_SystemMatrix* A) {
35 Paso_SystemMatrix_freePreconditioner(A);
36 }
37 /* call the iterative solver: */
38
39 void Paso_Solver(Paso_SystemMatrix* A,double* x,double* b,
40 Paso_Options* options,Paso_Performance* pp) {
41
42 double norm2_of_b,tol,tolerance,time_iter,net_time_start;
43 double *r=NULL,norm2_of_residual,last_norm2_of_residual,norm_max_of_b;
44 double norm2_of_b_local,norm_max_of_b_local,norm2_of_residual_local;
45 double norm_max_of_residual_local,norm_max_of_residual;
46 double last_norm_max_of_residual,*scaling;
47 #ifdef ESYS_MPI
48 double loc_norm;
49 #endif
50 dim_t i,totIter=0,cntIter,method;
51 bool_t finalizeIteration;
52 err_t errorCode=SOLVER_NO_ERROR;
53 dim_t numSol = Paso_SystemMatrix_getTotalNumCols(A);
54 dim_t numEqua = Paso_SystemMatrix_getTotalNumRows(A);
55 double blocktimer_precond, blocktimer_start = blocktimer_time();
56 double *x0=NULL;
57
58
59 Esys_resetError();
60 tolerance=options->tolerance;
61 if (tolerance < 100.* EPSILON) {
62 Esys_setError(VALUE_ERROR,"Paso_Solver: Tolerance is too small.");
63 }
64 if (tolerance >1.) {
65 Esys_setError(VALUE_ERROR,"Paso_Solver: Tolerance mut be less than one.");
66 }
67 method=Paso_Options_getSolver(options->method,PASO_PASO,options->symmetric,A->mpi_info);
68 /* check matrix type */
69 if ((A->type & MATRIX_FORMAT_CSC) || (A->type & MATRIX_FORMAT_OFFSET1) || (A->type & MATRIX_FORMAT_SYM) ) {
70 Esys_setError(TYPE_ERROR,"Paso_Solver: Iterative solver requires CSR format with unsymmetric storage scheme and index offset 0.");
71 }
72 if (A->col_block_size != A->row_block_size) {
73 Esys_setError(TYPE_ERROR,"Paso_Solver: Iterative solver requires row and column block sizes to be equal.");
74 }
75 if (Paso_SystemMatrix_getGlobalNumCols(A) != Paso_SystemMatrix_getGlobalNumRows(A)) {
76 Esys_setError(TYPE_ERROR,"Paso_Solver: Iterative solver requires a square matrix.");
77 return;
78 }
79 /*if (A->block_size != 1 && options->preconditioner==PASO_AMG) {
80 Esys_setError(TYPE_ERROR,"Paso_Solver: AMG on systems not supported yet.");
81 }
82 */
83 time_iter=Esys_timer();
84 /* this for testing only */
85 if (method==PASO_NONLINEAR_GMRES) {
86 Paso_Function* F=NULL;
87 F=Paso_Function_LinearSystem_alloc(A,b,options);
88 Paso_SystemMatrix_solvePreconditioner(A,x,b);
89 errorCode=Paso_Solver_NewtonGMRES(F,x,options,pp);
90 if (errorCode!=NO_ERROR) {
91 Esys_setError(SYSTEM_ERROR,"Paso_Solver_NewtonGMRES: an error has occured.");
92 }
93 Paso_Function_LinearSystem_free(F);
94 return;
95 }
96
97 /* ========================= */
98 Performance_startMonitor(pp,PERFORMANCE_ALL);
99 if (Esys_noError()) {
100 /* get normalization */
101 scaling=Paso_SystemMatrix_borrowNormalization(A);
102 if (Esys_noError()) {
103 /* get the norm of the right hand side */
104 norm2_of_b=0.;
105 norm_max_of_b=0.;
106 #pragma omp parallel private(norm2_of_b_local,norm_max_of_b_local)
107 {
108 norm2_of_b_local=0.;
109 norm_max_of_b_local=0.;
110 #pragma omp for private(i) schedule(static)
111 for (i = 0; i < numEqua ; ++i) {
112 norm2_of_b_local += b[i] * b[i];
113 norm_max_of_b_local = MAX(ABS(scaling[i]*b[i]),norm_max_of_b_local);
114 }
115 #pragma omp critical
116 {
117 norm2_of_b += norm2_of_b_local;
118 norm_max_of_b = MAX(norm_max_of_b_local,norm_max_of_b);
119 }
120 }
121 /* TODO: use one call */
122 #ifdef ESYS_MPI
123 {
124 loc_norm = norm2_of_b;
125 MPI_Allreduce(&loc_norm,&norm2_of_b, 1, MPI_DOUBLE, MPI_SUM, A->mpi_info->comm);
126 loc_norm = norm_max_of_b;
127 MPI_Allreduce(&loc_norm,&norm_max_of_b, 1, MPI_DOUBLE, MPI_MAX, A->mpi_info->comm);
128 }
129 #endif
130 norm2_of_b=sqrt(norm2_of_b);
131 /* if norm2_of_b==0 we are ready: x=0 */
132 if ( IS_NAN(norm2_of_b) || IS_NAN(norm_max_of_b) ) {
133 Esys_setError(VALUE_ERROR,"Paso_Solver: Matrix or right hand side contains undefined values.");
134 } else if (norm2_of_b <=0.) {
135
136 #pragma omp parallel for private(i) schedule(static)
137 for (i = 0; i < numSol; i++) x[i]=0.;
138 if (options->verbose) printf("right hand side is identical zero.\n");
139
140 } else {
141
142 if (options->verbose) {
143 printf("Paso_Solver: l2/lmax-norm of right hand side is %e/%e.\n",norm2_of_b,norm_max_of_b);
144 printf("Paso_Solver: l2/lmax-stopping criterion is %e/%e.\n",norm2_of_b*tolerance,norm_max_of_b*tolerance);
145 switch (method) {
146 case PASO_BICGSTAB:
147 printf("Paso_Solver: Iterative method is BiCGStab.\n");
148 break;
149 case PASO_PCG:
150 printf("Paso_Solver: Iterative method is PCG.\n");
151 break;
152 case PASO_TFQMR:
153 printf("Paso_Solver: Iterative method is TFQMR.\n");
154 break;
155 case PASO_MINRES:
156 printf("Paso_Solver: Iterative method is MINRES.\n");
157 break;
158 case PASO_PRES20:
159 printf("Paso_Solver: Iterative method is PRES20.\n");
160 break;
161 case PASO_GMRES:
162 if (options->restart>0) {
163 printf("Paso_Solver: Iterative method is GMRES(%d,%d)\n",options->truncation,options->restart);
164 } else {
165 printf("Paso_Solver: Iterative method is GMRES(%d)\n",options->truncation);
166 }
167 break;
168 }
169
170 }
171
172 /* construct the preconditioner */
173 blocktimer_precond = blocktimer_time();
174 Performance_startMonitor(pp,PERFORMANCE_PRECONDITIONER_INIT);
175 Paso_SystemMatrix_setPreconditioner(A,options);
176 Performance_stopMonitor(pp,PERFORMANCE_PRECONDITIONER_INIT);
177 blocktimer_increment("Paso_Solver_setPreconditioner()", blocktimer_precond);
178 options->set_up_time=Esys_timer()-time_iter;
179 if (! Esys_noError()) return;
180
181
182 /* get an initial guess by evaluating the preconditioner */
183 Paso_SystemMatrix_solvePreconditioner(A,x,b);
184
185 /* start the iteration process :*/
186 r=TMPMEMALLOC(numEqua,double);
187 x0=TMPMEMALLOC(numEqua,double);
188 Esys_checkPtr(r);
189 Esys_checkPtr(x0);
190 if (Esys_noError()) {
191
192 totIter = 1;
193 finalizeIteration = FALSE;
194 last_norm2_of_residual=norm2_of_b;
195 last_norm_max_of_residual=norm_max_of_b;
196 net_time_start=Esys_timer();
197
198 /* Loop */
199 while (! finalizeIteration) {
200 cntIter = options->iter_max - totIter;
201 finalizeIteration = TRUE;
202
203 /* Set initial residual. */
204 norm2_of_residual = 0;
205 norm_max_of_residual = 0;
206 #pragma omp parallel for private(i) schedule(static)
207 for (i = 0; i < numEqua; i++) r[i]=b[i];
208 Paso_SystemMatrix_MatrixVector_CSR_OFFSET0(DBLE(-1), A, x, DBLE(1), r);
209 #pragma omp parallel private(norm2_of_residual_local,norm_max_of_residual_local)
210 {
211 norm2_of_residual_local = 0;
212 norm_max_of_residual_local = 0;
213 #pragma omp for private(i) schedule(static)
214 for (i = 0; i < numEqua; i++) {
215 norm2_of_residual_local+= r[i] * r[i];
216 norm_max_of_residual_local=MAX(ABS(scaling[i]*r[i]),norm_max_of_residual_local);
217 }
218 #pragma omp critical
219 {
220 norm2_of_residual += norm2_of_residual_local;
221 norm_max_of_residual = MAX(norm_max_of_residual_local,norm_max_of_residual);
222 }
223 }
224 /* TODO: use one call */
225 #ifdef ESYS_MPI
226 loc_norm = norm2_of_residual;
227 MPI_Allreduce(&loc_norm,&norm2_of_residual, 1, MPI_DOUBLE, MPI_SUM, A->mpi_info->comm);
228 loc_norm = norm_max_of_residual;
229 MPI_Allreduce(&loc_norm,&norm_max_of_residual, 1, MPI_DOUBLE, MPI_MAX, A->mpi_info->comm);
230 #endif
231 norm2_of_residual =sqrt(norm2_of_residual);
232 options->residual_norm=norm2_of_residual;
233
234 if (options->verbose) printf("Paso_Solver: Step %5d: l2/lmax-norm of residual is %e/%e",totIter,norm2_of_residual,norm_max_of_residual);
235
236 if (totIter>1 && norm2_of_residual>=last_norm2_of_residual && norm_max_of_residual>=last_norm_max_of_residual) {
237
238 if (options->verbose) printf(" divergence!\n");
239 Esys_setError(DIVERGED, "Paso_Solver: No improvement during iteration. Iterative solver gives up.");
240
241 } else {
242
243 /* */
244 if (norm2_of_residual>tolerance*norm2_of_b || norm_max_of_residual>tolerance*norm_max_of_b ) {
245
246 tol=tolerance*MIN(norm2_of_b,0.1*norm2_of_residual/norm_max_of_residual*norm_max_of_b);
247 if (options->verbose) printf(" (new tolerance = %e).\n",tol);
248
249 last_norm2_of_residual=norm2_of_residual;
250 last_norm_max_of_residual=norm_max_of_residual;
251
252 /* call the solver */
253 switch (method) {
254 case PASO_BICGSTAB:
255 errorCode = Paso_Solver_BiCGStab(A, r, x, &cntIter, &tol, pp);
256 break;
257 case PASO_PCG:
258 errorCode = Paso_Solver_PCG(A, r, x, &cntIter, &tol, pp);
259 break;
260 case PASO_TFQMR:
261 tol=tolerance*norm2_of_residual/norm2_of_b;
262 errorCode = Paso_Solver_TFQMR(A, r, x0, &cntIter, &tol, pp);
263 #pragma omp for private(i) schedule(static)
264 for (i = 0; i < numEqua; i++) {
265 x[i]+= x0[i];
266 }
267 break;
268 case PASO_MINRES:
269 /* tol=tolerance*norm2_of_residual/norm2_of_b; */
270 errorCode = Paso_Solver_MINRES(A, r, x, &cntIter, &tol, pp);
271 break;
272 case PASO_PRES20:
273 errorCode = Paso_Solver_GMRES(A, r, x, &cntIter, &tol,5,20, pp);
274 break;
275 case PASO_GMRES:
276 errorCode = Paso_Solver_GMRES(A, r, x, &cntIter, &tol,options->truncation,options->restart, pp);
277 break;
278 }
279
280 totIter += cntIter;
281
282 /* error handling */
283 if (errorCode==SOLVER_NO_ERROR) {
284 finalizeIteration = FALSE;
285 } else if (errorCode==SOLVER_MAXITER_REACHED) {
286 Esys_setError(DIVERGED,"Paso_Solver: maximum number of iteration step reached.\nReturned solution does not fulfil stopping criterion.");
287 if (options->verbose) printf("Paso_Solver: Maximum number of iterations reached.!\n");
288 } else if (errorCode == SOLVER_INPUT_ERROR ) {
289 Esys_setError(SYSTEM_ERROR,"Paso_Solver: illegal dimension in iterative solver.");
290 if (options->verbose) printf("Paso_Solver: Internal error!\n");
291 } else if (errorCode == SOLVER_NEGATIVE_NORM_ERROR) {
292 Esys_setError(VALUE_ERROR,"Paso_Solver: negative energy norm (try other solver or preconditioner).");
293 if (options->verbose) printf("Paso_Solver: negative energy norm (try other solver or preconditioner)!\n");
294 } else if ( errorCode == SOLVER_BREAKDOWN ) {
295 if (cntIter <= 1) {
296 Esys_setError(ZERO_DIVISION_ERROR, "Paso_Solver: fatal break down in iterative solver.");
297 if (options->verbose) printf("Paso_Solver: Uncurable break down!\n");
298 } else {
299 if (options->verbose) printf("Paso_Solver: Breakdown at iter %d (residual = %e). Restarting ...\n", totIter, tol);
300 finalizeIteration = FALSE;
301 }
302 } else {
303 Esys_setError(SYSTEM_ERROR,"Paso_Solver:generic error in solver.");
304 if (options->verbose) printf("Paso_Solver: generic error in solver!\n");
305 }
306 } else {
307 if (options->verbose) printf(". convergence! \n");
308 options->converged=TRUE;
309 }
310 }
311 } /* while */
312 options->net_time=Esys_timer()-net_time_start;
313 }
314 MEMFREE(r);
315 MEMFREE(x0);
316 options->num_level=0;
317 options->num_inner_iter=0;
318 options->num_iter=totIter;
319 }
320 }
321 }
322 options->time=Esys_timer()-time_iter;
323 Performance_stopMonitor(pp,PERFORMANCE_ALL);
324 blocktimer_increment("Paso_Solver()", blocktimer_start);
325
326 }

Properties

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

  ViewVC Help
Powered by ViewVC 1.1.26