MueLu  Version of the Day
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
MueLu_GMRESSolver_def.hpp
Go to the documentation of this file.
1 // @HEADER
2 // *****************************************************************************
3 // MueLu: A package for multigrid based preconditioning
4 //
5 // Copyright 2012 NTESS and the MueLu contributors.
6 // SPDX-License-Identifier: BSD-3-Clause
7 // *****************************************************************************
8 // @HEADER
9 
10 #ifndef MUELU_GMRESSOLVER_DEF_HPP
11 #define MUELU_GMRESSOLVER_DEF_HPP
12 
13 #include <Teuchos_LAPACK.hpp>
14 
15 #include <Xpetra_MatrixFactory.hpp>
16 #include <Xpetra_MatrixFactory2.hpp>
17 #include <Xpetra_MatrixMatrix.hpp>
18 #include <Xpetra_IO.hpp>
19 
21 
22 #include "MueLu_Constraint.hpp"
23 #include "MueLu_Monitor.hpp"
24 #include "MueLu_Utilities.hpp"
25 
26 namespace MueLu {
27 
28 template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
30  : nIts_(Its) {}
31 
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];
37  v[i] = w1;
38  v[i + 1] = w2;
39  }
40 }
41 
42 template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
43 void GMRESSolver<Scalar, LocalOrdinal, GlobalOrdinal, Node>::Iterate(const Matrix& Aref, const Constraint& C, const Matrix& P0, RCP<Matrix>& finalP) const {
44  PrintMonitor m(*this, "GMRES iterations");
45 
46  finalP = MatrixFactory2::BuildCopy(rcpFromRef(P0));
47  if (nIts_ == 0)
48  return;
49 
51  "For now, solving Hessenberg system works only for a single iteration");
52 
54 
55  RCP<const Matrix> A = rcpFromRef(Aref);
56 
57  const auto rowMap = A->getRowMap();
59 
60  // FIXME: Don't know why, but in the MATLAB code we have D = I. Follow that for now.
61 #if 0
62  A->getLocalDiagCopy(*invDiag);
63  invDiag->reciprocal(*invDiag);
64 #else
65  invDiag->putScalar(one);
66 #endif
67 
68  Teuchos::FancyOStream& mmfancy = this->GetOStream(Statistics2);
69 
70  // Initial P0 would only be used for multiplication
71  RCP<Matrix> X = rcp_const_cast<Matrix>(rcpFromRef(P0)), tmpAP, newV;
72  std::vector<RCP<Matrix> > V(nIts_ + 1);
73 
74  // T is used only for projecting onto
75  RCP<CrsMatrix> T_ = CrsMatrixFactory::Build(C.GetPattern());
76  T_->fillComplete(P0.getDomainMap(), P0.getRangeMap());
77  RCP<Matrix> T = rcp(new CrsMatrixWrap(T_));
78 
79  SC rho;
80  {
81  // R_0 = -D^{-1}*A*X_0
82  // V_0 = R_0 / ||R_0||_F
83  tmpAP = MatrixMatrix::Multiply(*A, false, *X, false, mmfancy, true /*doFillComplete*/, true /*optimizeStorage*/);
84  C.Apply(*tmpAP, *T);
85 
86  V[0] = MatrixFactory2::BuildCopy(T);
87  V[0]->leftScale(*invDiag);
88 
89  rho = sqrt(Utilities::Frobenius(*V[0], *V[0]));
90 
91  V[0]->scale(-one / rho);
92  }
93 
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);
98  g[0] = rho;
99 
100 #define I(i, j) ((i) + (j) * (nIts_ + 1)) // column ordering
101  for (size_t i = 0; i < nIts_; i++) {
102  // V_{i+1} = D^{-1}*A*V_i
103  tmpAP = MatrixMatrix::Multiply(*A, false, *V[i], false, mmfancy, true /*doFillComplete*/, true /*optimizeStorage*/);
104  C.Apply(*tmpAP, *T);
105 
106  V[i + 1] = MatrixFactory2::BuildCopy(T);
107  V[i + 1]->leftScale(*invDiag);
108 
109  // Update Hessenberg matrix
110  for (size_t j = 0; j <= i; j++) {
111  h[I(j, i)] = Utilities::Frobenius(*V[i + 1], *V[j]);
112 
113  // V_{i+1} = V_{i+1} - h(j,i+1)*V_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());
118  V[i + 1].swap(newV);
119 #else
120  // FIXME: this does not work now. Fails with the following exception:
121  // what(): ../../packages/tpetra/core/ext/TpetraExt_MatrixMatrix_def.hpp:408:
122  //
123  // Throw number = 1
124  //
125  // Throw test that evaluated to true: B.isLocallyIndexed()
126  //
127  // TpetraExt::MatrixMatrix::Add(): ERROR, input matrix B must not be locally indexed
128  MatrixMatrix::TwoMatrixAdd(*V[j], false, -h[I(j, i)], *V[i + 1], one);
129 #endif
130  }
131  h[I(i + 1, i)] = sqrt(Utilities::Frobenius(*V[i + 1], *V[i + 1]));
132 
133  // NOTE: potentially we'll need some reorthogonalization code here
134  // The matching MATLAB code is
135  // normav = norm(v.num(k+1).matrix, 'fro');
136  // normav2 = h(k+1,k);
137  // if (reorth == -1 && normav + .001*normav2 == normav)
138  // for j = 1:k
139  // hr = v(:,j)'*v(:,k+1); % hr=v(:,k+1)'*v(:,j);
140  // h(j,k) = h(j,k)+hr;
141  // v(:,k+1) = v(:,k+1)-hr*v(:,j);
142  // end
143  // h(k+1,k) = norm(v(:,k+1));
144  // end
145 
146  // Check for nonsymmetric case
147  if (h[I(i + 1, i)] != zero) {
148  // Normalize V_i
149  V[i + 1]->scale(one / h[I(i + 1, i)]);
150  }
151 
152  if (i > 0)
153  givapp(&c[0], &s[0], &h[I(0, i)], i); // Due to column ordering &h[...] is a column
154 
155  SC nu = sqrt(h[I(i, i)] * h[I(i, i)] + h[I(i + 1, i)] * h[I(i + 1, i)]);
156  if (nu != zero) {
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;
161 
162  givapp(&c[i], &s[i], &g[i], 1);
163  }
164  }
165 
166  // Solve Hessenberg system
167  // y = solve(H, \rho e_1)
168  std::vector<SC> y(nIts_);
169  if (nIts_ == 1) {
170  y[0] = g[0] / h[I(0, 0)];
171  }
172 #undef I
173 
174  // Compute final
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());
180  finalP.swap(newV);
181 #else
182  MatrixMatrix::TwoMatrixAdd(*V[i], false, y[i], *finalP, one);
183 #endif
184  }
185 }
186 
187 } // namespace MueLu
188 
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
#define I(i, j)
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.
Scalar SC
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.