OGS 6.2.2-87-g988ee9c30.dirty.20200123122242
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 namespace MathLib
18 {
19 
30 namespace LinAlg
31 {
32 
33 // Matrix or Vector
34 
36 template<typename MatrixOrVector>
37 void copy(MatrixOrVector const& x, MatrixOrVector& y)
38 {
39  y = x;
40 }
41 
43 template<typename MatrixOrVector>
44 void scale(MatrixOrVector& x, double const a)
45 {
46  x *= a;
47 }
48 
50 template<typename MatrixOrVector>
51 void aypx(MatrixOrVector& y, double const a, MatrixOrVector const& x)
52 {
53  y = a*y + x;
54 }
55 
57 template<typename MatrixOrVector>
58 void axpy(MatrixOrVector& y, double const a, MatrixOrVector const& x)
59 {
60  y += a*x;
61 }
62 
64 template<typename MatrixOrVector>
65 void axpby(MatrixOrVector& y, double const a, double const b, MatrixOrVector const& x)
66 {
67  y = a*x + b*y;
68 }
69 
71 template<typename MatrixOrVector>
72 void componentwiseDivide(MatrixOrVector& w,
73  MatrixOrVector const& x, MatrixOrVector const& y);
74 
76 template<typename MatrixOrVector>
77 double norm1(MatrixOrVector const& x);
78 
80 template<typename MatrixOrVector>
81 double norm2(MatrixOrVector const& x);
82 
84 template<typename MatrixOrVector>
85 double normMax(MatrixOrVector const& x);
86 
87 template<typename MatrixOrVector>
88 double norm(MatrixOrVector const& x, MathLib::VecNormType type)
89 {
90  switch (type) {
92  return norm1(x);
94  return norm2(x);
96  return normMax(x);
97  default:
98  OGS_FATAL("Invalid norm type given.");
99  }
100 }
101 
102 template <typename Matrix>
103 void finalizeAssembly(Matrix& /*A*/);
104 
105 // Matrix and Vector
106 
113 template<typename Matrix, typename Vector>
114 void matMult(Matrix const& A, Vector const& x, Vector& y)
115 {
116  assert(&x != &y);
117  y = A*x;
118 }
119 
126 template<typename Matrix, typename Vector>
127 void matMultAdd(Matrix const& A, Vector const& v1, Vector const& v2, Vector& v3)
128 {
129  assert(&v1 != &v3);
130  v3 = v2 + A*v1;
131 }
132 
133 } // namespace LinAlg
134 } // namespace MathLib
135 
136 // Global PETScMatrix/PETScVector //////////////////////////////////////////
137 #ifdef USE_PETSC
138 
139 namespace MathLib {
140 
141 class PETScMatrix;
142 class PETScVector;
143 
144 namespace LinAlg
145 {
146 
147 // Vector
148 
151 void setLocalAccessibleVector(PETScVector const& x);
152 
153 void set(PETScVector& x, double const a);
154 
155 void copy(PETScVector const& x, PETScVector& y);
156 
157 void scale(PETScVector& x, double const a);
158 
159 // y = a*y + X
160 void aypx(PETScVector& y, double const a, PETScVector const& x);
161 
162 // y = a*x + y
163 void axpy(PETScVector& y, double const a, PETScVector const& x);
164 
165 // y = a*x + y
166 void axpby(PETScVector& y, double const a, double const b, PETScVector const& x);
167 
168 
169 // Matrix
170 
171 void copy(PETScMatrix const& A, PETScMatrix& B);
172 
173 // A = a*A
174 void scale(PETScMatrix& A, double const a);
175 
176 // Y = a*Y + X
177 void aypx(PETScMatrix& Y, double const a, PETScMatrix const& X);
178 
179 // Y = a*X + Y
180 void axpy(PETScMatrix& Y, double const a, PETScMatrix const& X);
181 
182 
183 // Matrix and Vector
184 
185 // v3 = A*v1 + v2
186 void matMult(PETScMatrix const& A, PETScVector const& x, PETScVector& y);
187 
188 // y = A*x
189 void matMultAdd(PETScMatrix const& A, PETScVector const& v1,
190  PETScVector const& v2, PETScVector& v3);
191 
194 
195 }} // namespaces
196 
197 
198 // Sparse global EigenMatrix/EigenVector //////////////////////////////////////////
199 #elif defined(OGS_USE_EIGEN)
200 
201 namespace MathLib {
202 
203 class EigenMatrix;
204 class EigenVector;
205 
206 namespace LinAlg
207 {
208 
209 // Vector
210 
219 void setLocalAccessibleVector(EigenVector const& x);
220 
221 void set(EigenVector& x, double const a);
222 
223 void copy(EigenVector const& x, EigenVector& y);
224 
225 void scale(EigenVector& x, double const a);
226 
227 // y = a*y + X
228 void aypx(EigenVector& y, double const a, EigenVector const& x);
229 
230 // y = a*x + y
231 void axpy(EigenVector& y, double const a, EigenVector const& x);
232 
233 // y = a*x + y
234 void axpby(EigenVector& y, double const a, double const b, EigenVector const& x);
235 
236 
237 // Matrix
238 
239 void copy(EigenMatrix const& A, EigenMatrix& B);
240 
241 // A = a*A
242 void scale(EigenMatrix& A, double const a);
243 
244 // Y = a*Y + X
245 void aypx(EigenMatrix& Y, double const a, EigenMatrix const& X);
246 
247 // Y = a*X + Y
248 void axpy(EigenMatrix& Y, double const a, EigenMatrix const& X);
249 
250 
251 // Matrix and Vector
252 
253 // y = A*x
254 void matMult(EigenMatrix const& A, EigenVector const& x, EigenVector& y);
255 
256 // v3 = A*v1 + v2
257 void matMultAdd(EigenMatrix const& A, EigenVector const& v1,
258  EigenVector const& v2, EigenVector& v3);
259 
262 
263 } // namespace LinAlg
264 
265 } // namespace MathLib
266 
267 #endif
Wrapper class for PETSc vector.
Definition: PETScVector.h:32
void aypx(MatrixOrVector &y, double const a, MatrixOrVector const &x)
Computes .
Definition: LinAlg.h:51
void matMult(Matrix const &A, Vector const &x, Vector &y)
Definition: LinAlg.h:114
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:21
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:65
double norm(MatrixOrVector const &x, MathLib::VecNormType type)
Definition: LinAlg.h:88
Global vector based on Eigen vector.
Definition: EigenVector.h:26
Wrapper class for PETSc matrix routines for matrix.
Definition: PETScMatrix.h:32
#define OGS_FATAL(fmt,...)
Definition: Error.h:64
void axpy(MatrixOrVector &y, double const a, MatrixOrVector const &x)
Computes .
Definition: LinAlg.h:58
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:37
std::array< double, 3 > Vector
Definition: VariableType.h:26
void matMultAdd(Matrix const &A, Vector const &v1, Vector const &v2, Vector &v3)
Definition: LinAlg.h:127
void scale(MatrixOrVector &x, double const a)
Scales x by a.
Definition: LinAlg.h:44
double norm1(MatrixOrVector const &x)
Computes the Manhattan norm of x.