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 |
MEMFREE(in); |
41 |
} |
42 |
} |
43 |
/* call the iterative solver: */ |
44 |
|
45 |
void Paso_Solver_setPreconditioner(Paso_SystemMatrix* A,Paso_Options* options) { |
46 |
Paso_Solver_Preconditioner* prec=NULL; |
47 |
if (A->solver==NULL) { |
48 |
/* allocate structure to hold preconditioner */ |
49 |
prec=MEMALLOC(1,Paso_Solver_Preconditioner); |
50 |
if (Paso_checkPtr(prec)) return; |
51 |
prec->type=UNKNOWN; |
52 |
prec->rilu=NULL; |
53 |
prec->ilu=NULL; |
54 |
prec->jacobi=NULL; |
55 |
prec->gs=NULL; |
56 |
A->solver=prec; |
57 |
switch (options->preconditioner) { |
58 |
default: |
59 |
case PASO_JACOBI: |
60 |
if (options->verbose) printf("Jacobi preconditioner is used.\n"); |
61 |
prec->jacobi=Paso_Solver_getJacobi(A->mainBlock); |
62 |
prec->type=PASO_JACOBI; |
63 |
break; |
64 |
case PASO_ILU0: |
65 |
if (options->verbose) printf("ILU preconditioner is used.\n"); |
66 |
prec->ilu=Paso_Solver_getILU(A->mainBlock,options->verbose); |
67 |
prec->type=PASO_ILU0; |
68 |
break; |
69 |
case PASO_RILU: |
70 |
if (options->verbose) printf("RILU preconditioner is used.\n"); |
71 |
prec->rilu=Paso_Solver_getRILU(A->mainBlock,options->verbose); |
72 |
prec->type=PASO_RILU; |
73 |
break; |
74 |
case PASO_GS: |
75 |
if (options->verbose) printf("Gauss-Seidel preconditioner is used.\n"); |
76 |
prec->gs=Paso_Solver_getGS(A->mainBlock,options->verbose); |
77 |
prec->gs->sweeps=options->sweeps; |
78 |
prec->type=PASO_GS; |
79 |
break; |
80 |
} |
81 |
if (! Paso_MPIInfo_noError(A->mpi_info ) ){ |
82 |
Paso_Preconditioner_free(prec); |
83 |
A->solver=NULL; |
84 |
} |
85 |
} |
86 |
} |
87 |
|
88 |
/* applies the preconditioner */ |
89 |
/* has to be called within a parallel reqion */ |
90 |
/* barrier synchronization is performed before the evaluation to make sure that the input vector is available */ |
91 |
void Paso_Solver_solvePreconditioner(Paso_SystemMatrix* A,double* x,double* b){ |
92 |
Paso_Solver_Preconditioner* prec=(Paso_Solver_Preconditioner*) A->solver; |
93 |
switch (prec->type) { |
94 |
default: |
95 |
case PASO_JACOBI: |
96 |
Paso_Solver_solveJacobi(prec->jacobi,x,b); |
97 |
break; |
98 |
case PASO_ILU0: |
99 |
Paso_Solver_solveILU(prec->ilu,x,b); |
100 |
break; |
101 |
case PASO_RILU: |
102 |
Paso_Solver_solveRILU(prec->rilu,x,b); |
103 |
break; |
104 |
case PASO_GS: |
105 |
Paso_Solver_solveGS(prec->gs,x,b); |
106 |
break; |
107 |
} |
108 |
} |