OGS
LinAlg.cpp
Go to the documentation of this file.
1
11#include "LinAlg.h"
12
13// TODO reorder LinAlg function signatures?
14
15// Global PETScMatrix/PETScVector //////////////////////////////////////////
16#ifdef USE_PETSC
17
20
21namespace MathLib
22{
23namespace LinAlg
24{
25// Vector
26
28{
30}
31
32void set(PETScVector& x, PetscScalar const a)
33{
34 VecSet(x.getRawVector(), a);
35}
36
37void copy(PETScVector const& x, PETScVector& y)
38{
39 if (!y.getRawVector())
40 y.shallowCopy(x);
41 VecCopy(x.getRawVector(), y.getRawVector());
42}
43
44void scale(PETScVector& x, PetscScalar const a)
45{
46 VecScale(x.getRawVector(), a);
47}
48
49// y = a*y + X
50void aypx(PETScVector& y, PetscScalar const a, PETScVector const& x)
51{
52 // TODO check sizes
53 VecAYPX(y.getRawVector(), a, x.getRawVector());
54}
55
56// y = a*x + y
57void axpy(PETScVector& y, PetscScalar const a, PETScVector const& x)
58{
59 // TODO check sizes
60 VecAXPY(y.getRawVector(), a, x.getRawVector());
61}
62
63// y = a*x + b*y
64void axpby(PETScVector& y, PetscScalar const a, PetscScalar const b,
65 PETScVector const& x)
66{
67 // TODO check sizes
68 VecAXPBY(y.getRawVector(), a, b, x.getRawVector());
69}
70
71// Explicit specialization
72// Computes w = x/y componentwise.
73template <>
75 PETScVector const& y)
76{
77 VecPointwiseDivide(w.getRawVector(), x.getRawVector(), y.getRawVector());
78}
79
80// Explicit specialization
81// Computes the Manhattan norm of x
82template <>
83double norm1(PETScVector const& x)
84{
85 PetscScalar norm = 0.;
86 VecNorm(x.getRawVector(), NORM_1, &norm);
87 return norm;
88}
89
90// Explicit specialization
91// Computes the Euclidean norm of x
92template <>
93double norm2(PETScVector const& x)
94{
95 PetscScalar norm = 0.;
96 VecNorm(x.getRawVector(), NORM_2, &norm);
97 return norm;
98}
99
100// Explicit specialization
101// Computes the Maximum norm of x
102template <>
103double normMax(PETScVector const& x)
104{
105 PetscScalar norm = 0.;
106 VecNorm(x.getRawVector(), NORM_INFINITY, &norm);
107 return norm;
108}
109
110// Matrix
111
112void copy(PETScMatrix const& A, PETScMatrix& B)
113{
114 B = A;
115}
116
117// A = a*A
118void scale(PETScMatrix& A, PetscScalar const a)
119{
120 MatScale(A.getRawMatrix(), a);
121}
122
123// Y = a*Y + X
124void aypx(PETScMatrix& Y, PetscScalar const a, PETScMatrix const& X)
125{
126 // TODO check sizes
127 // TODO sparsity pattern, currently they are assumed to be different (slow)
128 MatAYPX(Y.getRawMatrix(), a, X.getRawMatrix(), DIFFERENT_NONZERO_PATTERN);
129}
130
131// Y = a*X + Y
132void axpy(PETScMatrix& Y, PetscScalar const a, PETScMatrix const& X)
133{
134 // TODO check sizes
135 // TODO sparsity pattern, currently they are assumed to be different (slow)
136 MatAXPY(Y.getRawMatrix(), a, X.getRawMatrix(), DIFFERENT_NONZERO_PATTERN);
137}
138
139// Matrix and Vector
140
141// v3 = A*v1 + v2
142void matMult(PETScMatrix const& A, PETScVector const& x, PETScVector& y)
143{
144 // TODO check sizes
145 assert(&x != &y);
146 if (!y.getRawVector())
147 y.shallowCopy(x);
148 MatMult(A.getRawMatrix(), x.getRawVector(), y.getRawVector());
149}
150
151// v3 = A*v1 + v2
152void matMultAdd(PETScMatrix const& A, PETScVector const& v1,
153 PETScVector const& v2, PETScVector& v3)
154{
155 // TODO check sizes
156 assert(&v1 != &v3);
157 if (!v3.getRawVector())
158 v3.shallowCopy(v1);
159 MatMultAdd(A.getRawMatrix(), v1.getRawVector(), v2.getRawVector(),
160 v3.getRawVector());
161}
162
164{
165 A.finalizeAssembly(MAT_FINAL_ASSEMBLY);
166}
167
169{
171}
172
173} // namespace LinAlg
174} // namespace MathLib
175
176// Sparse global EigenMatrix/EigenVector
177// //////////////////////////////////////////
178#else
179
182
183namespace MathLib
184{
185namespace LinAlg
186{
187// Vector
188
189void setLocalAccessibleVector(EigenVector const& /*x*/) {}
190
191void set(EigenVector& x, double const a)
192{
193 x.getRawVector().setConstant(a);
194}
195
196void copy(EigenVector const& x, EigenVector& y)
197{
198 y = x;
199}
200
201void scale(EigenVector& x, double const a)
202{
203 x.getRawVector() *= a;
204}
205
206// y = a*y + X
207void aypx(EigenVector& y, double const a, EigenVector const& x)
208{
209 // TODO: does that break anything?
210 y.getRawVector() = a * y.getRawVector() + x.getRawVector();
211}
212
213// y = a*x + y
214void axpy(EigenVector& y, double const a, EigenVector const& x)
215{
216 // TODO: does that break anything?
217 y.getRawVector() += a * x.getRawVector();
218}
219
220// y = a*x + y
221void axpby(EigenVector& y, double const a, double const b, EigenVector const& x)
222{
223 // TODO: does that break anything?
224 y.getRawVector() = a * x.getRawVector() + b * y.getRawVector();
225}
226
227// Explicit specialization
228// Computes w = x/y componentwise.
229template <>
230void componentwiseDivide(EigenVector& w, EigenVector const& x,
231 EigenVector const& y)
232{
233 w.getRawVector().noalias() =
234 x.getRawVector().cwiseQuotient(y.getRawVector());
235}
236
237// Explicit specialization
238// Computes the Manhattan norm of x
239template <>
240double norm1(EigenVector const& x)
241{
242 return x.getRawVector().lpNorm<1>();
243}
244
245// Explicit specialization
246// Euclidean norm
247template <>
248double norm2(EigenVector const& x)
249{
250 return x.getRawVector().norm();
251}
252
253// Explicit specialization
254// Computes the Maximum norm of x
255template <>
256double normMax(EigenVector const& x)
257{
258 return x.getRawVector().lpNorm<Eigen::Infinity>();
259}
260
261// Matrix
262
263void copy(EigenMatrix const& A, EigenMatrix& B)
264{
265 B = A;
266}
267
268// A = a*A
269void scale(EigenMatrix& A, double const a)
270{
271 // TODO: does that break anything?
272 A.getRawMatrix() *= a;
273}
274
275// Y = a*Y + X
276void aypx(EigenMatrix& Y, double const a, EigenMatrix const& X)
277{
278 // TODO: does that break anything?
279 Y.getRawMatrix() = a * Y.getRawMatrix() + X.getRawMatrix();
280}
281
282// Y = a*X + Y
283void axpy(EigenMatrix& Y, double const a, EigenMatrix const& X)
284{
285 // TODO: does that break anything?
286 Y.getRawMatrix() = a * X.getRawMatrix() + Y.getRawMatrix();
287}
288
289// Matrix and Vector
290
291// v3 = A*v1 + v2
292void matMult(EigenMatrix const& A, EigenVector const& x, EigenVector& y)
293{
294 assert(&x != &y);
295 y.getRawVector() = A.getRawMatrix() * x.getRawVector();
296}
297
298// v3 = A*v1 + v2
299void matMultAdd(EigenMatrix const& A, EigenVector const& v1,
300 EigenVector const& v2, EigenVector& v3)
301{
302 assert(&v1 != &v3);
303 // TODO: does that break anything?
304 v3.getRawVector() =
305 v2.getRawVector() + A.getRawMatrix() * v1.getRawVector();
306}
307
308void finalizeAssembly(EigenMatrix& x)
309{
310 x.getRawMatrix().makeCompressed();
311}
312
313void finalizeAssembly(EigenVector& /*x*/) {}
314
315} // namespace LinAlg
316
317} // namespace MathLib
318
319#endif
Declaration of class PETScMatrix, which provides an interface to PETSc matrix routines.
Declaration of class PETScVector, which provides an interface to PETSc vector routines.
Wrapper class for PETSc matrix routines for matrix.
Definition: PETScMatrix.h:32
Mat & getRawMatrix()
Get matrix reference.
Definition: PETScMatrix.h:86
void finalizeAssembly(const MatAssemblyType asm_type=MAT_FINAL_ASSEMBLY)
Perform MPI collection of assembled entries in buffer.
Definition: PETScMatrix.h:67
Wrapper class for PETSc vector.
Definition: PETScVector.h:33
void finalizeAssembly()
Perform MPI collection of assembled entries in buffer.
void setLocalAccessibleVector() const
void shallowCopy(const PETScVector &v)
PETSc_Vec & getRawVector()
Exposes the underlying PETSc vector.
Definition: PETScVector.h:171
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