MueLu Version of the Day
Loading...
Searching...
No Matches
MueLu_SchurComplementFactory_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_SCHURCOMPLEMENTFACTORY_DEF_HPP_
47#define MUELU_SCHURCOMPLEMENTFACTORY_DEF_HPP_
48
49#include <Xpetra_BlockedCrsMatrix.hpp>
50#include <Xpetra_MultiVectorFactory.hpp>
51#include <Xpetra_VectorFactory.hpp>
52#include <Xpetra_MatrixFactory.hpp>
53#include <Xpetra_Matrix.hpp>
54#include <Xpetra_MatrixMatrix.hpp>
55#include <Xpetra_TripleMatrixMultiply.hpp>
56#include <Xpetra_CrsMatrixWrap.hpp>
57#include <Xpetra_BlockedCrsMatrix.hpp>
58#include <Xpetra_CrsMatrix.hpp>
59
60#include "MueLu_Level.hpp"
61#include "MueLu_Monitor.hpp"
62#include "MueLu_Utilities.hpp"
63#include "MueLu_SchurComplementFactory.hpp"
65
66namespace MueLu {
67
68 template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
70 RCP<ParameterList> validParamList = rcp(new ParameterList());
71
72 const SC one = Teuchos::ScalarTraits<SC>::one();
73
74 validParamList->set<RCP<const FactoryBase> >("A" , NoFactory::getRCP(), "Generating factory of the matrix A used for building Schur complement (must be a 2x2 blocked operator)");
75 validParamList->set<RCP<const FactoryBase> >("Ainv" , Teuchos::null, "Generating factory of the inverse matrix used in the Schur complement");
76
77 validParamList->set<SC> ("omega", one, "Scaling parameter in S = A(1,1) - 1/omega A(1,0) Ainv A(0,1)");
78
79 return validParamList;
80 }
81
82 template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
84 Input(currentLevel, "A");
85
86 // Get default or user-given inverse approximation factory
87 RCP<const FactoryBase> AinvFact = GetFactory("Ainv");
88 currentLevel.DeclareInput("Ainv", AinvFact.get(), this);
89 }
90
91 template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
93 FactoryMonitor m(*this, "Build", currentLevel);
94
95 RCP<Matrix> A = Get<RCP<Matrix> >(currentLevel, "A");
96 RCP<BlockedCrsMatrix> bA = rcp_dynamic_cast<BlockedCrsMatrix>(A);
97
98 TEUCHOS_TEST_FOR_EXCEPTION(bA.is_null(), Exceptions::BadCast,
99 "MueLu::SchurComplementFactory::Build: input matrix A is not of type BlockedCrsMatrix!");
100 TEUCHOS_TEST_FOR_EXCEPTION(bA->Rows() != 2 || bA->Cols() != 2, Exceptions::RuntimeError,
101 "MueLu::SchurComplementFactory::Build: input matrix A is a " << bA->Rows() << "x" << bA->Cols() << " block matrix. We expect a 2x2 blocked operator.");
102
103 // Calculate Schur Complement
104 RCP<Matrix> Ainv = currentLevel.Get<RCP<Matrix> >("Ainv", this->GetFactory("Ainv").get());
105 RCP<Matrix> S = ComputeSchurComplement(bA, Ainv);
106
107 GetOStream(Statistics1) << "S has " << S->getGlobalNumRows() << "x" << S->getGlobalNumCols() << " rows and columns." << std::endl;
108
109 // NOTE: "A" generated by this factory is actually the Schur complement
110 // matrix, but it is required as all smoothers expect "A"
111 Set(currentLevel, "A", S);
112 }
113
114 template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
115 RCP<Xpetra::Matrix<Scalar, LocalOrdinal, GlobalOrdinal, Node>>
117
118 using STS = Teuchos::ScalarTraits<SC>;
119 const SC zero = STS::zero(), one = STS::one();
120
121 RCP<Matrix> A01 = bA->getMatrix(0,1);
122 RCP<Matrix> A10 = bA->getMatrix(1,0);
123 RCP<Matrix> A11 = bA->getMatrix(1,1);
124
125 RCP<BlockedCrsMatrix> bA01 = Teuchos::rcp_dynamic_cast<BlockedCrsMatrix>(A01);
126 const bool isBlocked = (bA01 == Teuchos::null ? false : true);
127
128 const ParameterList& pL = GetParameterList();
129 const SC omega = pL.get<Scalar>("omega");
130
131 TEUCHOS_TEST_FOR_EXCEPTION(omega == zero, Exceptions::RuntimeError,
132 "MueLu::SchurComplementFactory::Build: Scaling parameter omega must not be zero to avoid division by zero.");
133
134 RCP<Matrix> S = Teuchos::null; // Schur complement
135 RCP<Matrix> D = Teuchos::null; // temporary result for A10*Ainv*A01
136
137 // only if the off-diagonal blocks A10 and A01 are non-zero we have to do the MM multiplication
138 if(A01.is_null() == false && A10.is_null() == false) {
139 // scale with -1/omega
140 Ainv->scale(Teuchos::as<Scalar>(-one/omega));
141
142 // build Schur complement operator
143 if (!isBlocked) {
144 RCP<ParameterList> myparams = rcp(new ParameterList);
145 myparams->set("compute global constants", true);
146
147 // -1/omega*Ainv*A01
148 TEUCHOS_TEST_FOR_EXCEPTION(A01->getRangeMap()->isSameAs(*(Ainv->getDomainMap())) == false, Exceptions::RuntimeError,
149 "MueLu::SchurComplementFactory::Build: RangeMap of A01 and domain map of Ainv are not the same.");
150 RCP<Matrix> C = MatrixMatrix::Multiply(*Ainv, false, *A01, false, GetOStream(Statistics2), true, true, std::string("SchurComplementFactory"), myparams);
151
152 // -1/omega*A10*Ainv*A01
153 TEUCHOS_TEST_FOR_EXCEPTION(A01->getRangeMap()->isSameAs(*(A10->getDomainMap())) == false, Exceptions::RuntimeError,
154 "MueLu::SchurComplementFactory::Build: RangeMap of A10 and domain map A01 are not the same.");
155 D = MatrixMatrix::Multiply(*A10, false, *C, false, GetOStream(Statistics2), true, true, std::string("SchurComplementFactory"), myparams);
156 }
157 else {
158 // nested blocking
159 auto bA10 = Teuchos::rcp_dynamic_cast<BlockedCrsMatrix>(A10);
160 auto bAinv = Teuchos::rcp_dynamic_cast<BlockedCrsMatrix>(Ainv);
161 TEUCHOS_TEST_FOR_EXCEPTION(bAinv == Teuchos::null, Exceptions::RuntimeError,
162 "MueLu::SchurComplementFactory::Build: Casting Ainv to BlockedCrsMatrix not possible.");
163
164 // -1/omega*bAinv*bA01
165 TEUCHOS_TEST_FOR_EXCEPTION(bA01->Rows() != bAinv->Cols(), Exceptions::RuntimeError,
166 "MueLu::SchurComplementFactory::Build: Block rows and cols of bA01 and bAinv are not compatible.");
167 RCP<BlockedCrsMatrix> C = MatrixMatrix::TwoMatrixMultiplyBlock(*bAinv, false, *bA01, false, GetOStream(Statistics2));
168
169 // -1/omega*A10*Ainv*A01
170 TEUCHOS_TEST_FOR_EXCEPTION(bA10->Rows() != bA01->Cols(), Exceptions::RuntimeError,
171 "MueLu::SchurComplementFactory::Build: Block rows and cols of bA10 and bA01 are not compatible.");
172 D = MatrixMatrix::TwoMatrixMultiplyBlock(*bA10, false, *C, false, GetOStream(Statistics2));
173 }
174 if (!A11.is_null()) {
175 MatrixMatrix::TwoMatrixAdd(*A11, false, one, *D, false, one, S, GetOStream(Statistics2));
176 S->fillComplete();
177
178 TEUCHOS_TEST_FOR_EXCEPTION(A11->getRangeMap()->isSameAs(*(S->getRangeMap())) == false, Exceptions::RuntimeError,
179 "MueLu::SchurComplementFactory::Build: RangeMap of A11 and S are not the same.");
180 TEUCHOS_TEST_FOR_EXCEPTION(A11->getDomainMap()->isSameAs(*(S->getDomainMap())) == false, Exceptions::RuntimeError,
181 "MueLu::SchurComplementFactory::Build: DomainMap of A11 and S are not the same.");
182 }
183 else {
184 S = MatrixFactory::BuildCopy(D);
185 }
186 }
187 else {
188 if (!A11.is_null()) {
189 S = MatrixFactory::BuildCopy(A11);
190 } else {
191 S = MatrixFactory::Build(A11->getRowMap(), 10 /*A11->getLocalMaxNumRowEntries()*/);
192 S->fillComplete(A11->getDomainMap(),A11->getRangeMap());
193 }
194 }
195
196 // Check whether Schur complement operator is a 1x1 block matrix.
197 // If so, unwrap it and return the CrsMatrix based Matrix object
198 // We need this, as single-block smoothers expect it this way.
199 // In case of Thyra GIDs we obtain a Schur complement operator in Thyra GIDs
200 // This may make some special handling in feeding the SchurComplement solver Apply routine
201 // necessary!
202 if (isBlocked) {
203 RCP<BlockedCrsMatrix> bS = Teuchos::rcp_dynamic_cast<BlockedCrsMatrix>(S);
204
205 if (bS != Teuchos::null && bS->Rows() == 1 && bS->Cols() == 1) {
206 RCP<Matrix> temp = bS->getCrsMatrix();
207 S.swap(temp);
208 }
209 }
210
211 return S;
212 }
213
214} // namespace MueLu
215
216#endif /* MUELU_SCHURCOMPLEMENTFACTORY_DEF_HPP_ */
MueLu::DefaultScalar Scalar
Exception indicating invalid cast attempted.
Exception throws to report errors in the internal logical of the program.
Timer to be used in factories. Similar to Monitor but with additional timers.
Class that holds all level-specific information.
Definition: MueLu_Level.hpp:99
void DeclareInput(const std::string &ename, const FactoryBase *factory, const FactoryBase *requestedBy=NoFactory::get())
Callback from FactoryBase::CallDeclareInput() and FactoryBase::DeclareInput()
T & Get(const std::string &ename, const FactoryBase *factory=NoFactory::get())
Get data without decrementing associated storage counter (i.e., read-only access)....
static const RCP< const NoFactory > getRCP()
Static Get() functions.
RCP< const ParameterList > GetValidParameterList() const
Return a const parameter list of valid parameters that setParameterList() will accept.
void DeclareInput(Level &currentLevel) const
Input.
void Build(Level &currentLevel) const
Build an object with this factory.
RCP< Matrix > ComputeSchurComplement(RCP< BlockedCrsMatrix > &bA, RCP< Matrix > &Ainv) const
Schur complement calculation method.
Namespace for MueLu classes and methods.
@ Statistics2
Print even more statistics.
@ Statistics1
Print more statistics.