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

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

Parent Directory Parent Directory | Revision Log Revision Log


Revision 3193 - (show annotations)
Tue Sep 21 06:56:44 2010 UTC (8 years, 11 months ago) by gross
File MIME type: text/plain
File size: 7775 byte(s)
more clean up in AMG
1
2 /*******************************************************
3 *
4 * Copyright (c) 2003-2010 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 "PasoUtil.h"
29 #include "Preconditioner.h"
30
31 /***********************************************************************************/
32
33 /* free space */
34
35 void Paso_Preconditioner_free(Paso_Preconditioner* in) {
36 if (in!=NULL) {
37 Paso_Preconditioner_Smoother_free(in->jacobi);
38 Paso_Preconditioner_Smoother_free(in->gs);
39 Paso_Preconditioner_LocalAMG_free(in->localamg);
40
41 /*********************************/
42 Paso_Solver_ILU_free(in->ilu);
43 Paso_Solver_RILU_free(in->rilu);
44 Paso_Solver_AMLI_free(in->amli);
45 Paso_Solver_AMLI_System_free(in->amliSystem);
46 /*********************************/
47
48 MEMFREE(in);
49 }
50 }
51
52 Paso_Preconditioner* Paso_Preconditioner_alloc(Paso_SystemMatrix* A,Paso_Options* options) {
53
54 Paso_Preconditioner* prec=NULL;
55 dim_t i;
56
57 prec=MEMALLOC(1,Paso_Preconditioner);
58
59 if (! Paso_checkPtr(prec)) {
60
61 prec->type=UNKNOWN;
62
63
64 prec->jacobi=NULL;
65 prec->gs=NULL;
66 prec->localamg=NULL;
67 /*********************************/
68 prec->rilu=NULL;
69 prec->ilu=NULL;
70
71 prec->amli=NULL;
72 prec->amliSystem=NULL;
73 /*********************************/
74
75 if (options->verbose && options->use_local_preconditioner) printf("Apply preconditioner locally only.\n");
76
77 switch (options->preconditioner) {
78 default:
79 case PASO_JACOBI:
80 if (options->verbose) printf("Jacobi(%d) preconditioner is used.\n",options->sweeps);
81 prec->jacobi=Paso_Preconditioner_Smoother_alloc(A, TRUE, options->use_local_preconditioner, options->verbose);
82 prec->type=PASO_JACOBI;
83 prec->sweeps=options->sweeps;
84 break;
85 case PASO_GS:
86 if (options->verbose) printf("Gauss-Seidel(%d) preconditioner is used.\n",options->sweeps);
87 prec->gs=Paso_Preconditioner_Smoother_alloc(A, FALSE, options->use_local_preconditioner, options->verbose);
88 prec->type=PASO_GS;
89 prec->sweeps=options->sweeps;
90 break;
91 case PASO_AMG:
92 if (options->verbose) printf("AMG preconditioner is used.\n");
93 prec->localamg=Paso_Preconditioner_LocalAMG_alloc(A->mainBlock,options->level_max,options);
94 Paso_MPIInfo_noError(A->mpi_info);
95 prec->type=PASO_AMG;
96 break;
97
98 /***************************************************************************************/
99 case PASO_ILU0:
100 if (options->verbose) printf("ILU preconditioner is used.\n");
101 prec->ilu=Paso_Solver_getILU(A->mainBlock,options->verbose);
102 prec->type=PASO_ILU0;
103 Paso_MPIInfo_noError(A->mpi_info);
104 break;
105 case PASO_RILU:
106 if (options->verbose) printf("RILU preconditioner is used.\n");
107 prec->rilu=Paso_Solver_getRILU(A->mainBlock,options->verbose);
108 Paso_MPIInfo_noError(A->mpi_info);
109 prec->type=PASO_RILU;
110 break;
111
112
113
114 case PASO_AMLI:
115
116 prec->amliSystem=MEMALLOC(1,Paso_Solver_AMLI_System);
117 if (! Paso_checkPtr(prec->amliSystem)) {
118
119 prec->amliSystem->block_size=A->row_block_size;
120
121 for (i=0;i<A->row_block_size;++i) {
122 prec->amliSystem->amliblock[i]=NULL;
123 prec->amliSystem->block[i]=NULL;
124 }
125
126 if (options->verbose) printf("AMLI preconditioner is used.\n");
127
128 /*For performace reasons we check if block_size is one. If yes, then we do not need to separate blocks.*/
129 if (A->row_block_size==1) {
130 prec->amli=Paso_Solver_getAMLI(A->mainBlock,options->level_max,options);
131 }
132 else {
133 for (i=0;i<A->row_block_size;++i) {
134 prec->amliSystem->block[i]=Paso_SparseMatrix_getBlock(A->mainBlock,i+1);
135 prec->amliSystem->amliblock[i]=Paso_Solver_getAMLI(prec->amliSystem->block[i],options->level_max,options);
136 }
137 }
138 prec->type=PASO_AMLI;
139 }
140 break;
141
142 }
143 }
144 if (! Paso_MPIInfo_noError(A->mpi_info ) ){
145 Paso_Preconditioner_free(prec);
146 return NULL;
147 } else {
148 return prec;
149 }
150 }
151
152 /* applies the preconditioner */
153 /* has to be called within a parallel reqion */
154 /* barrier synchronization is performed before the evaluation to make sure that the input vector is available */
155 void Paso_Preconditioner_solve(Paso_Preconditioner* prec, Paso_SystemMatrix* A,double* x,double* b){
156 dim_t i,j;
157 dim_t n=A->mainBlock->numRows;
158 switch (prec->type) {
159 default:
160 case PASO_JACOBI:
161 Paso_Preconditioner_Smoother_solve(A, prec->jacobi,x,b,prec->sweeps, FALSE);
162 break;
163 case PASO_GS:
164 Paso_Preconditioner_Smoother_solve(A, prec->gs,x,b,prec->sweeps, FALSE);
165 break;
166 case PASO_AMG:
167 Paso_Preconditioner_LocalAMG_solve(prec->localamg,x,b);
168 break;
169
170 /*=========================================================*/
171 case PASO_ILU0:
172 Paso_Solver_solveILU(A->mainBlock, prec->ilu,x,b);
173 break;
174 case PASO_RILU:
175 Paso_Solver_solveRILU(prec->rilu,x,b);
176 break;
177
178
179 case PASO_AMLI:
180
181 /*For performace reasons we check if block_size is one. If yes, then we do not need to do unnecessary copying.*/
182 if (A->row_block_size==1) {
183 Paso_Solver_solveAMLI(prec->amli,x,b);
184 }
185 else {
186 double **xx;
187 double **bb;
188 xx=MEMALLOC(A->row_block_size,double*);
189 bb=MEMALLOC(A->row_block_size,double*);
190 if (Paso_checkPtr(xx) || Paso_checkPtr(bb)) return;
191 for (i=0;i<A->row_block_size;i++) {
192 xx[i]=MEMALLOC(n,double);
193 bb[i]=MEMALLOC(n,double);
194 if (Paso_checkPtr(xx[i]) && Paso_checkPtr(bb[i])) return;
195 }
196
197 /*#pragma omp parallel for private(i,j) schedule(static)*/
198 for (i=0;i<n;i++) {
199 for (j=0;j<A->row_block_size;j++) {
200 bb[j][i]=b[A->row_block_size*i+j];
201 }
202 }
203
204
205 for (i=0;i<A->row_block_size;i++) {
206 Paso_Solver_solveAMLI(prec->amliSystem->amliblock[i],xx[i],bb[i]);
207 }
208
209 /*#pragma omp parallel for private(i,j) schedule(static)*/
210 for (i=0;i<n;i++) {
211 for (j=0;j<A->row_block_size;j++) {
212 x[A->row_block_size*i+j]=xx[j][i];
213 }
214 }
215
216 for (i=0;i<A->row_block_size;i++) {
217 MEMFREE(xx[i]);
218 MEMFREE(bb[i]);
219 }
220 MEMFREE(xx);
221 MEMFREE(bb);
222 }
223 break;
224 /*=========================================================*/
225 }
226
227 }

Properties

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

  ViewVC Help
Powered by ViewVC 1.1.26