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

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

Parent Directory Parent Directory | Revision Log Revision Log


Revision 2670 - (show annotations)
Thu Sep 17 06:17:57 2009 UTC (10 years, 1 month ago) by artak
File MIME type: text/plain
File size: 8503 byte(s)
Bug that made MPI tests to fail is fixed.
1
2 /*******************************************************
3 *
4 * Copyright (c) 2003-2009 by University of Queensland
5 * Earth Systems Science Computational Center (ESSCC)
6 * http://www.uq.edu.au/esscc
7 *
8 * Primary Business: Queensland, Australia
9 * Licensed under the Open Software License version 3.0
10 * http://www.opensource.org/licenses/osl-3.0.php
11 *
12 *******************************************************/
13
14
15 /**************************************************************/
16
17 /* Paso: SystemMatrix: sets-up the preconditioner */
18
19 /**************************************************************/
20
21 /* Copyrights by ACcESS Australia 2003/04 */
22 /* Author: Lutz Gross, l.gross@uq.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_Preconditioner_free(Paso_Solver_Preconditioner* in) {
35 if (in!=NULL) {
36 Paso_Solver_ILU_free(in->ilu);
37 Paso_Solver_RILU_free(in->rilu);
38 Paso_Solver_Jacobi_free(in->jacobi);
39 Paso_Solver_GS_free(in->gs);
40 Paso_Solver_AMG_free(in->amg);
41 Paso_Solver_AMG_System_free(in->amgSystem);
42 MEMFREE(in);
43 }
44 }
45 /* call the iterative solver: */
46
47 void Paso_Solver_setPreconditioner(Paso_SystemMatrix* A,Paso_Options* options) {
48
49 Paso_Solver_Preconditioner* prec=NULL;
50 dim_t i;
51 if (A->solver==NULL) {
52 /* allocate structure to hold preconditioner */
53 prec=MEMALLOC(1,Paso_Solver_Preconditioner);
54 prec->amgSystem=MEMALLOC(1,Paso_Solver_AMG_System);
55 if (Paso_checkPtr(prec)) return;
56 if (Paso_checkPtr(prec->amgSystem)) return;
57 prec->type=UNKNOWN;
58 prec->rilu=NULL;
59 prec->ilu=NULL;
60 prec->jacobi=NULL;
61 prec->gs=NULL;
62 prec->amg=NULL;
63
64 prec->amgSystem->block_size=A->row_block_size;
65 for (i=0;i<A->row_block_size;++i) {
66 prec->amgSystem->amgblock[i]=NULL;
67 prec->amgSystem->block[i]=NULL;
68 }
69
70
71 A->solver=prec;
72 switch (options->preconditioner) {
73 default:
74 case PASO_JACOBI:
75 if (options->verbose) printf("Jacobi preconditioner is used.\n");
76 prec->jacobi=Paso_Solver_getJacobi(A->mainBlock);
77 prec->type=PASO_JACOBI;
78 break;
79 case PASO_ILU0:
80 if (options->verbose) printf("ILU preconditioner is used.\n");
81 prec->ilu=Paso_Solver_getILU(A->mainBlock,options->verbose);
82 prec->type=PASO_ILU0;
83 break;
84 case PASO_RILU:
85 if (options->verbose) printf("RILU preconditioner is used.\n");
86 prec->rilu=Paso_Solver_getRILU(A->mainBlock,options->verbose);
87 prec->type=PASO_RILU;
88 break;
89 case PASO_GS:
90 if (options->verbose) printf("Gauss-Seidel preconditioner is used.\n");
91 prec->gs=Paso_Solver_getGS(A->mainBlock,options->verbose);
92 prec->gs->sweeps=options->sweeps;
93 prec->type=PASO_GS;
94 break;
95 case PASO_AMG:
96 if (options->verbose) printf("AMG preconditioner is used.\n");
97
98 /*For performace reasons we check if block_size is one. If yes, then we do not need to separate blocks.*/
99 if (A->row_block_size==1) {
100 prec->amg=Paso_Solver_getAMG(A->mainBlock,options->level_max,options);
101 }
102 else {
103 for (i=0;i<A->row_block_size;++i) {
104 prec->amgSystem->block[i]=Paso_SparseMatrix_getBlock(A->mainBlock,i+1);
105 prec->amgSystem->amgblock[i]=Paso_Solver_getAMG(prec->amgSystem->block[i],options->level_max,options);
106 }
107 }
108 prec->type=PASO_AMG;
109 break;
110
111 }
112 if (! Paso_MPIInfo_noError(A->mpi_info ) ){
113 Paso_Preconditioner_free(prec);
114 A->solver=NULL;
115 }
116 }
117 }
118
119 /* applies the preconditioner */
120 /* has to be called within a parallel reqion */
121 /* barrier synchronization is performed before the evaluation to make sure that the input vector is available */
122 void Paso_Solver_solvePreconditioner(Paso_SystemMatrix* A,double* x,double* b){
123 Paso_Solver_Preconditioner* prec=(Paso_Solver_Preconditioner*) A->solver;
124 dim_t i,j;
125 dim_t n=A->mainBlock->numRows;
126 double **xx;
127 double **bb;
128 xx=MEMALLOC(A->row_block_size,double*);
129 bb=MEMALLOC(A->row_block_size,double*);
130 if (Paso_checkPtr(xx) && Paso_checkPtr(bb)) return;
131
132 switch (prec->type) {
133 default:
134 case PASO_JACOBI:
135 Paso_Solver_solveJacobi(prec->jacobi,x,b);
136 break;
137 case PASO_ILU0:
138 Paso_Solver_solveILU(prec->ilu,x,b);
139 break;
140 case PASO_RILU:
141 Paso_Solver_solveRILU(prec->rilu,x,b);
142 break;
143 case PASO_GS:
144 /* Gauss-Seidel preconditioner P=U^{-1}DL^{-1} is used here with sweeps paramenter.
145 We want to solve x_new=x_old+P^{-1}(b-Ax_old). So for fisrt 3 we will have the following:
146 x_0=0
147 x_1=P^{-1}(b)
148
149 b_old=b
150
151 b_new=b_old+b
152 x_2=x_1+P^{-1}(b-Ax_1)=P^{-1}(b)+P^{-1}(b)-P^{-1}AP^{-1}(b))
153 =P^{-1}(2b-AP^{-1}b)=P^{-1}(b_new-AP^{-1}b_old)
154 b_old=b_new
155
156 b_new=b_old+b
157 x_3=....=P^{-1}(b_new-AP^{-1}b_old)
158 b_old=b_new
159
160 So for n- steps we will use loop, but we always make sure that every value calculated only once!
161 */
162
163 Paso_Solver_solveGS(prec->gs,x,b);
164 if (prec->gs->sweeps>1) {
165 double *bold=MEMALLOC(prec->gs->n*prec->gs->n_block,double);
166 double *bnew=MEMALLOC(prec->gs->n*prec->gs->n_block,double);
167
168 #pragma omp parallel for private(i) schedule(static)
169 for (i=0;i<prec->gs->n*prec->gs->n_block;++i) bold[i]=b[i];
170
171 while(prec->gs->sweeps>1) {
172 #pragma omp parallel for private(i) schedule(static)
173 for (i=0;i<prec->gs->n*prec->gs->n_block;++i) bnew[i]=bold[i]+b[i];
174 /* Compute the residual b=b-Ax*/
175 Paso_SystemMatrix_MatrixVector_CSR_OFFSET0(DBLE(-1), A, x, DBLE(1), bnew);
176 /* Go round again*/
177 Paso_Solver_solveGS(prec->gs,x,bnew);
178 #pragma omp parallel for private(i) schedule(static)
179 for (i=0;i<prec->gs->n*prec->gs->n_block;++i) bold[i]=bnew[i];
180 prec->gs->sweeps=prec->gs->sweeps-1;
181 }
182 MEMFREE(bold);
183 MEMFREE(bnew);
184 }
185 break;
186 case PASO_AMG:
187
188 /*For performace reasons we check if block_size is one. If yes, then we do not need to do unnecessary copying.*/
189 if (A->row_block_size==1) {
190 Paso_Solver_solveAMG(prec->amg,x,b);
191 }
192 else {
193
194
195 for (i=0;i<A->row_block_size;i++) {
196 xx[i]=MEMALLOC(n,double);
197 bb[i]=MEMALLOC(n,double);
198 if (Paso_checkPtr(xx[i]) && Paso_checkPtr(bb[i])) return;
199 }
200
201 #pragma omp parallel for private(i,j) schedule(static)
202 for (i=0;i<n;i++) {
203 for (j=0;j<A->row_block_size;j++) {
204 bb[j][i]=b[A->row_block_size*i+j];
205 }
206 }
207
208
209 for (i=0;i<A->row_block_size;i++) {
210 Paso_Solver_solveAMG(prec->amgSystem->amgblock[i],xx[i],bb[i]);
211 }
212
213 #pragma omp parallel for private(i,j) schedule(static)
214 for (i=0;i<n;i++) {
215 for (j=0;j<A->row_block_size;j++) {
216 x[A->row_block_size*i+j]=xx[j][i];
217 }
218 }
219
220 for (i=0;i<A->row_block_size;i++) {
221 MEMFREE(xx[i]);
222 MEMFREE(bb[i]);
223 }
224
225 }
226 break;
227 }
228 MEMFREE(xx);
229 MEMFREE(bb);
230 }

Properties

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

  ViewVC Help
Powered by ViewVC 1.1.26