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

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

Parent Directory Parent Directory | Revision Log Revision Log


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

Properties

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

  ViewVC Help
Powered by ViewVC 1.1.26