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
|
// Gmsh - Copyright (C) 1997-2020 C. Geuzaine, J.-F. Remacle
//
// See the LICENSE.txt file for license information. Please report all
// issues on https://gitlab.onelab.info/gmsh/gmsh/issues.
#ifndef LINEAR_SYSTEM_GMM_H
#define LINEAR_SYSTEM_GMM_H
#include <string>
#include "GmshConfig.h"
#include "GmshMessage.h"
#include "linearSystem.h"
#if defined(HAVE_GMM)
#undef BB // can be defined by FlGui.h, and clashes with gmm arg name
#include <gmm.h>
// Consider using linearSystemCSRGmm instead: assembly is much faster
template <class scalar> class linearSystemGmm : public linearSystem<scalar> {
protected:
std::vector<scalar> *_x; // nonLinearSystemGmm has to access to this vector
std::vector<scalar> *_b; // idem
gmm::row_matrix<gmm::wsvector<scalar> > *_a;
private:
std::string _method;
double _tol;
int _noisy;
public:
linearSystemGmm(const std::string &method = "gmres", double tol = 1e-8,
int noisy = 0)
: _x(0), _b(0), _a(0), _method(method), _tol(tol), _noisy(noisy) {}
virtual bool isAllocated() const { return _a != 0; }
virtual void allocate(int nbRows)
{
clear();
_a = new gmm::row_matrix<gmm::wsvector<scalar> >(nbRows, nbRows);
_b = new std::vector<scalar>(nbRows);
_x = new std::vector<scalar>(nbRows);
}
virtual ~linearSystemGmm() { clear(); }
virtual void clear()
{
if(_a) {
delete _a;
delete _b;
delete _x;
}
_a = 0;
}
virtual void addToMatrix(int row, int col, const scalar &val)
{
if(val != 0.0) (*_a)(row, col) += val;
}
virtual void getFromMatrix(int row, int col, scalar &val) const
{
val = (*_a)(row, col);
}
virtual void addToRightHandSide(int row, const scalar &val, int ith = 0)
{
if(val != 0.0) (*_b)[row] += val;
}
virtual void getFromRightHandSide(int row, scalar &val) const
{
val = (*_b)[row];
}
virtual void addToSolution(int row, const scalar &val)
{
if(val != 0.0) (*_x)[row] += val;
}
virtual void getFromSolution(int row, scalar &val) const { val = (*_x)[row]; }
virtual void zeroMatrix() { gmm::clear(*_a); }
virtual void zeroRightHandSide()
{
for(std::size_t i = 0; i < _b->size(); i++) (*_b)[i] = 0.;
}
virtual void zeroSolution()
{
for(std::size_t i = 0; i < _x->size(); i++) (*_x)[i] = 0.;
}
virtual double normInfRightHandSide() const
{
double nor = 0.;
double temp;
for(std::size_t i = 0; i < _b->size(); i++) {
temp = abs((*_b)[i]); // this is valid also for complex
// if(temp<0) temp = -temp;
if(nor < temp) nor = temp;
}
return nor;
}
void setPrec(double p) { _tol = p; }
void setNoisy(int n) { _noisy = n; }
void setGmres(int n) { _method = (n ? "gmres" : "cg"); }
virtual int systemSolve()
{
#if defined(HAVE_MUMPS)
if(_method == "mumps"){
gmm::MUMPS_solve(*_a, *_x, *_b);
return 1;
}
#else
//gmm::ilutp_precond<gmm::row_matrix<gmm::wsvector<scalar> > > P(*_a, 25, 0.);
gmm::ildltt_precond<gmm::row_matrix<gmm::wsvector<scalar> > > P(*_a, 30, 1.e-10);
//gmm::ilu_precond<gmm::row_matrix<gmm::wsvector<scalar> > > P(*_a);
gmm::iteration iter(_tol);
iter.set_noisy(_noisy);
if(_method == "gmres")
gmm::gmres(*_a, *_x, *_b, P, 100, iter);
else
gmm::cg(*_a, *_x, *_b, P, iter);
if(!iter.converged())
Msg::Warning("Iterative linear solver has not converged (res = %g)",
iter.get_res());
#endif
return 1;
}
};
#else
template <class scalar> class linearSystemGmm : public linearSystem<scalar> {
public:
linearSystemGmm(const std::string &method = "gmres", double tol = 1e-8,
int noisy = 0)
{
Msg::Error("Gmm++ is not available in this version of Gmsh");
}
virtual bool isAllocated() const { return false; }
virtual void allocate(int nbRows) {}
virtual void addToMatrix(int row, int col, const scalar &val) {}
virtual void getFromMatrix(int row, int col, scalar &val) const {}
virtual void addToRightHandSide(int row, const scalar &val, int ith = 0) {}
virtual void getFromRightHandSide(int row, scalar &val) const {}
virtual void addToSolution(int row, const scalar &val) {}
virtual void getFromSolution(int row, scalar &val) const {}
virtual void zeroMatrix() {}
virtual void zeroRightHandSide() {}
virtual void zeroSolution() {}
virtual int systemSolve() { return 0; }
virtual double normInfRightHandSide() const { return 0.; }
void setPrec(double p) {}
virtual void clear() {}
void setNoisy(int n) {}
void setGmres(int n) {}
};
#endif
#endif
|
