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

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

Parent Directory Parent Directory | Revision Log Revision Log


Revision 150 - (hide annotations)
Thu Sep 15 03:44:45 2005 UTC (14 years, 6 months ago) by jgs
Original Path: trunk/esys2/paso/src/Solvers/Solver.c
File MIME type: text/plain
File size: 10060 byte(s)
Merge of development branch dev-02 back to main trunk on 2005-09-15

1 jgs 150 /* $Id$ */
2    
3     /**************************************************************/
4    
5     /* Paso: SystemMatrix: controls iterative linear system solvers */
6    
7     /**************************************************************/
8    
9     /* Copyrights by ACcESS Australia 2003/04 */
10     /* Author: gross@access.edu.au */
11    
12     /**************************************************************/
13    
14     #include "Paso.h"
15     #include "SystemMatrix.h"
16     #include "Solver.h"
17    
18     #if PTR_OFFSET !=0 || INDEX_OFFSET!=0
19     #error Paso library usage requires PTR_OFFSET=0 and INDEX_OFFSET=0
20     #endif
21    
22    
23     /***********************************************************************************/
24    
25     /* free space */
26    
27     void Paso_Solver_free(Paso_SystemMatrix* A) {
28     Paso_Preconditioner_free(A->iterative);
29     A->iterative=NULL;
30     }
31     /* call the iterative solver: */
32    
33     void Paso_Solver(Paso_SystemMatrix* A,double* x,double* b,Paso_Options* options) {
34     double norm2_of_b,tol,tolerance,time_iter,time_prec,*r=NULL,norm2_of_residual,last_norm2_of_residual,norm_max_of_b,
35     norm2_of_b_local,norm_max_of_b_local,norm2_of_residual_local,norm_max_of_residual_local,norm_max_of_residual,
36     last_norm_max_of_residual,*scaling;
37     dim_t i,totIter,cntIter,method;
38     bool_t finalizeIteration;
39     err_t errorCode;
40     dim_t n_col = A->num_cols * A-> col_block_size;
41     dim_t n_row = A->num_rows * A-> row_block_size;
42    
43     tolerance=MAX(options->tolerance,EPSILON);
44     Paso_resetError();
45     method=Paso_Options_getSolver(options->method,PASO_PASO,options->symmetric);
46     /* check matrix type */
47     if (A->type!=CSR) {
48     Paso_setError(TYPE_ERROR,"Iterative solver requires CSR format.");
49     }
50     if (A->col_block_size != A->row_block_size) {
51     Paso_setError(TYPE_ERROR,"Iterative solver requires row and column block sizes to be equal.");
52     }
53     if (A->num_cols != A->num_rows) {
54     Paso_setError(TYPE_ERROR,"Iterative solver requires requires a square matrix.");
55     return;
56     }
57     if (Paso_noError()) {
58     /* get normalization */
59     scaling=Paso_SystemMatrix_borrowNormalization(A);
60     if (! Paso_noError()) return;
61     /* get the norm of the right hand side */
62     norm2_of_b=0.;
63     norm_max_of_b=0.;
64     #pragma omp parallel private(norm2_of_b_local,norm_max_of_b_local)
65     {
66     norm2_of_b_local=0.;
67     norm_max_of_b_local=0.;
68     #pragma omp for private(i) schedule(static)
69     for (i = 0; i < n_row ; ++i) {
70     norm2_of_b_local += b[i] * b[i];
71     norm_max_of_b_local = MAX(ABS(scaling[i]*b[i]),norm_max_of_b_local);
72     }
73     #pragma omp critical
74     {
75     norm2_of_b += norm2_of_b_local;
76     norm_max_of_b = MAX(norm_max_of_b_local,norm_max_of_b);
77     }
78     }
79     norm2_of_b=sqrt(norm2_of_b);
80    
81     /* if norm2_of_b==0 we are ready: x=0 */
82     if (norm2_of_b <=0.) {
83     #pragma omp parallel for private(i) schedule(static)
84     for (i = 0; i < n_col; i++) x[i]=0.;
85     if (options->verbose) printf("right hand side is identical zero.\n");
86     return;
87     }
88     if (options->verbose) {
89     printf("l2/lmax-norm of right hand side is %e/%e.\n",norm2_of_b,norm_max_of_b);
90     printf("l2/lmax-stopping criterion is %e/%e.\n",norm2_of_b*tolerance,norm_max_of_b*tolerance);
91     switch (method) {
92     case PASO_BICGSTAB:
93     printf("Iterative method is BiCGStab.\n");
94     break;
95     case PASO_PCG:
96     printf("Iterative method is PCG.\n");
97     break;
98     case PASO_PRES20:
99     printf("Iterative method is PRES20.\n");
100     break;
101     case PASO_GMRES:
102     if (options->restart>0) {
103     printf("Iterative method is GMRES(%d,%d)\n",options->truncation,options->restart);
104     } else {
105     printf("Iterative method is GMRES(%d)\n",options->truncation);
106     }
107     break;
108     }
109     }
110     /* construct the preconditioner */
111    
112     time_prec=Paso_timer();
113     Paso_Solver_setPreconditioner(A,options);
114     time_prec=Paso_timer()-time_prec;
115     if (! Paso_noError()) return;
116    
117     /* get an initial guess by evaluating the preconditioner */
118     time_iter=Paso_timer();
119     #pragma omp parallel
120     Paso_Solver_solvePreconditioner(A,x,b);
121     if (! Paso_noError()) return;
122     /* start the iteration process :*/
123     r=TMPMEMALLOC(n_row,double);
124     if (Paso_checkPtr(r)) return;
125     totIter = 0;
126     finalizeIteration = FALSE;
127     last_norm2_of_residual=norm2_of_b;
128     last_norm_max_of_residual=norm_max_of_b;
129     while (! finalizeIteration) {
130     cntIter = options->iter_max - totIter;
131     finalizeIteration = TRUE;
132     /* Set initial residual. */
133     norm2_of_residual = 0;
134     norm_max_of_residual = 0;
135     #pragma omp parallel private(norm2_of_residual_local,norm_max_of_residual_local)
136     {
137     #pragma omp for private(i) schedule(static)
138     for (i = 0; i < n_row; i++) r[i]=b[i];
139    
140     Paso_SystemMatrix_MatrixVector(DBLE(-1), A, x, DBLE(1), r);
141    
142     norm2_of_residual_local = 0;
143     norm_max_of_residual_local = 0;
144     #pragma omp for private(i) schedule(static)
145     for (i = 0; i < n_row; i++) {
146     norm2_of_residual_local+= r[i] * r[i];
147     norm_max_of_residual_local=MAX(ABS(scaling[i]*r[i]),norm_max_of_residual_local);
148     }
149     #pragma omp critical
150     {
151     norm2_of_residual += norm2_of_residual_local;
152     norm_max_of_residual = MAX(norm_max_of_residual_local,norm_max_of_residual);
153     }
154     }
155     norm2_of_residual =sqrt(norm2_of_residual);
156    
157     if (options->verbose) printf("Step %5d: l2/lmax-norm of residual is %e/%e",totIter,norm2_of_residual,norm_max_of_residual);
158     if (totIter>0 && norm2_of_residual>=last_norm2_of_residual && norm_max_of_residual>=last_norm_max_of_residual) {
159     if (options->verbose) printf(" divergence!\n");
160     Paso_setError(WARNING, "No improvement during iteration. Iterative solver gives up.");
161     } else {
162     /* */
163     if (norm2_of_residual>tolerance*norm2_of_b || norm_max_of_residual>tolerance*norm_max_of_b ) {
164     tol=tolerance*MIN(norm2_of_b,0.1*norm2_of_residual/norm_max_of_residual*norm_max_of_b);
165     if (options->verbose) printf(" (new tolerance = %e).\n",tol);
166     last_norm2_of_residual=norm2_of_residual;
167     last_norm_max_of_residual=norm_max_of_residual;
168     /* call the solver */
169     switch (method) {
170     case PASO_BICGSTAB:
171     errorCode = Paso_Solver_BiCGStab(A, r, x, &cntIter, &tol);
172     break;
173     case PASO_PCG:
174     errorCode = Paso_Solver_PCG(A, r, x, &cntIter, &tol);
175     break;
176     case PASO_PRES20:
177     errorCode = Paso_Solver_GMRES(A, r, x, &cntIter, &tol,5,20);
178     break;
179     case PASO_GMRES:
180     errorCode = Paso_Solver_GMRES(A, r, x, &cntIter, &tol,options->truncation,options->restart);
181     break;
182     }
183     totIter += cntIter;
184    
185     /* error handling */
186     if (errorCode==NO_ERROR) {
187     finalizeIteration = FALSE;
188     } else if (errorCode==SOLVER_MAXITER_REACHED) {
189     Paso_setError(WARNING,"maximum number of iteration step reached.\nReturned solution does not fulfil stopping criterion.");
190     } else if (errorCode == SOLVER_INPUT_ERROR ) {
191     Paso_setError(SYSTEM_ERROR,"illegal dimension in iterative solver.");
192     } else if ( errorCode == SOLVER_BREAKDOWN ) {
193     if (cntIter <= 1) {
194     Paso_setError(ZERO_DIVISION_ERROR, "fatal break down in iterative solver.");
195     } else {
196     if (options->verbose) printf("Breakdown at iter %d (residual = %e). Restarting ...\n", cntIter+totIter, tol);
197     finalizeIteration = FALSE;
198     }
199     } else if (errorCode == SOLVER_MEMORY_ERROR) {
200     Paso_setError(MEMORY_ERROR,"memory allocation failed.");
201     } else if (errorCode !=SOLVER_NO_ERROR ) {
202     Paso_setError(SYSTEM_ERROR,"unidentified error in iterative solver.");
203     }
204     } else {
205     if (options->verbose) printf(". convergence! \n");
206     }
207     }
208     } /* while */
209     MEMFREE(r);
210     time_iter=Paso_timer()-time_iter;
211     if (options->verbose) {
212     printf("timing: solver: %.4e sec\n",time_prec+time_iter);
213     printf("timing: preconditioner: %.4e sec\n",time_prec);
214     if (totIter>0) printf("timing: per iteration step: %.4e sec\n",time_iter/totIter);
215     }
216     }
217     }
218    
219     /*
220     * $Log$
221     * Revision 1.2 2005/09/15 03:44:40 jgs
222     * Merge of development branch dev-02 back to main trunk on 2005-09-15
223     *
224     * Revision 1.1.2.2 2005/09/07 00:59:09 gross
225     * some inconsistent renaming fixed to make the linking work.
226     *
227     * Revision 1.1.2.1 2005/09/05 06:29:49 gross
228     * These files have been extracted from finley to define a stand alone libray for iterative
229     * linear solvers on the ALTIX. main entry through Paso_solve. this version compiles but
230     * has not been tested yet.
231     *
232     *
233     */

Properties

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

  ViewVC Help
Powered by ViewVC 1.1.26