1 |
/******************************************************* |
2 |
* |
3 |
* Copyright (c) 2003-2008 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: gross@access.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; |
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,cntIter,method; |
52 |
bool_t finalizeIteration; |
53 |
err_t errorCode=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 |
|
58 |
Paso_resetError(); |
59 |
tolerance=options->tolerance; |
60 |
if (tolerance < 100.* EPSILON) { |
61 |
Paso_setError(VALUE_ERROR,"Paso_Solver: Tolerance is too small."); |
62 |
} |
63 |
if (tolerance >1.) { |
64 |
Paso_setError(VALUE_ERROR,"Paso_Solver: Tolerance mut be less than one."); |
65 |
} |
66 |
method=Paso_Options_getSolver(options->method,PASO_PASO,options->symmetric,A->mpi_info); |
67 |
/* check matrix type */ |
68 |
if ((A->type & MATRIX_FORMAT_CSC) || (A->type & MATRIX_FORMAT_OFFSET1) || (A->type & MATRIX_FORMAT_SYM) ) { |
69 |
Paso_setError(TYPE_ERROR,"Paso_Solver: Iterative solver requires CSR format with unsymmetric storage scheme and index offset 0."); |
70 |
} |
71 |
if (A->col_block_size != A->row_block_size) { |
72 |
Paso_setError(TYPE_ERROR,"Paso_Solver: Iterative solver requires row and column block sizes to be equal."); |
73 |
} |
74 |
if (Paso_SystemMatrix_getGlobalNumCols(A) != Paso_SystemMatrix_getGlobalNumRows(A)) { |
75 |
Paso_setError(TYPE_ERROR,"Paso_Solver: Iterative solver requires a square matrix."); |
76 |
return; |
77 |
} |
78 |
/* this for testing only */ |
79 |
if (method==PASO_NONLINEAR_GMRES) { |
80 |
Paso_Function* F=NULL; |
81 |
F=Paso_Function_LinearSystem_alloc(A,b,options); |
82 |
Paso_Solver_solvePreconditioner(A,x,b); |
83 |
errorCode=Paso_Solver_NewtonGMRES(F,x,options,pp); |
84 |
if (errorCode!=NO_ERROR) { |
85 |
Paso_setError(SYSTEM_ERROR,"Paso_Solver_NewtonGMRES: an error has occured."); |
86 |
} |
87 |
Paso_Function_LinearSystem_free(F); |
88 |
return; |
89 |
} |
90 |
/* ========================= */ |
91 |
Performance_startMonitor(pp,PERFORMANCE_ALL); |
92 |
if (Paso_noError()) { |
93 |
/* get normalization */ |
94 |
scaling=Paso_SystemMatrix_borrowNormalization(A); |
95 |
if (Paso_noError()) { |
96 |
/* get the norm of the right hand side */ |
97 |
norm2_of_b=0.; |
98 |
norm_max_of_b=0.; |
99 |
#pragma omp parallel private(norm2_of_b_local,norm_max_of_b_local) |
100 |
{ |
101 |
norm2_of_b_local=0.; |
102 |
norm_max_of_b_local=0.; |
103 |
#pragma omp for private(i) schedule(static) |
104 |
for (i = 0; i < numEqua ; ++i) { |
105 |
norm2_of_b_local += b[i] * b[i]; |
106 |
norm_max_of_b_local = MAX(ABS(scaling[i]*b[i]),norm_max_of_b_local); |
107 |
} |
108 |
#pragma omp critical |
109 |
{ |
110 |
norm2_of_b += norm2_of_b_local; |
111 |
norm_max_of_b = MAX(norm_max_of_b_local,norm_max_of_b); |
112 |
} |
113 |
} |
114 |
/* TODO: use one call */ |
115 |
#ifdef PASO_MPI |
116 |
{ |
117 |
loc_norm = norm2_of_b; |
118 |
MPI_Allreduce(&loc_norm,&norm2_of_b, 1, MPI_DOUBLE, MPI_SUM, A->mpi_info->comm); |
119 |
loc_norm = norm_max_of_b; |
120 |
MPI_Allreduce(&loc_norm,&norm_max_of_b, 1, MPI_DOUBLE, MPI_MAX, A->mpi_info->comm); |
121 |
} |
122 |
#endif |
123 |
norm2_of_b=sqrt(norm2_of_b); |
124 |
|
125 |
/* if norm2_of_b==0 we are ready: x=0 */ |
126 |
if (norm2_of_b <=0.) { |
127 |
#pragma omp parallel for private(i) schedule(static) |
128 |
for (i = 0; i < numSol; i++) x[i]=0.; |
129 |
if (options->verbose) printf("right hand side is identical zero.\n"); |
130 |
} else { |
131 |
if (options->verbose) { |
132 |
printf("Paso_Solver: l2/lmax-norm of right hand side is %e/%e.\n",norm2_of_b,norm_max_of_b); |
133 |
printf("Paso_Solver: l2/lmax-stopping criterion is %e/%e.\n",norm2_of_b*tolerance,norm_max_of_b*tolerance); |
134 |
switch (method) { |
135 |
case PASO_BICGSTAB: |
136 |
printf("Paso_Solver: Iterative method is BiCGStab.\n"); |
137 |
break; |
138 |
case PASO_PCG: |
139 |
printf("Paso_Solver: Iterative method is PCG.\n"); |
140 |
break; |
141 |
case PASO_TFQMR: |
142 |
printf("Paso_Solver: Iterative method is TFQMR.\n"); |
143 |
break; |
144 |
case PASO_MINRES: |
145 |
printf("Paso_Solver: Iterative method is MINRES.\n"); |
146 |
break; |
147 |
case PASO_PRES20: |
148 |
printf("Paso_Solver: Iterative method is PRES20.\n"); |
149 |
break; |
150 |
case PASO_GMRES: |
151 |
if (options->restart>0) { |
152 |
printf("Paso_Solver: Iterative method is GMRES(%d,%d)\n",options->truncation,options->restart); |
153 |
} else { |
154 |
printf("Paso_Solver: Iterative method is GMRES(%d)\n",options->truncation); |
155 |
} |
156 |
break; |
157 |
} |
158 |
} |
159 |
/* construct the preconditioner */ |
160 |
blocktimer_precond = blocktimer_time(); |
161 |
Performance_startMonitor(pp,PERFORMANCE_PRECONDITIONER_INIT); |
162 |
Paso_Solver_setPreconditioner(A,options); |
163 |
Performance_stopMonitor(pp,PERFORMANCE_PRECONDITIONER_INIT); |
164 |
blocktimer_increment("Paso_Solver_setPreconditioner()", blocktimer_precond); |
165 |
if (! Paso_noError()) return; |
166 |
time_iter=Paso_timer(); |
167 |
/* get an initial guess by evaluating the preconditioner */ |
168 |
Paso_Solver_solvePreconditioner(A,x,b); |
169 |
/* start the iteration process :*/ |
170 |
r=TMPMEMALLOC(numEqua,double); |
171 |
Paso_checkPtr(r); |
172 |
if (Paso_noError()) { |
173 |
totIter = 1; |
174 |
finalizeIteration = FALSE; |
175 |
last_norm2_of_residual=norm2_of_b; |
176 |
last_norm_max_of_residual=norm_max_of_b; |
177 |
/* Loop */ |
178 |
while (! finalizeIteration) { |
179 |
cntIter = options->iter_max - totIter; |
180 |
finalizeIteration = TRUE; |
181 |
/* Set initial residual. */ |
182 |
norm2_of_residual = 0; |
183 |
norm_max_of_residual = 0; |
184 |
#pragma omp parallel for private(i) schedule(static) |
185 |
for (i = 0; i < numEqua; i++) r[i]=b[i]; |
186 |
Paso_SystemMatrix_MatrixVector_CSR_OFFSET0(DBLE(-1), A, x, DBLE(1), r); |
187 |
#pragma omp parallel private(norm2_of_residual_local,norm_max_of_residual_local) |
188 |
{ |
189 |
norm2_of_residual_local = 0; |
190 |
norm_max_of_residual_local = 0; |
191 |
#pragma omp for private(i) schedule(static) |
192 |
for (i = 0; i < numEqua; i++) { |
193 |
norm2_of_residual_local+= r[i] * r[i]; |
194 |
norm_max_of_residual_local=MAX(ABS(scaling[i]*r[i]),norm_max_of_residual_local); |
195 |
} |
196 |
#pragma omp critical |
197 |
{ |
198 |
norm2_of_residual += norm2_of_residual_local; |
199 |
norm_max_of_residual = MAX(norm_max_of_residual_local,norm_max_of_residual); |
200 |
} |
201 |
} |
202 |
/* TODO: use one call */ |
203 |
#ifdef PASO_MPI |
204 |
loc_norm = norm2_of_residual; |
205 |
MPI_Allreduce(&loc_norm,&norm2_of_residual, 1, MPI_DOUBLE, MPI_SUM, A->mpi_info->comm); |
206 |
loc_norm = norm_max_of_residual; |
207 |
MPI_Allreduce(&loc_norm,&norm_max_of_residual, 1, MPI_DOUBLE, MPI_MAX, A->mpi_info->comm); |
208 |
#endif |
209 |
norm2_of_residual =sqrt(norm2_of_residual); |
210 |
|
211 |
if (options->verbose) printf("Paso_Solver: Step %5d: l2/lmax-norm of residual is %e/%e",totIter,norm2_of_residual,norm_max_of_residual); |
212 |
if (totIter>1 && norm2_of_residual>=last_norm2_of_residual && norm_max_of_residual>=last_norm_max_of_residual) { |
213 |
if (options->verbose) printf(" divergence!\n"); |
214 |
Paso_setError(WARNING, "Paso_Solver: No improvement during iteration. Iterative solver gives up."); |
215 |
} else { |
216 |
/* */ |
217 |
if (norm2_of_residual>tolerance*norm2_of_b || norm_max_of_residual>tolerance*norm_max_of_b ) { |
218 |
tol=tolerance*MIN(norm2_of_b,0.1*norm2_of_residual/norm_max_of_residual*norm_max_of_b); |
219 |
if (options->verbose) printf(" (new tolerance = %e).\n",tol); |
220 |
last_norm2_of_residual=norm2_of_residual; |
221 |
last_norm_max_of_residual=norm_max_of_residual; |
222 |
/* call the solver */ |
223 |
switch (method) { |
224 |
case PASO_BICGSTAB: |
225 |
errorCode = Paso_Solver_BiCGStab(A, r, x, &cntIter, &tol, pp); |
226 |
break; |
227 |
case PASO_PCG: |
228 |
errorCode = Paso_Solver_PCG(A, r, x, &cntIter, &tol, pp); |
229 |
break; |
230 |
case PASO_TFQMR: |
231 |
errorCode = Paso_Solver_TFQMR(A, r, x, &cntIter, &tol, pp); |
232 |
break; |
233 |
case PASO_MINRES: |
234 |
errorCode = Paso_Solver_MINRES(A, r, x, &cntIter, &tol, pp); |
235 |
break; |
236 |
case PASO_PRES20: |
237 |
errorCode = Paso_Solver_GMRES(A, r, x, &cntIter, &tol,5,20, pp); |
238 |
break; |
239 |
case PASO_GMRES: |
240 |
errorCode = Paso_Solver_GMRES(A, r, x, &cntIter, &tol,options->truncation,options->restart, pp); |
241 |
break; |
242 |
} |
243 |
totIter += cntIter; |
244 |
/* error handling */ |
245 |
if (errorCode==NO_ERROR) { |
246 |
finalizeIteration = FALSE; |
247 |
} else if (errorCode==SOLVER_MAXITER_REACHED) { |
248 |
Paso_setError(WARNING,"Paso_Solver: maximum number of iteration step reached.\nReturned solution does not fulfil stopping criterion."); |
249 |
if (options->verbose) printf("Paso_Solver: Maximum number of iterations reached.!\n"); |
250 |
} else if (errorCode == SOLVER_INPUT_ERROR ) { |
251 |
Paso_setError(SYSTEM_ERROR,"Paso_Solver: illegal dimension in iterative solver."); |
252 |
if (options->verbose) printf("Paso_Solver: Internal error!\n"); |
253 |
} else if ( errorCode == SOLVER_BREAKDOWN ) { |
254 |
if (cntIter <= 1) { |
255 |
Paso_setError(ZERO_DIVISION_ERROR, "Paso_Solver: fatal break down in iterative solver."); |
256 |
if (options->verbose) printf("Paso_Solver: Uncurable break down!\n"); |
257 |
} else { |
258 |
if (options->verbose) printf("Paso_Solver: Breakdown at iter %d (residual = %e). Restarting ...\n", cntIter+totIter, tol); |
259 |
finalizeIteration = FALSE; |
260 |
} |
261 |
} else { |
262 |
if (options->verbose) printf(". convergence! \n"); |
263 |
} |
264 |
} |
265 |
} /* while */ |
266 |
} |
267 |
MEMFREE(r); |
268 |
time_iter=Paso_timer()-time_iter; |
269 |
if (options->verbose) { |
270 |
printf("\ntiming: Paso_Solver: %.4e sec\n",time_iter); |
271 |
if (totIter>1) { |
272 |
if(totIter==options->iter_max) { |
273 |
printf("timing: Total MAX steps, time per iteration step: %.4e sec\n",time_iter/totIter); |
274 |
} else { |
275 |
printf("timing: Total %d steps, time per iteration step: %.4e sec\n",totIter,time_iter/totIter); |
276 |
} |
277 |
} |
278 |
else { |
279 |
printf("timing: Total 1 step, time per iteration step: %.4e sec\n",time_iter); |
280 |
} |
281 |
} |
282 |
} |
283 |
} |
284 |
} |
285 |
} |
286 |
Performance_stopMonitor(pp,PERFORMANCE_ALL); |
287 |
blocktimer_increment("Paso_Solver()", blocktimer_start); |
288 |
} |