MueLu  Version of the Day
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 
52 #include "MueLu_ReorderBlockAFactory_decl.hpp"
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 
66 namespace MueLu {
67 
68  template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
69  RCP<const ParameterList> ReorderBlockAFactory<Scalar, LocalOrdinal, GlobalOrdinal, Node>::GetValidParameterList() const {
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>
91  void ReorderBlockAFactory<Scalar, LocalOrdinal, GlobalOrdinal, Node>::DeclareInput(Level& currentLevel) const {
92  Input(currentLevel, "A");
93  }
94 
95  template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
96  void ReorderBlockAFactory<Scalar, LocalOrdinal, GlobalOrdinal, Node>::Build(Level& currentLevel) const {
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->getNodeNumRows() != bOp->getNodeNumRows(), Exceptions::RuntimeError, "Split operator not consistent with input operator (different number of node rows).");
136  TEUCHOS_TEST_FOR_EXCEPTION(Ain->getNodeNumEntries() != bOp->getNodeNumEntries(), 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 
MueLu::Runtime1
Description of what is happening (more verbose)
Definition: MueLu_VerbosityLevel.hpp:63
MueLu::ReorderBlockAFactory::Build
void Build(Level &currentLevel) const
Build an object with this factory.
Definition: MueLu_ReorderBlockAFactory_def.hpp:96
MueLu::FactoryMonitor
Timer to be used in factories. Similar to Monitor but with additional timers.
Definition: MueLu_Monitor.hpp:228
MueLu::Warnings0
Important warning messages (one line)
Definition: MueLu_VerbosityLevel.hpp:57
MueLu_Monitor.hpp
MueLu::Exceptions::BadCast
Exception indicating invalid cast attempted.
Definition: MueLu_Exceptions.hpp:57
MueLu
Namespace for MueLu classes and methods.
Definition: MueLu_BrickAggregationFactory_decl.hpp:76
MueLu::Debug
Print additional debugging information.
Definition: MueLu_VerbosityLevel.hpp:79
MueLu_Level.hpp
MueLu_ReorderBlockAFactory_decl.hpp
MueLu::Exceptions::RuntimeError
Exception throws to report errors in the internal logical of the program.
Definition: MueLu_Exceptions.hpp:70
MueLu::ReorderBlockAFactory::DeclareInput
void DeclareInput(Level &currentLevel) const
Specifies the data that this class needs, and the factories that generate that data.
Definition: MueLu_ReorderBlockAFactory_def.hpp:91
MueLu::Statistics1
Print more statistics.
Definition: MueLu_VerbosityLevel.hpp:71
MueLu::ReorderBlockAFactory::GetValidParameterList
RCP< const ParameterList > GetValidParameterList() const
Input.
Definition: MueLu_ReorderBlockAFactory_def.hpp:69
MueLu::NoFactory::get
static const NoFactory * get()
Definition: MueLu_NoFactory.cpp:59
MueLu::Level
Class that holds all level-specific information.
Definition: MueLu_Level.hpp:99
MueLu::NoFactory::getRCP
static const RCP< const NoFactory > getRCP()
Static Get() functions.
Definition: MueLu_NoFactory.hpp:90