10 #ifndef MUELU_GMRESSOLVER_DEF_HPP
11 #define MUELU_GMRESSOLVER_DEF_HPP
16 #include <Xpetra_MatrixFactory2.hpp>
22 #include "MueLu_Constraint.hpp"
24 #include "MueLu_Utilities.hpp"
28 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
32 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
34 for (
int i = 0; i < k; i++) {
35 SC w1 = c[i] * v[i] - s[i] * v[i + 1];
36 SC w2 = s[i] * v[i] + c[i] * v[i + 1];
42 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
46 finalP = MatrixFactory2::BuildCopy(rcpFromRef(P0));
51 "For now, solving Hessenberg system works only for a single iteration");
57 const auto rowMap = A->getRowMap();
62 A->getLocalDiagCopy(*invDiag);
63 invDiag->reciprocal(*invDiag);
65 invDiag->putScalar(one);
71 RCP<Matrix> X = rcp_const_cast<Matrix>(rcpFromRef(P0)), tmpAP, newV;
72 std::vector<RCP<Matrix> > V(nIts_ + 1);
76 T_->fillComplete(P0.getDomainMap(), P0.getRangeMap());
83 tmpAP = MatrixMatrix::Multiply(*A,
false, *X,
false, mmfancy,
true ,
true );
86 V[0] = MatrixFactory2::BuildCopy(T);
87 V[0]->leftScale(*invDiag);
91 V[0]->scale(-one / rho);
94 std::vector<SC> h((nIts_ + 1) * (nIts_ + 1));
95 std::vector<SC> c(nIts_ + 1, 0.0);
96 std::vector<SC> s(nIts_ + 1, 0.0);
97 std::vector<SC> g(nIts_ + 1, 0.0);
100 #define I(i, j) ((i) + (j) * (nIts_ + 1)) // column ordering
101 for (
size_t i = 0; i < nIts_; i++) {
103 tmpAP = MatrixMatrix::Multiply(*A,
false, *V[i],
false, mmfancy,
true ,
true );
106 V[i + 1] = MatrixFactory2::BuildCopy(T);
107 V[i + 1]->leftScale(*invDiag);
110 for (
size_t j = 0; j <= i; j++) {
114 #ifndef TWO_ARG_MATRIX_ADD
115 newV = Teuchos::null;
116 MatrixMatrix::TwoMatrixAdd(*V[j],
false, -h[
I(j, i)], *V[i + 1],
false, one, newV, mmfancy);
117 newV->fillComplete(V[i + 1]->getDomainMap(), V[i + 1]->getRangeMap());
128 MatrixMatrix::TwoMatrixAdd(*V[j],
false, -h[
I(j, i)], *V[i + 1], one);
147 if (h[
I(i + 1, i)] != zero) {
149 V[i + 1]->scale(one / h[
I(i + 1, i)]);
153 givapp(&c[0], &s[0], &h[
I(0, i)], i);
155 SC nu = sqrt(h[
I(i, i)] * h[
I(i, i)] + h[
I(i + 1, i)] * h[
I(i + 1, i)]);
157 c[i] = h[
I(i, i)] / nu;
158 s[i] = -h[
I(i + 1, i)] / nu;
159 h[
I(i, i)] = c[i] * h[
I(i, i)] - s[i] * h[
I(i + 1, i)];
160 h[
I(i + 1, i)] = zero;
162 givapp(&c[i], &s[i], &g[i], 1);
168 std::vector<SC> y(nIts_);
170 y[0] = g[0] / h[
I(0, 0)];
175 for (
size_t i = 0; i < nIts_; i++) {
176 #ifndef TWO_ARG_MATRIX_ADD
177 newV = Teuchos::null;
178 MatrixMatrix::TwoMatrixAdd(*V[i],
false, y[i], *finalP,
false, one, newV, mmfancy);
179 newV->fillComplete(finalP->getDomainMap(), finalP->getRangeMap());
182 MatrixMatrix::TwoMatrixAdd(*V[i],
false, y[i], *finalP, one);
189 #endif // ifndef MUELU_GMRESSOLVER_DECL_HPP
void givapp(SC *c, SC *s, SC *v, int k) const
#define TEUCHOS_TEST_FOR_EXCEPTION(throw_exception_test, Exception, msg)
Constraint space information for the potential prolongator.
void swap(RCP< T > &r_ptr)
RCP< const CrsGraph > GetPattern() const
Print even more statistics.
TEUCHOS_DEPRECATED RCP< T > rcp(T *p, Dealloc_T dealloc, bool owns_mem)
MueLu::DefaultScalar Scalar
static RCP< Vector > Build(const Teuchos::RCP< const Map > &map, bool zeroOut=true)
void Iterate(const Matrix &A, const Constraint &C, const Matrix &P0, RCP< Matrix > &P) const
Iterate.
void Apply(const Matrix &P, Matrix &Projected) const
Apply constraint.
Exception throws to report errors in the internal logical of the program.
static Scalar Frobenius(const Xpetra::Matrix< Scalar, LocalOrdinal, GlobalOrdinal, Node > &A, const Xpetra::Matrix< Scalar, LocalOrdinal, GlobalOrdinal, Node > &B)
Frobenius inner product of two matrices.