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

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

Parent Directory Parent Directory | Revision Log Revision Log


Revision 2548 - (hide annotations)
Mon Jul 20 06:20:06 2009 UTC (10 years, 3 months ago) by jfenwick
File MIME type: text/plain
File size: 5874 byte(s)
Updating copyright notices
1 ksteube 1312
2     /*******************************************************
3 ksteube 1811 *
4 jfenwick 2548 * Copyright (c) 2003-2009 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     /* Author: gross@access.edu.au */
23    
24     /**************************************************************/
25    
26 gross 700 #include "Paso.h"
27     #include "SystemMatrix.h"
28 jgs 150 #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 gross 430 Paso_Solver_RILU_free(in->rilu);
38 jgs 150 Paso_Solver_Jacobi_free(in->jacobi);
39 artak 1819 Paso_Solver_GS_free(in->gs);
40 artak 1842 Paso_Solver_AMG_free(in->amg);
41 jgs 150 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 gross 425 if (A->solver==NULL) {
49 jgs 150 /* allocate structure to hold preconditioner */
50     prec=MEMALLOC(1,Paso_Solver_Preconditioner);
51     if (Paso_checkPtr(prec)) return;
52     prec->type=UNKNOWN;
53 gross 430 prec->rilu=NULL;
54 jgs 150 prec->ilu=NULL;
55     prec->jacobi=NULL;
56 artak 1819 prec->gs=NULL;
57 artak 1842 prec->amg=NULL;
58 gross 425 A->solver=prec;
59 jgs 150 switch (options->preconditioner) {
60     default:
61     case PASO_JACOBI:
62     if (options->verbose) printf("Jacobi preconditioner is used.\n");
63 ksteube 1312 prec->jacobi=Paso_Solver_getJacobi(A->mainBlock);
64 jgs 150 prec->type=PASO_JACOBI;
65     break;
66     case PASO_ILU0:
67     if (options->verbose) printf("ILU preconditioner is used.\n");
68 ksteube 1312 prec->ilu=Paso_Solver_getILU(A->mainBlock,options->verbose);
69 jgs 150 prec->type=PASO_ILU0;
70     break;
71 gross 430 case PASO_RILU:
72     if (options->verbose) printf("RILU preconditioner is used.\n");
73 ksteube 1312 prec->rilu=Paso_Solver_getRILU(A->mainBlock,options->verbose);
74 gross 430 prec->type=PASO_RILU;
75     break;
76 artak 1819 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 artak 1823 prec->gs->sweeps=options->sweeps;
80 artak 1819 prec->type=PASO_GS;
81     break;
82 artak 1842 case PASO_AMG:
83     if (options->verbose) printf("AMG preconditioner is used.\n");
84 artak 2520 prec->amg=Paso_Solver_getAMG(A->mainBlock,options->level_max,options);
85 artak 1842 prec->type=PASO_AMG;
86     break;
87    
88 jgs 150 }
89 ksteube 1312 if (! Paso_MPIInfo_noError(A->mpi_info ) ){
90 jgs 150 Paso_Preconditioner_free(prec);
91 gross 425 A->solver=NULL;
92 jgs 150 }
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 gross 425 Paso_Solver_Preconditioner* prec=(Paso_Solver_Preconditioner*) A->solver;
101 artak 1842
102 jgs 150 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 gross 430 case PASO_RILU:
111     Paso_Solver_solveRILU(prec->rilu,x,b);
112     break;
113 artak 1819 case PASO_GS:
114 artak 1842 /* 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 artak 1819 Paso_Solver_solveGS(prec->gs,x,b);
134 artak 1842 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 artak 1819 break;
158 artak 1842 case PASO_AMG:
159     Paso_Solver_solveAMG(prec->amg,x,b);
160     break;
161 jgs 150 }
162 artak 1842
163 jgs 150 }

Properties

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

  ViewVC Help
Powered by ViewVC 1.1.26