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

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

Parent Directory Parent Directory | Revision Log Revision Log


Revision 155 - (show annotations)
Wed Nov 9 02:02:19 2005 UTC (13 years, 11 months ago) by jgs
File MIME type: text/plain
File size: 9964 byte(s)
move all directories from trunk/esys2 into trunk and remove esys2

1 /* $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 if (options->verbose) printf("Maximum number of iterations reached.!\n");
191 } else if (errorCode == SOLVER_INPUT_ERROR ) {
192 Paso_setError(SYSTEM_ERROR,"illegal dimension in iterative solver.");
193 if (options->verbose) printf("Internal error!\n");
194 } else if ( errorCode == SOLVER_BREAKDOWN ) {
195 if (cntIter <= 1) {
196 Paso_setError(ZERO_DIVISION_ERROR, "fatal break down in iterative solver.");
197 if (options->verbose) printf("Uncurable break down!\n");
198 } else {
199 if (options->verbose) printf("Breakdown at iter %d (residual = %e). Restarting ...\n", cntIter+totIter, tol);
200 finalizeIteration = FALSE;
201 }
202 } else if (errorCode == SOLVER_MEMORY_ERROR) {
203 Paso_setError(MEMORY_ERROR,"memory allocation failed.");
204 if (options->verbose) printf("Memory allocation failed!\n");
205 } else if (errorCode !=SOLVER_NO_ERROR ) {
206 Paso_setError(SYSTEM_ERROR,"unidentified error in iterative solver.");
207 if (options->verbose) printf("Unidentified error!\n");
208 }
209 } else {
210 if (options->verbose) printf(". convergence! \n");
211 }
212 }
213 } /* while */
214 MEMFREE(r);
215 time_iter=Paso_timer()-time_iter;
216 if (options->verbose) {
217 printf("timing: solver: %.4e sec\n",time_prec+time_iter);
218 printf("timing: preconditioner: %.4e sec\n",time_prec);
219 if (totIter>0) printf("timing: per iteration step: %.4e sec\n",time_iter/totIter);
220 }
221 }
222 }

Properties

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

  ViewVC Help
Powered by ViewVC 1.1.26