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

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

Parent Directory Parent Directory | Revision Log Revision Log


Revision 3159 - (hide annotations)
Mon Sep 6 06:59:31 2010 UTC (9 years ago) by gross
File MIME type: text/plain
File size: 7522 byte(s)
AMG is using now a simpler interface to the Smoothers.
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 3158 Paso_Preconditioner_Smoother_free(in->jacobi);
38     Paso_Preconditioner_Smoother_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 gross 3158 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 jgs 150 prec->type=PASO_JACOBI;
77 gross 3158 prec->sweeps=options->sweeps;
78 jgs 150 break;
79 gross 3094 case PASO_GS:
80 gross 3158 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 gross 3094 break;
85     /***************************************************************************************/
86 jgs 150 case PASO_ILU0:
87     if (options->verbose) printf("ILU preconditioner is used.\n");
88 ksteube 1312 prec->ilu=Paso_Solver_getILU(A->mainBlock,options->verbose);
89 jgs 150 prec->type=PASO_ILU0;
90 gross 3094 Paso_MPIInfo_noError(A->mpi_info);
91 jgs 150 break;
92 gross 430 case PASO_RILU:
93     if (options->verbose) printf("RILU preconditioner is used.\n");
94 ksteube 1312 prec->rilu=Paso_Solver_getRILU(A->mainBlock,options->verbose);
95 gross 3094 Paso_MPIInfo_noError(A->mpi_info);
96 gross 430 prec->type=PASO_RILU;
97     break;
98 gross 3094
99 artak 1842 case PASO_AMG:
100 artak 3013 if (options->verbose) printf("AMG preconditioner is used.\n");
101 gross 3094 prec->amg=Paso_Solver_getAMG(A->mainBlock,options->level_max,options);
102     Paso_MPIInfo_noError(A->mpi_info);
103 artak 1842 prec->type=PASO_AMG;
104     break;
105 artak 2992
106 artak 2760 case PASO_AMLI:
107 artak 3017
108     prec->amliSystem=MEMALLOC(1,Paso_Solver_AMLI_System);
109 gross 3120 if (! Paso_checkPtr(prec->amliSystem)) {
110 artak 3017
111     prec->amliSystem->block_size=A->row_block_size;
112    
113 gross 3005 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 artak 2760
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 gross 3120 }
132 artak 2760 break;
133 artak 1842
134 jgs 150 }
135     }
136 gross 3120 if (! Paso_MPIInfo_noError(A->mpi_info ) ){
137     Paso_Preconditioner_free(prec);
138     return NULL;
139     } else {
140     return prec;
141     }
142 jgs 150 }
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 gross 3120 void Paso_Preconditioner_solve(Paso_Solver_Preconditioner* prec, Paso_SystemMatrix* A,double* x,double* b){
148 artak 2662 dim_t i,j;
149 artak 2670 dim_t n=A->mainBlock->numRows;
150 jgs 150 switch (prec->type) {
151     default:
152     case PASO_JACOBI:
153 gross 3159 Paso_Preconditioner_Smoother_solve(A, prec->jacobi,x,b,prec->sweeps, FALSE);
154 gross 3125 break;
155 gross 3094 case PASO_GS:
156 gross 3159 Paso_Preconditioner_Smoother_solve(A, prec->gs,x,b,prec->sweeps, FALSE);
157 gross 3094 break;
158 jgs 150 case PASO_ILU0:
159 gross 3094 Paso_Solver_solveILU(A->mainBlock, prec->ilu,x,b);
160 jgs 150 break;
161 gross 430 case PASO_RILU:
162 gross 3094 Paso_Solver_solveRILU(prec->rilu,x,b);
163 gross 430 break;
164 gross 3094
165 artak 2662 case PASO_AMG:
166 artak 3013 Paso_Solver_solveAMG(prec->amg,x,b);
167 gross 3140 break;
168 artak 2760 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 gross 3094 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 artak 2760 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 gross 3094 MEMFREE(xx);
210     MEMFREE(bb);
211 artak 2760 }
212     break;
213 jgs 150 }
214 gross 3094
215 jgs 150 }

Properties

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

  ViewVC Help
Powered by ViewVC 1.1.26