MueLu  Version of the Day
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
MueLu_BlockedRAPFactory_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 #ifndef MUELU_BLOCKEDRAPFACTORY_DEF_HPP
47 #define MUELU_BLOCKEDRAPFACTORY_DEF_HPP
48 
50 #include <Xpetra_MatrixFactory.hpp>
51 #include <Xpetra_Matrix.hpp>
52 #include <Xpetra_MatrixMatrix.hpp>
53 #include <Xpetra_VectorFactory.hpp>
54 #include <Xpetra_Vector.hpp>
55 
57 
58 #include "MueLu_MasterList.hpp"
59 #include "MueLu_Monitor.hpp"
60 #include "MueLu_PerfUtils.hpp"
62 #include "MueLu_Utilities.hpp"
63 
64 namespace MueLu {
65 
66  template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
68  : checkAc_(false), repairZeroDiagonals_(false)
69  { }
70 
71  template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
73  RCP<ParameterList> validParamList = rcp(new ParameterList());
74 
75 #define SET_VALID_ENTRY(name) validParamList->setEntry(name, MasterList::getEntry(name))
76  SET_VALID_ENTRY("transpose: use implicit");
77 #undef SET_VALID_ENTRY
78  validParamList->set< RCP<const FactoryBase> >("A", null, "Generating factory of the matrix A used during the prolongator smoothing process");
79  validParamList->set< RCP<const FactoryBase> >("P", null, "Prolongator factory");
80  validParamList->set< RCP<const FactoryBase> >("R", null, "Restrictor factory");
81 
82  return validParamList;
83  }
84 
85  template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
87  const Teuchos::ParameterList& pL = GetParameterList();
88  if (pL.get<bool>("transpose: use implicit") == false)
89  Input(coarseLevel, "R");
90 
91  Input(fineLevel, "A");
92  Input(coarseLevel, "P");
93 
94  // call DeclareInput of all user-given transfer factories
95  for (std::vector<RCP<const FactoryBase> >::const_iterator it = transferFacts_.begin(); it != transferFacts_.end(); ++it)
96  (*it)->CallDeclareInput(coarseLevel);
97  }
98 
99  template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
100  void BlockedRAPFactory<Scalar, LocalOrdinal, GlobalOrdinal, Node>::Build(Level &fineLevel, Level &coarseLevel) const { //FIXME make fineLevel const!!
101  FactoryMonitor m(*this, "Computing Ac (block)", coarseLevel);
102 
103  const ParameterList& pL = GetParameterList();
104 
105  RCP<Matrix> A = Get< RCP<Matrix> >(fineLevel, "A");
106  RCP<Matrix> P = Get< RCP<Matrix> >(coarseLevel, "P");
107 
108 
109  RCP<BlockedCrsMatrix> bA = rcp_dynamic_cast<BlockedCrsMatrix>(A);
110  RCP<BlockedCrsMatrix> bP = rcp_dynamic_cast<BlockedCrsMatrix>(P);
111  TEUCHOS_TEST_FOR_EXCEPTION(bA.is_null() || bP.is_null(), Exceptions::BadCast, "Matrices A and P must be of type BlockedCrsMatrix.");
112 
113 
116  {
117  SubFactoryMonitor subM(*this, "MxM: A x P", coarseLevel);
118 
119  // Triple matrix product for BlockedCrsMatrixClass
120  TEUCHOS_TEST_FOR_EXCEPTION((bA->Cols() != bP->Rows()), Exceptions::BadCast,
121  "Block matrix dimensions do not match: "
122  "A is " << bA->Rows() << "x" << bA->Cols() <<
123  "P is " << bP->Rows() << "x" << bP->Cols());
124 
125  bAP = MatrixMatrix::TwoMatrixMultiplyBlock(*bA, false, *bP, false, GetOStream(Statistics2), true, true);
126  }
127 
128 
129  // If we do not modify matrix later, allow optimization of storage.
130  // This is necessary for new faster Epetra MM kernels.
131  bool doOptimizeStorage = !checkAc_;
132 
133  const bool doTranspose = true;
134  const bool doFillComplete = true;
135  if (pL.get<bool>("transpose: use implicit") == true) {
136  SubFactoryMonitor m2(*this, "MxM: P' x (AP) (implicit)", coarseLevel);
137  bAc = MatrixMatrix::TwoMatrixMultiplyBlock(*bP, doTranspose, *bAP, !doTranspose, GetOStream(Statistics2), doFillComplete, doOptimizeStorage);
138 
139  } else {
140  RCP<Matrix> R = Get< RCP<Matrix> >(coarseLevel, "R");
141  RCP<BlockedCrsMatrix> bR = rcp_dynamic_cast<BlockedCrsMatrix>(R);
142  TEUCHOS_TEST_FOR_EXCEPTION(bR.is_null(), Exceptions::BadCast, "Matrix R must be of type BlockedCrsMatrix.");
143 
144  TEUCHOS_TEST_FOR_EXCEPTION(bA->Rows() != bR->Cols(), Exceptions::BadCast,
145  "Block matrix dimensions do not match: "
146  "R is " << bR->Rows() << "x" << bR->Cols() <<
147  "A is " << bA->Rows() << "x" << bA->Cols());
148 
149  SubFactoryMonitor m2(*this, "MxM: R x (AP) (explicit)", coarseLevel);
150  bAc = MatrixMatrix::TwoMatrixMultiplyBlock(*bR, !doTranspose, *bAP, !doTranspose, GetOStream(Statistics2), doFillComplete, doOptimizeStorage);
151  }
152 
153 
154  if (checkAc_)
155  CheckMainDiagonal(bAc);
156 
157  GetOStream(Statistics1) << PerfUtils::PrintMatrixInfo(*bAc, "Ac (blocked)");
158 
159  // static int run = 1;
160  // RCP<CrsMatrixWrap> A11 = rcp(new CrsMatrixWrap(bAc->getMatrix(0,0)));
161  // Utils::Write(toString(run) + "_A_11.mm", *A11);
162  // if (!bAc->getMatrix(1,1).is_null()) {
163  // RCP<CrsMatrixWrap> A22 = rcp(new CrsMatrixWrap(bAc->getMatrix(1,1)));
164  // Utils::Write(toString(run) + "_A_22.mm", *A22);
165  // }
166  // RCP<CrsMatrixWrap> Am = rcp(new CrsMatrixWrap(bAc->Merge()));
167  // Utils::Write(toString(run) + "_A.mm", *Am);
168  // run++;
169 
170  Set<RCP <Matrix> >(coarseLevel, "A", bAc);
171 
172  if (transferFacts_.begin() != transferFacts_.end()) {
173  SubFactoryMonitor m1(*this, "Projections", coarseLevel);
174 
175  // call Build of all user-given transfer factories
176  for (std::vector<RCP<const FactoryBase> >::const_iterator it = transferFacts_.begin(); it != transferFacts_.end(); ++it) {
177  RCP<const FactoryBase> fac = *it;
178 
179  GetOStream(Runtime0) << "BlockRAPFactory: call transfer factory: " << fac->description() << std::endl;
180 
181  fac->CallBuild(coarseLevel);
182 
183  // AP (11/11/13): I am not sure exactly why we need to call Release, but we do need it to get rid
184  // of dangling data for CoordinatesTransferFactory
185  coarseLevel.Release(*fac);
186  }
187  }
188  }
189 
190 
191  template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
193  RCP<Matrix> c00 = bAc->getMatrix(0, 0);
194  RCP<Matrix> Aout = MatrixFactory::Build(c00->getRowMap(), c00->getGlobalMaxNumRowEntries(), Xpetra::StaticProfile);
195 
196  RCP<Vector> diagVec = VectorFactory::Build(c00->getRowMap());
197  c00->getLocalDiagCopy(*diagVec);
198  ArrayRCP<SC> diagVal = diagVec->getDataNonConst(0);
199 
200  // loop over local rows
201  for (size_t row = 0; row < c00->getNodeNumRows(); row++) {
202  // get global row id
203  GO grid = c00->getRowMap()->getGlobalElement(row); // global row id
204 
205  ArrayView<const LO> indices;
206  ArrayView<const SC> vals;
207  c00->getLocalRowView(row, indices, vals);
208 
209  // just copy all values in output
212 
213  // just copy values
214  for (size_t i = 0; i < as<size_t>(indices.size()); i++) {
215  GO gcid = c00->getColMap()->getGlobalElement(indices[i]); // LID -> GID (column)
216  indout [i] = gcid;
217  valout [i] = vals[i];
218  }
219 
220  Aout->insertGlobalValues(grid, indout.view(0, indout.size()), valout.view(0, valout.size()));
221  if (diagVal[row] == Teuchos::ScalarTraits<SC>::zero() && repairZeroDiagonals) {
222  // always overwrite diagonal entry
223  Aout->insertGlobalValues(grid, Teuchos::tuple<GO>(grid), Teuchos::tuple<SC>(1.0));
224  }
225  }
226 
227  Aout->fillComplete(c00->getDomainMap(), c00->getRangeMap());
228 
229  bAc->setMatrix(0, 0, Aout);
230  }
231 
232  template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
234  // check if it's a TwoLevelFactoryBase based transfer factory
235  TEUCHOS_TEST_FOR_EXCEPTION(rcp_dynamic_cast<const TwoLevelFactoryBase>(factory) == Teuchos::null, Exceptions::BadCast,
236  "Transfer factory is not derived from TwoLevelFactoryBase. This is very strange. "
237  "(Note: you can remove this exception if there's a good reason for)");
238  transferFacts_.push_back(factory);
239  }
240 
241 } //namespace MueLu
242 
243 #define MUELU_BLOCKEDRAPFACTORY_SHORT
244 #endif // MUELU_BLOCKEDRAPFACTORY_DEF_HPP
245 
246 // TODO add plausibility check
247 // TODO add CheckMainDiagonal for Blocked operator
248 // Avoid copying block matrix!
249 // create new empty Operator
Exception indicating invalid cast attempted.
virtual void CallBuild(Level &requestedLevel) const =0
GlobalOrdinal GO
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)
RCP< const ParameterList > GetValidParameterList() const
Return a const parameter list of valid parameters that setParameterList() will accept.
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.
size_type size() const
One-liner description of what is happening.
Print even more statistics.
void Build(Level &fineLevel, Level &coarseLevel) const
Build an object with this factory.
void DeclareInput(Level &fineLevel, Level &coarseLevel) 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
Timer to be used in factories. Similar to SubMonitor but adds a timer level by level.
void AddTransferFactory(const RCP< const FactoryBase > &factory)
Add transfer factory in the end of list of transfer factories in RepartitionAcFactory.
static std::string PrintMatrixInfo(const Matrix &A, const std::string &msgTag, RCP< const Teuchos::ParameterList > params=Teuchos::null)
static void CheckMainDiagonal(RCP< BlockedCrsMatrix > &bAc, bool repairZeroDiagonals=false)
#define SET_VALID_ENTRY(name)
virtual std::string description() const
Return a simple one-line description of this object.