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

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

Parent Directory Parent Directory | Revision Log Revision Log


Revision 682 - (hide annotations)
Mon Mar 27 02:43:09 2006 UTC (13 years ago) by robwdcock
Original Path: trunk/paso/src/Solvers/Solver.c
File MIME type: text/plain
File size: 11613 byte(s)
+ NEW BUILD SYSTEM

This commit contains the new build system with cross-platform support.
Most things work are before though you can have more control.

ENVIRONMENT settings have changed:
+ You no longer require LD_LIBRARY_PATH or PYTHONPATH to point to the
esysroot for building and testing performed via scons
+ ACcESS altix users: It is recommended you change your modules to load
the latest intel compiler and other libraries required by boost to match
the setup in svn (you can override). The correct modules are as follows

module load intel_cc.9.0.026
export
MODULEPATH=${MODULEPATH}:/data/raid2/toolspp4/modulefiles/gcc-3.3.6
module load boost/1.33.0/python-2.4.1
module load python/2.4.1
module load numarray/1.3.3


1 jgs 150 /* $Id$ */
2    
3 dhawcroft 631
4     /*
5     ********************************************************************************
6 dhawcroft 633 * Copyright 2006 by ACcESS MNRF *
7 dhawcroft 631 * *
8     * http://www.access.edu.au *
9     * Primary Business: Queensland, Australia *
10     * Licensed under the Open Software License version 3.0 *
11     * http://www.opensource.org/licenses/osl-3.0.php *
12     ********************************************************************************
13     */
14    
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 robwdcock 682 #include "../Paso.h"
27     #include "../SystemMatrix.h"
28 jgs 150 #include "Solver.h"
29    
30     /***********************************************************************************/
31    
32     /* free space */
33    
34     void Paso_Solver_free(Paso_SystemMatrix* A) {
35 gross 425 Paso_Preconditioner_free(A->solver);
36     A->solver=NULL;
37 jgs 150 }
38     /* call the iterative solver: */
39    
40 gross 584 void Paso_Solver(Paso_SystemMatrix* A,double* x,double* b,Paso_Options* options,Paso_Performance* pp) {
41     double norm2_of_b,tol,tolerance,time_iter,*r=NULL,norm2_of_residual,last_norm2_of_residual,norm_max_of_b,
42 jgs 150 norm2_of_b_local,norm_max_of_b_local,norm2_of_residual_local,norm_max_of_residual_local,norm_max_of_residual,
43     last_norm_max_of_residual,*scaling;
44 gross 415 dim_t i,totIter,cntIter,method;
45     bool_t finalizeIteration;
46     err_t errorCode;
47     dim_t n_col = A->num_cols * A-> col_block_size;
48     dim_t n_row = A->num_rows * A-> row_block_size;
49    
50     tolerance=MAX(options->tolerance,EPSILON);
51     Paso_resetError();
52     method=Paso_Options_getSolver(options->method,PASO_PASO,options->symmetric);
53     /* check matrix type */
54     if ((A->type & MATRIX_FORMAT_CSC) || (A->type & MATRIX_FORMAT_OFFSET1) || (A->type & MATRIX_FORMAT_SYM) ) {
55     Paso_setError(TYPE_ERROR,"Iterative solver requires CSR format with unsymmetric storage scheme and index offset 0.");
56 jgs 150 }
57 gross 415 if (A->col_block_size != A->row_block_size) {
58     Paso_setError(TYPE_ERROR,"Iterative solver requires row and column block sizes to be equal.");
59     }
60 jgs 150 if (A->num_cols != A->num_rows) {
61 gross 415 Paso_setError(TYPE_ERROR,"Iterative solver requires requires a square matrix.");
62     return;
63 jgs 150 }
64 gross 584 Performance_startMonitor(pp,PERFORMANCE_ALL);
65 jgs 150 if (Paso_noError()) {
66     /* get normalization */
67     scaling=Paso_SystemMatrix_borrowNormalization(A);
68 gross 584 if (Paso_noError()) {
69     /* get the norm of the right hand side */
70     norm2_of_b=0.;
71     norm_max_of_b=0.;
72     #pragma omp parallel private(norm2_of_b_local,norm_max_of_b_local)
73     {
74 jgs 150 norm2_of_b_local=0.;
75     norm_max_of_b_local=0.;
76     #pragma omp for private(i) schedule(static)
77     for (i = 0; i < n_row ; ++i) {
78     norm2_of_b_local += b[i] * b[i];
79     norm_max_of_b_local = MAX(ABS(scaling[i]*b[i]),norm_max_of_b_local);
80     }
81     #pragma omp critical
82     {
83     norm2_of_b += norm2_of_b_local;
84     norm_max_of_b = MAX(norm_max_of_b_local,norm_max_of_b);
85     }
86 gross 584 }
87     norm2_of_b=sqrt(norm2_of_b);
88 jgs 150
89 gross 584 /* if norm2_of_b==0 we are ready: x=0 */
90     if (norm2_of_b <=0.) {
91     #pragma omp parallel for private(i) schedule(static)
92     for (i = 0; i < n_col; i++) x[i]=0.;
93     if (options->verbose) printf("right hand side is identical zero.\n");
94     } else {
95     if (options->verbose) {
96     printf("l2/lmax-norm of right hand side is %e/%e.\n",norm2_of_b,norm_max_of_b);
97     printf("l2/lmax-stopping criterion is %e/%e.\n",norm2_of_b*tolerance,norm_max_of_b*tolerance);
98     switch (method) {
99     case PASO_BICGSTAB:
100     printf("Iterative method is BiCGStab.\n");
101     break;
102     case PASO_PCG:
103     printf("Iterative method is PCG.\n");
104     break;
105     case PASO_PRES20:
106     printf("Iterative method is PRES20.\n");
107     break;
108     case PASO_GMRES:
109     if (options->restart>0) {
110     printf("Iterative method is GMRES(%d,%d)\n",options->truncation,options->restart);
111     } else {
112     printf("Iterative method is GMRES(%d)\n",options->truncation);
113     }
114     break;
115     }
116     }
117     /* construct the preconditioner */
118    
119     Performance_startMonitor(pp,PERFORMANCE_PRECONDITIONER_INIT);
120     Paso_Solver_setPreconditioner(A,options);
121     Performance_stopMonitor(pp,PERFORMANCE_PRECONDITIONER_INIT);
122     if (! Paso_noError()) return;
123    
124     /* get an initial guess by evaluating the preconditioner */
125     time_iter=Paso_timer();
126     #pragma omp parallel
127     Paso_Solver_solvePreconditioner(A,x,b);
128     if (! Paso_noError()) return;
129     /* start the iteration process :*/
130     r=TMPMEMALLOC(n_row,double);
131     if (Paso_checkPtr(r)) return;
132     totIter = 0;
133     finalizeIteration = FALSE;
134     last_norm2_of_residual=norm2_of_b;
135     last_norm_max_of_residual=norm_max_of_b;
136     while (! finalizeIteration) {
137     cntIter = options->iter_max - totIter;
138     finalizeIteration = TRUE;
139     /* Set initial residual. */
140     norm2_of_residual = 0;
141     norm_max_of_residual = 0;
142     #pragma omp parallel private(norm2_of_residual_local,norm_max_of_residual_local)
143     {
144     #pragma omp for private(i) schedule(static)
145     for (i = 0; i < n_row; i++) r[i]=b[i];
146    
147     Paso_SystemMatrix_MatrixVector_CSR_OFFSET0(DBLE(-1), A, x, DBLE(1), r);
148    
149     norm2_of_residual_local = 0;
150     norm_max_of_residual_local = 0;
151     #pragma omp for private(i) schedule(static)
152     for (i = 0; i < n_row; i++) {
153     norm2_of_residual_local+= r[i] * r[i];
154     norm_max_of_residual_local=MAX(ABS(scaling[i]*r[i]),norm_max_of_residual_local);
155     }
156     #pragma omp critical
157     {
158     norm2_of_residual += norm2_of_residual_local;
159     norm_max_of_residual = MAX(norm_max_of_residual_local,norm_max_of_residual);
160     }
161     }
162     norm2_of_residual =sqrt(norm2_of_residual);
163    
164     if (options->verbose) printf("Step %5d: l2/lmax-norm of residual is %e/%e",totIter,norm2_of_residual,norm_max_of_residual);
165     if (totIter>0 && norm2_of_residual>=last_norm2_of_residual && norm_max_of_residual>=last_norm_max_of_residual) {
166     if (options->verbose) printf(" divergence!\n");
167     Paso_setError(WARNING, "No improvement during iteration. Iterative solver gives up.");
168     } else {
169     /* */
170     if (norm2_of_residual>tolerance*norm2_of_b || norm_max_of_residual>tolerance*norm_max_of_b ) {
171     tol=tolerance*MIN(norm2_of_b,0.1*norm2_of_residual/norm_max_of_residual*norm_max_of_b);
172     if (options->verbose) printf(" (new tolerance = %e).\n",tol);
173     last_norm2_of_residual=norm2_of_residual;
174     last_norm_max_of_residual=norm_max_of_residual;
175     /* call the solver */
176     switch (method) {
177     case PASO_BICGSTAB:
178     errorCode = Paso_Solver_BiCGStab(A, r, x, &cntIter, &tol,pp);
179     break;
180     case PASO_PCG:
181     errorCode = Paso_Solver_PCG(A, r, x, &cntIter, &tol,pp);
182     break;
183     case PASO_PRES20:
184     errorCode = Paso_Solver_GMRES(A, r, x, &cntIter, &tol,5,20,pp);
185     break;
186     case PASO_GMRES:
187     errorCode = Paso_Solver_GMRES(A, r, x, &cntIter, &tol,options->truncation,options->restart,pp);
188     break;
189     }
190     totIter += cntIter;
191    
192     /* error handling */
193     if (errorCode==NO_ERROR) {
194     finalizeIteration = FALSE;
195     } else if (errorCode==SOLVER_MAXITER_REACHED) {
196     Paso_setError(WARNING,"maximum number of iteration step reached.\nReturned solution does not fulfil stopping criterion.");
197     if (options->verbose) printf("Maximum number of iterations reached.!\n");
198     } else if (errorCode == SOLVER_INPUT_ERROR ) {
199     Paso_setError(SYSTEM_ERROR,"illegal dimension in iterative solver.");
200     if (options->verbose) printf("Internal error!\n");
201     } else if ( errorCode == SOLVER_BREAKDOWN ) {
202     if (cntIter <= 1) {
203     Paso_setError(ZERO_DIVISION_ERROR, "fatal break down in iterative solver.");
204     if (options->verbose) printf("Uncurable break down!\n");
205     } else {
206     if (options->verbose) printf("Breakdown at iter %d (residual = %e). Restarting ...\n", cntIter+totIter, tol);
207     finalizeIteration = FALSE;
208     }
209     } else if (errorCode == SOLVER_MEMORY_ERROR) {
210     Paso_setError(MEMORY_ERROR,"memory allocation failed.");
211     if (options->verbose) printf("Memory allocation failed!\n");
212     } else if (errorCode !=SOLVER_NO_ERROR ) {
213     Paso_setError(SYSTEM_ERROR,"unidentified error in iterative solver.");
214     if (options->verbose) printf("Unidentified error!\n");
215     }
216 jgs 150 } else {
217 gross 584 if (options->verbose) printf(". convergence! \n");
218 jgs 150 }
219 gross 584 }
220     } /* while */
221     MEMFREE(r);
222     time_iter=Paso_timer()-time_iter;
223     if (options->verbose) {
224     printf("timing: solver: %.4e sec\n",time_iter);
225     if (totIter>0) printf("timing: per iteration step: %.4e sec\n",time_iter/totIter);
226 jgs 150 }
227     }
228     }
229 gross 584 }
230     Performance_stopMonitor(pp,PERFORMANCE_ALL);
231 jgs 150 }

Properties

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

  ViewVC Help
Powered by ViewVC 1.1.26