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

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

Parent Directory Parent Directory | Revision Log Revision Log


Revision 1798 - (hide annotations)
Wed Sep 17 06:21:12 2008 UTC (10 years, 7 months ago) by gross
File MIME type: text/plain
File size: 13037 byte(s)
Fixes for the JacobeanFreeNewton scheme. Still needs to be tested under OPENMP but runs under MPI.


1 ksteube 1312
2 jgs 150 /* $Id$ */
3    
4 ksteube 1312 /*******************************************************
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 dhawcroft 631
16 jgs 150 /**************************************************************/
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 gross 700 #include "Paso.h"
28     #include "SystemMatrix.h"
29 jgs 150 #include "Solver.h"
30 ksteube 1312 #include "escript/blocktimer.h"
31 jgs 150
32     /***********************************************************************************/
33    
34     /* free space */
35    
36     void Paso_Solver_free(Paso_SystemMatrix* A) {
37 phornby 1628 Paso_Preconditioner_free(A->solver);
38     A->solver=NULL;
39 jgs 150 }
40     /* call the iterative solver: */
41    
42 phornby 1628 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 ksteube 1707 double blocktimer_precond, blocktimer_start = blocktimer_time();
59 gross 415
60 gross 1639 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 artak 1787 Paso_setError(TYPE_ERROR,"Iterative solver requires a square matrix.");
72 gross 1639 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 gross 1798 Paso_Solver_solvePreconditioner(A,x,b);
79 gross 1639 errorCode=Paso_Solver_NewtonGMRES(F,x,options,pp);
80 gross 1798 if (errorCode!=NO_ERROR) {
81 gross 1639 Paso_setError(SYSTEM_ERROR,"Paso_Solver_NewtonGMRES: an error has occured.");
82     }
83     Paso_Function_LinearSystem_free(F);
84 gross 1798 return;
85 gross 1639 }
86     /* ========================= */
87     Performance_startMonitor(pp,PERFORMANCE_ALL);
88     if (Paso_noError()) {
89     /* get normalization */
90     scaling=Paso_SystemMatrix_borrowNormalization(A);
91     if (Paso_noError()) {
92     /* get the norm of the right hand side */
93     norm2_of_b=0.;
94     norm_max_of_b=0.;
95     #pragma omp parallel private(norm2_of_b_local,norm_max_of_b_local)
96     {
97     norm2_of_b_local=0.;
98     norm_max_of_b_local=0.;
99     #pragma omp for private(i) schedule(static)
100     for (i = 0; i < numEqua ; ++i) {
101     norm2_of_b_local += b[i] * b[i];
102     norm_max_of_b_local = MAX(ABS(scaling[i]*b[i]),norm_max_of_b_local);
103     }
104     #pragma omp critical
105     {
106     norm2_of_b += norm2_of_b_local;
107     norm_max_of_b = MAX(norm_max_of_b_local,norm_max_of_b);
108     }
109     }
110     /* TODO: use one call */
111     #ifdef PASO_MPI
112     {
113     loc_norm = norm2_of_b;
114     MPI_Allreduce(&loc_norm,&norm2_of_b, 1, MPI_DOUBLE, MPI_SUM, A->mpi_info->comm);
115     loc_norm = norm_max_of_b;
116     MPI_Allreduce(&loc_norm,&norm_max_of_b, 1, MPI_DOUBLE, MPI_MAX, A->mpi_info->comm);
117     }
118     #endif
119     norm2_of_b=sqrt(norm2_of_b);
120 jgs 150
121 phornby 1628 /* if norm2_of_b==0 we are ready: x=0 */
122     if (norm2_of_b <=0.) {
123 gross 1798 #pragma omp parallel for private(i) schedule(static)
124 phornby 1628 for (i = 0; i < numSol; i++) x[i]=0.;
125     if (options->verbose) printf("right hand side is identical zero.\n");
126     } else {
127     if (options->verbose) {
128     printf("l2/lmax-norm of right hand side is %e/%e.\n",norm2_of_b,norm_max_of_b);
129     printf("l2/lmax-stopping criterion is %e/%e.\n",norm2_of_b*tolerance,norm_max_of_b*tolerance);
130     switch (method) {
131     case PASO_BICGSTAB:
132     printf("Iterative method is BiCGStab.\n");
133     break;
134     case PASO_PCG:
135     printf("Iterative method is PCG.\n");
136     break;
137 artak 1703 case PASO_TFQMR:
138     printf("Iterative method is TFQMR.\n");
139     break;
140 artak 1787 case PASO_MINRES:
141     printf("Iterative method is MINRES.\n");
142     break;
143 phornby 1628 case PASO_PRES20:
144     printf("Iterative method is PRES20.\n");
145     break;
146     case PASO_GMRES:
147     if (options->restart>0) {
148     printf("Iterative method is GMRES(%d,%d)\n",options->truncation,options->restart);
149     } else {
150     printf("Iterative method is GMRES(%d)\n",options->truncation);
151 gross 584 }
152 phornby 1628 break;
153     }
154     }
155     /* construct the preconditioner */
156 gross 584
157 ksteube 1707 blocktimer_precond = blocktimer_time();
158 phornby 1628 Performance_startMonitor(pp,PERFORMANCE_PRECONDITIONER_INIT);
159     Paso_Solver_setPreconditioner(A,options);
160     Performance_stopMonitor(pp,PERFORMANCE_PRECONDITIONER_INIT);
161 ksteube 1707 blocktimer_increment("Paso_Solver_setPreconditioner()", blocktimer_precond);
162 phornby 1628 if (! Paso_noError()) return;
163 ksteube 1312 time_iter=Paso_timer();
164 gross 584 /* get an initial guess by evaluating the preconditioner */
165 gross 1556 Paso_Solver_solvePreconditioner(A,x,b);
166 gross 584 /* start the iteration process :*/
167 ksteube 1312 r=TMPMEMALLOC(numEqua,double);
168     Paso_checkPtr(r);
169     if (Paso_noError()) {
170     totIter = 0;
171     finalizeIteration = FALSE;
172     last_norm2_of_residual=norm2_of_b;
173     last_norm_max_of_residual=norm_max_of_b;
174 artak 1787 /* Loop */
175 ksteube 1312 while (! finalizeIteration) {
176     cntIter = options->iter_max - totIter;
177     finalizeIteration = TRUE;
178     /* Set initial residual. */
179     norm2_of_residual = 0;
180     norm_max_of_residual = 0;
181 gross 1556 #pragma omp parallel for private(i) schedule(static)
182     for (i = 0; i < numEqua; i++) r[i]=b[i];
183     Paso_SystemMatrix_MatrixVector_CSR_OFFSET0(DBLE(-1), A, x, DBLE(1), r);
184    
185 ksteube 1312 #pragma omp parallel private(norm2_of_residual_local,norm_max_of_residual_local)
186     {
187     norm2_of_residual_local = 0;
188     norm_max_of_residual_local = 0;
189     #pragma omp for private(i) schedule(static)
190     for (i = 0; i < numEqua; i++) {
191     norm2_of_residual_local+= r[i] * r[i];
192     norm_max_of_residual_local=MAX(ABS(scaling[i]*r[i]),norm_max_of_residual_local);
193     }
194     #pragma omp critical
195     {
196     norm2_of_residual += norm2_of_residual_local;
197     norm_max_of_residual = MAX(norm_max_of_residual_local,norm_max_of_residual);
198     }
199 gross 584 }
200 ksteube 1312 /* TODO: use one call */
201     #ifdef PASO_MPI
202     loc_norm = norm2_of_residual;
203     MPI_Allreduce(&loc_norm,&norm2_of_residual, 1, MPI_DOUBLE, MPI_SUM, A->mpi_info->comm);
204     loc_norm = norm_max_of_residual;
205     MPI_Allreduce(&loc_norm,&norm_max_of_residual, 1, MPI_DOUBLE, MPI_MAX, A->mpi_info->comm);
206     #endif
207     norm2_of_residual =sqrt(norm2_of_residual);
208    
209 phornby 1628 if (options->verbose) printf("Step %5d: l2/lmax-norm of residual is %e/%e",totIter,norm2_of_residual,norm_max_of_residual);
210     if (totIter>0 && norm2_of_residual>=last_norm2_of_residual && norm_max_of_residual>=last_norm_max_of_residual) {
211     if (options->verbose) printf(" divergence!\n");
212     Paso_setError(WARNING, "No improvement during iteration. Iterative solver gives up.");
213     } else {
214     /* */
215     if (norm2_of_residual>tolerance*norm2_of_b || norm_max_of_residual>tolerance*norm_max_of_b ) {
216     tol=tolerance*MIN(norm2_of_b,0.1*norm2_of_residual/norm_max_of_residual*norm_max_of_b);
217     if (options->verbose) printf(" (new tolerance = %e).\n",tol);
218     last_norm2_of_residual=norm2_of_residual;
219     last_norm_max_of_residual=norm_max_of_residual;
220     /* call the solver */
221     switch (method) {
222     case PASO_BICGSTAB:
223     errorCode = Paso_Solver_BiCGStab(A, r, x, &cntIter, &tol, pp);
224     break;
225     case PASO_PCG:
226     errorCode = Paso_Solver_PCG(A, r, x, &cntIter, &tol, pp);
227     break;
228 artak 1703 case PASO_TFQMR:
229 artak 1787 errorCode = Paso_Solver_TFQMR(A, r, x, &cntIter, &tol, pp);
230 artak 1703 break;
231 artak 1787 case PASO_MINRES:
232     errorCode = Paso_Solver_MINRES(A, r, x, &cntIter, &tolerance, &tol, pp);
233     break;
234 phornby 1628 case PASO_PRES20:
235     errorCode = Paso_Solver_GMRES(A, r, x, &cntIter, &tol,5,20, pp);
236     break;
237     case PASO_GMRES:
238     errorCode = Paso_Solver_GMRES(A, r, x, &cntIter, &tol,options->truncation,options->restart, pp);
239     break;
240     }
241     totIter += cntIter;
242     /* error handling */
243     if (errorCode==NO_ERROR) {
244     finalizeIteration = FALSE;
245     } else if (errorCode==SOLVER_MAXITER_REACHED) {
246     Paso_setError(WARNING,"maximum number of iteration step reached.\nReturned solution does not fulfil stopping criterion.");
247     if (options->verbose) printf("Maximum number of iterations reached.!\n");
248     } else if (errorCode == SOLVER_INPUT_ERROR ) {
249     Paso_setError(SYSTEM_ERROR,"illegal dimension in iterative solver.");
250     if (options->verbose) printf("Internal error!\n");
251     } else if ( errorCode == SOLVER_BREAKDOWN ) {
252     if (cntIter <= 1) {
253     Paso_setError(ZERO_DIVISION_ERROR, "fatal break down in iterative solver.");
254     if (options->verbose) printf("Uncurable break down!\n");
255     } else {
256     if (options->verbose) printf("Breakdown at iter %d (residual = %e). Restarting ...\n", cntIter+totIter, tol);
257     finalizeIteration = FALSE;
258 ksteube 1312 }
259     } else {
260     if (options->verbose) printf(". convergence! \n");
261     }
262 jgs 150 }
263 ksteube 1312 } /* while */
264     }
265 gross 584 MEMFREE(r);
266     time_iter=Paso_timer()-time_iter;
267     if (options->verbose) {
268 gross 1639 printf("\ntiming: solver: %.4e sec\n",time_iter);
269 gross 584 if (totIter>0) printf("timing: per iteration step: %.4e sec\n",time_iter/totIter);
270 jgs 150 }
271     }
272     }
273 gross 584 }
274 phornby 1628 }
275     Performance_stopMonitor(pp,PERFORMANCE_ALL);
276     blocktimer_increment("Paso_Solver()", blocktimer_start);
277 jgs 150 }

Properties

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

  ViewVC Help
Powered by ViewVC 1.1.26