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 |
prec->type=UNKNOWN; |
57 |
prec->rilu=NULL; |
58 |
prec->ilu=NULL; |
59 |
prec->jacobi=NULL; |
60 |
prec->gs=NULL; |
61 |
prec->amg=NULL; |
62 |
|
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 |
A->solver=prec; |
71 |
switch (options->preconditioner) { |
72 |
default: |
73 |
case PASO_JACOBI: |
74 |
if (options->verbose) printf("Jacobi preconditioner is used.\n"); |
75 |
prec->jacobi=Paso_Solver_getJacobi(A->mainBlock); |
76 |
prec->type=PASO_JACOBI; |
77 |
break; |
78 |
case PASO_ILU0: |
79 |
if (options->verbose) printf("ILU preconditioner is used.\n"); |
80 |
prec->ilu=Paso_Solver_getILU(A->mainBlock,options->verbose); |
81 |
prec->type=PASO_ILU0; |
82 |
break; |
83 |
case PASO_RILU: |
84 |
if (options->verbose) printf("RILU preconditioner is used.\n"); |
85 |
prec->rilu=Paso_Solver_getRILU(A->mainBlock,options->verbose); |
86 |
prec->type=PASO_RILU; |
87 |
break; |
88 |
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 |
prec->gs->sweeps=options->sweeps; |
92 |
prec->type=PASO_GS; |
93 |
break; |
94 |
case PASO_AMG: |
95 |
if (options->verbose) printf("AMG preconditioner is used.\n"); |
96 |
|
97 |
/*For performace reasons we check if block_size is one. If yes, then we do not need to separate blocks.*/ |
98 |
if (A->row_block_size==1) { |
99 |
prec->amg=Paso_Solver_getAMG(A->mainBlock,options->level_max,options); |
100 |
} |
101 |
else { |
102 |
for (i=0;i<A->row_block_size;++i) { |
103 |
prec->amgSystem->block[i]=Paso_SparseMatrix_getBlock(A->mainBlock,i+1); |
104 |
prec->amgSystem->amgblock[i]=Paso_Solver_getAMG(prec->amgSystem->block[i],options->level_max,options); |
105 |
} |
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=Paso_SystemMatrix_getGlobalNumRows(A); |
126 |
double **xx; |
127 |
double **bb; |
128 |
|
129 |
switch (prec->type) { |
130 |
default: |
131 |
case PASO_JACOBI: |
132 |
Paso_Solver_solveJacobi(prec->jacobi,x,b); |
133 |
break; |
134 |
case PASO_ILU0: |
135 |
Paso_Solver_solveILU(prec->ilu,x,b); |
136 |
break; |
137 |
case PASO_RILU: |
138 |
Paso_Solver_solveRILU(prec->rilu,x,b); |
139 |
break; |
140 |
case PASO_GS: |
141 |
/* Gauss-Seidel preconditioner P=U^{-1}DL^{-1} is used here with sweeps paramenter. |
142 |
We want to solve x_new=x_old+P^{-1}(b-Ax_old). So for fisrt 3 we will have the following: |
143 |
x_0=0 |
144 |
x_1=P^{-1}(b) |
145 |
|
146 |
b_old=b |
147 |
|
148 |
b_new=b_old+b |
149 |
x_2=x_1+P^{-1}(b-Ax_1)=P^{-1}(b)+P^{-1}(b)-P^{-1}AP^{-1}(b)) |
150 |
=P^{-1}(2b-AP^{-1}b)=P^{-1}(b_new-AP^{-1}b_old) |
151 |
b_old=b_new |
152 |
|
153 |
b_new=b_old+b |
154 |
x_3=....=P^{-1}(b_new-AP^{-1}b_old) |
155 |
b_old=b_new |
156 |
|
157 |
So for n- steps we will use loop, but we always make sure that every value calculated only once! |
158 |
*/ |
159 |
|
160 |
Paso_Solver_solveGS(prec->gs,x,b); |
161 |
if (prec->gs->sweeps>1) { |
162 |
double *bold=MEMALLOC(prec->gs->n*prec->gs->n_block,double); |
163 |
double *bnew=MEMALLOC(prec->gs->n*prec->gs->n_block,double); |
164 |
|
165 |
#pragma omp parallel for private(i) schedule(static) |
166 |
for (i=0;i<prec->gs->n*prec->gs->n_block;++i) bold[i]=b[i]; |
167 |
|
168 |
while(prec->gs->sweeps>1) { |
169 |
#pragma omp parallel for private(i) schedule(static) |
170 |
for (i=0;i<prec->gs->n*prec->gs->n_block;++i) bnew[i]=bold[i]+b[i]; |
171 |
/* Compute the residual b=b-Ax*/ |
172 |
Paso_SystemMatrix_MatrixVector_CSR_OFFSET0(DBLE(-1), A, x, DBLE(1), bnew); |
173 |
/* Go round again*/ |
174 |
Paso_Solver_solveGS(prec->gs,x,bnew); |
175 |
#pragma omp parallel for private(i) schedule(static) |
176 |
for (i=0;i<prec->gs->n*prec->gs->n_block;++i) bold[i]=bnew[i]; |
177 |
prec->gs->sweeps=prec->gs->sweeps-1; |
178 |
} |
179 |
MEMFREE(bold); |
180 |
MEMFREE(bnew); |
181 |
} |
182 |
break; |
183 |
case PASO_AMG: |
184 |
|
185 |
/*For performace reasons we check if block_size is one. If yes, then we do not need to do unnecessary copying.*/ |
186 |
if (A->row_block_size==1) { |
187 |
Paso_Solver_solveAMG(prec->amg,x,b); |
188 |
} |
189 |
else { |
190 |
|
191 |
xx=MEMALLOC(A->row_block_size,double*); |
192 |
bb=MEMALLOC(A->row_block_size,double*); |
193 |
|
194 |
for (i=0;i<A->row_block_size;i++) { |
195 |
xx[i]=MEMALLOC(n,double); |
196 |
bb[i]=MEMALLOC(n,double); |
197 |
} |
198 |
|
199 |
#pragma omp parallel for private(i,j) schedule(static) |
200 |
for (i=0;i<n;i++) { |
201 |
for (j=0;j<A->row_block_size;j++) { |
202 |
bb[j][i]=b[A->row_block_size*i+j]; |
203 |
} |
204 |
} |
205 |
|
206 |
for (i=0;i<A->row_block_size;i++) { |
207 |
Paso_Solver_solveAMG(prec->amgSystem->amgblock[i],xx[i],bb[i]); |
208 |
} |
209 |
|
210 |
#pragma omp parallel for private(i,j) schedule(static) |
211 |
for (i=0;i<n;i++) { |
212 |
for (j=0;j<A->row_block_size;j++) { |
213 |
x[A->row_block_size*i+j]=xx[j][i]; |
214 |
} |
215 |
} |
216 |
|
217 |
for (i=0;i<A->row_block_size;i++) { |
218 |
MEMFREE(xx[i]); |
219 |
MEMFREE(bb[i]); |
220 |
} |
221 |
|
222 |
MEMFREE(xx); |
223 |
MEMFREE(bb); |
224 |
} |
225 |
|
226 |
break; |
227 |
} |
228 |
|
229 |
} |