MueLu  Version of the Day
MueLu_CGSolver_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_CGSOLVER_DEF_HPP
47 #define MUELU_CGSOLVER_DEF_HPP
48 
49 #include <Xpetra_MatrixFactory.hpp>
50 #include <Xpetra_MatrixMatrix.hpp>
51 
52 #include "MueLu_Utilities.hpp"
53 #include "MueLu_Constraint.hpp"
54 #include "MueLu_Monitor.hpp"
55 
56 
57 #include "MueLu_CGSolver.hpp"
58 
59 
60 
61 namespace MueLu {
62 
63  template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
65  : nIts_(Its)
66  { }
67 
68  template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
69  void CGSolver<Scalar, LocalOrdinal, GlobalOrdinal, Node>::Iterate(const Matrix& Aref, const Constraint& C, const Matrix& P0, RCP<Matrix>& finalP) const {
70  PrintMonitor m(*this, "CG iterations");
71 
72  if (nIts_ == 0) {
73  finalP = Xpetra::MatrixFactory2<Scalar, LocalOrdinal, GlobalOrdinal, Node>::BuildCopy(rcpFromRef(P0));
74  return;
75  }
76 
77  // Note: this function matrix notations follow Saad's "Iterative methods", ed. 2, pg. 246
78  // So, X is the unknown prolongator, P's are conjugate directions, Z's are preconditioned P's
79  RCP<const Matrix> A = rcpFromRef(Aref);
80 
81  RCP<Matrix> X, P, R, Z, AP;
82  RCP<Matrix> newX, tmpAP;
83 #ifndef TWO_ARG_MATRIX_ADD
84  RCP<Matrix> newR, newP;
85 #endif
86 
87  SC oldRZ, newRZ, alpha, beta, app;
88 
89  bool useTpetra = (A->getRowMap()->lib() == Xpetra::UseTpetra);
90 
91  Teuchos::FancyOStream& mmfancy = this->GetOStream(Statistics2);
92 
93  // T is used only for projecting onto
94  RCP<CrsMatrix> T_ = Xpetra::CrsMatrixFactory<Scalar, LocalOrdinal, GlobalOrdinal, Node>::Build(C.GetPattern());
95  T_->fillComplete(P0.getDomainMap(), P0.getRangeMap());
96  RCP<Matrix> T = rcp(new CrsMatrixWrap(T_));
97 
98  SC one = Teuchos::ScalarTraits<SC>::one();
99 
100  Teuchos::ArrayRCP<const SC> D = Utilities::GetMatrixDiagonal(*A);
101 
102  // Initial P0 would only be used for multiplication
103  X = rcp_const_cast<Matrix>(rcpFromRef(P0));
104 
105  bool doFillComplete = true;
106  // bool optimizeStorage = false;
107  bool optimizeStorage = true;
108 
109  tmpAP = MatrixMatrix::Multiply(*A, false, *X, false, mmfancy, doFillComplete, optimizeStorage);
110  C.Apply(*tmpAP, *T);
111 
112  // R_0 = -A*X_0
113  R = Xpetra::MatrixFactory2<Scalar, LocalOrdinal, GlobalOrdinal, Node>::BuildCopy(T);
114 #ifdef HAVE_MUELU_TPETRA
115 #ifdef HAVE_MUELU_TPETRA_INST_INT_INT
116  // TAW: Oct 16 2015: MueLu::Utilities returns the Tpetra::CrsMatrix object which would not be instantiated!
117  // Catching this in Op2NonConstTpetraCrs is not possible as this does not affect the return type
118  // Tpetra::CrsMatrix!
119  if (useTpetra)
120  Utilities::Op2NonConstTpetraCrs(R)->resumeFill();
121 #else
122  this->GetOStream(Warnings0) << "WARNING: MueLu_CGSolver: calling Xpetra::CrsMatrix::resumeFill instead of Tpetra::CrsMatrix::resumeFill. The results should be verified in this case." << std::endl;
123  R->resumeFill();
124 #endif
125 #endif
126  R->scale(-one);
127  if (useTpetra)
128  R->fillComplete(R->getDomainMap(), R->getRangeMap());
129 
130  // Z_0 = M^{-1}R_0
131  Z = Xpetra::MatrixFactory2<Scalar, LocalOrdinal, GlobalOrdinal, Node>::BuildCopy(R);
132  Utilities::MyOldScaleMatrix(*Z, D, true, true, false);
133 
134  // P_0 = Z_0
135  P = Xpetra::MatrixFactory2<Scalar, LocalOrdinal, GlobalOrdinal, Node>::BuildCopy(Z);
136 
137  oldRZ = Utilities::Frobenius(*R, *Z);
138 
139  for (size_t k = 0; k < nIts_; k++) {
140  // AP = constrain(A*P)
141  if (k == 0 || useTpetra) {
142  // Construct the MxM pattern from scratch
143  // This is done by default for Tpetra as the three argument version requires tmpAP
144  // to *not* be locally indexed which defeats the purpose
145  // TODO: need a three argument Tpetra version which allows reuse of already fill-completed matrix
146  tmpAP = MatrixMatrix::Multiply(*A, false, *P, false, mmfancy, doFillComplete, optimizeStorage);
147  } else {
148  // Reuse the MxM pattern
149  tmpAP = MatrixMatrix::Multiply(*A, false, *P, false, tmpAP, mmfancy, doFillComplete, optimizeStorage);
150  }
151  C.Apply(*tmpAP, *T);
152  AP = T;
153 
154  app = Utilities::Frobenius(*AP, *P);
155  if (Teuchos::ScalarTraits<SC>::magnitude(app) < Teuchos::ScalarTraits<SC>::sfmin()) {
156  // It happens, for instance, if P = 0
157  // For example, if we use TentativePFactory for both nonzero pattern and initial guess
158  // I think it might also happen because of numerical breakdown, but we don't test for that yet
159  if (k == 0)
160  X = Xpetra::MatrixFactory2<Scalar, LocalOrdinal, GlobalOrdinal, Node>::BuildCopy(rcpFromRef(P0));
161  break;
162  }
163 
164  // alpha = (R_k, Z_k)/(A*P_k, P_k)
165  alpha = oldRZ / app;
166  this->GetOStream(Runtime1,1) << "alpha = " << alpha << std::endl;
167 
168  // X_{k+1} = X_k + alpha*P_k
169 #ifndef TWO_ARG_MATRIX_ADD
170  newX = Teuchos::null;
171  Xpetra::MatrixMatrix<Scalar, LocalOrdinal, GlobalOrdinal, Node>::TwoMatrixAdd(*P, false, alpha, *X, false, Teuchos::ScalarTraits<Scalar>::one(), newX, mmfancy);
172  newX->fillComplete(P0.getDomainMap(), P0.getRangeMap());
173  X.swap(newX);
174 #else
175  MatrixMatrix::TwoMatrixAdd(*P, false, alpha, *X, one);
176 #endif
177 
178  if (k == nIts_ - 1)
179  break;
180 
181  // R_{k+1} = R_k - alpha*A*P_k
182 #ifndef TWO_ARG_MATRIX_ADD
183  newR = Teuchos::null;
184  Xpetra::MatrixMatrix<Scalar, LocalOrdinal, GlobalOrdinal, Node>::TwoMatrixAdd(*AP, false, -alpha, *R, false, Teuchos::ScalarTraits<Scalar>::one(), newR, mmfancy);
185  newR->fillComplete(P0.getDomainMap(), P0.getRangeMap());
186  R.swap(newR);
187 #else
188  MatrixMatrix::TwoMatrixAdd(*AP, false, -alpha, *R, one);
189 #endif
190 
191  // Z_{k+1} = M^{-1} R_{k+1}
192  Z = Xpetra::MatrixFactory2<Scalar, LocalOrdinal, GlobalOrdinal, Node>::BuildCopy(R);
193  Utilities::MyOldScaleMatrix(*Z, D, true, true, false);
194 
195  // beta = (R_{k+1}, Z_{k+1})/(R_k, Z_k)
196  newRZ = Utilities::Frobenius(*R, *Z);
197  beta = newRZ / oldRZ;
198 
199  // P_{k+1} = Z_{k+1} + beta*P_k
200 #ifndef TWO_ARG_MATRIX_ADD
201  newP = Teuchos::null;
202  MatrixMatrix::TwoMatrixAdd(*P, false, beta, *Z, false, Teuchos::ScalarTraits<Scalar>::one(), newP, mmfancy);
203  newP->fillComplete(P0.getDomainMap(), P0.getRangeMap());
204  P.swap(newP);
205 #else
206  MatrixMatrix::TwoMatrixAdd(*Z, false, one, *P, beta);
207 #endif
208 
209  oldRZ = newRZ;
210  }
211 
212  finalP = X;
213  }
214 
215 } // namespace MueLu
216 
217 #endif //ifndef MUELU_CGSOLVER_DECL_HPP
static Teuchos::ArrayRCP< Scalar > GetMatrixDiagonal(const Xpetra::Matrix< Scalar, LocalOrdinal, GlobalOrdinal, Node > &A)
Important warning messages (one line)
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)
Constraint space information for the potential prolongator.
Namespace for MueLu classes and methods.
RCP< const CrsGraph > GetPattern() const
Print even more statistics.
void Iterate(const Matrix &A, const Constraint &C, const Matrix &P0, RCP< Matrix > &P) const
Iterate.
static RCP< Tpetra::CrsMatrix< Scalar, LocalOrdinal, GlobalOrdinal, Node > > Op2NonConstTpetraCrs(RCP< Xpetra::Matrix< Scalar, LocalOrdinal, GlobalOrdinal, Node > > Op)
Teuchos::FancyOStream & GetOStream(MsgType type, int thisProcRankOnly=0) const
Get an output stream for outputting the input message type.
void Apply(const Matrix &P, Matrix &Projected) const
Apply constraint.
size_t nIts_
Number of performed iterations.
Description of what is happening (more verbose)
static Scalar Frobenius(const Xpetra::Matrix< Scalar, LocalOrdinal, GlobalOrdinal, Node > &A, const Xpetra::Matrix< Scalar, LocalOrdinal, GlobalOrdinal, Node > &B)