MueLu Version of the Day
Loading...
Searching...
No Matches
MueLu_GMRESSolver_def.hpp
Go to the documentation of this file.
1// @HEADER
2//
3// ***********************************************************************
4//
5// MueLu: A package for multigrid based preconditioning
6// Copyright 2012 Sandia Corporation
7//
8// Under the terms of Contract DE-AC04-94AL85000 with Sandia Corporation,
9// the U.S. Government retains certain rights in this software.
10//
11// Redistribution and use in source and binary forms, with or without
12// modification, are permitted provided that the following conditions are
13// met:
14//
15// 1. Redistributions of source code must retain the above copyright
16// notice, this list of conditions and the following disclaimer.
17//
18// 2. Redistributions in binary form must reproduce the above copyright
19// notice, this list of conditions and the following disclaimer in the
20// documentation and/or other materials provided with the distribution.
21//
22// 3. Neither the name of the Corporation nor the names of the
23// contributors may be used to endorse or promote products derived from
24// this software without specific prior written permission.
25//
26// THIS SOFTWARE IS PROVIDED BY SANDIA CORPORATION "AS IS" AND ANY
27// EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
28// IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR
29// PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL SANDIA CORPORATION OR THE
30// CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL,
31// EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO,
32// PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR
33// PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF
34// LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING
35// NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
36// SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
37//
38// Questions? Contact
39// Jonathan Hu (jhu@sandia.gov)
40// Andrey Prokopenko (aprokop@sandia.gov)
41// Ray Tuminaro (rstumin@sandia.gov)
42//
43// ***********************************************************************
44//
45// @HEADER
46#ifndef MUELU_GMRESSOLVER_DEF_HPP
47#define MUELU_GMRESSOLVER_DEF_HPP
48
49#include <Teuchos_LAPACK.hpp>
50
51#include <Xpetra_MatrixFactory.hpp>
52#include <Xpetra_MatrixMatrix.hpp>
53#include <Xpetra_IO.hpp>
54
55#include "MueLu_GMRESSolver.hpp"
56
57#include "MueLu_Constraint.hpp"
58#include "MueLu_Monitor.hpp"
59#include "MueLu_Utilities.hpp"
60
61
62namespace MueLu {
63
64 template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
68
69 template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
71 for (int i = 0; i < k; i++) {
72 SC w1 = c[i]*v[i] - s[i]*v[i+1];
73 SC w2 = s[i]*v[i] + c[i]*v[i+1];
74 v[i] = w1;
75 v[i+1] = w2;
76 }
77 }
78
79 template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
80 void GMRESSolver<Scalar, LocalOrdinal, GlobalOrdinal, Node>::Iterate(const Matrix& Aref, const Constraint& C, const Matrix& P0, RCP<Matrix>& finalP) const {
81 PrintMonitor m(*this, "GMRES iterations");
82
83 finalP = MatrixFactory2::BuildCopy(rcpFromRef(P0));
84 if (nIts_ == 0)
85 return;
86
87 TEUCHOS_TEST_FOR_EXCEPTION(nIts_ > 1, Exceptions::RuntimeError,
88 "For now, solving Hessenberg system works only for a single iteration");
89
90 SC one = Teuchos::ScalarTraits<SC>::one(), zero = Teuchos::ScalarTraits<SC>::zero();
91
92 RCP<const Matrix> A = rcpFromRef(Aref);
93 //bool useTpetra = (A->getRowMap()->lib() == Xpetra::UseTpetra);
94
95 // FIXME: Don't know why, but in the MATLAB code we have D = I. Follow that for now.
96#if 0
97 ArrayRCP<const SC> D = Utilities::GetMatrixDiagonal(*A);
98#else
99 ArrayRCP<const SC> D(A->getLocalNumRows(), one);
100#endif
101
102 Teuchos::FancyOStream& mmfancy = this->GetOStream(Statistics2);
103
104 // Initial P0 would only be used for multiplication
105 RCP<Matrix> X = rcp_const_cast<Matrix>(rcpFromRef(P0)), tmpAP, newV;
106 std::vector<RCP<Matrix> > V(nIts_+1);
107
108 // T is used only for projecting onto
109 RCP<CrsMatrix> T_ = CrsMatrixFactory::Build(C.GetPattern());
110 T_->fillComplete(P0.getDomainMap(), P0.getRangeMap());
111 RCP<Matrix> T = rcp(new CrsMatrixWrap(T_));
112
113 SC rho;
114 {
115 // R_0 = -D^{-1}*A*X_0
116 // V_0 = R_0 / ||R_0||_F
117 tmpAP = MatrixMatrix::Multiply(*A, false, *X, false, mmfancy, true/*doFillComplete*/, true/*optimizeStorage*/);
118 C.Apply(*tmpAP, *T);
119
120 V[0] = MatrixFactory2::BuildCopy(T);
121 Utilities::MyOldScaleMatrix(*V[0], D, true/*doInverse*/, true/*doFillComplete*/, false/*doOptimizeStorage*/);
122
123 rho = sqrt(Utilities::Frobenius(*V[0], *V[0]));
124
125 V[0]->scale(-one/rho);
126 }
127
128 std::vector<SC> h((nIts_+1) * (nIts_+1));
129 std::vector<SC> c(nIts_+1, 0.0);
130 std::vector<SC> s(nIts_+1, 0.0);
131 std::vector<SC> g(nIts_+1, 0.0);
132 g[0] = rho;
133
134#define I(i,j) ((i) + (j)*(nIts_+1)) // column ordering
135 for (size_t i = 0; i < nIts_; i++) {
136 // V_{i+1} = D^{-1}*A*V_i
137 tmpAP = MatrixMatrix::Multiply(*A, false, *V[i], false, mmfancy, true/*doFillComplete*/, true/*optimizeStorage*/);
138 C.Apply(*tmpAP, *T);
139
140 V[i+1] = MatrixFactory2::BuildCopy(T);
141 Utilities::MyOldScaleMatrix(*V[i+1], D, true/*doInverse*/, true/*doFillComplete*/, false/*doOptimizeStorage*/);
142
143 // Update Hessenberg matrix
144 for (size_t j = 0; j <= i; j++) {
145 h[I(j,i)] = Utilities::Frobenius(*V[i+1], *V[j]);
146
147 // V_{i+1} = V_{i+1} - h(j,i+1)*V_j
148#ifndef TWO_ARG_MATRIX_ADD
149 newV = Teuchos::null;
150 MatrixMatrix::TwoMatrixAdd(*V[j], false, -h[I(j,i)], *V[i+1], false, one, newV, mmfancy);
151 newV->fillComplete(V[i+1]->getDomainMap(), V[i+1]->getRangeMap());
152 V[i+1].swap(newV);
153#else
154 // FIXME: this does not work now. Fails with the following exception:
155 // what(): ../../packages/tpetra/core/ext/TpetraExt_MatrixMatrix_def.hpp:408:
156 //
157 // Throw number = 1
158 //
159 // Throw test that evaluated to true: B.isLocallyIndexed()
160 //
161 // TpetraExt::MatrixMatrix::Add(): ERROR, input matrix B must not be locally indexed
162 MatrixMatrix::TwoMatrixAdd(*V[j], false, -h[I(j,i)], *V[i+1], one);
163#endif
164 }
165 h[I(i+1,i)] = sqrt(Utilities::Frobenius(*V[i+1], *V[i+1]));
166
167 // NOTE: potentially we'll need some reorthogonalization code here
168 // The matching MATLAB code is
169 // normav = norm(v.num(k+1).matrix, 'fro');
170 // normav2 = h(k+1,k);
171 // if (reorth == -1 && normav + .001*normav2 == normav)
172 // for j = 1:k
173 // hr = v(:,j)'*v(:,k+1); % hr=v(:,k+1)'*v(:,j);
174 // h(j,k) = h(j,k)+hr;
175 // v(:,k+1) = v(:,k+1)-hr*v(:,j);
176 // end
177 // h(k+1,k) = norm(v(:,k+1));
178 // end
179
180 // Check for nonsymmetric case
181 if (h[I(i+1,i)] != zero) {
182 // Normalize V_i
183 V[i+1]->scale(one/h[I(i+1,i)]);
184 }
185
186 if (i > 0)
187 givapp(&c[0], &s[0], &h[I(0,i)], i); // Due to column ordering &h[...] is a column
188
189 SC nu = sqrt(h[I(i,i)]*h[I(i,i)] + h[I(i+1,i)]*h[I(i+1,i)]);
190 if (nu != zero) {
191 c[i] = h[I(i, i)] / nu;
192 s[i] = -h[I(i+1,i)] / nu;
193 h[I(i,i)] = c[i] * h[I(i,i)] - s[i] * h[I(i+1,i)];
194 h[I(i+1,i)] = zero;
195
196 givapp(&c[i], &s[i], &g[i], 1);
197 }
198 }
199
200 // Solve Hessenberg system
201 // y = solve(H, \rho e_1)
202 std::vector<SC> y(nIts_);
203 if (nIts_ == 1) {
204 y[0] = g[0] / h[I(0,0)];
205 }
206#undef I
207
208 // Compute final
209 for (size_t i = 0; i < nIts_; i++) {
210#ifndef TWO_ARG_MATRIX_ADD
211 newV = Teuchos::null;
212 MatrixMatrix::TwoMatrixAdd(*V[i], false, y[i], *finalP, false, one, newV, mmfancy);
213 newV->fillComplete(finalP->getDomainMap(), finalP->getRangeMap());
214 finalP.swap(newV);
215#else
216 MatrixMatrix::TwoMatrixAdd(*V[i], false, y[i], *finalP, one);
217#endif
218 }
219 }
220
221} // namespace MueLu
222
223#endif //ifndef MUELU_GMRESSOLVER_DECL_HPP
#define I(i, j)
MueLu::DefaultScalar Scalar
Constraint space information for the potential prolongator.
RCP< const CrsGraph > GetPattern() const
void Apply(const Matrix &P, Matrix &Projected) const
Apply constraint.
Exception throws to report errors in the internal logical of the program.
void givapp(SC *c, SC *s, SC *v, int k) const
void Iterate(const Matrix &A, const Constraint &C, const Matrix &P0, RCP< Matrix > &P) const
Iterate.
static Teuchos::ArrayRCP< Scalar > GetMatrixDiagonal(const Xpetra::Matrix< Scalar, LocalOrdinal, GlobalOrdinal, Node > &A)
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)
static Scalar Frobenius(const Xpetra::Matrix< Scalar, LocalOrdinal, GlobalOrdinal, Node > &A, const Xpetra::Matrix< Scalar, LocalOrdinal, GlobalOrdinal, Node > &B)
Namespace for MueLu classes and methods.
@ Statistics2
Print even more statistics.