MueLu Version of the Day
Loading...
Searching...
No Matches
MueLu_IndefBlockedDiagonalSmoother_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 * MueLu_IndefBlockedDiagonalSmoother_def.hpp
48 *
49 * Created on: 13 May 2014
50 * Author: wiesner
51 */
52
53#ifndef MUELU_INDEFBLOCKEDDIAGONALSMOOTHER_DEF_HPP_
54#define MUELU_INDEFBLOCKEDDIAGONALSMOOTHER_DEF_HPP_
55
56#include "Teuchos_ArrayViewDecl.hpp"
57#include "Teuchos_ScalarTraits.hpp"
58
59#include "MueLu_ConfigDefs.hpp"
60
61#include <Xpetra_Matrix.hpp>
62#include <Xpetra_CrsMatrixWrap.hpp>
63#include <Xpetra_BlockedCrsMatrix.hpp>
64#include <Xpetra_MultiVectorFactory.hpp>
65#include <Xpetra_VectorFactory.hpp>
66#include <Xpetra_ReorderedBlockedCrsMatrix.hpp>
67
69#include "MueLu_Level.hpp"
70#include "MueLu_Utilities.hpp"
71#include "MueLu_Monitor.hpp"
72#include "MueLu_HierarchyUtils.hpp"
74#include "MueLu_SubBlockAFactory.hpp"
75
76// include files for default FactoryManager
77#include "MueLu_SchurComplementFactory.hpp"
78#include "MueLu_DirectSolver.hpp"
79#include "MueLu_SmootherFactory.hpp"
80#include "MueLu_FactoryManager.hpp"
81
82namespace MueLu {
83
84 template <class Scalar,class LocalOrdinal, class GlobalOrdinal, class Node>
89
90 template <class Scalar,class LocalOrdinal, class GlobalOrdinal, class Node>
92
93 template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
95 RCP<ParameterList> validParamList = rcp(new ParameterList());
96
97 validParamList->set< RCP<const FactoryBase> >("A", Teuchos::null, "Generating factory of the matrix A (must be a 2x2 block matrix)");
98 validParamList->set< Scalar > ("Damping factor", 1.0, "Damping/Scaling factor");
99 validParamList->set< LocalOrdinal > ("Sweeps", 1, "Number of SIMPLE sweeps (default = 1)");
100 //validParamList->set< bool > ("UseSIMPLEC", false, "Use SIMPLEC instead of SIMPLE (default = false)");
101
102 return validParamList;
103 }
104
106 template <class Scalar,class LocalOrdinal, class GlobalOrdinal, class Node>
108 TEUCHOS_TEST_FOR_EXCEPTION(pos < 0, Exceptions::RuntimeError, "MueLu::IndefBlockedDiagonalSmoother::AddFactoryManager: parameter \'pos\' must not be negative! error.");
109
110 size_t myPos = Teuchos::as<size_t>(pos);
111
112 if (myPos < FactManager_.size()) {
113 // replace existing entris in FactManager_ vector
114 FactManager_.at(myPos) = FactManager;
115 } else if( myPos == FactManager_.size()) {
116 // add new Factory manager in the end of the vector
117 FactManager_.push_back(FactManager);
118 } else { // if(myPos > FactManager_.size())
119 RCP<Teuchos::FancyOStream> out = Teuchos::fancyOStream(Teuchos::rcpFromRef(std::cout));
120 *out << "Warning: cannot add new FactoryManager at proper position " << pos << ". The FactoryManager is just appended to the end. Check this!" << std::endl;
121
122 // add new Factory manager in the end of the vector
123 FactManager_.push_back(FactManager);
124 }
125 }
126
127 template <class Scalar,class LocalOrdinal, class GlobalOrdinal, class Node>
129 currentLevel.DeclareInput("A",this->GetFactory("A").get());
130
131 TEUCHOS_TEST_FOR_EXCEPTION(FactManager_.size() != 2, Exceptions::RuntimeError,"MueLu::IndefBlockedDiagonalSmoother::DeclareInput: You have to declare two FactoryManagers with a \"Smoother\" object: One for predicting the primary variable and one for the SchurComplement system. The smoother for the SchurComplement system needs a SchurComplementFactory as input for variable \"A\"!");
132
133 // loop over all factory managers for the subblocks of blocked operator A
134 std::vector<Teuchos::RCP<const FactoryManagerBase> >::const_iterator it;
135 for(it = FactManager_.begin(); it!=FactManager_.end(); ++it) {
136 SetFactoryManager currentSFM (rcpFromRef(currentLevel), *it);
137
138 // request "Smoother" for current subblock row.
139 currentLevel.DeclareInput("PreSmoother",(*it)->GetFactory("Smoother").get());
140 }
141
142 }
143
144 template <class Scalar,class LocalOrdinal, class GlobalOrdinal, class Node>
146 FactoryMonitor m(*this, "Setup for indefinite blocked diagonal smoother", currentLevel);
147
148 if (SmootherPrototype::IsSetup() == true)
149 this->GetOStream(Warnings0) << "MueLu::IndefBlockedDiagonalSmoother::Setup(): Setup() has already been called";
150
151 // extract blocked operator A from current level
152 A_ = Factory::Get<RCP<Matrix> > (currentLevel, "A");
153
154 RCP<BlockedCrsMatrix> bA = Teuchos::rcp_dynamic_cast<BlockedCrsMatrix>(A_);
155 TEUCHOS_TEST_FOR_EXCEPTION(bA == Teuchos::null, Exceptions::BadCast, "MueLu::IndefBlockedDiagonalSmoother::Setup: input matrix A is not of type BlockedCrsMatrix! error.");
156
157 // store map extractors
158 rangeMapExtractor_ = bA->getRangeMapExtractor();
159 domainMapExtractor_ = bA->getDomainMapExtractor();
160
161 // Store the blocks in local member variables
162 Teuchos::RCP<Matrix> A00 = bA->getMatrix(0, 0);
163 Teuchos::RCP<Matrix> A01 = bA->getMatrix(0, 1);
164 Teuchos::RCP<Matrix> A10 = bA->getMatrix(1, 0);
165 Teuchos::RCP<Matrix> A11 = bA->getMatrix(1, 1);
166
167 F_ = A00;
168 Z_ = A11;
169
170 /*const ParameterList & pL = Factory::GetParameterList();
171 bool bSIMPLEC = pL.get<bool>("UseSIMPLEC");
172
173 // Create the inverse of the diagonal of F
174 RCP<Vector> diagFVector = VectorFactory::Build(F_->getRowMap());
175 if(!bSIMPLEC) {
176 F_->getLocalDiagCopy(*diagFVector); // extract diagonal of F
177 diagFVector->reciprocal(*diagFVector); // build reciprocal
178 } else {
179 const RCP<const Map> rowmap = F_->getRowMap();
180 size_t locSize = rowmap->getLocalNumElements();
181 Teuchos::ArrayRCP<SC> diag = diagFVector->getDataNonConst(0);
182 Teuchos::ArrayView<const LO> cols;
183 Teuchos::ArrayView<const SC> vals;
184 for (size_t i=0; i<locSize; ++i) { // loop over rows
185 F_->getLocalRowView(i,cols,vals);
186 Scalar absRowSum = Teuchos::ScalarTraits<Scalar>::zero();
187 for (LO j=0; j<cols.size(); ++j) { // loop over cols
188 absRowSum += Teuchos::ScalarTraits<Scalar>::magnitude(vals[j]);
189 }
190 diag[i] = absRowSum;
191 }
192 diagFVector->reciprocal(*diagFVector); // build reciprocal
193 }
194 diagFinv_ = diagFVector;*/
195
196 // Set the Smoother
197 // carefully switch to the SubFactoryManagers (defined by the users)
198 {
199 RCP<const FactoryManagerBase> velpredictFactManager = FactManager_.at(0);
200 SetFactoryManager currentSFM (rcpFromRef(currentLevel), velpredictFactManager);
201 velPredictSmoo_ = currentLevel.Get< RCP<SmootherBase> >("PreSmoother",velpredictFactManager->GetFactory("Smoother").get());
202 }
203 {
204 RCP<const FactoryManagerBase> schurFactManager = FactManager_.at(1);
205 SetFactoryManager currentSFM (rcpFromRef(currentLevel), schurFactManager);
206 schurCompSmoo_ = currentLevel.Get< RCP<SmootherBase> >("PreSmoother", schurFactManager->GetFactory("Smoother").get());
207 }
209 }
210
211 template <class Scalar,class LocalOrdinal, class GlobalOrdinal, class Node>
212 void IndefBlockedDiagonalSmoother<Scalar, LocalOrdinal, GlobalOrdinal, Node>::Apply(MultiVector& X, const MultiVector& B, bool /* InitialGuessIsZero */) const
213 {
214 TEUCHOS_TEST_FOR_EXCEPTION(SmootherPrototype::IsSetup() == false, Exceptions::RuntimeError, "MueLu::IndefBlockedDiagonalSmoother::Apply(): Setup() has not been called");
215
216 Teuchos::RCP<Teuchos::FancyOStream> fos = Teuchos::getFancyOStream(Teuchos::rcpFromRef(std::cout));
217
218 SC zero = Teuchos::ScalarTraits<SC>::zero(), one = Teuchos::ScalarTraits<SC>::one();
219
220 // extract parameters from internal parameter list
221 const ParameterList & pL = Factory::GetParameterList();
222 LocalOrdinal nSweeps = pL.get<LocalOrdinal>("Sweeps");
223 Scalar omega = pL.get<Scalar>("Damping factor");
224
225 bool bRangeThyraMode = rangeMapExtractor_->getThyraMode();
226 bool bDomainThyraMode = domainMapExtractor_->getThyraMode();
227
228 // wrap current solution vector in RCP
229 RCP<MultiVector> rcpX = Teuchos::rcpFromRef(X);
230 RCP<const MultiVector> rcpB = Teuchos::rcpFromRef(B);
231
232 // make sure that both rcpX and rcpB are BlockedMultiVector objects
233 bool bCopyResultX = false;
234 RCP<BlockedCrsMatrix> bA = Teuchos::rcp_dynamic_cast<BlockedCrsMatrix>(A_);
235 MUELU_TEST_FOR_EXCEPTION(bA.is_null() == true, Exceptions::RuntimeError, "MueLu::BlockedGaussSeidelSmoother::Apply(): A_ must be a BlockedCrsMatrix");
236 RCP<BlockedMultiVector> bX = Teuchos::rcp_dynamic_cast<BlockedMultiVector>(rcpX);
237 RCP<const BlockedMultiVector> bB = Teuchos::rcp_dynamic_cast<const BlockedMultiVector>(rcpB);
238
239 if(bX.is_null() == true) {
240 RCP<MultiVector> test = Teuchos::rcp(new BlockedMultiVector(bA->getBlockedDomainMap(),rcpX));
241 rcpX.swap(test);
242 bCopyResultX = true;
243 }
244
245 if(bB.is_null() == true) {
246 RCP<const MultiVector> test = Teuchos::rcp(new BlockedMultiVector(bA->getBlockedRangeMap(),rcpB));
247 rcpB.swap(test);
248 }
249
250 // we now can guarantee that X and B are blocked multi vectors
251 bX = Teuchos::rcp_dynamic_cast<BlockedMultiVector>(rcpX);
252 bB = Teuchos::rcp_dynamic_cast<const BlockedMultiVector>(rcpB);
253
254 // check the type of operator
255 RCP<Xpetra::ReorderedBlockedCrsMatrix<Scalar,LocalOrdinal,GlobalOrdinal,Node> > rbA = Teuchos::rcp_dynamic_cast<Xpetra::ReorderedBlockedCrsMatrix<Scalar,LocalOrdinal,GlobalOrdinal,Node> >(bA);
256 if(rbA.is_null() == false) {
257 // A is a ReorderedBlockedCrsMatrix
258 Teuchos::RCP<const Xpetra::BlockReorderManager > brm = rbA->getBlockReorderManager();
259
260 // check type of X vector
261 if(bX->getBlockedMap()->getNumMaps() != bA->getDomainMapExtractor()->NumMaps()) {
262 // X is a blocked multi vector but incompatible to the reordered blocked operator A
263 Teuchos::RCP<MultiVector> test =
264 buildReorderedBlockedMultiVector(brm, bX);
265 rcpX.swap(test);
266 }
267 if(bB->getBlockedMap()->getNumMaps() != bA->getRangeMapExtractor()->NumMaps()) {
268 // B is a blocked multi vector but incompatible to the reordered blocked operator A
269 Teuchos::RCP<const MultiVector> test =
270 buildReorderedBlockedMultiVector(brm, bB);
271 rcpB.swap(test);
272 }
273 }
274
275 // Throughout the rest of the algorithm rcpX and rcpB are used for solution vector and RHS
276
277 // create residual vector
278 // contains current residual of current solution X with rhs B
279 RCP<MultiVector> residual = MultiVectorFactory::Build(rcpB->getMap(), rcpB->getNumVectors());
280 RCP<BlockedMultiVector> bresidual = Teuchos::rcp_dynamic_cast<BlockedMultiVector>(residual);
281 Teuchos::RCP<MultiVector> r1 = bresidual->getMultiVector(0,bRangeThyraMode);
282 Teuchos::RCP<MultiVector> r2 = bresidual->getMultiVector(1,bRangeThyraMode);
283
284 // helper vector 1
285 RCP<MultiVector> xtilde = MultiVectorFactory::Build(rcpX->getMap(), rcpX->getNumVectors());
286 RCP<BlockedMultiVector> bxtilde = Teuchos::rcp_dynamic_cast<BlockedMultiVector>(xtilde);
287 RCP<MultiVector> xtilde1 = bxtilde->getMultiVector(0,bDomainThyraMode);
288 RCP<MultiVector> xtilde2 = bxtilde->getMultiVector(1,bDomainThyraMode);
289
290 // incrementally improve solution vector X
291 for (LocalOrdinal run = 0; run < nSweeps; ++run) {
292 // 1) calculate current residual
293 residual->update(one,*rcpB,zero); // residual = B
294 A_->apply(*rcpX, *residual, Teuchos::NO_TRANS, -one, one);
295
296 // split residual vector
297 //Teuchos::RCP<MultiVector> r1 = rangeMapExtractor_->ExtractVector(residual, 0, bRangeThyraMode);
298 //Teuchos::RCP<MultiVector> r2 = rangeMapExtractor_->ExtractVector(residual, 1, bRangeThyraMode);
299
300 // 2) solve F * \Delta \tilde{x}_1 = r_1
301 // start with zero guess \Delta \tilde{x}_1
302 //RCP<MultiVector> xtilde1 = MultiVectorFactory::Build(r1->getMap(),rcpX->getNumVectors(),true);
303 //RCP<MultiVector> xtilde2 = MultiVectorFactory::Build(r2->getMap(),rcpX->getNumVectors(),true);
304 bxtilde->putScalar(zero);
305 velPredictSmoo_->Apply(*xtilde1,*r1);
306
307 // 3) solve SchurComp equation
308 // start with zero guess \Delta \tilde{x}_2
309 schurCompSmoo_->Apply(*xtilde2,*r2);
310#if 1
311 // 4) update solution vector
312 rcpX->update(omega,*bxtilde,one);
313#else
314 Teuchos::RCP<MultiVector> x1 = domainMapExtractor_->ExtractVector(rcpX, 0, bDomainThyraMode);
315 Teuchos::RCP<MultiVector> x2 = domainMapExtractor_->ExtractVector(rcpX, 1, bDomainThyraMode);
316
317 // 5) update solution vector with increments xhat1 and xhat2
318 // rescale increment for x2 with omega_
319 x1->update(omega,*xtilde1,one); // x1 = x1_old + omega xtilde1
320 x2->update(omega,*xtilde2,one); // x2 = x2_old + omega xtilde2
321
322 // write back solution in global vector X
323 domainMapExtractor_->InsertVector(x1, 0, rcpX, bDomainThyraMode);
324 domainMapExtractor_->InsertVector(x2, 1, rcpX, bDomainThyraMode);
325#endif
326 }
327
328 if (bCopyResultX == true) {
329 RCP<MultiVector> Xmerged = bX->Merge();
330 X.update(one, *Xmerged, zero);
331 }
332 }
333
334 template <class Scalar,class LocalOrdinal, class GlobalOrdinal, class Node>
335 RCP<MueLu::SmootherPrototype<Scalar, LocalOrdinal, GlobalOrdinal, Node> >
339
340 template <class Scalar,class LocalOrdinal, class GlobalOrdinal, class Node>
342 std::ostringstream out;
344 out << "{type = " << type_ << "}";
345 return out.str();
346 }
347
348 template <class Scalar,class LocalOrdinal, class GlobalOrdinal, class Node>
351
352 if (verbLevel & Parameters0) {
353 out0 << "Prec. type: " << type_ << /*" Sweeps: " << nSweeps_ << " damping: " << omega_ <<*/ std::endl;
354 }
355
356 if (verbLevel & Debug) {
357 out0 << "IsSetup: " << Teuchos::toString(SmootherPrototype::IsSetup()) << std::endl;
358 }
359 }
360 template <class Scalar,class LocalOrdinal, class GlobalOrdinal, class Node>
362 // FIXME: This is a placeholder
363 return Teuchos::OrdinalTraits<size_t>::invalid();
364 }
365
366
367
368} // namespace MueLu
369
370#endif /* MUELU_INDEFBLOCKEDDIAGONALSMOOTHER_DEF_HPP_ */
#define MUELU_DESCRIBE
Helper macro for implementing Describable::describe() for BaseClass objects.
#define MUELU_TEST_FOR_EXCEPTION(throw_exception_test, Exception, msg)
MueLu::DefaultLocalOrdinal LocalOrdinal
MueLu::DefaultScalar Scalar
virtual std::string description() const
Return a simple one-line description of this object.
Exception indicating invalid cast attempted.
Exception throws to report errors in the internal logical of the program.
Timer to be used in factories. Similar to Monitor but with additional timers.
Cheap Blocked diagonal smoother for indefinite 2x2 block matrices.
std::string description() const
Return a simple one-line description of this object.
size_t getNodeSmootherComplexity() const
Get a rough estimate of cost per iteration.
void AddFactoryManager(RCP< const FactoryManagerBase > FactManager, int pos=0)
Add a factory manager for BraessSarazin internal SchurComplement handling.
void print(Teuchos::FancyOStream &out, const VerbLevel verbLevel=Default) const
Print the object with some verbosity level to an FancyOStream object.
void Apply(MultiVector &X, const MultiVector &B, bool InitialGuessIsZero=false) const
Apply the Braess Sarazin smoother.
RCP< const ParameterList > GetValidParameterList() const
Input.
Class that holds all level-specific information.
void DeclareInput(const std::string &ename, const FactoryBase *factory, const FactoryBase *requestedBy=NoFactory::get())
Callback from FactoryBase::CallDeclareInput() and FactoryBase::DeclareInput()
T & Get(const std::string &ename, const FactoryBase *factory=NoFactory::get())
Get data without decrementing associated storage counter (i.e., read-only access)....
virtual const Teuchos::ParameterList & GetParameterList() const =0
An exception safe way to call the method 'Level::SetFactoryManager()'.
bool IsSetup() const
Get the state of a smoother prototype.
Namespace for MueLu classes and methods.
@ Warnings0
Important warning messages (one line)
@ Debug
Print additional debugging information.
@ Parameters0
Print class parameters.