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 |
#ifndef INC_SOLVER |
16 |
#define INC_SOLVER |
17 |
|
18 |
#include "SystemMatrix.h" |
19 |
#include "performance.h" |
20 |
#include "Functions.h" |
21 |
|
22 |
#define PASO_TRACE |
23 |
/* error codes used in the solver */ |
24 |
#define SOLVER_NO_ERROR 0 |
25 |
#define SOLVER_MAXITER_REACHED 1 |
26 |
#define SOLVER_INPUT_ERROR -1 |
27 |
#define SOLVER_MEMORY_ERROR -9 |
28 |
#define SOLVER_BREAKDOWN -10 |
29 |
#define SOLVER_NEGATIVE_NORM_ERROR -11 |
30 |
|
31 |
|
32 |
static double ONE=1.; |
33 |
static double ZERO=0.; |
34 |
static double TOLERANCE_FOR_SCALARS=0.; |
35 |
|
36 |
/* jacobi preconditioner */ |
37 |
|
38 |
typedef struct Paso_Solver_Jacobi { |
39 |
dim_t n_block; |
40 |
dim_t n; |
41 |
double* values; |
42 |
index_t* pivot; |
43 |
} Paso_Solver_Jacobi; |
44 |
|
45 |
|
46 |
/* ILU preconditioner */ |
47 |
struct Paso_Solver_ILU { |
48 |
dim_t n_block; |
49 |
dim_t n; |
50 |
index_t num_colors; |
51 |
index_t* colorOf; |
52 |
index_t* main_iptr; |
53 |
double* factors; |
54 |
Paso_Pattern* pattern; |
55 |
}; |
56 |
typedef struct Paso_Solver_ILU Paso_Solver_ILU; |
57 |
|
58 |
/* GS preconditioner */ |
59 |
struct Paso_Solver_GS { |
60 |
dim_t n_block; |
61 |
dim_t n; |
62 |
index_t num_colors; |
63 |
index_t* colorOf; |
64 |
index_t* main_iptr; |
65 |
double* diag; |
66 |
Paso_SparseMatrix * factors; |
67 |
Paso_Pattern* pattern; |
68 |
dim_t sweeps; |
69 |
double* x_old; |
70 |
}; |
71 |
typedef struct Paso_Solver_GS Paso_Solver_GS; |
72 |
|
73 |
/* RILU preconditioner */ |
74 |
struct Paso_Solver_RILU { |
75 |
dim_t n; |
76 |
dim_t n_block; |
77 |
dim_t n_F; |
78 |
dim_t n_C; |
79 |
double* inv_A_FF; |
80 |
index_t* A_FF_pivot; |
81 |
Paso_SparseMatrix * A_FC; |
82 |
Paso_SparseMatrix * A_CF; |
83 |
index_t* rows_in_F; |
84 |
index_t* rows_in_C; |
85 |
index_t* mask_F; |
86 |
index_t* mask_C; |
87 |
double* x_F; |
88 |
double* b_F; |
89 |
double* x_C; |
90 |
double* b_C; |
91 |
struct Paso_Solver_RILU * RILU_of_Schur; |
92 |
}; |
93 |
typedef struct Paso_Solver_RILU Paso_Solver_RILU; |
94 |
|
95 |
/* AMG preconditioner */ |
96 |
struct Paso_Solver_AMG { |
97 |
dim_t n; |
98 |
dim_t level; |
99 |
dim_t n_block; |
100 |
dim_t n_F; |
101 |
dim_t n_C; |
102 |
double* inv_A_FF; |
103 |
index_t* A_FF_pivot; |
104 |
Paso_SparseMatrix * A_FC; |
105 |
Paso_SparseMatrix * A_CF; |
106 |
index_t* rows_in_F; |
107 |
index_t* rows_in_C; |
108 |
index_t* mask_F; |
109 |
index_t* mask_C; |
110 |
double* x_F; |
111 |
double* b_F; |
112 |
double* x_C; |
113 |
double* b_C; |
114 |
Paso_SparseMatrix * A; |
115 |
Paso_Solver_GS* GS; |
116 |
struct Paso_Solver_AMG * AMG_of_Schur; |
117 |
}; |
118 |
typedef struct Paso_Solver_AMG Paso_Solver_AMG; |
119 |
|
120 |
|
121 |
/* general preconditioner interface */ |
122 |
|
123 |
typedef struct Paso_Solver_Preconditioner { |
124 |
dim_t type; |
125 |
/* jacobi preconditioner */ |
126 |
Paso_Solver_Jacobi* jacobi; |
127 |
/* ilu preconditioner */ |
128 |
Paso_Solver_ILU* ilu; |
129 |
/* rilu preconditioner */ |
130 |
Paso_Solver_RILU* rilu; |
131 |
/* Gauss-Seidel preconditioner */ |
132 |
Paso_Solver_GS* gs; |
133 |
/* amg preconditioner */ |
134 |
Paso_Solver_AMG* amg; |
135 |
|
136 |
} Paso_Solver_Preconditioner; |
137 |
|
138 |
void Paso_Solver(Paso_SystemMatrix*,double*,double*,Paso_Options*,Paso_Performance* pp); |
139 |
void Paso_Solver_free(Paso_SystemMatrix*); |
140 |
err_t Paso_Solver_BiCGStab( Paso_SystemMatrix * A, double* B, double * X, dim_t *iter, double * tolerance, Paso_Performance* pp); |
141 |
err_t Paso_Solver_PCG( Paso_SystemMatrix * A, double* B, double * X, dim_t *iter, double * tolerance, Paso_Performance* pp); |
142 |
err_t Paso_Solver_TFQMR( Paso_SystemMatrix * A, double* B, double * X, dim_t *iter, double * tolerance, Paso_Performance* pp); |
143 |
err_t Paso_Solver_MINRES( Paso_SystemMatrix * A, double* B, double * X, dim_t *iter, double * tolerance, Paso_Performance* pp); |
144 |
err_t Paso_Solver_GMRES(Paso_SystemMatrix * A, double * r, double * x, dim_t *num_iter, double * tolerance,dim_t length_of_recursion,dim_t restart, Paso_Performance* pp); |
145 |
void Paso_Preconditioner_free(Paso_Solver_Preconditioner*); |
146 |
void Paso_Solver_setPreconditioner(Paso_SystemMatrix* A,Paso_Options* options); |
147 |
void Paso_Solver_solvePreconditioner(Paso_SystemMatrix* A,double*,double*); |
148 |
void Paso_Solver_applyBlockDiagonalMatrix(dim_t n_block,dim_t n,double* D,index_t* pivot,double* x,double* b); |
149 |
|
150 |
void Paso_Solver_ILU_free(Paso_Solver_ILU * in); |
151 |
Paso_Solver_ILU* Paso_Solver_getILU(Paso_SparseMatrix * A_p,bool_t verbose); |
152 |
void Paso_Solver_solveILU(Paso_Solver_ILU * ilu, double * x, double * b); |
153 |
|
154 |
void Paso_Solver_GS_free(Paso_Solver_GS * in); |
155 |
Paso_Solver_GS* Paso_Solver_getGS(Paso_SparseMatrix * A_p,bool_t verbose); |
156 |
void Paso_Solver_solveGS1(Paso_Solver_GS * gs, double * x, double * b); |
157 |
void Paso_Solver_solveGS(Paso_Solver_GS * gs, double * x, double * b); |
158 |
|
159 |
void Paso_Solver_RILU_free(Paso_Solver_RILU * in); |
160 |
Paso_Solver_RILU* Paso_Solver_getRILU(Paso_SparseMatrix * A_p,bool_t verbose); |
161 |
void Paso_Solver_solveRILU(Paso_Solver_RILU * rilu, double * x, double * b); |
162 |
|
163 |
void Paso_Solver_AMG_free(Paso_Solver_AMG * in); |
164 |
Paso_Solver_AMG* Paso_Solver_getAMG(Paso_SparseMatrix * A_p,bool_t verbose,dim_t level); |
165 |
void Paso_Solver_solveAMG(Paso_Solver_AMG * amg, double * x, double * b); |
166 |
|
167 |
void Paso_Solver_updateIncompleteSchurComplement(Paso_SparseMatrix* A_CC, Paso_SparseMatrix *A_CF,double* invA_FF,index_t* A_FF_pivot, Paso_SparseMatrix *A_FC); |
168 |
Paso_Solver_Jacobi* Paso_Solver_getJacobi(Paso_SparseMatrix * A_p); |
169 |
void Paso_Solver_solveJacobi(Paso_Solver_Jacobi * prec, double * x, double * b); |
170 |
void Paso_Solver_Jacobi_free(Paso_Solver_Jacobi * in); |
171 |
|
172 |
err_t Paso_Solver_GMRES2(Paso_Function * F, const double* f0, const double* x0, double * x, dim_t *iter, double* tolerance, Paso_Performance* pp); |
173 |
err_t Paso_Solver_NewtonGMRES(Paso_Function *F, double *x, Paso_Options* options, Paso_Performance* pp); |
174 |
|
175 |
Paso_Function * Paso_Function_LinearSystem_alloc(Paso_SystemMatrix* A, double* b, Paso_Options* options); |
176 |
err_t Paso_Function_LinearSystem_call(Paso_Function * F,double* value, const double* arg, Paso_Performance *pp); |
177 |
void Paso_Function_LinearSystem_free(Paso_Function * F); |
178 |
err_t Paso_Function_LinearSystem_setInitialGuess(Paso_SystemMatrix* A, double* x, Paso_Performance *pp); |
179 |
|
180 |
#endif /* #ifndef INC_SOLVER */ |