/[escript]/branches/doubleplusgood/paso/src/Solver.cpp
ViewVC logotype

Contents of /branches/doubleplusgood/paso/src/Solver.cpp

Parent Directory Parent Directory | Revision Log Revision Log


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

Properties

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

  ViewVC Help
Powered by ViewVC 1.1.26