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

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

Parent Directory Parent Directory | Revision Log Revision Log


Revision 2705 - (hide annotations)
Thu Oct 1 23:21:10 2009 UTC (10 years ago) by artak
File MIME type: text/plain
File size: 8511 byte(s)
Some openMP codes in AMG systems are removed for testing.
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 artak 2670 if (Paso_checkPtr(prec->amgSystem)) return;
57 jgs 150 prec->type=UNKNOWN;
58 gross 430 prec->rilu=NULL;
59 jgs 150 prec->ilu=NULL;
60     prec->jacobi=NULL;
61 artak 1819 prec->gs=NULL;
62 artak 1842 prec->amg=NULL;
63 artak 2662
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 gross 425 A->solver=prec;
72 jgs 150 switch (options->preconditioner) {
73     default:
74     case PASO_JACOBI:
75     if (options->verbose) printf("Jacobi preconditioner is used.\n");
76 ksteube 1312 prec->jacobi=Paso_Solver_getJacobi(A->mainBlock);
77 jgs 150 prec->type=PASO_JACOBI;
78     break;
79     case PASO_ILU0:
80     if (options->verbose) printf("ILU preconditioner is used.\n");
81 ksteube 1312 prec->ilu=Paso_Solver_getILU(A->mainBlock,options->verbose);
82 jgs 150 prec->type=PASO_ILU0;
83     break;
84 gross 430 case PASO_RILU:
85     if (options->verbose) printf("RILU preconditioner is used.\n");
86 ksteube 1312 prec->rilu=Paso_Solver_getRILU(A->mainBlock,options->verbose);
87 gross 430 prec->type=PASO_RILU;
88     break;
89 artak 1819 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 artak 1823 prec->gs->sweeps=options->sweeps;
93 artak 1819 prec->type=PASO_GS;
94     break;
95 artak 1842 case PASO_AMG:
96     if (options->verbose) printf("AMG preconditioner is used.\n");
97 artak 2659
98 artak 2662 /*For performace reasons we check if block_size is one. If yes, then we do not need to separate blocks.*/
99 artak 2659 if (A->row_block_size==1) {
100     prec->amg=Paso_Solver_getAMG(A->mainBlock,options->level_max,options);
101     }
102 artak 2662 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 artak 2659 }
108 artak 1842 prec->type=PASO_AMG;
109     break;
110    
111 jgs 150 }
112 ksteube 1312 if (! Paso_MPIInfo_noError(A->mpi_info ) ){
113 artak 2659 Paso_Preconditioner_free(prec);
114     A->solver=NULL;
115 jgs 150 }
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 gross 425 Paso_Solver_Preconditioner* prec=(Paso_Solver_Preconditioner*) A->solver;
124 artak 2662 dim_t i,j;
125 artak 2670 dim_t n=A->mainBlock->numRows;
126 artak 2662 double **xx;
127     double **bb;
128 artak 2670 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 artak 2662
132 jgs 150 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 gross 430 case PASO_RILU:
141     Paso_Solver_solveRILU(prec->rilu,x,b);
142     break;
143 artak 2662 case PASO_GS:
144 artak 1842 /* 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 artak 1819 Paso_Solver_solveGS(prec->gs,x,b);
164 artak 1842 if (prec->gs->sweeps>1) {
165 artak 2662 double *bold=MEMALLOC(prec->gs->n*prec->gs->n_block,double);
166     double *bnew=MEMALLOC(prec->gs->n*prec->gs->n_block,double);
167 artak 1842
168 artak 2662 #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 artak 1842
171 artak 2662 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 artak 1842 }
185 artak 1819 break;
186 artak 2662 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 artak 2659 if (A->row_block_size==1) {
190     Paso_Solver_solveAMG(prec->amg,x,b);
191     }
192 artak 2662 else {
193    
194 artak 2659
195 artak 2670 for (i=0;i<A->row_block_size;i++) {
196 artak 2662 xx[i]=MEMALLOC(n,double);
197     bb[i]=MEMALLOC(n,double);
198 artak 2670 if (Paso_checkPtr(xx[i]) && Paso_checkPtr(bb[i])) return;
199 artak 2659 }
200 artak 2670
201 artak 2705 /*#pragma omp parallel for private(i,j) schedule(static)*/
202 artak 2662 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 artak 2660 }
207 artak 2670
208    
209 artak 2662 for (i=0;i<A->row_block_size;i++) {
210     Paso_Solver_solveAMG(prec->amgSystem->amgblock[i],xx[i],bb[i]);
211 artak 2659 }
212 artak 2670
213 artak 2705 /*#pragma omp parallel for private(i,j) schedule(static)*/
214 artak 2662 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 artak 2670
220 artak 2662 for (i=0;i<A->row_block_size;i++) {
221     MEMFREE(xx[i]);
222     MEMFREE(bb[i]);
223     }
224 artak 2670
225 artak 2659 }
226     break;
227 jgs 150 }
228 artak 2670 MEMFREE(xx);
229     MEMFREE(bb);
230 jgs 150 }

Properties

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

  ViewVC Help
Powered by ViewVC 1.1.26