/[escript]/trunk/esys2/finley/src/finleyC/Solvers/Solver.c
ViewVC logotype

Contents of /trunk/esys2/finley/src/finleyC/Solvers/Solver.c

Parent Directory Parent Directory | Revision Log Revision Log


Revision 97 - (show annotations)
Tue Dec 14 05:39:33 2004 UTC (15 years, 1 month ago) by jgs
File MIME type: text/plain
File size: 9225 byte(s)
*** empty log message ***

1 /* $Id$ */
2
3 /**************************************************************/
4
5 /* Finley: 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 "Finley.h"
15 #include "System.h"
16 #include "Solver.h"
17 #include "Common.h"
18 #if ITERATIVE_SOLVER == NO_LIB
19 #if PTR_OFFSET !=0 || INDEX_OFFSET!=0
20 #error Finley library usage requires PTR_OFFSET=0 and INDEX_OFFSET=0
21 #endif
22 #endif
23
24
25 /***********************************************************************************/
26
27 /* free space */
28
29 void Finley_Solver_free(Finley_SystemMatrix* A) {
30 #if ITERATIVE_SOLVER == NO_LIB
31 Finley_Preconditioner_free(A->iterative);
32 A->iterative=NULL;
33 #endif
34 }
35 /* call the iterative solver: */
36
37 void Finley_Solver(Finley_SystemMatrix* A,double* x,double* b,Finley_SolverOptions* options) {
38 #if ITERATIVE_SOLVER == NO_LIB
39 double norm2OfB,tol,tolerance,time_iter,time_prec,*r=NULL,norm_of_residual,last_norm_of_residual;
40 int i,totIter,cntIter,finalizeIteration,errorCode,method;
41 int n_col = A->num_cols * A-> col_block_size;
42 int n_row = A->num_rows * A-> row_block_size;
43
44 tolerance=MAX(options->tolerance,EPSILON);
45 Finley_ErrorCode=NO_ERROR;
46 /* check matrix type */
47 if (A->type!=CSR) {
48 Finley_ErrorCode = TYPE_ERROR;
49 sprintf(Finley_ErrorMsg, "Iterative solver requires CSR format.");
50 return;
51 }
52 if (A->col_block_size != A->row_block_size) {
53 Finley_ErrorCode = TYPE_ERROR;
54 sprintf(Finley_ErrorMsg, "preconditioner requires row and column block sizes to be equal.");
55 return;
56 }
57 if (A->num_cols != A->num_rows) {
58 Finley_ErrorCode = TYPE_ERROR;
59 sprintf(Finley_ErrorMsg, "This method requires a square matrix.");
60 return;
61 }
62 /* select iterative method */
63 switch (options->method) {
64 default:
65 if (options->symmetric) {
66 method=ESCRIPT_PCG;
67 } else {
68 method=ESCRIPT_BICGSTAB;
69 }
70 break;
71 case ESCRIPT_BICGSTAB:
72 method=ESCRIPT_BICGSTAB;
73 break;
74 case ESCRIPT_PCG:
75 method=ESCRIPT_PCG;
76 break;
77 case ESCRIPT_PRES20:
78 method=ESCRIPT_PRES20;
79 break;
80 case ESCRIPT_GMRES:
81 method=ESCRIPT_GMRES;
82 break;
83 }
84 if (options->verbose) {
85 switch (method) {
86 case ESCRIPT_BICGSTAB:
87 printf("Iterative method is BiCGStab.\n");
88 break;
89 case ESCRIPT_PCG:
90 printf("Iterative method is PCG.\n");
91 break;
92 case ESCRIPT_PRES20:
93 printf("Iterative method is PRES20.\n");
94 break;
95 case ESCRIPT_GMRES:
96 if (options->restart>0) {
97 printf("Iterative method is GMRES(%d,%d)\n",options->truncation,options->restart);
98 } else {
99 printf("Iterative method is GMRES(%d)\n",options->truncation);
100 }
101 break;
102
103 }
104 }
105 /* get the norm of the right hand side */
106 norm2OfB=0;
107 #pragma omp parallel for private(i) reduction(+:norm2OfB) schedule(static)
108 for (i = 0; i < n_row ; i++) norm2OfB += b[i] * b[i];
109 norm2OfB=sqrt(norm2OfB);
110
111 /* if norm2OfB==0 we are ready: x=0 */
112
113 if (norm2OfB <=0.) {
114 #pragma omp parallel for private(i) schedule(static)
115 for (i = 0; i < n_col; i++) x[i]=0.;
116 if (options->verbose) printf("right hand side is identical zero.\n");
117 return;
118 }
119
120
121 /* construct the preconditioner */
122
123 time_prec=Finley_timer();
124 Finley_Solver_setPreconditioner(A,options);
125 time_prec=Finley_timer()-time_prec;
126 if (Finley_ErrorCode!=NO_ERROR) return;
127
128 /* get an initial guess by evaluating the preconditioner */
129 time_iter=Finley_timer();
130 #pragma omp parallel
131 Finley_Solver_solvePreconditioner(A,x,b);
132 if (Finley_ErrorCode!=NO_ERROR) return;
133 /* start the iteration process :*/
134 r=TMPMEMALLOC(n_row,double);
135 if (Finley_checkPtr(r)) return;
136 totIter = 0;
137 finalizeIteration = FALSE;
138 norm_of_residual=norm2OfB;
139 while (! finalizeIteration) {
140 cntIter = options->iter_max - totIter;
141 finalizeIteration = TRUE;
142 /* Set initial residual. */
143 norm_of_residual = 0;
144 #pragma omp parallel
145 {
146 #pragma omp for private(i) schedule(static)
147 for (i = 0; i < n_row; i++) r[i]=b[i];
148 Finley_RawScaledSystemMatrixVector(DBLE(-1), A, x, DBLE(1), r);
149 #pragma omp for private(i) reduction(+:norm_of_residual) schedule(static)
150 for (i = 0; i < n_row; i++) norm_of_residual+= r[i] * r[i];
151 }
152 norm_of_residual =sqrt(norm_of_residual);
153 if (totIter>0 && norm_of_residual>=last_norm_of_residual) {
154 Finley_ErrorCode=WARNING;
155 sprintf(Finley_ErrorMsg, "No improvement during iteration.\nIterative solver gives up.");
156 } else {
157 /* */
158 if (norm_of_residual>tolerance*norm2OfB) {
159 if (totIter>0) {
160 tol=MAX(tolerance*tol/norm_of_residual,EPSILON*10.)*norm2OfB;
161 if (options->verbose) printf("new stopping criterion is %e.\n",tol);
162 } else {
163 tol=tolerance*norm2OfB;
164 if (options->verbose) printf("Stopping criterion is %e.\n",tol);
165 }
166 last_norm_of_residual=norm_of_residual;
167 /* call the solver */
168 switch (method) {
169 case ESCRIPT_BICGSTAB:
170 errorCode = Finley_Solver_BiCGStab(A, r, x, &cntIter, &tol);
171 break;
172 case ESCRIPT_PCG:
173 errorCode = Finley_Solver_PCG(A, r, x, &cntIter, &tol);
174 break;
175 case ESCRIPT_PRES20:
176 errorCode = Finley_Solver_GMRES(A, r, x, &cntIter, &tol,5,20);
177 break;
178 case ESCRIPT_GMRES:
179 errorCode = Finley_Solver_GMRES(A, r, x, &cntIter, &tol,options->truncation,options->restart);
180 break;
181 }
182 totIter += cntIter;
183
184 /* error handling */
185 if (errorCode==NO_ERROR) {
186 finalizeIteration = FALSE;
187 } else if (errorCode==SOLVER_MAXITER_REACHED) {
188 Finley_ErrorCode=WARNING;
189 sprintf(Finley_ErrorMsg,"maximum number of iteration step reached.\nReturned solution does not fulfil stopping criterion.");
190 } else if (errorCode == SOLVER_INPUT_ERROR ) {
191 Finley_ErrorCode=SYSTEM_ERROR;
192 sprintf(Finley_ErrorMsg,"illegal dimension in iterative solver.");
193 } else if ( errorCode == SOLVER_BREAKDOWN ) {
194 if (cntIter <= 1) {
195 Finley_ErrorCode =ZERO_DIVISION_ERROR;
196 sprintf(Finley_ErrorMsg, "fatal break down in iterative solver.");
197 } else {
198 if (options->verbose) printf("Breakdown at iter %d (residual = %e). Restarting ...\n", cntIter+totIter, tol);
199 finalizeIteration = FALSE;
200 }
201 } else if (errorCode == SOLVER_MEMORY_ERROR) {
202 Finley_ErrorCode=MEMORY_ERROR;
203 sprintf(Finley_ErrorMsg,"memory allocation failed.");
204 } else if (errorCode !=SOLVER_NO_ERROR ) {
205 Finley_ErrorCode=SYSTEM_ERROR;
206 sprintf(Finley_ErrorMsg,"unidentified error in iterative solver.");
207 }
208 }
209 }
210 } /* while */
211 MEMFREE(r);
212 time_iter=Finley_timer()-time_iter;
213 if (options->verbose) {
214 printf("Iterative solver reached tolerance %.5e after %d steps.\n",norm_of_residual,totIter);
215 printf("timing: solver: %.4e sec\n",time_prec+time_iter);
216 printf("timing: preconditioner: %.4e sec\n",time_prec);
217 if (totIter>0) printf("timing: per iteration step: %.4e sec\n",time_iter/totIter);
218 }
219 #else
220 Finley_ErrorCode=SYSTEM_ERROR;
221 sprintf(Finley_ErrorMsg,"No native Finley solver available.");
222 #endif
223 }
224
225 /*
226 * $Log$
227 * Revision 1.2 2004/12/14 05:39:32 jgs
228 * *** empty log message ***
229 *
230 * Revision 1.1.1.1.2.4 2004/12/07 10:12:06 gross
231 * GMRES added
232 *
233 * Revision 1.1.1.1.2.3 2004/11/24 01:37:17 gross
234 * some changes dealing with the integer overflow in memory allocation. Finley solves 4M unknowns now
235 *
236 * Revision 1.1.1.1.2.2 2004/11/15 06:05:27 gross
237 * poisson solver added
238 *
239 * Revision 1.1.1.1.2.1 2004/11/12 06:58:20 gross
240 * a lot of changes to get the linearPDE class running: most important change is that there is no matrix format exposed to the user anymore. the format is chosen by the Domain according to the solver and symmetry
241 *
242 * Revision 1.1.1.1 2004/10/26 06:53:57 jgs
243 * initial import of project esys2
244 *
245 * Revision 1.1 2004/07/02 04:21:14 gross
246 * Finley C code has been included
247 *
248 *
249 */

Properties

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

  ViewVC Help
Powered by ViewVC 1.1.26