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
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);
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
54}
55
56// y = a*x + y
57void axpy(PETScVector& y, PetscScalar const a, PETScVector const& x)
58{
59 // TODO check sizes
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.
73// \note that VecPointwiseDivide avoids to divide by values that are
74// identically zero such as
75// for (int i=0; i<n; i++)
76// {
77// w[i] = y[i] == 0.0 ? 0.0 : x[i] / y[i];
78// }
79//
80template <>
86
87// Explicit specialization
88// Computes the Manhattan norm of x
89template <>
90double norm1(PETScVector const& x)
91{
92 PetscScalar norm = 0.;
94 return norm;
95}
96
97// Explicit specialization
98// Computes the Euclidean norm of x
99template <>
100double norm2(PETScVector const& x)
101{
102 PetscScalar norm = 0.;
104 return norm;
105}
106
107// Explicit specialization
108// Computes the Maximum norm of x
109template <>
110double normMax(PETScVector const& x)
111{
112 PetscScalar norm = 0.;
114 return norm;
115}
116
117// Matrix
118
119void copy(PETScMatrix const& A, PETScMatrix& B)
120{
121 B = A;
122}
123
124// A = a*A
126{
127 MatScale(A.getRawMatrix(), a);
128}
129
130// Y = a*Y + X
131void aypx(PETScMatrix& Y, PetscScalar const a, PETScMatrix const& X)
132{
133 // TODO check sizes
134 // TODO sparsity pattern, currently they are assumed to be different (slow)
136}
137
138// Y = a*X + Y
139void axpy(PETScMatrix& Y, PetscScalar const a, PETScMatrix const& X)
140{
141 // TODO check sizes
142 // TODO sparsity pattern, currently they are assumed to be different (slow)
144}
145
146// Matrix and Vector
147
148// v3 = A*v1 + v2
149void matMult(PETScMatrix const& A, PETScVector const& x, PETScVector& y)
150{
151 // TODO check sizes
152 assert(&x != &y);
153 if (!y.getRawVector())
154 y.shallowCopy(x);
156}
157
158// v3 = A*v1 + v2
159void matMultAdd(PETScMatrix const& A, PETScVector const& v1,
160 PETScVector const& v2, PETScVector& v3)
161{
162 // TODO check sizes
163 assert(&v1 != &v3);
164 if (!v3.getRawVector())
165 v3.shallowCopy(v1);
166 MatMultAdd(A.getRawMatrix(), v1.getRawVector(), v2.getRawVector(),
167 v3.getRawVector());
168}
169
170void linearSysNormalize(PETScMatrix const& /*A*/, PETScMatrix& /*new_A*/,
171 PETScVector const& /*b*/, PETScVector& /*new_b*/)
172{
173 // The following block is deactivated, because there is no tests yet for the
174 // normalization operation in PETSc. This will be a task for later.
175 /*
176 assert(&A != &new_A);
177 assert(&b != &new_b);
178
179 PetscInt n_rows(0);
180 PetscInt n_cols(0);
181 MatGetSize(A.getRawMatrix(), &n_rows, &n_cols);
182 // only when A matrix is not square
183 if (n_rows != n_cols)
184 {
185 // new_b = A^T * b
186 MatMultTranspose(A.getRawMatrix(), b.getRawVector(),
187 new_b.getRawVector());
188 // new_A = A^T * A
189 MatTranspose(A.getRawMatrix(), MAT_INITIAL_MATRIX,
190 &(new_A.getRawMatrix()));
191 }
192 */
193 OGS_FATAL(
194 "Normalization operation is not implemented yet for PETSc library! "
195 "Program terminated.");
196}
197
202
207
208} // namespace LinAlg
209} // namespace MathLib
210
211// Sparse global EigenMatrix/EigenVector
212// //////////////////////////////////////////
213#else
214
217
218namespace MathLib
219{
220namespace LinAlg
221{
222// Vector
223
224void setLocalAccessibleVector(EigenVector const& /*x*/) {}
225
226void set(EigenVector& x, double const a)
227{
228 x.getRawVector().setConstant(a);
229}
230
231void copy(EigenVector const& x, EigenVector& y)
232{
233 y = x;
234}
235
236void scale(EigenVector& x, double const a)
237{
238 x.getRawVector() *= a;
239}
240
241// y = a*y + X
242void aypx(EigenVector& y, double const a, EigenVector const& x)
243{
244 // TODO: does that break anything?
245 y.getRawVector() = a * y.getRawVector() + x.getRawVector();
246}
247
248// y = a*x + y
249void axpy(EigenVector& y, double const a, EigenVector const& x)
250{
251 // TODO: does that break anything?
252 y.getRawVector() += a * x.getRawVector();
253}
254
255// y = a*x + b*y
256void axpby(EigenVector& y, double const a, double const b, EigenVector const& x)
257{
258 // TODO: does that break anything?
259 y.getRawVector() = a * x.getRawVector() + b * y.getRawVector();
260}
261
262// Explicit specialization
263// Computes w = x/y componentwise.
264template <>
265void componentwiseDivide(EigenVector& w, EigenVector const& x,
266 EigenVector const& y)
267{
268 w.getRawVector().noalias() = x.getRawVector().binaryExpr(
269 y.getRawVector(),
270 [](auto const x, auto const y) { return y == 0 ? 0.0 : x / y; });
271}
272
273// Explicit specialization
274// Computes the Manhattan norm of x
275template <>
276double norm1(EigenVector const& x)
277{
278 return x.getRawVector().lpNorm<1>();
279}
280
281// Explicit specialization
282// Euclidean norm
283template <>
284double norm2(EigenVector const& x)
285{
286 return x.getRawVector().norm();
287}
288
289// Explicit specialization
290// Computes the Maximum norm of x
291template <>
292double normMax(EigenVector const& x)
293{
294 return x.getRawVector().lpNorm<Eigen::Infinity>();
295}
296
297// Matrix
298
299void copy(EigenMatrix const& A, EigenMatrix& B)
300{
301 B = A;
302}
303
304// A = a*A
305void scale(EigenMatrix& A, double const a)
306{
307 // TODO: does that break anything?
308 A.getRawMatrix() *= a;
309}
310
311// Y = a*Y + X
312void aypx(EigenMatrix& Y, double const a, EigenMatrix const& X)
313{
314 // TODO: does that break anything?
315 Y.getRawMatrix() = a * Y.getRawMatrix() + X.getRawMatrix();
316}
317
318// Y = a*X + Y
319void axpy(EigenMatrix& Y, double const a, EigenMatrix const& X)
320{
321 // TODO: does that break anything?
322 Y.getRawMatrix() = a * X.getRawMatrix() + Y.getRawMatrix();
323}
324
325// Matrix and Vector
326
327// v3 = A*v1 + v2
328void matMult(EigenMatrix const& A, EigenVector const& x, EigenVector& y)
329{
330 assert(&x != &y);
331 y.getRawVector() = A.getRawMatrix() * x.getRawVector();
332}
333
334// v3 = A*v1 + v2
335void matMultAdd(EigenMatrix const& A, EigenVector const& v1,
336 EigenVector const& v2, EigenVector& v3)
337{
338 assert(&v1 != &v3);
339 // TODO: does that break anything?
340 v3.getRawVector() =
341 v2.getRawVector() + A.getRawMatrix() * v1.getRawVector();
342}
343
344void linearSysNormalize(EigenMatrix const& A, EigenMatrix& new_A,
345 EigenVector const& b, EigenVector& new_b)
346{
347 // make sure that new_A and new_b are not the same memory
348 assert(&A != &new_A);
349 assert(&b != &new_b);
350
351 if (A.getRawMatrix().rows() == A.getRawMatrix().cols())
352 {
353 WARN(
354 "The number of rows and columns are the same for the LHS matrix."
355 "Are you sure you still need to normalize the LHS matrix and RHS "
356 "vector? ");
357 }
358
359 new_b.getRawVector() = A.getRawMatrix().transpose() * b.getRawVector();
360 new_A.getRawMatrix() = A.getRawMatrix().transpose() * A.getRawMatrix();
361}
362
363void finalizeAssembly(EigenMatrix& x)
364{
365 x.getRawMatrix().makeCompressed();
366}
367
368void finalizeAssembly(EigenVector& /*x*/) {}
369
370} // namespace LinAlg
371
372} // namespace MathLib
373
374#endif
#define OGS_FATAL(...)
Definition Error.h:26
void WARN(fmt::format_string< Args... > fmt, Args &&... args)
Definition Logging.h:40
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.
double norm(MatrixOrVector const &x, MathLib::VecNormType type)
Definition LinAlg.h:96
void linearSysNormalize(PETScMatrix const &, PETScMatrix &, PETScVector const &, PETScVector &)
Definition LinAlg.cpp:170
void finalizeAssembly(PETScMatrix &A)
Definition LinAlg.cpp:198
double norm1(PETScVector const &x)
Definition LinAlg.cpp:90
double normMax(PETScVector const &x)
Definition LinAlg.cpp:110
void copy(PETScVector const &x, PETScVector &y)
Definition LinAlg.cpp:37
void componentwiseDivide(PETScVector &w, PETScVector const &x, PETScVector const &y)
Definition LinAlg.cpp:81
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:149
void matMultAdd(PETScMatrix const &A, PETScVector const &v1, PETScVector const &v2, PETScVector &v3)
Definition LinAlg.cpp:159
double norm2(PETScVector const &x)
Definition LinAlg.cpp:100
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