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

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

Parent Directory Parent Directory | Revision Log Revision Log


Revision 2475 - (show annotations)
Wed Jun 17 01:48:46 2009 UTC (10 years, 3 months ago) by artak
Original Path: trunk/paso/src/Solver_preconditioner.c
File MIME type: text/plain
File size: 5940 byte(s)
AMG now takes into account new SolverOptions class, namely one can select coarsening algorithms from python. Not all options are included into AMG, this will be implemented later on.
1
2 /*******************************************************
3 *
4 * Copyright (c) 2003-2008 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: gross@access.edu.au */
23
24 /**************************************************************/
25
26 #include "Paso.h"
27 #include "SystemMatrix.h"
28 #include "Solver.h"
29
30 /***********************************************************************************/
31
32 /* free space */
33
34 void Paso_Preconditioner_free(Paso_Solver_Preconditioner* in) {
35 if (in!=NULL) {
36 Paso_Solver_ILU_free(in->ilu);
37 Paso_Solver_RILU_free(in->rilu);
38 Paso_Solver_Jacobi_free(in->jacobi);
39 Paso_Solver_GS_free(in->gs);
40 Paso_Solver_AMG_free(in->amg);
41 MEMFREE(in);
42 }
43 }
44 /* call the iterative solver: */
45
46 void Paso_Solver_setPreconditioner(Paso_SystemMatrix* A,Paso_Options* options) {
47 Paso_Solver_Preconditioner* prec=NULL;
48 if (A->solver==NULL) {
49 /* allocate structure to hold preconditioner */
50 prec=MEMALLOC(1,Paso_Solver_Preconditioner);
51 if (Paso_checkPtr(prec)) return;
52 prec->type=UNKNOWN;
53 prec->rilu=NULL;
54 prec->ilu=NULL;
55 prec->jacobi=NULL;
56 prec->gs=NULL;
57 prec->amg=NULL;
58 A->solver=prec;
59 switch (options->preconditioner) {
60 default:
61 case PASO_JACOBI:
62 if (options->verbose) printf("Jacobi preconditioner is used.\n");
63 prec->jacobi=Paso_Solver_getJacobi(A->mainBlock);
64 prec->type=PASO_JACOBI;
65 break;
66 case PASO_ILU0:
67 if (options->verbose) printf("ILU preconditioner is used.\n");
68 prec->ilu=Paso_Solver_getILU(A->mainBlock,options->verbose);
69 prec->type=PASO_ILU0;
70 break;
71 case PASO_RILU:
72 if (options->verbose) printf("RILU preconditioner is used.\n");
73 prec->rilu=Paso_Solver_getRILU(A->mainBlock,options->verbose);
74 prec->type=PASO_RILU;
75 break;
76 case PASO_GS:
77 if (options->verbose) printf("Gauss-Seidel preconditioner is used.\n");
78 prec->gs=Paso_Solver_getGS(A->mainBlock,options->verbose);
79 prec->gs->sweeps=options->sweeps;
80 prec->type=PASO_GS;
81 break;
82 case PASO_AMG:
83 if (options->verbose) printf("AMG preconditioner is used.\n");
84 prec->amg=Paso_Solver_getAMG(A->mainBlock,options->verbose,options->level_max,options->coarsening_threshold,options->coarsening_method);
85 prec->type=PASO_AMG;
86 break;
87
88 }
89 if (! Paso_MPIInfo_noError(A->mpi_info ) ){
90 Paso_Preconditioner_free(prec);
91 A->solver=NULL;
92 }
93 }
94 }
95
96 /* applies the preconditioner */
97 /* has to be called within a parallel reqion */
98 /* barrier synchronization is performed before the evaluation to make sure that the input vector is available */
99 void Paso_Solver_solvePreconditioner(Paso_SystemMatrix* A,double* x,double* b){
100 Paso_Solver_Preconditioner* prec=(Paso_Solver_Preconditioner*) A->solver;
101
102 switch (prec->type) {
103 default:
104 case PASO_JACOBI:
105 Paso_Solver_solveJacobi(prec->jacobi,x,b);
106 break;
107 case PASO_ILU0:
108 Paso_Solver_solveILU(prec->ilu,x,b);
109 break;
110 case PASO_RILU:
111 Paso_Solver_solveRILU(prec->rilu,x,b);
112 break;
113 case PASO_GS:
114 /* Gauss-Seidel preconditioner P=U^{-1}DL^{-1} is used here with sweeps paramenter.
115 We want to solve x_new=x_old+P^{-1}(b-Ax_old). So for fisrt 3 we will have the following:
116 x_0=0
117 x_1=P^{-1}(b)
118
119 b_old=b
120
121 b_new=b_old+b
122 x_2=x_1+P^{-1}(b-Ax_1)=P^{-1}(b)+P^{-1}(b)-P^{-1}AP^{-1}(b))
123 =P^{-1}(2b-AP^{-1}b)=P^{-1}(b_new-AP^{-1}b_old)
124 b_old=b_new
125
126 b_new=b_old+b
127 x_3=....=P^{-1}(b_new-AP^{-1}b_old)
128 b_old=b_new
129
130 So for n- steps we will use loop, but we always make sure that every value calculated only once!
131 */
132
133 Paso_Solver_solveGS(prec->gs,x,b);
134 if (prec->gs->sweeps>1) {
135 double *bold=MEMALLOC(prec->gs->n*prec->gs->n_block,double);
136 double *bnew=MEMALLOC(prec->gs->n*prec->gs->n_block,double);
137 dim_t i;
138 #pragma omp parallel for private(i) schedule(static)
139 for (i=0;i<prec->gs->n*prec->gs->n_block;++i) bold[i]=b[i];
140
141 while(prec->gs->sweeps>1) {
142 #pragma omp parallel for private(i) schedule(static)
143 for (i=0;i<prec->gs->n*prec->gs->n_block;++i) bnew[i]=bold[i]+b[i];
144 /* Compute the residual b=b-Ax*/
145 Paso_SystemMatrix_MatrixVector_CSR_OFFSET0(DBLE(-1), A, x, DBLE(1), bnew);
146 /* Go round again*/
147 Paso_Solver_solveGS(prec->gs,x,bnew);
148 #pragma omp parallel for private(i) schedule(static)
149 for (i=0;i<prec->gs->n*prec->gs->n_block;++i) bold[i]=bnew[i];
150 prec->gs->sweeps=prec->gs->sweeps-1;
151 }
152
153 MEMFREE(bold);
154 MEMFREE(bnew);
155
156 }
157 break;
158 case PASO_AMG:
159 Paso_Solver_solveAMG(prec->amg,x,b);
160 break;
161 }
162
163 }

Properties

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

  ViewVC Help
Powered by ViewVC 1.1.26