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

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

Parent Directory Parent Directory | Revision Log Revision Log


Revision 2315 - (show annotations)
Wed Mar 18 00:38:48 2009 UTC (10 years, 8 months ago) by gross
File MIME type: text/plain
File size: 13927 byte(s)
fixes for MPI
1
2 /*******************************************************
3 *
4 * Copyright (c) 2003-2008 by University of Queensland
5 * Earth Systems Science Computational Center (ESSCC)
6 * http://www.uq.edu.au/esscc
7 *
8 * Primary Business: Queensland, Australia
9 * Licensed under the Open Software License version 3.0
10 * http://www.opensource.org/licenses/osl-3.0.php
11 *
12 *******************************************************/
13
14
15 /**************************************************************/
16
17 /* Paso: SystemMatrix: controls iterative linear system solvers */
18
19 /**************************************************************/
20
21 /* Copyrights by ACcESS Australia 2003/04 */
22 /* Author: gross@access.edu.au */
23
24 /**************************************************************/
25
26 #include "Paso.h"
27 #include "SystemMatrix.h"
28 #include "Solver.h"
29 #include "esysUtils/blocktimer.h"
30
31 /***********************************************************************************/
32
33 /* free space */
34
35 void Paso_Solver_free(Paso_SystemMatrix* A) {
36 Paso_Preconditioner_free(A->solver);
37 A->solver=NULL;
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;
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,*scaling;
49 #ifdef PASO_MPI
50 double loc_norm;
51 #endif
52 dim_t i,totIter,cntIter,method;
53 bool_t finalizeIteration;
54 err_t errorCode=NO_ERROR;
55 dim_t numSol = Paso_SystemMatrix_getTotalNumCols(A);
56 dim_t numEqua = Paso_SystemMatrix_getTotalNumRows(A);
57 double blocktimer_precond, blocktimer_start = blocktimer_time();
58
59 Paso_resetError();
60 tolerance=options->tolerance;
61 if (tolerance < 100.* EPSILON) {
62 Paso_setError(VALUE_ERROR,"Paso_Solver: Tolerance is too small.");
63 }
64 if (tolerance >1.) {
65 Paso_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 Paso_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 Paso_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 Paso_setError(TYPE_ERROR,"Paso_Solver: Iterative solver requires a square matrix.");
77 return;
78 }
79 /* this for testing only */
80 if (method==PASO_NONLINEAR_GMRES) {
81 Paso_Function* F=NULL;
82 F=Paso_Function_LinearSystem_alloc(A,b,options);
83 Paso_Solver_solvePreconditioner(A,x,b);
84 errorCode=Paso_Solver_NewtonGMRES(F,x,options,pp);
85 if (errorCode!=NO_ERROR) {
86 Paso_setError(SYSTEM_ERROR,"Paso_Solver_NewtonGMRES: an error has occured.");
87 }
88 Paso_Function_LinearSystem_free(F);
89 return;
90 }
91 /* ========================= */
92 Performance_startMonitor(pp,PERFORMANCE_ALL);
93 if (Paso_noError()) {
94 /* get normalization */
95 scaling=Paso_SystemMatrix_borrowNormalization(A);
96 if (Paso_noError()) {
97 /* get the norm of the right hand side */
98 norm2_of_b=0.;
99 norm_max_of_b=0.;
100 #pragma omp parallel private(norm2_of_b_local,norm_max_of_b_local)
101 {
102 norm2_of_b_local=0.;
103 norm_max_of_b_local=0.;
104 #pragma omp for private(i) schedule(static)
105 for (i = 0; i < numEqua ; ++i) {
106 norm2_of_b_local += b[i] * b[i];
107 norm_max_of_b_local = MAX(ABS(scaling[i]*b[i]),norm_max_of_b_local);
108 }
109 #pragma omp critical
110 {
111 norm2_of_b += norm2_of_b_local;
112 norm_max_of_b = MAX(norm_max_of_b_local,norm_max_of_b);
113 }
114 }
115 /* TODO: use one call */
116 #ifdef PASO_MPI
117 {
118 loc_norm = norm2_of_b;
119 MPI_Allreduce(&loc_norm,&norm2_of_b, 1, MPI_DOUBLE, MPI_SUM, A->mpi_info->comm);
120 loc_norm = norm_max_of_b;
121 MPI_Allreduce(&loc_norm,&norm_max_of_b, 1, MPI_DOUBLE, MPI_MAX, A->mpi_info->comm);
122 }
123 #endif
124 norm2_of_b=sqrt(norm2_of_b);
125
126 /* if norm2_of_b==0 we are ready: x=0 */
127 if (norm2_of_b <=0.) {
128 #pragma omp parallel for private(i) schedule(static)
129 for (i = 0; i < numSol; i++) x[i]=0.;
130 if (options->verbose) printf("right hand side is identical zero.\n");
131 } else {
132 if (options->verbose) {
133 printf("Paso_Solver: l2/lmax-norm of right hand side is %e/%e.\n",norm2_of_b,norm_max_of_b);
134 printf("Paso_Solver: l2/lmax-stopping criterion is %e/%e.\n",norm2_of_b*tolerance,norm_max_of_b*tolerance);
135 switch (method) {
136 case PASO_BICGSTAB:
137 printf("Paso_Solver: Iterative method is BiCGStab.\n");
138 break;
139 case PASO_PCG:
140 printf("Paso_Solver: Iterative method is PCG.\n");
141 break;
142 case PASO_TFQMR:
143 printf("Paso_Solver: Iterative method is TFQMR.\n");
144 break;
145 case PASO_MINRES:
146 printf("Paso_Solver: Iterative method is MINRES.\n");
147 break;
148 case PASO_PRES20:
149 printf("Paso_Solver: Iterative method is PRES20.\n");
150 break;
151 case PASO_GMRES:
152 if (options->restart>0) {
153 printf("Paso_Solver: Iterative method is GMRES(%d,%d)\n",options->truncation,options->restart);
154 } else {
155 printf("Paso_Solver: Iterative method is GMRES(%d)\n",options->truncation);
156 }
157 break;
158 }
159 }
160 /* construct the preconditioner */
161 blocktimer_precond = blocktimer_time();
162 Performance_startMonitor(pp,PERFORMANCE_PRECONDITIONER_INIT);
163 Paso_Solver_setPreconditioner(A,options);
164 Performance_stopMonitor(pp,PERFORMANCE_PRECONDITIONER_INIT);
165 blocktimer_increment("Paso_Solver_setPreconditioner()", blocktimer_precond);
166 if (! Paso_noError()) return;
167 time_iter=Paso_timer();
168 /* get an initial guess by evaluating the preconditioner */
169 Paso_Solver_solvePreconditioner(A,x,b);
170 /* start the iteration process :*/
171 r=TMPMEMALLOC(numEqua,double);
172 Paso_checkPtr(r);
173 if (Paso_noError()) {
174 totIter = 1;
175 finalizeIteration = FALSE;
176 last_norm2_of_residual=norm2_of_b;
177 last_norm_max_of_residual=norm_max_of_b;
178 /* Loop */
179 while (! finalizeIteration) {
180 cntIter = options->iter_max - totIter;
181 finalizeIteration = TRUE;
182 /* Set initial residual. */
183 norm2_of_residual = 0;
184 norm_max_of_residual = 0;
185 #pragma omp parallel for private(i) schedule(static)
186 for (i = 0; i < numEqua; i++) r[i]=b[i];
187 Paso_SystemMatrix_MatrixVector_CSR_OFFSET0(DBLE(-1), A, x, DBLE(1), r);
188 #pragma omp parallel private(norm2_of_residual_local,norm_max_of_residual_local)
189 {
190 norm2_of_residual_local = 0;
191 norm_max_of_residual_local = 0;
192 #pragma omp for private(i) schedule(static)
193 for (i = 0; i < numEqua; i++) {
194 norm2_of_residual_local+= r[i] * r[i];
195 norm_max_of_residual_local=MAX(ABS(scaling[i]*r[i]),norm_max_of_residual_local);
196 }
197 #pragma omp critical
198 {
199 norm2_of_residual += norm2_of_residual_local;
200 norm_max_of_residual = MAX(norm_max_of_residual_local,norm_max_of_residual);
201 }
202 }
203 /* TODO: use one call */
204 #ifdef PASO_MPI
205 loc_norm = norm2_of_residual;
206 MPI_Allreduce(&loc_norm,&norm2_of_residual, 1, MPI_DOUBLE, MPI_SUM, A->mpi_info->comm);
207 loc_norm = norm_max_of_residual;
208 MPI_Allreduce(&loc_norm,&norm_max_of_residual, 1, MPI_DOUBLE, MPI_MAX, A->mpi_info->comm);
209 #endif
210 norm2_of_residual =sqrt(norm2_of_residual);
211
212 if (options->verbose) printf("Paso_Solver: Step %5d: l2/lmax-norm of residual is %e/%e",totIter,norm2_of_residual,norm_max_of_residual);
213 if (totIter>1 && norm2_of_residual>=last_norm2_of_residual && norm_max_of_residual>=last_norm_max_of_residual) {
214 if (options->verbose) printf(" divergence!\n");
215 Paso_setError(WARNING, "Paso_Solver: No improvement during iteration. Iterative solver gives up.");
216 } else {
217 /* */
218 if (norm2_of_residual>tolerance*norm2_of_b || norm_max_of_residual>tolerance*norm_max_of_b ) {
219 tol=tolerance*MIN(norm2_of_b,0.1*norm2_of_residual/norm_max_of_residual*norm_max_of_b);
220 if (options->verbose) printf(" (new tolerance = %e).\n",tol);
221 last_norm2_of_residual=norm2_of_residual;
222 last_norm_max_of_residual=norm_max_of_residual;
223 /* call the solver */
224 switch (method) {
225 case PASO_BICGSTAB:
226 errorCode = Paso_Solver_BiCGStab(A, r, x, &cntIter, &tol, pp);
227 break;
228 case PASO_PCG:
229 errorCode = Paso_Solver_PCG(A, r, x, &cntIter, &tol, pp);
230 break;
231 case PASO_TFQMR:
232 errorCode = Paso_Solver_TFQMR(A, r, x, &cntIter, &tol, pp);
233 break;
234 case PASO_MINRES:
235 errorCode = Paso_Solver_MINRES(A, r, x, &cntIter, &tol, pp);
236 break;
237 case PASO_PRES20:
238 errorCode = Paso_Solver_GMRES(A, r, x, &cntIter, &tol,5,20, pp);
239 break;
240 case PASO_GMRES:
241 errorCode = Paso_Solver_GMRES(A, r, x, &cntIter, &tol,options->truncation,options->restart, pp);
242 break;
243 }
244 totIter += cntIter;
245 /* error handling */
246 if (errorCode==NO_ERROR) {
247 finalizeIteration = FALSE;
248 } else if (errorCode==SOLVER_MAXITER_REACHED) {
249 Paso_setError(WARNING,"Paso_Solver: maximum number of iteration step reached.\nReturned solution does not fulfil stopping criterion.");
250 if (options->verbose) printf("Paso_Solver: Maximum number of iterations reached.!\n");
251 } else if (errorCode == SOLVER_INPUT_ERROR ) {
252 Paso_setError(SYSTEM_ERROR,"Paso_Solver: illegal dimension in iterative solver.");
253 if (options->verbose) printf("Paso_Solver: Internal error!\n");
254 } else if ( errorCode == SOLVER_BREAKDOWN ) {
255 if (cntIter <= 1) {
256 Paso_setError(ZERO_DIVISION_ERROR, "Paso_Solver: fatal break down in iterative solver.");
257 if (options->verbose) printf("Paso_Solver: Uncurable break down!\n");
258 } else {
259 if (options->verbose) printf("Paso_Solver: Breakdown at iter %d (residual = %e). Restarting ...\n", cntIter+totIter, tol);
260 finalizeIteration = FALSE;
261 }
262 } else {
263 if (options->verbose) printf(". convergence! \n");
264 }
265 }
266 } /* while */
267 }
268 MEMFREE(r);
269 time_iter=Paso_timer()-time_iter;
270 if (options->verbose) {
271 printf("\ntiming: Paso_Solver: %.4e sec\n",time_iter);
272 if (totIter>1) {
273 if(totIter==options->iter_max) {
274 printf("timing: Total MAX steps, time per iteration step: %.4e sec\n",time_iter/totIter);
275 } else {
276 printf("timing: Total %d steps, time per iteration step: %.4e sec\n",totIter,time_iter/totIter);
277 }
278 }
279 else {
280 printf("timing: Total 1 step, time per iteration step: %.4e sec\n",time_iter);
281 }
282 }
283 }
284 }
285 }
286 }
287 Performance_stopMonitor(pp,PERFORMANCE_ALL);
288 blocktimer_increment("Paso_Solver()", blocktimer_start);
289 }

Properties

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

  ViewVC Help
Powered by ViewVC 1.1.26