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

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

Parent Directory Parent Directory | Revision Log Revision Log


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

Properties

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

  ViewVC Help
Powered by ViewVC 1.1.26