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

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

Parent Directory Parent Directory | Revision Log Revision Log


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

Properties

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

  ViewVC Help
Powered by ViewVC 1.1.26