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

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

Parent Directory Parent Directory | Revision Log Revision Log


Revision 1312 - (show annotations)
Mon Sep 24 06:18:44 2007 UTC (11 years, 7 months ago) by ksteube
File MIME type: text/plain
File size: 12938 byte(s)
The MPI branch is hereby closed. All future work should be in trunk.

Previously in revision 1295 I merged the latest changes to trunk into trunk-mpi-branch.
In this revision I copied all files from trunk-mpi-branch over the corresponding
trunk files. I did not use 'svn merge', it was a copy.

1
2 /* $Id$ */
3
4 /*******************************************************
5 *
6 * Copyright 2003-2007 by ACceSS MNRF
7 * Copyright 2007 by University of Queensland
8 *
9 * http://esscc.uq.edu.au
10 * Primary Business: Queensland, Australia
11 * Licensed under the Open Software License version 3.0
12 * http://www.opensource.org/licenses/osl-3.0.php
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: gross@access.edu.au */
24
25 /**************************************************************/
26
27 #include "Paso.h"
28 #include "SystemMatrix.h"
29 #include "Solver.h"
30 #include "escript/blocktimer.h"
31
32 /***********************************************************************************/
33
34 /* free space */
35
36 void Paso_Solver_free(Paso_SystemMatrix* A) {
37 Paso_Preconditioner_free(A->solver);
38 A->solver=NULL;
39 }
40 /* call the iterative solver: */
41
42 void Paso_Solver(Paso_SystemMatrix* A,double* x,double* b,Paso_Options* options,Paso_Performance* pp) {
43 double norm2_of_b,tol,tolerance,time_iter,*r=NULL,norm2_of_residual,last_norm2_of_residual,norm_max_of_b,
44 norm2_of_b_local,norm_max_of_b_local,norm2_of_residual_local,norm_max_of_residual_local,norm_max_of_residual,
45 last_norm_max_of_residual,*scaling, loc_norm;
46 dim_t i,totIter,cntIter,method;
47 bool_t finalizeIteration;
48 err_t errorCode;
49 dim_t numSol = Paso_SystemMatrix_getTotalNumCols(A);
50 dim_t numEqua = Paso_SystemMatrix_getTotalNumRows(A);
51 double blocktimer_start = blocktimer_time();
52
53 tolerance=MAX(options->tolerance,EPSILON);
54 Paso_resetError();
55 method=Paso_Options_getSolver(options->method,PASO_PASO,options->symmetric);
56 /* check matrix type */
57 if ((A->type & MATRIX_FORMAT_CSC) || (A->type & MATRIX_FORMAT_OFFSET1) || (A->type & MATRIX_FORMAT_SYM) ) {
58 Paso_setError(TYPE_ERROR,"Iterative solver requires CSR format with unsymmetric storage scheme and index offset 0.");
59 }
60 if (A->col_block_size != A->row_block_size) {
61 Paso_setError(TYPE_ERROR,"Iterative solver requires row and column block sizes to be equal.");
62 }
63 if (Paso_SystemMatrix_getGlobalNumCols(A) != Paso_SystemMatrix_getGlobalNumRows(A)) {
64 Paso_setError(TYPE_ERROR,"Iterative solver requires requires a square matrix.");
65 return;
66 }
67 Performance_startMonitor(pp,PERFORMANCE_ALL);
68 if (Paso_noError()) {
69 /* get normalization */
70 scaling=Paso_SystemMatrix_borrowNormalization(A);
71 if (Paso_noError()) {
72 /* get the norm of the right hand side */
73 norm2_of_b=0.;
74 norm_max_of_b=0.;
75 #pragma omp parallel private(norm2_of_b_local,norm_max_of_b_local)
76 {
77 norm2_of_b_local=0.;
78 norm_max_of_b_local=0.;
79 #pragma omp for private(i) schedule(static)
80 for (i = 0; i < numEqua ; ++i) {
81 norm2_of_b_local += b[i] * b[i];
82 norm_max_of_b_local = MAX(ABS(scaling[i]*b[i]),norm_max_of_b_local);
83 }
84 #pragma omp critical
85 {
86 norm2_of_b += norm2_of_b_local;
87 norm_max_of_b = MAX(norm_max_of_b_local,norm_max_of_b);
88 }
89 }
90 /* TODO: use one call */
91 #ifdef PASO_MPI
92 {
93 loc_norm = norm2_of_b;
94 MPI_Allreduce(&loc_norm,&norm2_of_b, 1, MPI_DOUBLE, MPI_SUM, A->mpi_info->comm);
95 loc_norm = norm_max_of_b;
96 MPI_Allreduce(&loc_norm,&norm_max_of_b, 1, MPI_DOUBLE, MPI_MAX, A->mpi_info->comm);
97 }
98 #endif
99 norm2_of_b=sqrt(norm2_of_b);
100
101 /* if norm2_of_b==0 we are ready: x=0 */
102 if (norm2_of_b <=0.) {
103 #pragma omp parallel for private(i) schedule(static)
104 for (i = 0; i < numSol; i++) x[i]=0.;
105 if (options->verbose) printf("right hand side is identical zero.\n");
106 } else {
107 if (options->verbose) {
108 printf("l2/lmax-norm of right hand side is %e/%e.\n",norm2_of_b,norm_max_of_b);
109 printf("l2/lmax-stopping criterion is %e/%e.\n",norm2_of_b*tolerance,norm_max_of_b*tolerance);
110 switch (method) {
111 case PASO_BICGSTAB:
112 printf("Iterative method is BiCGStab.\n");
113 break;
114 case PASO_PCG:
115 printf("Iterative method is PCG.\n");
116 break;
117 case PASO_PRES20:
118 printf("Iterative method is PRES20.\n");
119 break;
120 case PASO_GMRES:
121 if (options->restart>0) {
122 printf("Iterative method is GMRES(%d,%d)\n",options->truncation,options->restart);
123 } else {
124 printf("Iterative method is GMRES(%d)\n",options->truncation);
125 }
126 break;
127 }
128 }
129 /* construct the preconditioner */
130
131 Performance_startMonitor(pp,PERFORMANCE_PRECONDITIONER_INIT);
132 Paso_Solver_setPreconditioner(A,options);
133 Performance_stopMonitor(pp,PERFORMANCE_PRECONDITIONER_INIT);
134 if (! Paso_noError()) return;
135
136 time_iter=Paso_timer();
137 Paso_SystemMatrix_allocBuffer(A);
138 /* get an initial guess by evaluating the preconditioner */
139 #pragma omp parallel
140 {
141 Paso_Solver_solvePreconditioner(A,x,b);
142 }
143 /* start the iteration process :*/
144 r=TMPMEMALLOC(numEqua,double);
145 Paso_checkPtr(r);
146 if (Paso_noError()) {
147 totIter = 0;
148 finalizeIteration = FALSE;
149 last_norm2_of_residual=norm2_of_b;
150 last_norm_max_of_residual=norm_max_of_b;
151 while (! finalizeIteration) {
152 cntIter = options->iter_max - totIter;
153 finalizeIteration = TRUE;
154 /* Set initial residual. */
155 norm2_of_residual = 0;
156 norm_max_of_residual = 0;
157 #pragma omp parallel private(norm2_of_residual_local,norm_max_of_residual_local)
158 {
159 #pragma omp for private(i) schedule(static)
160 for (i = 0; i < numEqua; i++) {
161 r[i]=b[i];
162 }
163
164 Paso_SystemMatrix_MatrixVector_CSR_OFFSET0(DBLE(-1), A, x, DBLE(1), r);
165
166 norm2_of_residual_local = 0;
167 norm_max_of_residual_local = 0;
168 #pragma omp for private(i) schedule(static)
169 for (i = 0; i < numEqua; i++) {
170 norm2_of_residual_local+= r[i] * r[i];
171 norm_max_of_residual_local=MAX(ABS(scaling[i]*r[i]),norm_max_of_residual_local);
172 }
173 #pragma omp critical
174 {
175 norm2_of_residual += norm2_of_residual_local;
176 norm_max_of_residual = MAX(norm_max_of_residual_local,norm_max_of_residual);
177 }
178 }
179 /* TODO: use one call */
180 #ifdef PASO_MPI
181 {
182 loc_norm = norm2_of_residual;
183 MPI_Allreduce(&loc_norm,&norm2_of_residual, 1, MPI_DOUBLE, MPI_SUM, A->mpi_info->comm);
184 loc_norm = norm_max_of_residual;
185 MPI_Allreduce(&loc_norm,&norm_max_of_residual, 1, MPI_DOUBLE, MPI_MAX, A->mpi_info->comm);
186 }
187 #endif
188 norm2_of_residual =sqrt(norm2_of_residual);
189
190 if (options->verbose) printf("Step %5d: l2/lmax-norm of residual is %e/%e",totIter,norm2_of_residual,norm_max_of_residual);
191 if (totIter>0 && norm2_of_residual>=last_norm2_of_residual && norm_max_of_residual>=last_norm_max_of_residual) {
192 if (options->verbose) printf(" divergence!\n");
193 Paso_setError(WARNING, "No improvement during iteration. Iterative solver gives up.");
194 } else {
195 /* */
196 if (norm2_of_residual>tolerance*norm2_of_b || norm_max_of_residual>tolerance*norm_max_of_b ) {
197 tol=tolerance*MIN(norm2_of_b,0.1*norm2_of_residual/norm_max_of_residual*norm_max_of_b);
198 if (options->verbose) printf(" (new tolerance = %e).\n",tol);
199 last_norm2_of_residual=norm2_of_residual;
200 last_norm_max_of_residual=norm_max_of_residual;
201 /* call the solver */
202 switch (method) {
203 case PASO_BICGSTAB:
204 errorCode = Paso_Solver_BiCGStab(A, r, x, &cntIter, &tol, pp);
205 break;
206 case PASO_PCG:
207 errorCode = Paso_Solver_PCG(A, r, x, &cntIter, &tol, pp);
208 break;
209 case PASO_PRES20:
210 errorCode = Paso_Solver_GMRES(A, r, x, &cntIter, &tol,5,20, pp);
211 break;
212 case PASO_GMRES:
213 errorCode = Paso_Solver_GMRES(A, r, x, &cntIter, &tol,options->truncation,options->restart, pp);
214 break;
215 }
216 totIter += cntIter;
217
218 /* error handling */
219 if (errorCode==NO_ERROR) {
220 finalizeIteration = FALSE;
221 } else if (errorCode==SOLVER_MAXITER_REACHED) {
222 Paso_setError(WARNING,"maximum number of iteration step reached.\nReturned solution does not fulfil stopping criterion.");
223 if (options->verbose) printf("Maximum number of iterations reached.!\n");
224 } else if (errorCode == SOLVER_INPUT_ERROR ) {
225 Paso_setError(SYSTEM_ERROR,"illegal dimension in iterative solver.");
226 if (options->verbose) printf("Internal error!\n");
227 } else if ( errorCode == SOLVER_BREAKDOWN ) {
228 if (cntIter <= 1) {
229 Paso_setError(ZERO_DIVISION_ERROR, "fatal break down in iterative solver.");
230 if (options->verbose) printf("Uncurable break down!\n");
231 } else {
232 if (options->verbose) printf("Breakdown at iter %d (residual = %e). Restarting ...\n", cntIter+totIter, tol);
233 finalizeIteration = FALSE;
234 }
235 } else if (errorCode == SOLVER_MEMORY_ERROR) {
236 Paso_setError(MEMORY_ERROR,"memory allocation failed.");
237 if (options->verbose) printf("Memory allocation failed!\n");
238 } else if (errorCode !=SOLVER_NO_ERROR ) {
239 Paso_setError(SYSTEM_ERROR,"unidentified error in iterative solver.");
240 if (options->verbose) printf("Unidentified error!\n");
241 }
242 } else {
243 if (options->verbose) printf(". convergence! \n");
244 }
245 }
246 } /* while */
247 }
248 MEMFREE(r);
249 Paso_SystemMatrix_freeBuffer(A);
250 time_iter=Paso_timer()-time_iter;
251 if (options->verbose) {
252 printf("timing: solver: %.4e sec\n",time_iter);
253 if (totIter>0) printf("timing: per iteration step: %.4e sec\n",time_iter/totIter);
254 }
255 }
256 }
257 }
258 Performance_stopMonitor(pp,PERFORMANCE_ALL);
259 blocktimer_increment("Paso_Solver()", blocktimer_start);
260 }

Properties

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

  ViewVC Help
Powered by ViewVC 1.1.26