OGS 6.1.0-1699-ge946d4c5f
LinAlg.h
Go to the documentation of this file.
1 
10 #pragma once
11 
12 #include <cassert>
13 #include "BaseLib/Error.h"
14 #include "LinAlgEnums.h"
15 
16 namespace MathLib
17 {
18 
29 namespace LinAlg
30 {
31 
32 // Matrix or Vector
33 
35 template<typename MatrixOrVector>
36 void copy(MatrixOrVector const& x, MatrixOrVector& y)
37 {
38  y = x;
39 }
40 
42 template<typename MatrixOrVector>
43 void scale(MatrixOrVector& x, double const a)
44 {
45  x *= a;
46 }
47 
49 template<typename MatrixOrVector>
50 void aypx(MatrixOrVector& y, double const a, MatrixOrVector const& x)
51 {
52  y = a*y + x;
53 }
54 
56 template<typename MatrixOrVector>
57 void axpy(MatrixOrVector& y, double const a, MatrixOrVector const& x)
58 {
59  y += a*x;
60 }
61 
63 template<typename MatrixOrVector>
64 void axpby(MatrixOrVector& y, double const a, double const b, MatrixOrVector const& x)
65 {
66  y = a*x + b*y;
67 }
68 
70 template<typename MatrixOrVector>
71 void componentwiseDivide(MatrixOrVector& w,
72  MatrixOrVector const& x, MatrixOrVector const& y);
73 
75 template<typename MatrixOrVector>
76 double norm1(MatrixOrVector const& x);
77 
79 template<typename MatrixOrVector>
80 double norm2(MatrixOrVector const& x);
81 
83 template<typename MatrixOrVector>
84 double normMax(MatrixOrVector const& x);
85 
86 template<typename MatrixOrVector>
87 double norm(MatrixOrVector const& x, MathLib::VecNormType type)
88 {
89  switch (type) {
91  return norm1(x);
93  return norm2(x);
95  return normMax(x);
96  default:
97  OGS_FATAL("Invalid norm type given.");
98  }
99 }
100 
101 template <typename Matrix>
102 void finalizeAssembly(Matrix& /*A*/);
103 
104 // Matrix and Vector
105 
112 template<typename Matrix, typename Vector>
113 void matMult(Matrix const& A, Vector const& x, Vector& y)
114 {
115  assert(&x != &y);
116  y = A*x;
117 }
118 
125 template<typename Matrix, typename Vector>
126 void matMultAdd(Matrix const& A, Vector const& v1, Vector const& v2, Vector& v3)
127 {
128  assert(&v1 != &v3);
129  v3 = v2 + A*v1;
130 }
131 
132 }} // namespaces
133 
134 
135 // Global PETScMatrix/PETScVector //////////////////////////////////////////
136 #ifdef USE_PETSC
137 
138 namespace MathLib {
139 
140 class PETScMatrix;
141 class PETScVector;
142 
143 namespace LinAlg
144 {
145 
146 // Vector
147 
150 void setLocalAccessibleVector(PETScVector const& x);
151 
152 void set(PETScVector& x, double const a);
153 
154 void copy(PETScVector const& x, PETScVector& y);
155 
156 void scale(PETScVector& x, double const a);
157 
158 // y = a*y + X
159 void aypx(PETScVector& y, double const a, PETScVector const& x);
160 
161 // y = a*x + y
162 void axpy(PETScVector& y, double const a, PETScVector const& x);
163 
164 // y = a*x + y
165 void axpby(PETScVector& y, double const a, double const b, PETScVector const& x);
166 
167 
168 // Matrix
169 
170 void copy(PETScMatrix const& A, PETScMatrix& B);
171 
172 // A = a*A
173 void scale(PETScMatrix& A, double const a);
174 
175 // Y = a*Y + X
176 void aypx(PETScMatrix& Y, double const a, PETScMatrix const& X);
177 
178 // Y = a*X + Y
179 void axpy(PETScMatrix& Y, double const a, PETScMatrix const& X);
180 
181 
182 // Matrix and Vector
183 
184 // v3 = A*v1 + v2
185 void matMult(PETScMatrix const& A, PETScVector const& x, PETScVector& y);
186 
187 // y = A*x
188 void matMultAdd(PETScMatrix const& A, PETScVector const& v1,
189  PETScVector const& v2, PETScVector& v3);
190 
193 
194 }} // namespaces
195 
196 
197 // Sparse global EigenMatrix/EigenVector //////////////////////////////////////////
198 #elif defined(OGS_USE_EIGEN)
199 
200 namespace MathLib {
201 
202 class EigenMatrix;
203 class EigenVector;
204 
205 namespace LinAlg
206 {
207 
208 // Vector
209 
218 void setLocalAccessibleVector(EigenVector const& x);
219 
220 void set(EigenVector& x, double const a);
221 
222 void copy(EigenVector const& x, EigenVector& y);
223 
224 void scale(EigenVector& x, double const a);
225 
226 // y = a*y + X
227 void aypx(EigenVector& y, double const a, EigenVector const& x);
228 
229 // y = a*x + y
230 void axpy(EigenVector& y, double const a, EigenVector const& x);
231 
232 // y = a*x + y
233 void axpby(EigenVector& y, double const a, double const b, EigenVector const& x);
234 
235 
236 // Matrix
237 
238 void copy(EigenMatrix const& A, EigenMatrix& B);
239 
240 // A = a*A
241 void scale(EigenMatrix& A, double const a);
242 
243 // Y = a*Y + X
244 void aypx(EigenMatrix& Y, double const a, EigenMatrix const& X);
245 
246 // Y = a*X + Y
247 void axpy(EigenMatrix& Y, double const a, EigenMatrix const& X);
248 
249 
250 // Matrix and Vector
251 
252 // y = A*x
253 void matMult(EigenMatrix const& A, EigenVector const& x, EigenVector& y);
254 
255 // v3 = A*v1 + v2
256 void matMultAdd(EigenMatrix const& A, EigenVector const& v1,
257  EigenVector const& v2, EigenVector& v3);
258 
261 
262 } // namespace LinAlg
263 
264 } // namespace MathLib
265 
266 #endif
Wrapper class for PETSc vector.
Definition: PETScVector.h:32
void aypx(MatrixOrVector &y, double const a, MatrixOrVector const &x)
Computes .
Definition: LinAlg.h:50
void matMult(Matrix const &A, Vector const &x, Vector &y)
Definition: LinAlg.h:113
double normMax(MatrixOrVector const &x)
Computes the maximum norm of x.
void finalizeAssembly(Matrix &)
VecNormType
Norm type. Not declared as class type in order to use the members as integers.
Definition: LinAlgEnums.h:20
void componentwiseDivide(MatrixOrVector &w, MatrixOrVector const &x, MatrixOrVector const &y)
Computes componentwise.
void axpby(MatrixOrVector &y, double const a, double const b, MatrixOrVector const &x)
Computes .
Definition: LinAlg.h:64
double norm(MatrixOrVector const &x, MathLib::VecNormType type)
Definition: LinAlg.h:87
Global vector based on Eigen vector.
Definition: EigenVector.h:25
Wrapper class for PETSc matrix routines for matrix.
Definition: PETScMatrix.h:32
#define OGS_FATAL(fmt,...)
Definition: Error.h:71
void axpy(MatrixOrVector &y, double const a, MatrixOrVector const &x)
Computes .
Definition: LinAlg.h:57
double norm2(MatrixOrVector const &x)
Computes the Euclidean norm of x.
void copy(MatrixOrVector const &x, MatrixOrVector &y)
Copies x to y.
Definition: LinAlg.h:36
std::array< double, 3 > Vector
Definition: VariableType.h:27
void matMultAdd(Matrix const &A, Vector const &v1, Vector const &v2, Vector &v3)
Definition: LinAlg.h:126
void scale(MatrixOrVector &x, double const a)
Scales x by a.
Definition: LinAlg.h:43
double norm1(MatrixOrVector const &x)
Computes the Manhattan norm of x.