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

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

Parent Directory Parent Directory | Revision Log Revision Log


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

Properties

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

  ViewVC Help
Powered by ViewVC 1.1.26