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