OGS
LinAlg.h
Go to the documentation of this file.
1
11#pragma once
12
13#include <cassert>
14#include "BaseLib/Error.h"
15#include "LinAlgEnums.h"
16
17#ifdef USE_PETSC
18#include <petscsystypes.h>
19#endif
20
21namespace MathLib
22{
23
34namespace LinAlg
35{
36
37// Matrix or Vector
38
40template<typename MatrixOrVector>
41void copy(MatrixOrVector const& x, MatrixOrVector& y)
42{
43 y = x;
44}
45
47template <typename MatrixOrVector>
48void scale(MatrixOrVector& x, double const a)
49{
50 x *= a;
51}
52
54template <typename MatrixOrVector>
55void aypx(MatrixOrVector& y, double const a, MatrixOrVector const& x)
56{
57 y = a*y + x;
58}
59
61template <typename MatrixOrVector>
62void axpy(MatrixOrVector& y, double const a, MatrixOrVector const& x)
63{
64 y += a*x;
65}
66
68template <typename MatrixOrVector>
69void axpby(MatrixOrVector& y, double const a, double const b,
70 MatrixOrVector const& x)
71{
72 y = a*x + b*y;
73}
74
76template<typename MatrixOrVector>
77void componentwiseDivide(MatrixOrVector& w,
78 MatrixOrVector const& x, MatrixOrVector const& y);
79
81template<typename MatrixOrVector>
82double norm1(MatrixOrVector const& x);
83
85template<typename MatrixOrVector>
86double norm2(MatrixOrVector const& x);
87
89template<typename MatrixOrVector>
90double normMax(MatrixOrVector const& x);
91
92template<typename MatrixOrVector>
93double norm(MatrixOrVector const& x, MathLib::VecNormType type)
94{
95 switch (type) {
97 return norm1(x);
99 return norm2(x);
101 return normMax(x);
102 default:
103 OGS_FATAL("Invalid norm type given.");
104 }
105}
106
107template <typename Matrix>
108void finalizeAssembly(Matrix& /*A*/);
109
110// Matrix and Vector
111
118template<typename Matrix, typename Vector>
119void matMult(Matrix const& A, Vector const& x, Vector& y)
120{
121 assert(&x != &y);
122 y = A*x;
123}
124
131template<typename Matrix, typename Vector>
132void matMultAdd(Matrix const& A, Vector const& v1, Vector const& v2, Vector& v3)
133{
134 assert(&v1 != &v3);
135 v3 = v2 + A*v1;
136}
137
138} // namespace LinAlg
139} // namespace MathLib
140
141// Global PETScMatrix/PETScVector //////////////////////////////////////////
142#ifdef USE_PETSC
143
144namespace MathLib {
145
146class PETScMatrix;
147class PETScVector;
148
149namespace LinAlg
150{
151
152// Vector
153
156void setLocalAccessibleVector(PETScVector const& x);
157
158void set(PETScVector& x, PetscScalar const a);
159
160void copy(PETScVector const& x, PETScVector& y);
161
162void scale(PETScVector& x, PetscScalar const a);
163
164// y = a*y + X
165void aypx(PETScVector& y, PetscScalar const a, PETScVector const& x);
166
167// y = a*x + y
168void axpy(PETScVector& y, PetscScalar const a, PETScVector const& x);
169
170// y = a*x + y
171void axpby(PETScVector& y, PetscScalar const a, PetscScalar const b,
172 PETScVector const& x);
173
174// Matrix
175
176void copy(PETScMatrix const& A, PETScMatrix& B);
177
178// A = a*A
179void scale(PETScMatrix& A, PetscScalar const a);
180
181// Y = a*Y + X
182void aypx(PETScMatrix& Y, PetscScalar const a, PETScMatrix const& X);
183
184// Y = a*X + Y
185void axpy(PETScMatrix& Y, PetscScalar const a, PETScMatrix const& X);
186
187// Matrix and Vector
188
189// v3 = A*v1 + v2
190void matMult(PETScMatrix const& A, PETScVector const& x, PETScVector& y);
191
192// y = A*x
193void matMultAdd(PETScMatrix const& A, PETScVector const& v1,
194 PETScVector const& v2, PETScVector& v3);
195
196void finalizeAssembly(PETScMatrix& A);
197void finalizeAssembly(PETScVector& x);
198
199}} // namespaces
200
201
202// Sparse global EigenMatrix/EigenVector //////////////////////////////////////////
203#else
204
205namespace MathLib {
206
207class EigenMatrix;
208class EigenVector;
209
210namespace LinAlg
211{
212
213// Vector
214
223void setLocalAccessibleVector(EigenVector const& x);
224
225void set(EigenVector& x, double const a);
226
227void copy(EigenVector const& x, EigenVector& y);
228
229void scale(EigenVector& x, double const a);
230
231// y = a*y + X
232void aypx(EigenVector& y, double const a, EigenVector const& x);
233
234// y = a*x + y
235void axpy(EigenVector& y, double const a, EigenVector const& x);
236
237// y = a*x + y
238void axpby(EigenVector& y, double const a, double const b, EigenVector const& x);
239
240
241// Matrix
242
243void copy(EigenMatrix const& A, EigenMatrix& B);
244
245// A = a*A
246void scale(EigenMatrix& A, double const a);
247
248// Y = a*Y + X
249void aypx(EigenMatrix& Y, double const a, EigenMatrix const& X);
250
251// Y = a*X + Y
252void axpy(EigenMatrix& Y, double const a, EigenMatrix const& X);
253
254
255// Matrix and Vector
256
257// y = A*x
258void matMult(EigenMatrix const& A, EigenVector const& x, EigenVector& y);
259
260// v3 = A*v1 + v2
261void matMultAdd(EigenMatrix const& A, EigenVector const& v1,
262 EigenVector const& v2, EigenVector& v3);
263
264void finalizeAssembly(EigenMatrix& x);
265void finalizeAssembly(EigenVector& A);
266
267} // namespace LinAlg
268
269} // namespace MathLib
270
271#endif
#define OGS_FATAL(...)
Definition: Error.h:26
std::array< double, 3 > Vector
Definition: VariableType.h:28
double norm(MatrixOrVector const &x, MathLib::VecNormType type)
Definition: LinAlg.h:93
void finalizeAssembly(PETScMatrix &A)
Definition: LinAlg.cpp:163
double norm1(PETScVector const &x)
Definition: LinAlg.cpp:83
double normMax(PETScVector const &x)
Definition: LinAlg.cpp:103
void copy(PETScVector const &x, PETScVector &y)
Definition: LinAlg.cpp:37
void componentwiseDivide(PETScVector &w, PETScVector const &x, PETScVector const &y)
Definition: LinAlg.cpp:74
void setLocalAccessibleVector(PETScVector const &x)
Definition: LinAlg.cpp:27
void set(PETScVector &x, PetscScalar const a)
Definition: LinAlg.cpp:32
void matMult(PETScMatrix const &A, PETScVector const &x, PETScVector &y)
Definition: LinAlg.cpp:142
void matMultAdd(PETScMatrix const &A, PETScVector const &v1, PETScVector const &v2, PETScVector &v3)
Definition: LinAlg.cpp:152
double norm2(PETScVector const &x)
Definition: LinAlg.cpp:93
void scale(PETScVector &x, PetscScalar const a)
Definition: LinAlg.cpp:44
void aypx(PETScVector &y, PetscScalar const a, PETScVector const &x)
Definition: LinAlg.cpp:50
void axpy(PETScVector &y, PetscScalar const a, PETScVector const &x)
Definition: LinAlg.cpp:57
void axpby(PETScVector &y, PetscScalar const a, PetscScalar const b, PETScVector const &x)
Definition: LinAlg.cpp:64
VecNormType
Norm type. Not declared as class type in order to use the members as integers.
Definition: LinAlgEnums.h:22