MueLu Version of the Day
Loading...
Searching...
No Matches
MueLu_ReorderBlockAFactory_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// Tobias Wiesner (tawiesn@sandia.gov)
43//
44// ***********************************************************************
45//
46// @HEADER
47
48#ifndef MUELU_REORDERBLOCKAFACTORY_DEF_HPP_
49#define MUELU_REORDERBLOCKAFACTORY_DEF_HPP_
50
51
53
54#include <Xpetra_BlockReorderManager.hpp>
55#include <Xpetra_BlockedCrsMatrix.hpp>
56#include <Xpetra_ReorderedBlockedCrsMatrix.hpp>
57#include <Xpetra_MapExtractor.hpp>
58#include <Xpetra_MapExtractorFactory.hpp>
59#include <Xpetra_Matrix.hpp>
60#include <Xpetra_MatrixUtils.hpp>
61//#include <Xpetra_StridedMapFactory.hpp>
62
63#include "MueLu_Level.hpp"
64#include "MueLu_Monitor.hpp"
65
66namespace MueLu {
67
68 template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
70 RCP<ParameterList> validParamList = rcp(new ParameterList());
71
72 validParamList->set< RCP<const FactoryBase> >("A", MueLu::NoFactory::getRCP(), "Generating factory for A.");
73
74 validParamList->set< std::string > ("Reorder Type", "", "String describing the reordering of blocks");
75
76 // TODO not very elegant.
77 validParamList->set< RCP<const FactoryBase> >("Map1", Teuchos::null, "Generating factory of the fine level map associated with the (1,1) block in your n x n block matrix.");
78 validParamList->set< RCP<const FactoryBase> >("Map2", Teuchos::null, "Generating factory of the fine level map associated with the (2,2) block in your n x n block matrix.");
79 validParamList->set< RCP<const FactoryBase> >("Map3", Teuchos::null, "Generating factory of the fine level map associated with the (3,3) block in your n x n block matrix.");
80 validParamList->set< RCP<const FactoryBase> >("Map4", Teuchos::null, "Generating factory of the fine level map associated with the (4,4) block in your n x n block matrix.");
81 validParamList->set< RCP<const FactoryBase> >("Map5", Teuchos::null, "Generating factory of the fine level map associated with the (5,5) block in your n x n block matrix.");
82 validParamList->set< RCP<const FactoryBase> >("Map6", Teuchos::null, "Generating factory of the fine level map associated with the (6,6) block in your n x n block matrix.");
83 validParamList->set< RCP<const FactoryBase> >("Map7", Teuchos::null, "Generating factory of the fine level map associated with the (7,7) block in your n x n block matrix.");
84 validParamList->set< RCP<const FactoryBase> >("Map8", Teuchos::null, "Generating factory of the fine level map associated with the (8,8) block in your n x n block matrix.");
85 validParamList->set< RCP<const FactoryBase> >("Map9", Teuchos::null, "Generating factory of the fine level map associated with the (9,9) block in your n x n block matrix.");
86
87 return validParamList;
88 }
89
90 template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
92 Input(currentLevel, "A");
93 }
94
95 template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
97 FactoryMonitor m(*this, "ReorderBlockA factory", currentLevel);
98
99 const ParameterList& pL = GetParameterList();
100 std::string reorderStr = pL.get<std::string>("Reorder Type");
101
102 RCP<Matrix> Ain = Get<RCP<Matrix> >(currentLevel, "A");
103
104 RCP<BlockedCrsMatrix> A = rcp_dynamic_cast<BlockedCrsMatrix>(Ain);
105
106 // special case: we get a single block CrsMatrix object on the finest level and
107 // split it into a nxn blocked operator
108 if (A == Teuchos::null && currentLevel.GetLevelID() == 0) {
109
110 GetOStream(Warnings0) << "Split input matrix (Warning: this is a rather expensive operation)" << std::endl;
111
112 std::vector<Teuchos::RCP<const Map> > xmaps;
113
114 for(int it = 1; it < 10; it++) {
115 std::stringstream ss;
116 ss << "Map" << it;
117 if (currentLevel.IsAvailable(ss.str(), NoFactory::get())) {
118 RCP<const Map> submap = currentLevel.Get< RCP<const Map> >(ss.str(), NoFactory::get());
119 GetOStream(Runtime1) << "Use user-given submap #" << it << ": length dimension=" << submap->getGlobalNumElements() << std::endl;
120 xmaps.push_back(submap);
121 }
122 }
123
124 bool bThyraMode = false; // no support for Thyra mode (yet)
125 RCP<const MapExtractor> map_extractor = Xpetra::MapExtractorFactory<Scalar, LocalOrdinal, GlobalOrdinal, Node>::Build(Ain->getRowMap(),xmaps,bThyraMode);
126
127 // split null space vectors
128 //RCP<MultiVector> nullspace1 = map_extractor->ExtractVector(nullspace,0);
129 //RCP<MultiVector> nullspace2 = map_extractor->ExtractVector(nullspace,1);
130
131 Teuchos::RCP<Xpetra::BlockedCrsMatrix<Scalar,LocalOrdinal,GlobalOrdinal,Node> > bOp =
132 Xpetra::MatrixUtils<Scalar, LocalOrdinal, GlobalOrdinal, Node>::SplitMatrix(*Ain,map_extractor,map_extractor,Teuchos::null,bThyraMode);
133
134 TEUCHOS_TEST_FOR_EXCEPTION(Ain->getGlobalNumRows() != bOp->getGlobalNumRows(), Exceptions::RuntimeError, "Split operator not consistent with input operator (different number of rows).");
135 TEUCHOS_TEST_FOR_EXCEPTION(Ain->getLocalNumRows() != bOp->getLocalNumRows(), Exceptions::RuntimeError, "Split operator not consistent with input operator (different number of node rows).");
136 TEUCHOS_TEST_FOR_EXCEPTION(Ain->getLocalNumEntries() != bOp->getLocalNumEntries(), Exceptions::RuntimeError, "Split operator not consistent with input operator (different number of local entries).");
137 TEUCHOS_TEST_FOR_EXCEPTION(Ain->getGlobalNumEntries() != bOp->getGlobalNumEntries(), Exceptions::RuntimeError, "Split operator not consistent with input operator (different number of global entries).");
138
139 A = bOp;
140 }
141
142 // we have a blocked operator as input
143 TEUCHOS_TEST_FOR_EXCEPTION(A.is_null(), Exceptions::BadCast, "Input matrix A is not a BlockedCrsMatrix.");
144 GetOStream(Statistics1) << "Got a " << A->Rows() << "x" << A->Cols() << " blocked operator as input" << std::endl;
145
146 Teuchos::RCP<const Xpetra::BlockReorderManager> brm = Xpetra::blockedReorderFromString(reorderStr);
147 GetOStream(Debug) << "Reordering A using " << brm->toString() << std::endl;
148
149 Teuchos::RCP<const ReorderedBlockedCrsMatrix> brop =
150 Teuchos::rcp_dynamic_cast<const ReorderedBlockedCrsMatrix>(Xpetra::buildReorderedBlockedCrsMatrix<Scalar,LocalOrdinal,GlobalOrdinal,Node>(brm, A));
151
152 TEUCHOS_TEST_FOR_EXCEPTION(brop.is_null(), Exceptions::RuntimeError, "Block reordering of " << A->Rows() << "x" << A->Cols() << " blocked operator failed.");
153
154 GetOStream(Statistics1) << "Reordering A using " << brm->toString() << " block gives a " << brop->Rows() << "x" << brop->Cols() << " blocked operators" << std::endl;
155 GetOStream(Debug) << "Reordered operator has " << brop->getRangeMap()->getGlobalNumElements() << " rows and " << brop->getDomainMap()->getGlobalNumElements() << " columns" << std::endl;
156 GetOStream(Debug) << "Reordered operator: Use of Thyra style gids = " << brop->getRangeMapExtractor()->getThyraMode() << std::endl;
157
158 // get rid of const (we expect non-const operators stored in Level)
159 Teuchos::RCP<ReorderedBlockedCrsMatrix> bret =
160 Teuchos::rcp_const_cast<ReorderedBlockedCrsMatrix>(brop);
161
162 // strided maps
163 // blocked operators do not have strided maps (information could be misleading?)
164 //Op->CreateView("stridedMaps", srangeMap, sdomainMap);
165
166 currentLevel.Set("A", Teuchos::rcp_dynamic_cast<Matrix>(bret), this);
167 }
168
169} // namespace MueLu
170
171#endif /* MUELU_REORDERBLOCKAFACTORY_DEF_HPP_ */
172
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
bool IsAvailable(const std::string &ename, const FactoryBase *factory=NoFactory::get()) const
Test whether a need's value has been saved.
int GetLevelID() const
Return level number.
Definition: MueLu_Level.cpp:76
T & Get(const std::string &ename, const FactoryBase *factory=NoFactory::get())
Get data without decrementing associated storage counter (i.e., read-only access)....
void Set(const std::string &ename, const T &entry, const FactoryBase *factory=NoFactory::get())
static const RCP< const NoFactory > getRCP()
Static Get() functions.
static const NoFactory * get()
RCP< const ParameterList > GetValidParameterList() const
Input.
void Build(Level &currentLevel) const
Build an object with this factory.
void DeclareInput(Level &currentLevel) const
Specifies the data that this class needs, and the factories that generate that data.
Namespace for MueLu classes and methods.
@ Warnings0
Important warning messages (one line)
@ Debug
Print additional debugging information.
@ Statistics1
Print more statistics.
@ Runtime1
Description of what is happening (more verbose)