MueLu Version of the Day
Loading...
Searching...
No Matches
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
61namespace 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);
80 ArrayRCP<const SC> D = Utilities::GetMatrixDiagonal(*A);
81 bool useTpetra = (A->getRowMap()->lib() == Xpetra::UseTpetra);
82
83 Teuchos::FancyOStream& mmfancy = this->GetOStream(Statistics2);
84
85 SC one = Teuchos::ScalarTraits<SC>::one();
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
107 R = Xpetra::MatrixFactory2<Scalar, LocalOrdinal, GlobalOrdinal, Node>::BuildCopy(T);
108
109 R->scale(-one);
110
111 // Z_0 = M^{-1}R_0
112 Z = Xpetra::MatrixFactory2<Scalar, LocalOrdinal, GlobalOrdinal, Node>::BuildCopy(R);
113 Utilities::MyOldScaleMatrix(*Z, D, true, true, false);
114
115 // P_0 = Z_0
116 P = Xpetra::MatrixFactory2<Scalar, LocalOrdinal, GlobalOrdinal, Node>::BuildCopy(Z);
117
118 oldRZ = Utilities::Frobenius(*R, *Z);
119
120 for (size_t i = 0; i < nIts_; i++) {
121 // AP = constrain(A*P)
122 if (i == 0 || useTpetra) {
123 // Construct the MxM pattern from scratch
124 // This is done by default for Tpetra as the three argument version requires tmpAP
125 // to *not* be locally indexed which defeats the purpose
126 // TODO: need a three argument Tpetra version which allows reuse of already fill-completed matrix
127 tmpAP = MatrixMatrix::Multiply(*A, false, *P, false, mmfancy, true/*doFillComplete*/, true/*optimizeStorage*/);
128 } else {
129 // Reuse the MxM pattern
130 tmpAP = MatrixMatrix::Multiply(*A, false, *P, false, tmpAP, mmfancy, true/*doFillComplete*/, true/*optimizeStorage*/);
131 }
132 C.Apply(*tmpAP, *T);
133 AP = T;
134
135 app = Utilities::Frobenius(*AP, *P);
136 if (Teuchos::ScalarTraits<SC>::magnitude(app) < Teuchos::ScalarTraits<SC>::sfmin()) {
137 // It happens, for instance, if P = 0
138 // For example, if we use TentativePFactory for both nonzero pattern and initial guess
139 // I think it might also happen because of numerical breakdown, but we don't test for that yet
140 if (i == 0)
141 X = MatrixFactory2::BuildCopy(rcpFromRef(P0));
142 break;
143 }
144
145 // alpha = (R_i, Z_i)/(A*P_i, P_i)
146 alpha = oldRZ / app;
147 this->GetOStream(Runtime1,1) << "alpha = " << alpha << std::endl;
148
149 // X_{i+1} = X_i + alpha*P_i
150#ifndef TWO_ARG_MATRIX_ADD
151 newX = Teuchos::null;
152 MatrixMatrix::TwoMatrixAdd(*P, false, alpha, *X, false, one, newX, mmfancy);
153 newX->fillComplete(P0.getDomainMap(), P0.getRangeMap());
154 X.swap(newX);
155#else
156 MatrixMatrix::TwoMatrixAdd(*P, false, alpha, *X, one);
157#endif
158
159 if (i == nIts_ - 1)
160 break;
161
162 // R_{i+1} = R_i - alpha*A*P_i
163#ifndef TWO_ARG_MATRIX_ADD
164 newR = Teuchos::null;
165 MatrixMatrix::TwoMatrixAdd(*AP, false, -alpha, *R, false, one, newR, mmfancy);
166 newR->fillComplete(P0.getDomainMap(), P0.getRangeMap());
167 R.swap(newR);
168#else
169 MatrixMatrix::TwoMatrixAdd(*AP, false, -alpha, *R, one);
170#endif
171
172 // Z_{i+1} = M^{-1} R_{i+1}
173 Z = MatrixFactory2::BuildCopy(R);
174 Utilities::MyOldScaleMatrix(*Z, D, true, true, false);
175
176 // beta = (R_{i+1}, Z_{i+1})/(R_i, Z_i)
177 newRZ = Utilities::Frobenius(*R, *Z);
178 beta = newRZ / oldRZ;
179
180 // P_{i+1} = Z_{i+1} + beta*P_i
181#ifndef TWO_ARG_MATRIX_ADD
182 newP = Teuchos::null;
183 MatrixMatrix::TwoMatrixAdd(*P, false, beta, *Z, false, one, newP, mmfancy);
184 newP->fillComplete(P0.getDomainMap(), P0.getRangeMap());
185 P.swap(newP);
186#else
187 MatrixMatrix::TwoMatrixAdd(*Z, false, one, *P, beta);
188#endif
189
190 oldRZ = newRZ;
191 }
192
193 finalP = X;
194 }
195
196} // namespace MueLu
197
198#endif //ifndef MUELU_CGSOLVER_DECL_HPP
void Iterate(const Matrix &A, const Constraint &C, const Matrix &P0, RCP< Matrix > &P) const
Iterate.
Constraint space information for the potential prolongator.
RCP< const CrsGraph > GetPattern() const
void Apply(const Matrix &P, Matrix &Projected) const
Apply constraint.
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.
@ Runtime1
Description of what is happening (more verbose)