MueLu Version of the Day
Loading...
Searching...
No Matches
MueLu_BlockedDirectSolver_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/*
47 * MueLu_BlockedDirectSolver_def.hpp
48 *
49 * Created on: 09.02.2014
50 * Author: tobias
51 */
52
53#ifndef MUELU_BLOCKEDDIRECTSOLVER_DEF_HPP_
54#define MUELU_BLOCKEDDIRECTSOLVER_DEF_HPP_
55
56#include "Teuchos_ArrayViewDecl.hpp"
57#include "Teuchos_ScalarTraits.hpp"
58
59#include "MueLu_ConfigDefs.hpp"
60
61#include <Xpetra_Matrix.hpp>
62#include <Xpetra_BlockedCrsMatrix.hpp>
63#include <Xpetra_MultiVectorFactory.hpp>
64
65#include "MueLu_BlockedDirectSolver.hpp"
66#include "MueLu_MergedBlockedMatrixFactory.hpp"
67#include "MueLu_Level.hpp"
68#include "MueLu_Utilities.hpp"
69#include "MueLu_Monitor.hpp"
70#include "MueLu_HierarchyUtils.hpp"
72#include "MueLu_DirectSolver.hpp"
73
74namespace MueLu {
75
76 template <class Scalar,class LocalOrdinal, class GlobalOrdinal, class Node>
78 : type_("blocked direct solver")
79 {
80 MergedAFact_ = Teuchos::rcp(new MergedBlockedMatrixFactory());
81
82 Teuchos::ParameterList params;
83 s_ = Teuchos::rcp(new DirectSolver("", params));
84 }
85
86 template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
88 RCP<ParameterList> validParamList = rcp(new ParameterList());
89
90 validParamList->set< RCP<const FactoryBase> >("A", null, "Generating factory of the matrix A");
91
92 return validParamList;
93 }
94
95 template <class Scalar,class LocalOrdinal, class GlobalOrdinal, class Node>
97 // Note that we have a nested smoother/solver object (of type DirectSolver), so we have to declare the dependencies by hand
98 // call DeclareInput by hand, since this->Input(currentLevel, "A") would not properly free A in the release mode (dependencies)
99 // We need the blocked version of A as input for the MergedAFact_
100 currentLevel.DeclareInput("A",this->GetFactory("A").get());
101
102 // syncronize input factory for "A" and nested representation for "A"
103 MergedAFact_->SetFactory("A", this->GetFactory("A"));
104
105 // declare input factories for nested direct solver
106 s_->SetFactory("A",MergedAFact_);
107 s_->DeclareInput(currentLevel);
108 }
109
110 template <class Scalar,class LocalOrdinal, class GlobalOrdinal, class Node>
112 RCP<Teuchos::FancyOStream> out = Teuchos::fancyOStream(Teuchos::rcpFromRef(std::cout));
113
114 FactoryMonitor m(*this, "Setup BlockedDirectSolver", currentLevel);
115 if (this->IsSetup() == true)
116 this->GetOStream(Warnings0) << "MueLu::BlockedDirectSolver::Setup(): Setup() has already been called";
117
118 // extract blocked operator A from current level
119 A_ = Factory::Get< RCP<Matrix> >(currentLevel, "A"); // A needed for extracting map extractors
120 RCP<BlockedCrsMatrix> bA = Teuchos::rcp_dynamic_cast<BlockedCrsMatrix>(A_);
121 TEUCHOS_TEST_FOR_EXCEPTION(bA.is_null(), Exceptions::BadCast,
122 "MueLu::BlockedDirectSolver::Build: input matrix A is not of type BlockedCrsMatrix.");
123
124 s_->Setup(currentLevel);
125
126 this->IsSetup(true);
127 }
128
129 template <class Scalar,class LocalOrdinal, class GlobalOrdinal, class Node>
130 void BlockedDirectSolver<Scalar, LocalOrdinal, GlobalOrdinal, Node>::Apply(MultiVector &X, const MultiVector& B, bool InitialGuessIsZero) const {
131 TEUCHOS_TEST_FOR_EXCEPTION(this->IsSetup() == false, Exceptions::RuntimeError,
132 "MueLu::BlockedDirectSolver::Apply(): Setup() has not been called");
133
134 RCP<MultiVector> rcpX = Teuchos::rcpFromRef(X);
135 RCP<const MultiVector> rcpB = Teuchos::rcpFromRef(B);
136 RCP<BlockedMultiVector> bX = Teuchos::rcp_dynamic_cast<BlockedMultiVector>(rcpX);
137 RCP<const BlockedMultiVector> bB = Teuchos::rcp_dynamic_cast<const BlockedMultiVector>(rcpB);
138
139#ifdef HAVE_MUELU_DEBUG
140 RCP<BlockedCrsMatrix> bA = Teuchos::rcp_dynamic_cast<BlockedCrsMatrix>(A_);
141 if(bB.is_null() == false) {
142 //TEUCHOS_TEST_FOR_EXCEPTION(A_->getRangeMap()->isSameAs(*(B.getMap())) == false, Exceptions::RuntimeError, "MueLu::BlockedDirectSolver::Apply(): The map of RHS vector B is not the same as range map of the blocked operator A. Please check the map of B and A.");
143 } else {
144 TEUCHOS_TEST_FOR_EXCEPTION(bA->getFullRangeMap()->isSameAs(*(B.getMap())) == false, Exceptions::RuntimeError, "MueLu::BlockedDirectSolver::Apply(): The map of RHS vector B is not the same as range map of the blocked operator A. Please check the map of B and A.");
145 }
146 if(bX.is_null() == false) {
147 //TEUCHOS_TEST_FOR_EXCEPTION(A_->getDomainMap()->isSameAs(*(X.getMap())) == false, Exceptions::RuntimeError, "MueLu::BlockedDirectSolver::Apply(): The map of the solution vector X is not the same as domain map of the blocked operator A. Please check the map of X and A.");
148 } else {
149 TEUCHOS_TEST_FOR_EXCEPTION(bA->getFullDomainMap()->isSameAs(*(X.getMap())) == false, Exceptions::RuntimeError, "MueLu::BlockedDirectSolver::Apply(): The map of the solution vector X is not the same as domain map of the blocked operator A. Please check the map of X and A.");
150 }
151#endif
152
153 if (bB.is_null() == true && bX.is_null() == true) {
154 // standard case (neither B nor X are blocked)
155 s_->Apply(X, B, InitialGuessIsZero);
156 } else if(bB.is_null() == false && bX.is_null() == false) {
157 // both B and X are blocked
158 RCP<MultiVector> mergedX = bX->Merge();
159 RCP<const MultiVector> mergedB = bB->Merge();
160 s_->Apply(*mergedX, *mergedB, InitialGuessIsZero);
161 RCP<MultiVector> xx = Teuchos::rcp(new BlockedMultiVector(bX->getBlockedMap(),mergedX));
162 SC zero = Teuchos::ScalarTraits<SC>::zero(), one = Teuchos::ScalarTraits<SC>::one();
163 X.update(one,*xx,zero);
164 } else if (bB.is_null() == true && bX.is_null() == false) {
165 // the solution vector is blocked
166 RCP<MultiVector> mergedX = bX->Merge();
167 s_->Apply(*mergedX, B, InitialGuessIsZero);
168 RCP<MultiVector> xx = Teuchos::rcp(new BlockedMultiVector(bX->getBlockedMap(),mergedX));
169 SC zero = Teuchos::ScalarTraits<SC>::zero(), one = Teuchos::ScalarTraits<SC>::one();
170 X.update(one,*xx,zero);
171 } else if (bB.is_null() == false && bX.is_null() == true) {
172 // only the RHS vector is blocked
173 RCP<const MultiVector> mergedB = bB->Merge();
174 s_->Apply(X, *mergedB, InitialGuessIsZero);
175 }
176 }
177
178 template <class Scalar,class LocalOrdinal, class GlobalOrdinal, class Node>
179 RCP<MueLu::SmootherPrototype<Scalar, LocalOrdinal, GlobalOrdinal, Node> >
183
184 template <class Scalar,class LocalOrdinal, class GlobalOrdinal, class Node>
186 std::ostringstream out;
188 out << "{type = " << type_ << "}";
189 return out.str();
190 }
191
192 template <class Scalar,class LocalOrdinal, class GlobalOrdinal, class Node>
193 void BlockedDirectSolver<Scalar, LocalOrdinal, GlobalOrdinal, Node>::print(Teuchos::FancyOStream &out, const VerbLevel verbLevel) const {
195
196 if (verbLevel & Parameters0)
197 out0 << "Prec. type: " << type_ << std::endl;
198
199 if (verbLevel & Debug)
200 out0 << "IsSetup: " << Teuchos::toString(SmootherPrototype::IsSetup()) << std::endl;
201 }
202
203 template <class Scalar,class LocalOrdinal, class GlobalOrdinal, class Node>
205 // FIXME: This is a placeholder
206 return Teuchos::OrdinalTraits<size_t>::invalid();
207 }
208
209
210} // namespace MueLu
211
212
213
214#endif /* MUELU_BLOCKEDDIRECTSOLVER_DEF_HPP_ */
#define MUELU_DESCRIBE
Helper macro for implementing Describable::describe() for BaseClass objects.
direct solver for nxn blocked matrices
void Setup(Level &currentLevel)
Setup routine Call the underlaying Setup routine of the nested direct solver once the input block mat...
RCP< const ParameterList > GetValidParameterList() const
Input.
void DeclareInput(Level &currentLevel) const
Input.
RCP< DirectSolver > s_
Direct solver.
void Apply(MultiVector &X, const MultiVector &B, bool InitialGuessIsZero=false) const
Apply the direct solver. Solves the linear system AX=B using the constructed solver.
std::string description() const
Return a simple one-line description of this object.
void print(Teuchos::FancyOStream &out, const VerbLevel verbLevel=Default) const
Print the object with some verbosity level to an FancyOStream object.
RCP< MergedBlockedMatrixFactory > MergedAFact_
Factory to generate merged block matrix.
size_t getNodeSmootherComplexity() const
Get a rough estimate of cost per iteration.
RCP< SmootherPrototype > Copy() const
virtual std::string description() const
Return a simple one-line description of this object.
Class that encapsulates direct solvers. Autoselection of AmesosSmoother or Amesos2Smoother according ...
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.
void DeclareInput(const std::string &ename, const FactoryBase *factory, const FactoryBase *requestedBy=NoFactory::get())
Callback from FactoryBase::CallDeclareInput() and FactoryBase::DeclareInput()
bool IsSetup() const
Get the state of a smoother prototype.
Namespace for MueLu classes and methods.
@ Warnings0
Important warning messages (one line)
@ Debug
Print additional debugging information.
@ Parameters0
Print class parameters.