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

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

Parent Directory Parent Directory | Revision Log Revision Log


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

Properties

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

  ViewVC Help
Powered by ViewVC 1.1.26