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

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

Parent Directory Parent Directory | Revision Log Revision Log


Revision 2240 - (hide annotations)
Wed Feb 4 01:17:06 2009 UTC (12 years, 2 months ago) by artak
File MIME type: text/plain
File size: 13915 byte(s)
Some statistical print outs are added in solver and loadMM to Paso_SparseMatrix method is added
1 ksteube 1312
2     /*******************************************************
3 ksteube 1811 *
4     * Copyright (c) 2003-2008 by University of Queensland
5     * Earth Systems Science Computational Center (ESSCC)
6     * http://www.uq.edu.au/esscc
7     *
8     * Primary Business: Queensland, Australia
9     * Licensed under the Open Software License version 3.0
10     * http://www.opensource.org/licenses/osl-3.0.php
11     *
12     *******************************************************/
13 dhawcroft 631
14 ksteube 1811
15 jgs 150 /**************************************************************/
16    
17     /* Paso: SystemMatrix: controls iterative linear system solvers */
18    
19     /**************************************************************/
20    
21     /* Copyrights by ACcESS Australia 2003/04 */
22     /* Author: gross@access.edu.au */
23    
24     /**************************************************************/
25    
26 gross 700 #include "Paso.h"
27     #include "SystemMatrix.h"
28 jgs 150 #include "Solver.h"
29 phornby 2078 #include "esysUtils/blocktimer.h"
30 jgs 150
31     /***********************************************************************************/
32    
33     /* free space */
34    
35     void Paso_Solver_free(Paso_SystemMatrix* A) {
36 phornby 1628 Paso_Preconditioner_free(A->solver);
37     A->solver=NULL;
38 jgs 150 }
39     /* call the iterative solver: */
40    
41 phornby 1628 void Paso_Solver(Paso_SystemMatrix* A,double* x,double* b,
42     Paso_Options* options,Paso_Performance* pp) {
43    
44     double norm2_of_b,tol,tolerance,time_iter;
45     double *r=NULL,norm2_of_residual,last_norm2_of_residual,norm_max_of_b;
46     double norm2_of_b_local,norm_max_of_b_local,norm2_of_residual_local;
47     double norm_max_of_residual_local,norm_max_of_residual;
48     double last_norm_max_of_residual,*scaling;
49     #ifdef PASO_MPI
50     double loc_norm;
51     #endif
52     dim_t i,totIter,cntIter,method;
53     bool_t finalizeIteration;
54 jfenwick 1981 err_t errorCode=NO_ERROR;
55 phornby 1628 dim_t numSol = Paso_SystemMatrix_getTotalNumCols(A);
56     dim_t numEqua = Paso_SystemMatrix_getTotalNumRows(A);
57 ksteube 1707 double blocktimer_precond, blocktimer_start = blocktimer_time();
58 gross 415
59 gross 1639 Paso_resetError();
60 gross 2156 tolerance=options->tolerance;
61     if (tolerance < 100.* EPSILON) {
62     Paso_setError(VALUE_ERROR,"Paso_Solver: Tolerance is too small.");
63     }
64     if (tolerance >1.) {
65     Paso_setError(VALUE_ERROR,"Paso_Solver: Tolerance mut be less than one.");
66     }
67 gross 1639 method=Paso_Options_getSolver(options->method,PASO_PASO,options->symmetric);
68     /* check matrix type */
69     if ((A->type & MATRIX_FORMAT_CSC) || (A->type & MATRIX_FORMAT_OFFSET1) || (A->type & MATRIX_FORMAT_SYM) ) {
70 gross 2156 Paso_setError(TYPE_ERROR,"Paso_Solver: Iterative solver requires CSR format with unsymmetric storage scheme and index offset 0.");
71 gross 1639 }
72     if (A->col_block_size != A->row_block_size) {
73 gross 2156 Paso_setError(TYPE_ERROR,"Paso_Solver: Iterative solver requires row and column block sizes to be equal.");
74 gross 1639 }
75     if (Paso_SystemMatrix_getGlobalNumCols(A) != Paso_SystemMatrix_getGlobalNumRows(A)) {
76 gross 2156 Paso_setError(TYPE_ERROR,"Paso_Solver: Iterative solver requires a square matrix.");
77 gross 1639 return;
78     }
79     /* this for testing only */
80     if (method==PASO_NONLINEAR_GMRES) {
81     Paso_Function* F=NULL;
82     F=Paso_Function_LinearSystem_alloc(A,b,options);
83 gross 1798 Paso_Solver_solvePreconditioner(A,x,b);
84 gross 1639 errorCode=Paso_Solver_NewtonGMRES(F,x,options,pp);
85 gross 1798 if (errorCode!=NO_ERROR) {
86 gross 1639 Paso_setError(SYSTEM_ERROR,"Paso_Solver_NewtonGMRES: an error has occured.");
87     }
88     Paso_Function_LinearSystem_free(F);
89 gross 1798 return;
90 gross 1639 }
91     /* ========================= */
92     Performance_startMonitor(pp,PERFORMANCE_ALL);
93     if (Paso_noError()) {
94     /* get normalization */
95     scaling=Paso_SystemMatrix_borrowNormalization(A);
96     if (Paso_noError()) {
97     /* get the norm of the right hand side */
98     norm2_of_b=0.;
99     norm_max_of_b=0.;
100     #pragma omp parallel private(norm2_of_b_local,norm_max_of_b_local)
101     {
102     norm2_of_b_local=0.;
103     norm_max_of_b_local=0.;
104     #pragma omp for private(i) schedule(static)
105     for (i = 0; i < numEqua ; ++i) {
106     norm2_of_b_local += b[i] * b[i];
107     norm_max_of_b_local = MAX(ABS(scaling[i]*b[i]),norm_max_of_b_local);
108     }
109     #pragma omp critical
110     {
111     norm2_of_b += norm2_of_b_local;
112     norm_max_of_b = MAX(norm_max_of_b_local,norm_max_of_b);
113     }
114     }
115     /* TODO: use one call */
116     #ifdef PASO_MPI
117     {
118     loc_norm = norm2_of_b;
119     MPI_Allreduce(&loc_norm,&norm2_of_b, 1, MPI_DOUBLE, MPI_SUM, A->mpi_info->comm);
120     loc_norm = norm_max_of_b;
121     MPI_Allreduce(&loc_norm,&norm_max_of_b, 1, MPI_DOUBLE, MPI_MAX, A->mpi_info->comm);
122     }
123     #endif
124     norm2_of_b=sqrt(norm2_of_b);
125 jgs 150
126 phornby 1628 /* if norm2_of_b==0 we are ready: x=0 */
127     if (norm2_of_b <=0.) {
128 gross 1798 #pragma omp parallel for private(i) schedule(static)
129 phornby 1628 for (i = 0; i < numSol; i++) x[i]=0.;
130     if (options->verbose) printf("right hand side is identical zero.\n");
131     } else {
132     if (options->verbose) {
133 gross 2156 printf("Paso_Solver: l2/lmax-norm of right hand side is %e/%e.\n",norm2_of_b,norm_max_of_b);
134     printf("Paso_Solver: l2/lmax-stopping criterion is %e/%e.\n",norm2_of_b*tolerance,norm_max_of_b*tolerance);
135 phornby 1628 switch (method) {
136     case PASO_BICGSTAB:
137 gross 2156 printf("Paso_Solver: Iterative method is BiCGStab.\n");
138 phornby 1628 break;
139     case PASO_PCG:
140 gross 2156 printf("Paso_Solver: Iterative method is PCG.\n");
141 phornby 1628 break;
142 artak 1703 case PASO_TFQMR:
143 gross 2156 printf("Paso_Solver: Iterative method is TFQMR.\n");
144 artak 1703 break;
145 artak 1787 case PASO_MINRES:
146 gross 2156 printf("Paso_Solver: Iterative method is MINRES.\n");
147 artak 1787 break;
148 phornby 1628 case PASO_PRES20:
149 gross 2156 printf("Paso_Solver: Iterative method is PRES20.\n");
150 phornby 1628 break;
151     case PASO_GMRES:
152     if (options->restart>0) {
153 gross 2156 printf("Paso_Solver: Iterative method is GMRES(%d,%d)\n",options->truncation,options->restart);
154 phornby 1628 } else {
155 gross 2156 printf("Paso_Solver: Iterative method is GMRES(%d)\n",options->truncation);
156 gross 584 }
157 phornby 1628 break;
158     }
159     }
160     /* construct the preconditioner */
161 ksteube 1707 blocktimer_precond = blocktimer_time();
162 phornby 1628 Performance_startMonitor(pp,PERFORMANCE_PRECONDITIONER_INIT);
163     Paso_Solver_setPreconditioner(A,options);
164     Performance_stopMonitor(pp,PERFORMANCE_PRECONDITIONER_INIT);
165 ksteube 1707 blocktimer_increment("Paso_Solver_setPreconditioner()", blocktimer_precond);
166 phornby 1628 if (! Paso_noError()) return;
167 ksteube 1312 time_iter=Paso_timer();
168 gross 584 /* get an initial guess by evaluating the preconditioner */
169 gross 1556 Paso_Solver_solvePreconditioner(A,x,b);
170 gross 584 /* start the iteration process :*/
171 ksteube 1312 r=TMPMEMALLOC(numEqua,double);
172     Paso_checkPtr(r);
173     if (Paso_noError()) {
174     totIter = 0;
175     finalizeIteration = FALSE;
176     last_norm2_of_residual=norm2_of_b;
177     last_norm_max_of_residual=norm_max_of_b;
178 artak 1787 /* Loop */
179 ksteube 1312 while (! finalizeIteration) {
180     cntIter = options->iter_max - totIter;
181     finalizeIteration = TRUE;
182     /* Set initial residual. */
183     norm2_of_residual = 0;
184     norm_max_of_residual = 0;
185 gross 1556 #pragma omp parallel for private(i) schedule(static)
186     for (i = 0; i < numEqua; i++) r[i]=b[i];
187     Paso_SystemMatrix_MatrixVector_CSR_OFFSET0(DBLE(-1), A, x, DBLE(1), r);
188 ksteube 1312 #pragma omp parallel private(norm2_of_residual_local,norm_max_of_residual_local)
189     {
190     norm2_of_residual_local = 0;
191     norm_max_of_residual_local = 0;
192     #pragma omp for private(i) schedule(static)
193     for (i = 0; i < numEqua; i++) {
194     norm2_of_residual_local+= r[i] * r[i];
195     norm_max_of_residual_local=MAX(ABS(scaling[i]*r[i]),norm_max_of_residual_local);
196     }
197     #pragma omp critical
198     {
199     norm2_of_residual += norm2_of_residual_local;
200     norm_max_of_residual = MAX(norm_max_of_residual_local,norm_max_of_residual);
201     }
202 gross 584 }
203 ksteube 1312 /* TODO: use one call */
204     #ifdef PASO_MPI
205     loc_norm = norm2_of_residual;
206     MPI_Allreduce(&loc_norm,&norm2_of_residual, 1, MPI_DOUBLE, MPI_SUM, A->mpi_info->comm);
207     loc_norm = norm_max_of_residual;
208     MPI_Allreduce(&loc_norm,&norm_max_of_residual, 1, MPI_DOUBLE, MPI_MAX, A->mpi_info->comm);
209     #endif
210     norm2_of_residual =sqrt(norm2_of_residual);
211    
212 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);
213 phornby 1628 if (totIter>0 && norm2_of_residual>=last_norm2_of_residual && norm_max_of_residual>=last_norm_max_of_residual) {
214     if (options->verbose) printf(" divergence!\n");
215 gross 2156 Paso_setError(WARNING, "Paso_Solver: No improvement during iteration. Iterative solver gives up.");
216 phornby 1628 } else {
217     /* */
218     if (norm2_of_residual>tolerance*norm2_of_b || norm_max_of_residual>tolerance*norm_max_of_b ) {
219     tol=tolerance*MIN(norm2_of_b,0.1*norm2_of_residual/norm_max_of_residual*norm_max_of_b);
220     if (options->verbose) printf(" (new tolerance = %e).\n",tol);
221     last_norm2_of_residual=norm2_of_residual;
222     last_norm_max_of_residual=norm_max_of_residual;
223     /* call the solver */
224     switch (method) {
225     case PASO_BICGSTAB:
226     errorCode = Paso_Solver_BiCGStab(A, r, x, &cntIter, &tol, pp);
227     break;
228     case PASO_PCG:
229 artak 1819 errorCode = Paso_Solver_PCG(A, r, x, &cntIter, &tol, pp);
230 phornby 1628 break;
231 artak 1703 case PASO_TFQMR:
232 artak 1787 errorCode = Paso_Solver_TFQMR(A, r, x, &cntIter, &tol, pp);
233 artak 1703 break;
234 artak 1787 case PASO_MINRES:
235 artak 1862 errorCode = Paso_Solver_MINRES(A, r, x, &cntIter, &tol, pp);
236 artak 1787 break;
237 phornby 1628 case PASO_PRES20:
238     errorCode = Paso_Solver_GMRES(A, r, x, &cntIter, &tol,5,20, pp);
239     break;
240     case PASO_GMRES:
241     errorCode = Paso_Solver_GMRES(A, r, x, &cntIter, &tol,options->truncation,options->restart, pp);
242     break;
243     }
244     totIter += cntIter;
245     /* error handling */
246     if (errorCode==NO_ERROR) {
247     finalizeIteration = FALSE;
248     } else if (errorCode==SOLVER_MAXITER_REACHED) {
249 gross 2156 Paso_setError(WARNING,"Paso_Solver: maximum number of iteration step reached.\nReturned solution does not fulfil stopping criterion.");
250     if (options->verbose) printf("Paso_Solver: Maximum number of iterations reached.!\n");
251 phornby 1628 } else if (errorCode == SOLVER_INPUT_ERROR ) {
252 gross 2156 Paso_setError(SYSTEM_ERROR,"Paso_Solver: illegal dimension in iterative solver.");
253     if (options->verbose) printf("Paso_Solver: Internal error!\n");
254 phornby 1628 } else if ( errorCode == SOLVER_BREAKDOWN ) {
255     if (cntIter <= 1) {
256 gross 2156 Paso_setError(ZERO_DIVISION_ERROR, "Paso_Solver: fatal break down in iterative solver.");
257     if (options->verbose) printf("Paso_Solver: Uncurable break down!\n");
258 phornby 1628 } else {
259 gross 2156 if (options->verbose) printf("Paso_Solver: Breakdown at iter %d (residual = %e). Restarting ...\n", cntIter+totIter, tol);
260 phornby 1628 finalizeIteration = FALSE;
261 ksteube 1312 }
262     } else {
263     if (options->verbose) printf(". convergence! \n");
264     }
265 jgs 150 }
266 ksteube 1312 } /* while */
267     }
268 gross 584 MEMFREE(r);
269     time_iter=Paso_timer()-time_iter;
270     if (options->verbose) {
271 gross 2156 printf("\ntiming: Paso_Solver: %.4e sec\n",time_iter);
272 artak 2240 if (totIter>0) {
273     if(totIter==options->iter_max) {
274     printf("timing: Total MAX steps, time per iteration step: %.4e sec\n",time_iter/totIter);
275     } else {
276     printf("timing: Total %d steps, time per iteration step: %.4e sec\n",totIter,time_iter/totIter);
277     }
278     }
279     else {
280     printf("timing: Total 0 step, time per iteration step: %.4e sec\n",time_iter);
281     }
282 jgs 150 }
283     }
284     }
285 gross 584 }
286 phornby 1628 }
287     Performance_stopMonitor(pp,PERFORMANCE_ALL);
288     blocktimer_increment("Paso_Solver()", blocktimer_start);
289 jgs 150 }

Properties

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

  ViewVC Help
Powered by ViewVC 1.1.26