94 TEUCHOS_TEST_FOR_EXCEPTION(pos < 0,
Exceptions::RuntimeError,
"MueLu::BlockedJacobiSmoother::AddFactoryManager: parameter \'pos\' must not be negative! error.");
96 size_t myPos = Teuchos::as<size_t>(pos);
98 if (myPos < FactManager_.size()) {
100 FactManager_.at(myPos) = FactManager;
101 }
else if( myPos == FactManager_.size()) {
103 FactManager_.push_back(FactManager);
105 RCP<Teuchos::FancyOStream> out = Teuchos::fancyOStream(Teuchos::rcpFromRef(std::cout));
106 *out <<
"Warning: cannot add new FactoryManager at proper position " << pos <<
". The FactoryManager is just appended to the end. Check this!" << std::endl;
109 FactManager_.push_back(FactManager);
135 RCP<Teuchos::FancyOStream> out = Teuchos::fancyOStream(Teuchos::rcpFromRef(std::cout));
137 FactoryMonitor m(*
this,
"Setup blocked Jacobi Smoother", currentLevel);
141 A_ = Factory::Get< RCP<Matrix> >(currentLevel,
"A");
142 RCP<BlockedCrsMatrix> bA = Teuchos::rcp_dynamic_cast<BlockedCrsMatrix>(A_);
143 TEUCHOS_TEST_FOR_EXCEPTION(bA==Teuchos::null,
Exceptions::BadCast,
"MueLu::BlockedPFactory::Build: input matrix A is not of type BlockedCrsMatrix! error.");
146 TEUCHOS_TEST_FOR_EXCEPTION(bA->Rows() != FactManager_.size(),
Exceptions::RuntimeError,
"MueLu::BlockedJacobiSmoother::Setup: number of block rows of A is " << bA->Rows() <<
" and does not match number of SubFactoryManagers " << FactManager_.size() <<
". error.");
147 TEUCHOS_TEST_FOR_EXCEPTION(bA->Cols() != FactManager_.size(),
Exceptions::RuntimeError,
"MueLu::BlockedJacobiSmoother::Setup: number of block cols of A is " << bA->Cols() <<
" and does not match number of SubFactoryManagers " << FactManager_.size() <<
". error.");
150 rangeMapExtractor_ = bA->getRangeMapExtractor();
151 domainMapExtractor_ = bA->getDomainMapExtractor();
154 std::vector<Teuchos::RCP<const FactoryManagerBase> >::const_iterator it;
155 for(it = FactManager_.begin(); it!=FactManager_.end(); ++it) {
159 RCP<const SmootherBase> Smoo = currentLevel.
Get< RCP<SmootherBase> >(
"PreSmoother",(*it)->GetFactory(
"Smoother").get());
160 Inverse_.push_back(Smoo);
163 RCP<Matrix> Aii = currentLevel.
Get< RCP<Matrix> >(
"A",(*it)->GetFactory(
"A").get());
164 bIsBlockedOperator_.push_back(Teuchos::rcp_dynamic_cast<BlockedCrsMatrix>(Aii)!=Teuchos::null);
177 RCP<MultiVector> rcpDebugX = Teuchos::rcpFromRef(X);
178 RCP<const MultiVector> rcpDebugB = Teuchos::rcpFromRef(B);
179 RCP<BlockedMultiVector> rcpBDebugX = Teuchos::rcp_dynamic_cast<BlockedMultiVector>(rcpDebugX);
180 RCP<const BlockedMultiVector> rcpBDebugB = Teuchos::rcp_dynamic_cast<const BlockedMultiVector>(rcpDebugB);
182 if(rcpBDebugB.is_null() ==
false) {
189 if(rcpBDebugX.is_null() ==
false) {
199 SC zero = Teuchos::ScalarTraits<SC>::zero(), one = Teuchos::ScalarTraits<SC>::one();
202 RCP<MultiVector> rcpX = Teuchos::rcpFromRef(X);
203 RCP<const MultiVector> rcpB = Teuchos::rcpFromRef(B);
206 bool bCopyResultX =
false;
207 RCP<BlockedCrsMatrix> bA = Teuchos::rcp_dynamic_cast<BlockedCrsMatrix>(A_);
209 RCP<BlockedMultiVector> bX = Teuchos::rcp_dynamic_cast<BlockedMultiVector>(rcpX);
210 RCP<const BlockedMultiVector> bB = Teuchos::rcp_dynamic_cast<const BlockedMultiVector>(rcpB);
212 if(bX.is_null() ==
true) {
213 RCP<MultiVector> test = Teuchos::rcp(
new BlockedMultiVector(bA->getBlockedDomainMap(),rcpX));
218 if(bB.is_null() ==
true) {
219 RCP<const MultiVector> test = Teuchos::rcp(
new BlockedMultiVector(bA->getBlockedRangeMap(),rcpB));
224 bX = Teuchos::rcp_dynamic_cast<BlockedMultiVector>(rcpX);
225 bB = Teuchos::rcp_dynamic_cast<const BlockedMultiVector>(rcpB);
228 RCP<Xpetra::ReorderedBlockedCrsMatrix<Scalar,LocalOrdinal,GlobalOrdinal,Node> > rbA = Teuchos::rcp_dynamic_cast<Xpetra::ReorderedBlockedCrsMatrix<Scalar,LocalOrdinal,GlobalOrdinal,Node> >(bA);
229 if(rbA.is_null() ==
false) {
231 Teuchos::RCP<const Xpetra::BlockReorderManager > brm = rbA->getBlockReorderManager();
234 if(bX->getBlockedMap()->getNumMaps() != bA->getDomainMapExtractor()->NumMaps()) {
236 Teuchos::RCP<MultiVector> test =
237 buildReorderedBlockedMultiVector(brm, bX);
240 if(bB->getBlockedMap()->getNumMaps() != bA->getRangeMapExtractor()->NumMaps()) {
242 Teuchos::RCP<const MultiVector> test =
243 buildReorderedBlockedMultiVector(brm, bB);
288 RCP<MultiVector> residual = MultiVectorFactory::Build(rcpB->getMap(), rcpB->getNumVectors());
289 RCP<MultiVector> tempres = MultiVectorFactory::Build(rcpB->getMap(), rcpB->getNumVectors());
298 residual->update(1.0,*rcpB,0.0);
299 if(InitialGuessIsZero ==
false || run > 0){
300 bA->apply(*rcpX, *residual, Teuchos::NO_TRANS, -1.0, 1.0);
303 for(
size_t i = 0; i<Inverse_.size(); i++) {
306 bool bRangeThyraMode = rangeMapExtractor_->getThyraMode();
307 bool bDomainThyraMode = domainMapExtractor_->getThyraMode();
308 Teuchos::RCP<MultiVector> Xi = domainMapExtractor_->ExtractVector(rcpX, i, bDomainThyraMode);
309 Teuchos::RCP<MultiVector> ri = rangeMapExtractor_->ExtractVector(residual, i, bRangeThyraMode);
310 Teuchos::RCP<MultiVector> tXi = domainMapExtractor_->getVector(i, X.getNumVectors(), bDomainThyraMode);
313 Inverse_.at(i)->Apply(*tXi, *ri,
false);
316 if( InitialGuessIsZero && run == 0 ){
317 Xi->update(omega,*tXi,0.0);
319 Xi->update(omega,*tXi,1.0);
325 if (bCopyResultX ==
true) {
326 RCP<MultiVector> Xmerged = bX->Merge();
327 X.update(one, *Xmerged, zero);