MueLu  Version of the Day
MueLu_ZeroSubBlockAFactory_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 #ifndef MUELU_ZEROSUBBLOCKAFACTORY_DEF_HPP_
48 #define MUELU_ZEROSUBBLOCKAFACTORY_DEF_HPP_
49 
50 
52 
53 #include <Xpetra_BlockedCrsMatrix.hpp>
54 #include <Xpetra_Matrix.hpp>
55 #include <Xpetra_MatrixFactory.hpp>
56 #include <Xpetra_Vector.hpp>
57 #include <Xpetra_VectorFactory.hpp>
58 
59 #include "MueLu_Level.hpp"
60 #include "MueLu_Monitor.hpp"
61 
62 namespace MueLu {
63 
64  template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
66  RCP<ParameterList> validParamList = rcp(new ParameterList());
67 
68  validParamList->set<RCP<const FactoryBase>>("A", MueLu::NoFactory::getRCP(), "Generating factory for A.");
69 
70  validParamList->set<int>("block row", 0, "Block row of subblock in block matrix A");
71  validParamList->set<int>("block col", 0, "Block column of subblock in block matrix A");
72 
73  return validParamList;
74  }
75 
76  template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
78  Input(currentLevel, "A");
79  }
80 
81  template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
83  FactoryMonitor m(*this, "Build", currentLevel);
84 
85  const ParameterList& pL = GetParameterList();
86  const size_t row = Teuchos::as<size_t>(pL.get<int>("block row"));
87  const size_t col = Teuchos::as<size_t>(pL.get<int>("block col"));
88 
89  // Get the blocked input matrix
90  RCP<Matrix> Ain = currentLevel.Get<RCP<Matrix>>("A", this->GetFactory("A").get());
91  RCP<BlockedCrsMatrix> A = rcp_dynamic_cast<BlockedCrsMatrix>(Ain);
92 
93  // Perform some basic sanity checks
94  TEUCHOS_TEST_FOR_EXCEPTION(A.is_null(), Exceptions::BadCast, "Input matrix A is not a BlockedCrsMatrix.");
95  TEUCHOS_TEST_FOR_EXCEPTION(row > A->Rows(), Exceptions::RuntimeError, "row [" << row << "] > A.Rows() [" << A->Rows() << "].");
96  TEUCHOS_TEST_FOR_EXCEPTION(col > A->Cols(), Exceptions::RuntimeError, "col [" << col << "] > A.Cols() [" << A->Cols() << "].");
97 
98  // strided maps for range and domain map of sub matrix
99  RCP<const StridedMap> stridedRangeMap = Teuchos::null;
100  RCP<const StridedMap> stridedDomainMap = Teuchos::null;
101 
102  // extract map information from MapExtractor
103  RCP<const MapExtractor> rangeMapExtractor = A->getRangeMapExtractor();
104  RCP<const MapExtractor> domainMapExtractor = A->getDomainMapExtractor();
105 
106  // In case that both user-specified and internal striding information from the submaps
107  // does not contain valid striding information, try to extract it from the global maps
108  // in the map extractor.
109  if (stridedRangeMap.is_null()) {
110  TEUCHOS_TEST_FOR_EXCEPTION(rangeMapExtractor->getMap(row).is_null(), Exceptions::BadCast,
111  "Range map extractor contains non-strided maps in block row " << row << ". This should not be.");
112  stridedRangeMap = rcp_dynamic_cast<const StridedMap>(rangeMapExtractor->getMap(row));
113  }
114 
115  if (stridedDomainMap.is_null()) {
116  TEUCHOS_TEST_FOR_EXCEPTION(domainMapExtractor->getMap(row).is_null(), Exceptions::BadCast,
117  "Domain map extractor contains non-strided maps in block row " << row << ". This should not be.");
118  stridedDomainMap = rcp_dynamic_cast<const StridedMap>(domainMapExtractor->getMap(col));
119  }
120 
121  TEUCHOS_TEST_FOR_EXCEPTION(stridedRangeMap.is_null(), Exceptions::BadCast, "rangeMap " << row << " is not a strided map.");
122  TEUCHOS_TEST_FOR_EXCEPTION(stridedDomainMap.is_null(), Exceptions::BadCast, "domainMap " << col << " is not a strided map.");
123 
124  RCP<Matrix> Op = MatrixFactory::Build(stridedRangeMap, stridedDomainMap, static_cast<size_t>(0));
125  TEUCHOS_ASSERT(!Op.is_null());
126 
127  Op->fillComplete(stridedDomainMap, stridedRangeMap);
128  TEUCHOS_ASSERT(Op->isFillComplete());
129 
130  GetOStream(Statistics1) << "A(" << row << "," << col << ") is a single block and has strided maps:"
131  << "\n range map fixed block size = " << stridedRangeMap ->getFixedBlockSize() << ", strided block id = " << stridedRangeMap ->getStridedBlockId()
132  << "\n domain map fixed block size = " << stridedDomainMap->getFixedBlockSize() << ", strided block id = " << stridedDomainMap->getStridedBlockId() << std::endl;
133  GetOStream(Statistics2) << "A(" << row << "," << col << ") has " << Op->getGlobalNumRows() << "x" << Op->getGlobalNumCols() << " rows and columns." << std::endl;
134 
135  if (Op->IsView("stridedMaps") == true)
136  Op->RemoveView("stridedMaps");
137  Op->CreateView("stridedMaps", stridedRangeMap, stridedDomainMap);
138 
139  currentLevel.Set("A", Op, this);
140  }
141 
142 } // namespace MueLu
143 
144 #endif /* MUELU_ZEROSUBBLOCKAFACTORY_DEF_HPP_ */
Exception indicating invalid cast attempted.
T & Get(const std::string &ename, const FactoryBase *factory=NoFactory::get())
Get data without decrementing associated storage counter (i.e., read-only access). Usage: Level->Get< RCP<Matrix> >("A", factory) if factory == NULL => use default factory.
void DeclareInput(Level &currentLevel) const override
Specifies the data that this class needs, and the factories that generate that data.
void Build(Level &currentLevel) const override
Build a zero sub-block object with this factory.
Timer to be used in factories. Similar to Monitor but with additional timers.
Print more statistics.
Namespace for MueLu classes and methods.
Print even more statistics.
Class that holds all level-specific information.
Definition: MueLu_Level.hpp:99
RCP< const ParameterList > GetValidParameterList() const override
Input.
void Set(const std::string &ename, const T &entry, const FactoryBase *factory=NoFactory::get())
Exception throws to report errors in the internal logical of the program.
static const RCP< const NoFactory > getRCP()
Static Get() functions.