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

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

Parent Directory Parent Directory | Revision Log Revision Log


Revision 2520 - (show annotations)
Mon Jul 6 03:13:11 2009 UTC (10 years, 2 months ago) by artak
Original Path: trunk/paso/src/Solver_preconditioner.c
File MIME type: text/plain
File size: 5874 byte(s)
AMG now takes as an input Paso_options for selecting coarsenening algorithm and treshold variable.
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->level_max,options);
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