46 #ifndef MUELU_GMRESSOLVER_DEF_HPP
47 #define MUELU_GMRESSOLVER_DEF_HPP
57 #include "MueLu_Constraint.hpp"
59 #include "MueLu_Utilities.hpp"
63 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
67 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
69 for (
int i = 0; i < k; i++) {
70 SC w1 = c[i] * v[i] - s[i] * v[i + 1];
71 SC w2 = s[i] * v[i] + c[i] * v[i + 1];
77 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
81 finalP = MatrixFactory2::BuildCopy(rcpFromRef(P0));
86 "For now, solving Hessenberg system works only for a single iteration");
103 RCP<Matrix> X = rcp_const_cast<Matrix>(rcpFromRef(P0)), tmpAP, newV;
104 std::vector<RCP<Matrix> > V(nIts_ + 1);
108 T_->fillComplete(P0.getDomainMap(), P0.getRangeMap());
115 tmpAP = MatrixMatrix::Multiply(*A,
false, *X,
false, mmfancy,
true ,
true );
118 V[0] = MatrixFactory2::BuildCopy(T);
123 V[0]->scale(-one / rho);
126 std::vector<SC> h((nIts_ + 1) * (nIts_ + 1));
127 std::vector<SC> c(nIts_ + 1, 0.0);
128 std::vector<SC> s(nIts_ + 1, 0.0);
129 std::vector<SC> g(nIts_ + 1, 0.0);
132 #define I(i, j) ((i) + (j) * (nIts_ + 1)) // column ordering
133 for (
size_t i = 0; i < nIts_; i++) {
135 tmpAP = MatrixMatrix::Multiply(*A,
false, *V[i],
false, mmfancy,
true ,
true );
138 V[i + 1] = MatrixFactory2::BuildCopy(T);
142 for (
size_t j = 0; j <= i; j++) {
146 #ifndef TWO_ARG_MATRIX_ADD
147 newV = Teuchos::null;
148 MatrixMatrix::TwoMatrixAdd(*V[j],
false, -h[
I(j, i)], *V[i + 1],
false, one, newV, mmfancy);
149 newV->fillComplete(V[i + 1]->getDomainMap(), V[i + 1]->getRangeMap());
160 MatrixMatrix::TwoMatrixAdd(*V[j],
false, -h[
I(j, i)], *V[i + 1], one);
179 if (h[
I(i + 1, i)] != zero) {
181 V[i + 1]->scale(one / h[
I(i + 1, i)]);
185 givapp(&c[0], &s[0], &h[
I(0, i)], i);
187 SC nu = sqrt(h[
I(i, i)] * h[
I(i, i)] + h[
I(i + 1, i)] * h[
I(i + 1, i)]);
189 c[i] = h[
I(i, i)] / nu;
190 s[i] = -h[
I(i + 1, i)] / nu;
191 h[
I(i, i)] = c[i] * h[
I(i, i)] - s[i] * h[
I(i + 1, i)];
192 h[
I(i + 1, i)] = zero;
194 givapp(&c[i], &s[i], &g[i], 1);
200 std::vector<SC> y(nIts_);
202 y[0] = g[0] / h[
I(0, 0)];
207 for (
size_t i = 0; i < nIts_; i++) {
208 #ifndef TWO_ARG_MATRIX_ADD
209 newV = Teuchos::null;
210 MatrixMatrix::TwoMatrixAdd(*V[i],
false, y[i], *finalP,
false, one, newV, mmfancy);
211 newV->fillComplete(finalP->getDomainMap(), finalP->getRangeMap());
214 MatrixMatrix::TwoMatrixAdd(*V[i],
false, y[i], *finalP, one);
221 #endif // ifndef MUELU_GMRESSOLVER_DECL_HPP
void givapp(SC *c, SC *s, SC *v, int k) const
static void MyOldScaleMatrix(Xpetra::Matrix< Scalar, LocalOrdinal, GlobalOrdinal, Node > &Op, const Teuchos::ArrayRCP< const Scalar > &scalingVector, bool doInverse=true, bool doFillComplete=true, bool doOptimizeStorage=true)
#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 Teuchos::ArrayRCP< Scalar > GetMatrixDiagonal_arcp(const Matrix &A)
Extract Matrix Diagonal.
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.