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

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

Parent Directory Parent Directory | Revision Log Revision Log


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

Properties

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

  ViewVC Help
Powered by ViewVC 1.1.26