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
|
// 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_FULL_H
#define LINEAR_SYSTEM_FULL_H
// Interface to a linear system with a FULL matrix
#include "GmshMessage.h"
#include "linearSystem.h"
#include "fullMatrix.h"
#include <stdlib.h>
#include <set>
template <class scalar> class linearSystemFull : public linearSystem<scalar> {
private:
fullMatrix<scalar> *_a;
fullVector<scalar> *_b, *_x;
public:
linearSystemFull() : _a(0), _b(0), _x(0) {}
virtual bool isAllocated() const { return _a != 0; }
virtual void allocate(int nbRows)
{
clear();
_a = new fullMatrix<scalar>(nbRows, nbRows);
_b = new fullVector<scalar>(nbRows);
_x = new fullVector<scalar>(nbRows);
}
virtual ~linearSystemFull() { 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 addToSolution(int row, const scalar &val)
{
if(val != 0.0) (*_x)(row) += val;
}
virtual void getFromRightHandSide(int row, scalar &val) const
{
val = (*_b)(row);
}
virtual void getFromSolution(int row, scalar &val) const { val = (*_x)(row); }
virtual void zeroMatrix() { _a->setAll(0.); }
virtual void zeroRightHandSide() { _b->setAll(0.); }
virtual void zeroSolution() { _x->setAll(0.); }
virtual double normInfRightHandSide() const
{
double nor = 0.;
double temp;
for(int i = 0; i < _b->size(); i++) {
temp = (*_b)(i);
if(temp < 0) temp = -temp;
if(nor < temp) nor = temp;
}
return nor;
}
virtual int systemSolve()
{
if(_b->size()) _a->luSolve(*_b, *_x);
// _x->print("X in solve");
return 1;
}
};
#endif
|
