MueLu  Version of the Day
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
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 
66 namespace 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 
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 
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)
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 
Important warning messages (one line)
Exception indicating invalid cast attempted.
T & get(const std::string &name, T def_value)
ParameterList & set(std::string const &name, T const &value, std::string const &docString="", RCP< const ParameterEntryValidator > const &validator=null)
Timer to be used in factories. Similar to Monitor but with additional timers.
#define TEUCHOS_TEST_FOR_EXCEPTION(throw_exception_test, Exception, msg)
Print more statistics.
Print additional debugging information.
static const NoFactory * get()
RCP< const ParameterList > GetValidParameterList() const
Input.
TEUCHOS_DEPRECATED RCP< T > rcp(T *p, Dealloc_T dealloc, bool owns_mem)
Class that holds all level-specific information.
Definition: MueLu_Level.hpp:99
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.
Exception throws to report errors in the internal logical of the program.
Description of what is happening (more verbose)
static const RCP< const NoFactory > getRCP()
Static Get() functions.
bool is_null() const