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

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

Parent Directory Parent Directory | Revision Log Revision Log


Revision 700 - (hide annotations)
Thu Apr 6 00:13:40 2006 UTC (13 years, 5 months ago) by gross
File MIME type: text/plain
File size: 11607 byte(s)
A few changes in the build mechanism and the file structure so scons can build release tar files:

  * paso/src/Solver has been moved to paso/src 
  * all test_.py are now run_.py files and are assumed to be passing python tests. they can run by 
    scons py_tests and are part of the release test set
  * escript/py_src/test_ are moved to escript/test/python and are installed in to the build directory 
    (rather then the PYTHONPATH).
  * all py files in test/python which don't start with run_ or test_ are now 'local_py_tests'. they are installed i
    by not run automatically.
  * CppUnitTest is now treated as a escript module (against previous decisions).
  * scons realse builds nor tar/zip files with relvant source code (src and tests in seperate files)

the python tests don't pass yet due to path problems.


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 gross 700 #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