G+Smo  25.01.0
Geometry + Simulation Modules
 
Loading...
Searching...
No Matches
gsSimplePreconditioners.hpp
Go to the documentation of this file.
1
14namespace gismo
15{
16
17namespace internal
18{
19
20template<typename T>
21void gaussSeidelSweep(const gsSparseMatrix<T> & A, gsMatrix<T>& x, const gsMatrix<T>& f)
22{
23 GISMO_ASSERT( A.rows() == x.rows() && x.rows() == f.rows() && A.cols() == A.rows() && x.cols() == f.cols(),
24 "Dimensions do not match.");
25
26 GISMO_ASSERT( f.cols() == 1, "This operator is only implemented for a single right-hand side." );
27
28 // A is supposed to be symmetric, so it doesn't matter if it's stored in row- or column-major order
29 for (index_t i = 0; i < A.outerSize(); ++i)
30 {
31 T diag = 0;
32 T sum = 0;
33
34 for (typename gsSparseMatrix<T>::InnerIterator it(A,i); it; ++it)
35 {
36 sum += it.value() * x( it.index() ); // compute A.x
37 if (it.index() == i)
38 diag = it.value();
39 }
40
41 x(i) += (f(i) - sum) / diag;
42 }
43}
44
45template<typename T>
46void reverseGaussSeidelSweep(const gsSparseMatrix<T> & A, gsMatrix<T>& x, const gsMatrix<T>& f)
47{
48 GISMO_ASSERT( A.rows() == x.rows() && x.rows() == f.rows() && A.cols() == A.rows() && x.cols() == f.cols(),
49 "Dimensions do not match.");
50
51 GISMO_ASSERT( f.cols() == 1, "This operator is only implemented for a single right-hand side." );
52
53 // A is supposed to be symmetric, so it doesn't matter if it's stored in row- or column-major order
54 for (index_t i = A.outerSize() - 1; i >= 0; --i)
55 {
56 T diag = 0;
57 T sum = 0;
58
59 for (typename gsSparseMatrix<T>::InnerIterator it(A,i); it; ++it)
60 {
61 sum += it.value() * x( it.index() ); // compute A.x
62 if (it.index() == i)
63 diag = it.value();
64 }
65
66 x(i) += (f(i) - sum) / diag;
67 }
68}
69
70} // namespace internal
71
72} // namespace gismo
#define index_t
Definition gsConfig.h:32
#define GISMO_ASSERT(cond, message)
Definition gsDebug.h:89
The G+Smo namespace, containing all definitions for the library.