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

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

Parent Directory Parent Directory | Revision Log Revision Log


Revision 3283 - (show annotations)
Mon Oct 18 22:39:28 2010 UTC (10 years, 4 months ago) by gross
File MIME type: text/plain
File size: 5802 byte(s)
AMG reengineered, BUG is Smoother fixed.


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 Paso_Preconditioner_LocalSmoother_free(in->localamgsubstitute);
41 /*********************************/
42 Paso_Solver_ILU_free(in->ilu);
43 Paso_Solver_RILU_free(in->rilu);
44 /*********************************/
45
46 MEMFREE(in);
47 }
48 }
49
50 Paso_Preconditioner* Paso_Preconditioner_alloc(Paso_SystemMatrix* A,Paso_Options* options) {
51
52 Paso_Preconditioner* prec=NULL;
53
54 prec=MEMALLOC(1,Paso_Preconditioner);
55
56 if (! Esys_checkPtr(prec)) {
57
58 prec->type=UNKNOWN;
59
60
61 prec->jacobi=NULL;
62 prec->gs=NULL;
63 prec->localamg=NULL;
64 prec->localamgsubstitute=NULL;
65
66 /*********************************/
67 prec->rilu=NULL;
68 prec->ilu=NULL;
69 /*********************************/
70
71 if (options->verbose && options->use_local_preconditioner) printf("Paso: Apply preconditioner locally only.\n");
72
73 switch (options->preconditioner) {
74 default:
75 case PASO_JACOBI:
76 if (options->verbose) printf("Paso: Jacobi(%d) preconditioner is used.\n",options->sweeps);
77 prec->jacobi=Paso_Preconditioner_Smoother_alloc(A, TRUE, options->use_local_preconditioner, options->verbose);
78 prec->type=PASO_JACOBI;
79 prec->sweeps=options->sweeps;
80 break;
81 case PASO_GS:
82 if (options->verbose) printf("Paso: Gauss-Seidel(%d) preconditioner is used.\n",options->sweeps);
83 prec->gs=Paso_Preconditioner_Smoother_alloc(A, FALSE, options->use_local_preconditioner, options->verbose);
84 prec->type=PASO_GS;
85 prec->sweeps=options->sweeps;
86 break;
87 case PASO_AMLI:
88 case PASO_AMG:
89 if (options->verbose) printf("Paso: AMG preconditioner is used.\n");
90 prec->localamg=Paso_Preconditioner_LocalAMG_alloc(A->mainBlock,1,options);
91 prec->sweeps=options->sweeps;
92 /* if NULL is returned (and no error) no AMG has been constructed because the system is too small or not big enough
93 we now use the Smoother as a preconditioner */
94 if ( (Esys_noError()) && (prec->localamg == NULL) ) {
95 if (options->verbose) {
96 if (options->smoother == PASO_JACOBI) {
97 printf("Paso: Jacobi(%d) preconditioner is used.\n",prec->sweeps);
98 } else {
99 printf("Paso: Gauss-Seidel(%d) preconditioner is used.\n",prec->sweeps);
100 }
101 }
102 prec->localamgsubstitute=Paso_Preconditioner_LocalSmoother_alloc(A->mainBlock, (options->smoother == PASO_JACOBI), options->verbose);
103 } else {
104 if (options->verbose) printf("Paso: AMG preconditioner is used.\n");
105 }
106 prec->type=PASO_AMG;
107 Esys_MPIInfo_noError(A->mpi_info);
108 break;
109
110 /***************************************************************************************/
111 case PASO_ILU0:
112 if (options->verbose) printf("Paso: ILU preconditioner is used.\n");
113 prec->ilu=Paso_Solver_getILU(A->mainBlock,options->verbose);
114 prec->type=PASO_ILU0;
115 Esys_MPIInfo_noError(A->mpi_info);
116 break;
117 case PASO_RILU:
118 if (options->verbose) printf("Paso: RILU preconditioner is used.\n");
119 prec->rilu=Paso_Solver_getRILU(A->mainBlock,options->verbose);
120 Esys_MPIInfo_noError(A->mpi_info);
121 prec->type=PASO_RILU;
122 break;
123 }
124 }
125 if (! Esys_MPIInfo_noError(A->mpi_info ) ){
126 Paso_Preconditioner_free(prec);
127 return NULL;
128 } else {
129 return prec;
130 }
131 }
132
133 /* applies the preconditioner */
134 /* has to be called within a parallel reqion */
135 /* barrier synchronization is performed before the evaluation to make sure that the input vector is available */
136 void Paso_Preconditioner_solve(Paso_Preconditioner* prec, Paso_SystemMatrix* A,double* x,double* b){
137
138 switch (prec->type) {
139 default:
140 case PASO_JACOBI:
141 Paso_Preconditioner_Smoother_solve(A, prec->jacobi,x,b,prec->sweeps, FALSE);
142 break;
143 case PASO_GS:
144 Paso_Preconditioner_Smoother_solve(A, prec->gs,x,b,prec->sweeps, FALSE);
145 break;
146 case PASO_AMG:
147 if (prec->localamg == NULL) {
148 Paso_Preconditioner_LocalSmoother_solve(A->mainBlock, prec->localamgsubstitute,x,b,prec->sweeps, FALSE);
149 } else {
150 Paso_Preconditioner_LocalAMG_solve(A->mainBlock, prec->localamg,x,b);
151 }
152 break;
153
154 /*=========================================================*/
155 case PASO_ILU0:
156 Paso_Solver_solveILU(A->mainBlock, prec->ilu,x,b);
157 break;
158 case PASO_RILU:
159 Paso_Solver_solveRILU(prec->rilu,x,b);
160 break;
161
162 }
163
164 }

Properties

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

  ViewVC Help
Powered by ViewVC 1.1.26