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

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

Parent Directory Parent Directory | Revision Log Revision Log


Revision 2987 - (show annotations)
Tue Mar 16 01:32:43 2010 UTC (9 years, 10 months ago) by gross
File MIME type: text/plain
File size: 14662 byte(s)
FCT solver rewritten
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_Preconditioner_free(A->solver);
36 A->solver=NULL;
37 }
38 /* call the iterative solver: */
39
40 void Paso_Solver(Paso_SystemMatrix* A,double* x,double* b,
41 Paso_Options* options,Paso_Performance* pp) {
42
43 double norm2_of_b,tol,tolerance,time_iter,net_time_start;
44 double *r=NULL,norm2_of_residual,last_norm2_of_residual,norm_max_of_b;
45 double norm2_of_b_local,norm_max_of_b_local,norm2_of_residual_local;
46 double norm_max_of_residual_local,norm_max_of_residual;
47 double last_norm_max_of_residual,*scaling;
48 #ifdef PASO_MPI
49 double loc_norm;
50 #endif
51 dim_t i,totIter=0,cntIter,method;
52 bool_t finalizeIteration;
53 err_t errorCode=SOLVER_NO_ERROR;
54 dim_t numSol = Paso_SystemMatrix_getTotalNumCols(A);
55 dim_t numEqua = Paso_SystemMatrix_getTotalNumRows(A);
56 double blocktimer_precond, blocktimer_start = blocktimer_time();
57 double *x0=NULL;
58
59
60 Paso_resetError();
61 tolerance=options->tolerance;
62 if (tolerance < 100.* EPSILON) {
63 Paso_setError(VALUE_ERROR,"Paso_Solver: Tolerance is too small.");
64 }
65 if (tolerance >1.) {
66 Paso_setError(VALUE_ERROR,"Paso_Solver: Tolerance mut be less than one.");
67 }
68 method=Paso_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) || (A->type & MATRIX_FORMAT_SYM) ) {
71 Paso_setError(TYPE_ERROR,"Paso_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 Paso_setError(TYPE_ERROR,"Paso_Solver: Iterative solver requires row and column block sizes to be equal.");
75 }
76 if (Paso_SystemMatrix_getGlobalNumCols(A) != Paso_SystemMatrix_getGlobalNumRows(A)) {
77 Paso_setError(TYPE_ERROR,"Paso_Solver: Iterative solver requires a square matrix.");
78 return;
79 }
80 /*if (A->block_size != 1 && options->preconditioner==PASO_AMG) {
81 Paso_setError(TYPE_ERROR,"Paso_Solver: AMG on systems not supported yet.");
82 }
83 */
84 time_iter=Paso_timer();
85 /* this for testing only */
86 if (method==PASO_NONLINEAR_GMRES) {
87 Paso_Function* F=NULL;
88 F=Paso_Function_LinearSystem_alloc(A,b,options);
89 Paso_Solver_solvePreconditioner(A,x,b);
90 errorCode=Paso_Solver_NewtonGMRES(F,x,options,pp);
91 if (errorCode!=NO_ERROR) {
92 Paso_setError(SYSTEM_ERROR,"Paso_Solver_NewtonGMRES: an error has occured.");
93 }
94 Paso_Function_LinearSystem_free(F);
95 return;
96 }
97
98 /* ========================= */
99 Performance_startMonitor(pp,PERFORMANCE_ALL);
100 if (Paso_noError()) {
101 /* get normalization */
102 scaling=Paso_SystemMatrix_borrowNormalization(A);
103 if (Paso_noError()) {
104 /* get the norm of the right hand side */
105 norm2_of_b=0.;
106 norm_max_of_b=0.;
107 #pragma omp parallel private(norm2_of_b_local,norm_max_of_b_local)
108 {
109 norm2_of_b_local=0.;
110 norm_max_of_b_local=0.;
111 #pragma omp for private(i) schedule(static)
112 for (i = 0; i < numEqua ; ++i) {
113 norm2_of_b_local += b[i] * b[i];
114 norm_max_of_b_local = MAX(ABS(scaling[i]*b[i]),norm_max_of_b_local);
115 }
116 #pragma omp critical
117 {
118 norm2_of_b += norm2_of_b_local;
119 norm_max_of_b = MAX(norm_max_of_b_local,norm_max_of_b);
120 }
121 }
122 /* TODO: use one call */
123 #ifdef PASO_MPI
124 {
125 loc_norm = norm2_of_b;
126 MPI_Allreduce(&loc_norm,&norm2_of_b, 1, MPI_DOUBLE, MPI_SUM, A->mpi_info->comm);
127 loc_norm = norm_max_of_b;
128 MPI_Allreduce(&loc_norm,&norm_max_of_b, 1, MPI_DOUBLE, MPI_MAX, A->mpi_info->comm);
129 }
130 #endif
131 norm2_of_b=sqrt(norm2_of_b);
132 /* if norm2_of_b==0 we are ready: x=0 */
133 if ( IS_NAN(norm2_of_b) || IS_NAN(norm_max_of_b) ) {
134 Paso_setError(VALUE_ERROR,"Paso_Solver: Matrix or right hand side contains undefined values.");
135 } else if (norm2_of_b <=0.) {
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 } else {
140 if (options->verbose) {
141 printf("Paso_Solver: l2/lmax-norm of right hand side is %e/%e.\n",norm2_of_b,norm_max_of_b);
142 printf("Paso_Solver: l2/lmax-stopping criterion is %e/%e.\n",norm2_of_b*tolerance,norm_max_of_b*tolerance);
143 switch (method) {
144 case PASO_BICGSTAB:
145 printf("Paso_Solver: Iterative method is BiCGStab.\n");
146 break;
147 case PASO_PCG:
148 printf("Paso_Solver: Iterative method is PCG.\n");
149 break;
150 case PASO_TFQMR:
151 printf("Paso_Solver: Iterative method is TFQMR.\n");
152 break;
153 case PASO_MINRES:
154 printf("Paso_Solver: Iterative method is MINRES.\n");
155 break;
156 case PASO_PRES20:
157 printf("Paso_Solver: Iterative method is PRES20.\n");
158 break;
159 case PASO_GMRES:
160 if (options->restart>0) {
161 printf("Paso_Solver: Iterative method is GMRES(%d,%d)\n",options->truncation,options->restart);
162 } else {
163 printf("Paso_Solver: Iterative method is GMRES(%d)\n",options->truncation);
164 }
165 break;
166 }
167 }
168 /* construct the preconditioner */
169 blocktimer_precond = blocktimer_time();
170 Performance_startMonitor(pp,PERFORMANCE_PRECONDITIONER_INIT);
171 Paso_Solver_setPreconditioner(A,options);
172 Performance_stopMonitor(pp,PERFORMANCE_PRECONDITIONER_INIT);
173 blocktimer_increment("Paso_Solver_setPreconditioner()", blocktimer_precond);
174 options->set_up_time=Paso_timer()-time_iter;
175 if (! Paso_noError()) return;
176 /* get an initial guess by evaluating the preconditioner */
177 Paso_Solver_solvePreconditioner(A,x,b);
178 /* start the iteration process :*/
179 r=TMPMEMALLOC(numEqua,double);
180 x0=TMPMEMALLOC(numEqua,double);
181 Paso_checkPtr(r);
182 if (Paso_noError()) {
183 totIter = 1;
184 finalizeIteration = FALSE;
185 last_norm2_of_residual=norm2_of_b;
186 last_norm_max_of_residual=norm_max_of_b;
187 net_time_start=Paso_timer();
188 /* Loop */
189 while (! finalizeIteration) {
190 cntIter = options->iter_max - totIter;
191 finalizeIteration = TRUE;
192 /* Set initial residual. */
193 norm2_of_residual = 0;
194 norm_max_of_residual = 0;
195 #pragma omp parallel for private(i) schedule(static)
196 for (i = 0; i < numEqua; i++) r[i]=b[i];
197 Paso_SystemMatrix_MatrixVector_CSR_OFFSET0(DBLE(-1), A, x, DBLE(1), r);
198 #pragma omp parallel private(norm2_of_residual_local,norm_max_of_residual_local)
199 {
200 norm2_of_residual_local = 0;
201 norm_max_of_residual_local = 0;
202 #pragma omp for private(i) schedule(static)
203 for (i = 0; i < numEqua; i++) {
204 norm2_of_residual_local+= r[i] * r[i];
205 norm_max_of_residual_local=MAX(ABS(scaling[i]*r[i]),norm_max_of_residual_local);
206 }
207 #pragma omp critical
208 {
209 norm2_of_residual += norm2_of_residual_local;
210 norm_max_of_residual = MAX(norm_max_of_residual_local,norm_max_of_residual);
211 }
212 }
213 /* TODO: use one call */
214 #ifdef PASO_MPI
215 loc_norm = norm2_of_residual;
216 MPI_Allreduce(&loc_norm,&norm2_of_residual, 1, MPI_DOUBLE, MPI_SUM, A->mpi_info->comm);
217 loc_norm = norm_max_of_residual;
218 MPI_Allreduce(&loc_norm,&norm_max_of_residual, 1, MPI_DOUBLE, MPI_MAX, A->mpi_info->comm);
219 #endif
220 norm2_of_residual =sqrt(norm2_of_residual);
221 options->residual_norm=norm2_of_residual;
222
223 if (options->verbose) printf("Paso_Solver: Step %5d: l2/lmax-norm of residual is %e/%e",totIter,norm2_of_residual,norm_max_of_residual);
224 if (totIter>1 && norm2_of_residual>=last_norm2_of_residual && norm_max_of_residual>=last_norm_max_of_residual) {
225 if (options->verbose) printf(" divergence!\n");
226 Paso_setError(DIVERGED, "Paso_Solver: No improvement during iteration. Iterative solver gives up.");
227 } else {
228 /* */
229 if (norm2_of_residual>tolerance*norm2_of_b || norm_max_of_residual>tolerance*norm_max_of_b ) {
230 tol=tolerance*MIN(norm2_of_b,0.1*norm2_of_residual/norm_max_of_residual*norm_max_of_b);
231 if (options->verbose) printf(" (new tolerance = %e).\n",tol);
232 last_norm2_of_residual=norm2_of_residual;
233 last_norm_max_of_residual=norm_max_of_residual;
234 /* call the solver */
235 switch (method) {
236 case PASO_BICGSTAB:
237 errorCode = Paso_Solver_BiCGStab(A, r, x, &cntIter, &tol, pp);
238 break;
239 case PASO_PCG:
240 errorCode = Paso_Solver_PCG(A, r, x, &cntIter, &tol, pp);
241 break;
242 case PASO_TFQMR:
243 tol=tolerance*norm2_of_residual/norm2_of_b;
244 errorCode = Paso_Solver_TFQMR(A, r, x0, &cntIter, &tol, pp);
245 #pragma omp for private(i) schedule(static)
246 for (i = 0; i < numEqua; i++) {
247 x[i]+= x0[i];
248 }
249 break;
250 case PASO_MINRES:
251 tol=tolerance*norm2_of_residual/norm2_of_b;
252 errorCode = Paso_Solver_MINRES(A, r, x0, &cntIter, &tol, pp);
253 #pragma omp for private(i) schedule(static)
254 for (i = 0; i < numEqua; i++) {
255 x[i]+= x0[i];
256 }
257 break;
258 case PASO_PRES20:
259 errorCode = Paso_Solver_GMRES(A, r, x, &cntIter, &tol,5,20, pp);
260 break;
261 case PASO_GMRES:
262 errorCode = Paso_Solver_GMRES(A, r, x, &cntIter, &tol,options->truncation,options->restart, pp);
263 break;
264 }
265 totIter += cntIter;
266 /* error handling */
267 if (errorCode==SOLVER_NO_ERROR) {
268 finalizeIteration = FALSE;
269 } else if (errorCode==SOLVER_MAXITER_REACHED) {
270 Paso_setError(DIVERGED,"Paso_Solver: maximum number of iteration step reached.\nReturned solution does not fulfil stopping criterion.");
271 if (options->verbose) printf("Paso_Solver: Maximum number of iterations reached.!\n");
272 } else if (errorCode == SOLVER_INPUT_ERROR ) {
273 Paso_setError(SYSTEM_ERROR,"Paso_Solver: illegal dimension in iterative solver.");
274 if (options->verbose) printf("Paso_Solver: Internal error!\n");
275 } else if ( errorCode == SOLVER_BREAKDOWN ) {
276 if (cntIter <= 1) {
277 Paso_setError(ZERO_DIVISION_ERROR, "Paso_Solver: fatal break down in iterative solver.");
278 if (options->verbose) printf("Paso_Solver: Uncurable break down!\n");
279 } else {
280 if (options->verbose) printf("Paso_Solver: Breakdown at iter %d (residual = %e). Restarting ...\n", cntIter+totIter, tol);
281 finalizeIteration = FALSE;
282 }
283 }
284 } else {
285 if (options->verbose) printf(". convergence! \n");
286 options->converged=TRUE;
287 }
288 }
289 } /* while */
290 options->net_time=Paso_timer()-net_time_start;
291 }
292 MEMFREE(r);
293 MEMFREE(x0);
294 options->num_level=0;
295 options->num_inner_iter=0;
296 options->num_iter=totIter;
297 }
298 }
299 }
300 options->time=Paso_timer()-time_iter;
301 Performance_stopMonitor(pp,PERFORMANCE_ALL);
302 blocktimer_increment("Paso_Solver()", blocktimer_start);
303
304 }

Properties

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

  ViewVC Help
Powered by ViewVC 1.1.26