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

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

Parent Directory Parent Directory | Revision Log Revision Log


Revision 2358 - (hide annotations)
Wed Apr 1 22:25:24 2009 UTC (10 years ago) by gross
File MIME type: text/plain
File size: 13926 byte(s)
mpi is not really ready yet.
1 ksteube 1312 /*******************************************************
2 ksteube 1811 *
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 dhawcroft 631
13 ksteube 1811
14 jgs 150 /**************************************************************/
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 gross 700 #include "Paso.h"
26     #include "SystemMatrix.h"
27 jgs 150 #include "Solver.h"
28 phornby 2078 #include "esysUtils/blocktimer.h"
29 jgs 150
30     /***********************************************************************************/
31    
32     /* free space */
33    
34     void Paso_Solver_free(Paso_SystemMatrix* A) {
35 phornby 1628 Paso_Preconditioner_free(A->solver);
36     A->solver=NULL;
37 jgs 150 }
38     /* call the iterative solver: */
39    
40 phornby 1628 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 jfenwick 1981 err_t errorCode=NO_ERROR;
54 phornby 1628 dim_t numSol = Paso_SystemMatrix_getTotalNumCols(A);
55     dim_t numEqua = Paso_SystemMatrix_getTotalNumRows(A);
56 ksteube 1707 double blocktimer_precond, blocktimer_start = blocktimer_time();
57 gross 415
58 gross 1639 Paso_resetError();
59 gross 2156 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 gross 2315 method=Paso_Options_getSolver(options->method,PASO_PASO,options->symmetric,A->mpi_info);
67 gross 1639 /* check matrix type */
68     if ((A->type & MATRIX_FORMAT_CSC) || (A->type & MATRIX_FORMAT_OFFSET1) || (A->type & MATRIX_FORMAT_SYM) ) {
69 gross 2156 Paso_setError(TYPE_ERROR,"Paso_Solver: Iterative solver requires CSR format with unsymmetric storage scheme and index offset 0.");
70 gross 1639 }
71     if (A->col_block_size != A->row_block_size) {
72 gross 2156 Paso_setError(TYPE_ERROR,"Paso_Solver: Iterative solver requires row and column block sizes to be equal.");
73 gross 1639 }
74     if (Paso_SystemMatrix_getGlobalNumCols(A) != Paso_SystemMatrix_getGlobalNumRows(A)) {
75 gross 2156 Paso_setError(TYPE_ERROR,"Paso_Solver: Iterative solver requires a square matrix.");
76 gross 1639 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 gross 1798 Paso_Solver_solvePreconditioner(A,x,b);
83 gross 1639 errorCode=Paso_Solver_NewtonGMRES(F,x,options,pp);
84 gross 1798 if (errorCode!=NO_ERROR) {
85 gross 1639 Paso_setError(SYSTEM_ERROR,"Paso_Solver_NewtonGMRES: an error has occured.");
86     }
87     Paso_Function_LinearSystem_free(F);
88 gross 1798 return;
89 gross 1639 }
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 jgs 150
125 phornby 1628 /* if norm2_of_b==0 we are ready: x=0 */
126     if (norm2_of_b <=0.) {
127 gross 1798 #pragma omp parallel for private(i) schedule(static)
128 phornby 1628 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 gross 2156 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 phornby 1628 switch (method) {
135     case PASO_BICGSTAB:
136 gross 2156 printf("Paso_Solver: Iterative method is BiCGStab.\n");
137 phornby 1628 break;
138     case PASO_PCG:
139 gross 2156 printf("Paso_Solver: Iterative method is PCG.\n");
140 phornby 1628 break;
141 artak 1703 case PASO_TFQMR:
142 gross 2156 printf("Paso_Solver: Iterative method is TFQMR.\n");
143 artak 1703 break;
144 artak 1787 case PASO_MINRES:
145 gross 2156 printf("Paso_Solver: Iterative method is MINRES.\n");
146 artak 1787 break;
147 phornby 1628 case PASO_PRES20:
148 gross 2156 printf("Paso_Solver: Iterative method is PRES20.\n");
149 phornby 1628 break;
150     case PASO_GMRES:
151     if (options->restart>0) {
152 gross 2156 printf("Paso_Solver: Iterative method is GMRES(%d,%d)\n",options->truncation,options->restart);
153 phornby 1628 } else {
154 gross 2156 printf("Paso_Solver: Iterative method is GMRES(%d)\n",options->truncation);
155 gross 584 }
156 phornby 1628 break;
157     }
158     }
159     /* construct the preconditioner */
160 ksteube 1707 blocktimer_precond = blocktimer_time();
161 phornby 1628 Performance_startMonitor(pp,PERFORMANCE_PRECONDITIONER_INIT);
162     Paso_Solver_setPreconditioner(A,options);
163     Performance_stopMonitor(pp,PERFORMANCE_PRECONDITIONER_INIT);
164 ksteube 1707 blocktimer_increment("Paso_Solver_setPreconditioner()", blocktimer_precond);
165 phornby 1628 if (! Paso_noError()) return;
166 ksteube 1312 time_iter=Paso_timer();
167 gross 584 /* get an initial guess by evaluating the preconditioner */
168 gross 1556 Paso_Solver_solvePreconditioner(A,x,b);
169 gross 584 /* start the iteration process :*/
170 ksteube 1312 r=TMPMEMALLOC(numEqua,double);
171     Paso_checkPtr(r);
172     if (Paso_noError()) {
173 artak 2280 totIter = 1;
174 ksteube 1312 finalizeIteration = FALSE;
175     last_norm2_of_residual=norm2_of_b;
176     last_norm_max_of_residual=norm_max_of_b;
177 artak 1787 /* Loop */
178 ksteube 1312 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 gross 1556 #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 ksteube 1312 #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 gross 584 }
202 ksteube 1312 /* 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 gross 2156 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 artak 2280 if (totIter>1 && norm2_of_residual>=last_norm2_of_residual && norm_max_of_residual>=last_norm_max_of_residual) {
213 phornby 1628 if (options->verbose) printf(" divergence!\n");
214 gross 2156 Paso_setError(WARNING, "Paso_Solver: No improvement during iteration. Iterative solver gives up.");
215 phornby 1628 } 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 artak 1819 errorCode = Paso_Solver_PCG(A, r, x, &cntIter, &tol, pp);
229 phornby 1628 break;
230 artak 1703 case PASO_TFQMR:
231 artak 1787 errorCode = Paso_Solver_TFQMR(A, r, x, &cntIter, &tol, pp);
232 artak 1703 break;
233 artak 1787 case PASO_MINRES:
234 artak 1862 errorCode = Paso_Solver_MINRES(A, r, x, &cntIter, &tol, pp);
235 artak 1787 break;
236 phornby 1628 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 gross 2156 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 phornby 1628 } else if (errorCode == SOLVER_INPUT_ERROR ) {
251 gross 2156 Paso_setError(SYSTEM_ERROR,"Paso_Solver: illegal dimension in iterative solver.");
252     if (options->verbose) printf("Paso_Solver: Internal error!\n");
253 phornby 1628 } else if ( errorCode == SOLVER_BREAKDOWN ) {
254     if (cntIter <= 1) {
255 gross 2156 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 phornby 1628 } else {
258 gross 2156 if (options->verbose) printf("Paso_Solver: Breakdown at iter %d (residual = %e). Restarting ...\n", cntIter+totIter, tol);
259 phornby 1628 finalizeIteration = FALSE;
260 ksteube 1312 }
261     } else {
262     if (options->verbose) printf(". convergence! \n");
263     }
264 jgs 150 }
265 ksteube 1312 } /* while */
266     }
267 gross 584 MEMFREE(r);
268     time_iter=Paso_timer()-time_iter;
269     if (options->verbose) {
270 gross 2156 printf("\ntiming: Paso_Solver: %.4e sec\n",time_iter);
271 artak 2280 if (totIter>1) {
272 artak 2240 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 artak 2280 printf("timing: Total 1 step, time per iteration step: %.4e sec\n",time_iter);
280 artak 2240 }
281 jgs 150 }
282     }
283     }
284 gross 584 }
285 phornby 1628 }
286     Performance_stopMonitor(pp,PERFORMANCE_ALL);
287     blocktimer_increment("Paso_Solver()", blocktimer_start);
288 jgs 150 }

Properties

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

  ViewVC Help
Powered by ViewVC 1.1.26