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

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

Parent Directory Parent Directory | Revision Log Revision Log


Revision 682 - (show 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 /* $Id$ */
2
3
4 /*
5 ********************************************************************************
6 * Copyright 2006 by ACcESS MNRF *
7 * *
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 /**************************************************************/
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 #include "../Paso.h"
27 #include "../SystemMatrix.h"
28 #include "Solver.h"
29
30 /***********************************************************************************/
31
32 /* free space */
33
34 void Paso_Solver_free(Paso_SystemMatrix* A) {
35 Paso_Preconditioner_free(A->solver);
36 A->solver=NULL;
37 }
38 /* call the iterative solver: */
39
40 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 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 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 }
57 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 if (A->num_cols != A->num_rows) {
61 Paso_setError(TYPE_ERROR,"Iterative solver requires requires a square matrix.");
62 return;
63 }
64 Performance_startMonitor(pp,PERFORMANCE_ALL);
65 if (Paso_noError()) {
66 /* get normalization */
67 scaling=Paso_SystemMatrix_borrowNormalization(A);
68 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 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 }
87 norm2_of_b=sqrt(norm2_of_b);
88
89 /* 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 } else {
217 if (options->verbose) printf(". convergence! \n");
218 }
219 }
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 }
227 }
228 }
229 }
230 Performance_stopMonitor(pp,PERFORMANCE_ALL);
231 }

Properties

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

  ViewVC Help
Powered by ViewVC 1.1.26