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

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

Parent Directory Parent Directory | Revision Log Revision Log


Revision 5265 - (show annotations)
Mon Nov 17 04:59:25 2014 UTC (5 years, 4 months ago) by jfenwick
File size: 15787 byte(s)
Adding missing include which clang noticed.
1 /*****************************************************************************
2 *
3 * Copyright (c) 2003-2014 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 2012-2013 by School of Earth Sciences
12 * Development from 2014 by Centre for Geoscience Computing (GeoComp)
13 *
14 *****************************************************************************/
15
16
17 /****************************************************************************/
18
19 /* Paso: SystemMatrix: controls iterative linear system solvers */
20
21 /****************************************************************************/
22
23 /* Copyrights by ACcESS Australia 2003/04 */
24 /* Author: Lutz Gross, l.gross@uq.edu.au */
25
26 /****************************************************************************/
27 #include <iostream>
28 #include "Paso.h"
29 #include "SystemMatrix.h"
30 #include "Solver.h"
31 #include "esysUtils/blocktimer.h"
32
33 namespace paso {
34
35 void Solver_free(SystemMatrix* A)
36 {
37 A->freePreconditioner();
38 }
39
40 /// calls the iterative solver
41 void Solver(SystemMatrix_ptr A, double* x, double* b, Options* options,
42 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 finalizeIteration;
54 err_t errorCode=SOLVER_NO_ERROR;
55 const dim_t numSol = A->getTotalNumCols();
56 const dim_t numEqua = A->getTotalNumRows();
57 double blocktimer_precond, blocktimer_start = blocktimer_time();
58 double *x0=NULL;
59
60 Esys_resetError();
61 tolerance=options->tolerance;
62 if (tolerance < 100.* EPSILON) {
63 Esys_setError(VALUE_ERROR,"Solver: Tolerance is too small.");
64 }
65 if (tolerance >1.) {
66 Esys_setError(VALUE_ERROR,"Solver: Tolerance must be less than one.");
67 }
68 method=Options::getSolver(options->method, PASO_PASO, options->symmetric, A->mpi_info);
69 /* check matrix type */
70 if ((A->type & MATRIX_FORMAT_CSC) || (A->type & MATRIX_FORMAT_OFFSET1) ) {
71 Esys_setError(TYPE_ERROR,"Solver: Iterative solver requires CSR format with unsymmetric storage scheme and index offset 0.");
72 }
73 if (A->col_block_size != A->row_block_size) {
74 Esys_setError(TYPE_ERROR,"Solver: Iterative solver requires row and column block sizes to be equal.");
75 }
76 if (A->getGlobalNumCols() != A->getGlobalNumRows()) {
77 Esys_setError(TYPE_ERROR,"Solver: Iterative solver requires a square matrix.");
78 return;
79 }
80 time_iter=Esys_timer();
81 /* this for testing only */
82 if (method==PASO_NONLINEAR_GMRES) {
83 LinearSystem* F = new LinearSystem(A, b, options);
84 A->solvePreconditioner(x, b);
85 errorCode = Solver_NewtonGMRES(F, x, options, pp);
86 if (errorCode != NO_ERROR) {
87 Esys_setError(SYSTEM_ERROR,"Solver_NewtonGMRES: an error has occurred.");
88 }
89 delete F;
90 return;
91 }
92
93 r = new double[numEqua];
94 x0 = new double[numEqua];
95 A->balance();
96 options->num_level=0;
97 options->num_inner_iter=0;
98
99 /* ========================= */
100 if (Esys_noError()) {
101 Performance_startMonitor(pp, PERFORMANCE_ALL);
102 A->applyBalance(r, b, true);
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 += r[i] * r[i];
113 norm_max_of_b_local = MAX(ABS(r[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 #ifdef ESYS_MPI
122 /* TODO: use one call */
123 loc_norm = norm2_of_b;
124 MPI_Allreduce(&loc_norm,&norm2_of_b, 1, MPI_DOUBLE, MPI_SUM, A->mpi_info->comm);
125 loc_norm = norm_max_of_b;
126 MPI_Allreduce(&loc_norm,&norm_max_of_b, 1, MPI_DOUBLE, MPI_MAX, A->mpi_info->comm);
127 #endif
128 norm2_of_b=sqrt(norm2_of_b);
129 /* if norm2_of_b==0 we are ready: x=0 */
130 if (IS_NAN(norm2_of_b) || IS_NAN(norm_max_of_b)) {
131 Esys_setError(VALUE_ERROR, "Solver: Matrix or right hand side contains undefined values.");
132 } else if (norm2_of_b <= 0.) {
133 #pragma omp parallel for private(i) schedule(static)
134 for (i = 0; i < numSol; i++) x[i]=0.;
135 if (options->verbose)
136 std::cout << "right hand side is identical to zero." << std::endl;
137 } else {
138 if (options->verbose) {
139 std::cout << "Solver: l2/lmax-norm of right hand side is "
140 << norm2_of_b << "/" << norm_max_of_b << "." << std::endl
141 << "Solver: l2/lmax-stopping criterion is "
142 << norm2_of_b*tolerance << "/" << norm_max_of_b*tolerance
143 << "." << std::endl;
144 switch (method) {
145 case PASO_BICGSTAB:
146 std::cout << "Solver: Iterative method is BiCGStab.\n";
147 break;
148 case PASO_PCG:
149 std::cout << "Solver: Iterative method is PCG.\n";
150 break;
151 case PASO_TFQMR:
152 std::cout << "Solver: Iterative method is TFQMR.\n";
153 break;
154 case PASO_MINRES:
155 std::cout << "Solver: Iterative method is MINRES.\n";
156 break;
157 case PASO_PRES20:
158 std::cout << "Solver: Iterative method is PRES20.\n";
159 break;
160 case PASO_GMRES:
161 if (options->restart > 0) {
162 std::cout << "Solver: Iterative method is GMRES("
163 << options->truncation << ","
164 << options->restart << ")." << std::endl;
165 } else {
166 std::cout << "Solver: Iterative method is GMRES("
167 << options->truncation << ")." << std::endl;
168 }
169 break;
170 }
171 }
172
173 // construct the preconditioner
174 blocktimer_precond = blocktimer_time();
175 Performance_startMonitor(pp, PERFORMANCE_PRECONDITIONER_INIT);
176 A->setPreconditioner(options);
177 Performance_stopMonitor(pp, PERFORMANCE_PRECONDITIONER_INIT);
178 blocktimer_increment("Solver_setPreconditioner()", blocktimer_precond);
179 options->set_up_time=Esys_timer()-time_iter;
180 if (Esys_noError()) {
181 // get an initial guess by evaluating the preconditioner
182 A->solvePreconditioner(x, r);
183
184 totIter = 1;
185 finalizeIteration = false;
186 last_norm2_of_residual=norm2_of_b;
187 last_norm_max_of_residual=norm_max_of_b;
188 net_time_start=Esys_timer();
189
190 // main loop
191 while (!finalizeIteration) {
192 cntIter = options->iter_max - totIter;
193 finalizeIteration = true;
194
195 // Set initial residual
196 if (totIter > 1) {
197 // in the first iteration r = balance * b already
198 A->applyBalance(r, b, true);
199 }
200
201 SystemMatrix_MatrixVector_CSR_OFFSET0(-1., A, x, 1., r);
202 norm2_of_residual = 0;
203 norm_max_of_residual = 0;
204 #pragma omp parallel private(norm2_of_residual_local,norm_max_of_residual_local)
205 {
206 norm2_of_residual_local = 0;
207 norm_max_of_residual_local = 0;
208 #pragma omp for private(i) schedule(static)
209 for (i = 0; i < numEqua; i++) {
210 norm2_of_residual_local+= r[i] * r[i];
211 norm_max_of_residual_local=MAX(ABS(r[i]),norm_max_of_residual_local);
212 }
213 #pragma omp critical
214 {
215 norm2_of_residual += norm2_of_residual_local;
216 norm_max_of_residual = MAX(norm_max_of_residual_local,norm_max_of_residual);
217 }
218 }
219 #ifdef ESYS_MPI
220 // TODO: use one call
221 loc_norm = norm2_of_residual;
222 MPI_Allreduce(&loc_norm,&norm2_of_residual, 1, MPI_DOUBLE, MPI_SUM, A->mpi_info->comm);
223 loc_norm = norm_max_of_residual;
224 MPI_Allreduce(&loc_norm,&norm_max_of_residual, 1, MPI_DOUBLE, MPI_MAX, A->mpi_info->comm);
225 #endif
226 norm2_of_residual =sqrt(norm2_of_residual);
227 options->residual_norm=norm2_of_residual;
228
229 if (options->verbose)
230 std::cout << "Solver: Step " << totIter
231 << ": l2/lmax-norm of residual is "
232 << norm2_of_residual << "/" << norm_max_of_residual;
233
234 if (totIter > 1 &&
235 norm2_of_residual >= last_norm2_of_residual &&
236 norm_max_of_residual >= last_norm_max_of_residual) {
237
238 if (options->verbose) std::cout << " divergence!\n";
239 Esys_setError(DIVERGED, "Solver: No improvement during iteration. Iterative solver gives up.");
240
241 } else {
242 if (norm2_of_residual>tolerance*norm2_of_b ||
243 norm_max_of_residual>tolerance*norm_max_of_b ) {
244
245 tol=tolerance*MIN(norm2_of_b,0.1*norm2_of_residual/norm_max_of_residual*norm_max_of_b);
246 if (options->verbose)
247 std::cout << " (new tolerance = " << tol << ").\n";
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 = Solver_BiCGStab(A, r, x, &cntIter, &tol, pp);
256 break;
257 case PASO_PCG:
258 errorCode = 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 = 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 = Solver_MINRES(A, r, x, &cntIter, &tol, pp);
271 break;
272 case PASO_PRES20:
273 errorCode = Solver_GMRES(A, r, x, &cntIter, &tol, 5, 20, pp);
274 break;
275 case PASO_GMRES:
276 errorCode = 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, "Solver: maximum number of iteration steps reached.\nReturned solution does not fulfil stopping criterion.");
287 if (options->verbose)
288 std::cout << "Solver: Maximum number of "
289 "iterations reached." << std::endl;
290 } else if (errorCode == SOLVER_INPUT_ERROR) {
291 Esys_setError(SYSTEM_ERROR, "Solver: illegal dimension in iterative solver.");
292 if (options->verbose)
293 std::cout << "Solver: Internal error!\n";
294 } else if (errorCode == SOLVER_NEGATIVE_NORM_ERROR) {
295 Esys_setError(VALUE_ERROR, "Solver: negative energy norm (try other solver or preconditioner).");
296 if (options->verbose)
297 std::cout << "Solver: negative energy norm"
298 " (try other solver or preconditioner)!\n";
299 } else if (errorCode == SOLVER_BREAKDOWN) {
300 if (cntIter <= 1) {
301 Esys_setError(ZERO_DIVISION_ERROR, "Solver: fatal break down in iterative solver.");
302 if (options->verbose)
303 std::cout << "Solver: Uncurable break "
304 "down!" << std::endl;
305 } else {
306 if (options->verbose)
307 std::cout << "Solver: Breakdown at iter "
308 << totIter << " (residual = "
309 << tol << "). Restarting ...\n";
310 finalizeIteration = false;
311 }
312 } else {
313 Esys_setError(SYSTEM_ERROR, "Solver: Generic error in solver.");
314 if (options->verbose)
315 std::cout << "Solver: Generic error in solver!\n";
316 }
317 } else {
318 if (options->verbose)
319 std::cout << " convergence!" << std::endl;
320 options->converged = true;
321 }
322 }
323 } // while
324 options->net_time = Esys_timer()-net_time_start;
325 }
326 options->num_iter = totIter;
327 A->applyBalanceInPlace(x, false);
328 }
329 }
330 delete[] r;
331 delete[] x0;
332 options->time = Esys_timer()-time_iter;
333 Performance_stopMonitor(pp, PERFORMANCE_ALL);
334 blocktimer_increment("Solver()", blocktimer_start);
335 }
336
337 } // namespace paso
338

Properties

Name Value
svn:eol-style native
svn:keywords Author Date Id Revision
svn:mergeinfo /branches/amg_from_3530/paso/src/Solver.cpp:3531-3826 /branches/lapack2681/paso/src/Solver.cpp:2682-2741 /branches/pasowrap/paso/src/Solver.cpp:3661-3674 /branches/py3_attempt2/paso/src/Solver.cpp:3871-3891 /branches/restext/paso/src/Solver.cpp:2610-2624 /branches/ripleygmg_from_3668/paso/src/Solver.cpp:3669-3791 /branches/stage3.0/paso/src/Solver.cpp:2569-2590 /branches/symbolic_from_3470/paso/src/Solver.cpp:3471-3974 /branches/symbolic_from_3470/ripley/test/python/paso/src/Solver.cpp:3517-3974 /release/3.0/paso/src/Solver.cpp:2591-2601 /trunk/paso/src/Solver.cpp:4257-4344 /trunk/ripley/test/python/paso/src/Solver.cpp:3480-3515

  ViewVC Help
Powered by ViewVC 1.1.26