forked from cycheung/gmsh
-
Notifications
You must be signed in to change notification settings - Fork 0
/
Copy pathlinearSystemPETSc.h
156 lines (146 loc) · 5.18 KB
/
linearSystemPETSc.h
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
// Gmsh - Copyright (C) 1997-2015 C. Geuzaine, J.-F. Remacle
//
// See the LICENSE.txt file for license information. Please report all
// bugs and problems to the public mailing list <[email protected]>.
#ifndef _LINEAR_SYSTEM_PETSC_H_
#define _LINEAR_SYSTEM_PETSC_H_
// Interface to PETSc 3.x
// Options for PETSc can be provided on the command line, or in the
// file ~/.petscrc. By default, we use the following options (GMRES
// with ILU(6) preconditioner and RCMK renumbering):
//
// -ksp_type gmres
// -pc_type ilu
// -pc_factor_levels 6
// -pc_factor_mat_ordering rcm
// -ksp_rtol 1.e-8
//
// Other useful options include:
//
// -ksp_gmres_restart 100
// -ksp_monitor
//
// To use a direct solver (a sparse lu) instead of an iterative
// solver, use
//
// -ksp_type preonly
// -pc_type lu
//
// When PETSc is compiled with external direct solvers (UMFPACK,
// SuperLU, etc.), they can be selected like this:
//
// -pc_factor_mat_solver_package umfpack
#include "GmshConfig.h"
#include "GmshMessage.h"
#include "linearSystem.h"
#include "sparsityPattern.h"
#include "fullMatrix.h"
#include <string.h>
#include <vector>
#if defined(HAVE_PETSC)
#ifndef SWIG
#include "petsc.h"
#include "petscksp.h"
#else
typedef struct _p_Mat* Mat;
typedef struct _p_Vec* Vec;
typedef struct _p_KSP* KSP;
#endif
//support old PETSc version, try to avoid using PETSC_VERSION in other places
#if PETSC_VERSION_RELEASE != 0 && (PETSC_VERSION_MAJOR < 3 || (PETSC_VERSION_MAJOR == 3 && PETSC_VERSION_MINOR < 2))
#define KSPDestroy(k) KSPDestroy(*(k))
#define MatDestroy(m) MatDestroy(*(m))
#define VecDestroy(v) VecDestroy(*(v))
#define PetscViewerDestroy(v) PetscViewerDestroy(*(v))
#define PetscBool PetscTruth
#define PetscOptionsGetBool PetscOptionsGetTruth
#endif
// The petsc3.5 change log says:
// "KSPSetOperators() no longer has the MatStructure argument. The Mat objects now track that information themselves.
// Use KPS/PCSetReusePreconditioner() to prevent the recomputation of the preconditioner if the operator changed
// in the way that SAME_PRECONDITIONER did with KSPSetOperators()"
// So I guess this should be called with PETSC_TRUE as second argument only for SAME_PRECONDITIONER
// and false otherwise (i.e. for SAME_NONZERO_PATTERN, DIFFRENT_NONZERO_PATTERN) but it is a guess...
#if (PETSC_VERSION_MAJOR < 3 || (PETSC_VERSION_MAJOR == 3 && PETSC_VERSION_MINOR < 5))
#define KSPSetOperators(_ksp, _a, _b, SAME_PRECONDITIONER) KSPSetOperators(_ksp, _a, _b, SAME_PRECONDITIONER)
#else
# define SAME_PRECONDITIONER 999
# define KSPSetOperators(_ksp, _a, _b, OPT_PRECON) (KSPSetOperators(_ksp, _a, _b), KSPSetReusePreconditioner(_ksp,PetscBool(OPT_PRECON == SAME_PRECONDITIONER)))
#endif
template <class scalar>
class linearSystemPETSc : public linearSystem<scalar> {
protected:
bool _isAllocated, _kspAllocated, _entriesPreAllocated;
bool _matrixChangedSinceLastSolve;
bool _valuesNotAssembled; //cannot use MatAssembled since MatAssembled return false for an empty matrix
Mat _a;
Vec _b, _x;
KSP _ksp;
int _localRowStart, _localRowEnd, _localSize, _globalSize;
sparsityPattern _sparsity;
void _kspCreate();
void _assembleMatrixIfNeeded();
#ifndef SWIG
MPI_Comm _comm;
#endif
int _getBlockSizeFromParameters() const;
public:
~linearSystemPETSc();
void insertInSparsityPattern (int i, int j);
bool isAllocated() const { return _isAllocated; }
void preAllocateEntries();
virtual void allocate(int nbRows);
void print();
void clear();
void getFromMatrix(int row, int col, scalar &val) const;
void addToRightHandSide(int row, const scalar &val);
void getFromRightHandSide(int row, scalar &val) const;
double normInfRightHandSide() const;
int getNumKspIteration() const;
void addToMatrix(int row, int col, const scalar &val);
void addToSolution(int row, const scalar &val);
void getFromSolution(int row, scalar &val) const;
void zeroMatrix();
void zeroRightHandSide();
void zeroSolution();
void printMatlab(const char *filename) const;
virtual int systemSolve();
Mat getMatrix(){ return _a; }
virtual int matMult();
//std::vector<scalar> getData();
//std::vector<int> getRowPointers();
//std::vector<int> getColumnsIndices();
#ifndef SWIG
linearSystemPETSc(MPI_Comm com);
MPI_Comm& getComm() {return _comm;}
#endif
linearSystemPETSc();
};
#else
template <class scalar>
class linearSystemPETSc : public linearSystem<scalar> {
public :
linearSystemPETSc()
{
Msg::Error("PETSc is not available in this version of Gmsh");
}
bool isAllocated() const { return false; }
void allocate(int nbRows) {}
void clear(){}
void addToMatrix(int row, int col, const scalar &val) {}
void getFromMatrix(int row, int col, scalar &val) const {}
void addToRightHandSide(int row, const scalar &val) {}
void addToSolution(int row, const scalar &val) {}
void getFromRightHandSide(int row, scalar &val) const {}
void getFromSolution(int row, scalar &val) const {}
int getNumKspIteration() const {return 0;};
void zeroMatrix() {}
void zeroRightHandSide() {}
void zeroSolution() {}
void printMatlab(const char *filename) const{};
virtual int systemSolve() { return 0; }
double normInfRightHandSide() const{return 0;}
virtual int matMult() { return 0; }
};
#endif
#endif