MueLu  Version of the Day
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
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  // Note: this function matrix notations follow Saad's "Iterative methods", ed. 2, pg. 246
71  // So, X is the unknown prolongator, P's are conjugate directions, Z's are preconditioned P's
72  PrintMonitor m(*this, "CG iterations");
73 
74  if (nIts_ == 0) {
75  finalP = MatrixFactory2::BuildCopy(rcpFromRef(P0));
76  return;
77  }
78 
79  RCP<const Matrix> A = rcpFromRef(Aref);
81  bool useTpetra = (A->getRowMap()->lib() == Xpetra::UseTpetra);
82 
83  Teuchos::FancyOStream& mmfancy = this->GetOStream(Statistics2);
84 
86 
87  RCP<Matrix> X, P, R, Z, AP;
88  RCP<Matrix> newX, tmpAP;
89 #ifndef TWO_ARG_MATRIX_ADD
90  RCP<Matrix> newR, newP;
91 #endif
92 
93  SC oldRZ, newRZ, alpha, beta, app;
94 
95  // T is used only for projecting onto
96  RCP<CrsMatrix> T_ = CrsMatrixFactory::Build(C.GetPattern());
97  T_->fillComplete(P0.getDomainMap(), P0.getRangeMap());
98  RCP<Matrix> T = rcp(new CrsMatrixWrap(T_));
99 
100  // Initial P0 would only be used for multiplication
101  X = rcp_const_cast<Matrix>(rcpFromRef(P0));
102 
103  tmpAP = MatrixMatrix::Multiply(*A, false, *X, false, mmfancy, true/*doFillComplete*/, true/*optimizeStorage*/);
104  C.Apply(*tmpAP, *T);
105 
106  // R_0 = -A*X_0
108 
109  R->resumeFill();
110  R->scale(-one);
111  R->fillComplete(R->getDomainMap(), R->getRangeMap());
112 
113  // Z_0 = M^{-1}R_0
115  Utilities::MyOldScaleMatrix(*Z, D, true, true, false);
116 
117  // P_0 = Z_0
119 
120  oldRZ = Utilities::Frobenius(*R, *Z);
121 
122  for (size_t i = 0; i < nIts_; i++) {
123  // AP = constrain(A*P)
124  if (i == 0 || useTpetra) {
125  // Construct the MxM pattern from scratch
126  // This is done by default for Tpetra as the three argument version requires tmpAP
127  // to *not* be locally indexed which defeats the purpose
128  // TODO: need a three argument Tpetra version which allows reuse of already fill-completed matrix
129  tmpAP = MatrixMatrix::Multiply(*A, false, *P, false, mmfancy, true/*doFillComplete*/, true/*optimizeStorage*/);
130  } else {
131  // Reuse the MxM pattern
132  tmpAP = MatrixMatrix::Multiply(*A, false, *P, false, tmpAP, mmfancy, true/*doFillComplete*/, true/*optimizeStorage*/);
133  }
134  C.Apply(*tmpAP, *T);
135  AP = T;
136 
137  app = Utilities::Frobenius(*AP, *P);
139  // It happens, for instance, if P = 0
140  // For example, if we use TentativePFactory for both nonzero pattern and initial guess
141  // I think it might also happen because of numerical breakdown, but we don't test for that yet
142  if (i == 0)
143  X = MatrixFactory2::BuildCopy(rcpFromRef(P0));
144  break;
145  }
146 
147  // alpha = (R_i, Z_i)/(A*P_i, P_i)
148  alpha = oldRZ / app;
149  this->GetOStream(Runtime1,1) << "alpha = " << alpha << std::endl;
150 
151  // X_{i+1} = X_i + alpha*P_i
152 #ifndef TWO_ARG_MATRIX_ADD
153  newX = Teuchos::null;
154  MatrixMatrix::TwoMatrixAdd(*P, false, alpha, *X, false, one, newX, mmfancy);
155  newX->fillComplete(P0.getDomainMap(), P0.getRangeMap());
156  X.swap(newX);
157 #else
158  MatrixMatrix::TwoMatrixAdd(*P, false, alpha, *X, one);
159 #endif
160 
161  if (i == nIts_ - 1)
162  break;
163 
164  // R_{i+1} = R_i - alpha*A*P_i
165 #ifndef TWO_ARG_MATRIX_ADD
166  newR = Teuchos::null;
167  MatrixMatrix::TwoMatrixAdd(*AP, false, -alpha, *R, false, one, newR, mmfancy);
168  newR->fillComplete(P0.getDomainMap(), P0.getRangeMap());
169  R.swap(newR);
170 #else
171  MatrixMatrix::TwoMatrixAdd(*AP, false, -alpha, *R, one);
172 #endif
173 
174  // Z_{i+1} = M^{-1} R_{i+1}
175  Z = MatrixFactory2::BuildCopy(R);
176  Utilities::MyOldScaleMatrix(*Z, D, true, true, false);
177 
178  // beta = (R_{i+1}, Z_{i+1})/(R_i, Z_i)
179  newRZ = Utilities::Frobenius(*R, *Z);
180  beta = newRZ / oldRZ;
181 
182  // P_{i+1} = Z_{i+1} + beta*P_i
183 #ifndef TWO_ARG_MATRIX_ADD
184  newP = Teuchos::null;
185  MatrixMatrix::TwoMatrixAdd(*P, false, beta, *Z, false, one, newP, mmfancy);
186  newP->fillComplete(P0.getDomainMap(), P0.getRangeMap());
187  P.swap(newP);
188 #else
189  MatrixMatrix::TwoMatrixAdd(*Z, false, one, *P, beta);
190 #endif
191 
192  oldRZ = newRZ;
193  }
194 
195  finalP = X;
196  }
197 
198 } // namespace MueLu
199 
200 #endif //ifndef MUELU_CGSOLVER_DECL_HPP
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)
Constraint space information for the potential prolongator.
void swap(RCP< T > &r_ptr)
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.
TEUCHOS_DEPRECATED RCP< T > rcp(T *p, Dealloc_T dealloc, bool owns_mem)
static RCP< Xpetra::Matrix< Scalar, LocalOrdinal, GlobalOrdinal, Node > > BuildCopy(const RCP< const Xpetra::Matrix< Scalar, LocalOrdinal, GlobalOrdinal, Node > > A)
void Apply(const Matrix &P, Matrix &Projected) const
Apply constraint.
Scalar SC
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)