From bf9d36b96cadc5691a4d1a94a23347923c0396b7 Mon Sep 17 00:00:00 2001 From: Vinh Dang Date: Tue, 20 Jun 2023 14:46:36 -0700 Subject: [PATCH 01/16] Remove one conj in rot --- packages/teuchos/numerics/src/Teuchos_BLAS.hpp | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/packages/teuchos/numerics/src/Teuchos_BLAS.hpp b/packages/teuchos/numerics/src/Teuchos_BLAS.hpp index 9aab75d92013..79bcc9234288 100644 --- a/packages/teuchos/numerics/src/Teuchos_BLAS.hpp +++ b/packages/teuchos/numerics/src/Teuchos_BLAS.hpp @@ -531,7 +531,7 @@ namespace Teuchos if (n <= 0) return; if (incx==1 && incy==1) { for(OrdinalType i=0; i Date: Fri, 23 Jun 2023 10:05:47 +0200 Subject: [PATCH 02/16] Modify MultiVector processing If we enter the BlockSmoothers with a reordered/nested matrix, things go wrong in the Apply method if it is called with plain multi vectors. Thus changing the initial processing of the vectors in the apply method. --- .../MueLu_BlockedGaussSeidelSmoother_def.hpp | 96 ++++++++---------- .../MueLu_BlockedJacobiSmoother_def.hpp | 98 ++++++++----------- .../MueLu_BraessSarazinSmoother_def.hpp | 58 ++++++++--- ...MueLu_IndefBlockedDiagonalSmoother_def.hpp | 58 ++++++++--- .../MueLu_SimpleSmoother_def.hpp | 58 ++++++++--- .../MueLu_UzawaSmoother_def.hpp | 58 ++++++++--- .../Xpetra_ReorderedBlockedCrsMatrix.hpp | 3 + 7 files changed, 256 insertions(+), 173 deletions(-) diff --git a/packages/muelu/src/Smoothers/BlockedSmoothers/MueLu_BlockedGaussSeidelSmoother_def.hpp b/packages/muelu/src/Smoothers/BlockedSmoothers/MueLu_BlockedGaussSeidelSmoother_def.hpp index 986dc5d27919..fea0bee9d9ff 100644 --- a/packages/muelu/src/Smoothers/BlockedSmoothers/MueLu_BlockedGaussSeidelSmoother_def.hpp +++ b/packages/muelu/src/Smoothers/BlockedSmoothers/MueLu_BlockedGaussSeidelSmoother_def.hpp @@ -201,85 +201,73 @@ namespace MueLu { // make sure that both rcpX and rcpB are BlockedMultiVector objects bool bCopyResultX = false; + bool bReorderX = false; + bool bReorderB = false; RCP bA = Teuchos::rcp_dynamic_cast(A_); MUELU_TEST_FOR_EXCEPTION(bA.is_null() == true, Exceptions::RuntimeError, "MueLu::BlockedGaussSeidelSmoother::Apply(): A_ must be a BlockedCrsMatrix"); RCP bX = Teuchos::rcp_dynamic_cast(rcpX); RCP bB = Teuchos::rcp_dynamic_cast(rcpB); - if(bX.is_null() == true) { - RCP test = Teuchos::rcp(new BlockedMultiVector(bA->getBlockedDomainMap(),rcpX)); - rcpX.swap(test); - bCopyResultX = true; + // check the type of operator + RCP > rbA = + Teuchos::rcp_dynamic_cast >(bA); + + if(rbA.is_null() == false) { + // A is a ReorderedBlockedCrsMatrix and thus uses nested maps, retrieve BlockedCrsMatrix and use plain blocked + // map for the construction of the blocked multivectors + + // check type of X vector + if (bX.is_null() == true) { + RCP test = Teuchos::rcp(new BlockedMultiVector(rbA->getBlockedCrsMatrix()->getBlockedDomainMap(), rcpX)); + rcpX.swap(test); + bCopyResultX = true; + bReorderX = true; + } + + // check type of B vector + if (bB.is_null() == true) { + RCP test = Teuchos::rcp(new BlockedMultiVector(rbA->getBlockedCrsMatrix()->getBlockedRangeMap(), rcpB)); + rcpB.swap(test); + bReorderB = true; + } } + else { + // A is a BlockedCrsMatrix and uses a plain blocked map + if (bX.is_null() == true) { + RCP test = Teuchos::rcp(new BlockedMultiVector(bA->getBlockedDomainMap(), rcpX)); + rcpX.swap(test); + bCopyResultX = true; + } - if(bB.is_null() == true) { - RCP test = Teuchos::rcp(new BlockedMultiVector(bA->getBlockedRangeMap(),rcpB)); - rcpB.swap(test); + if(bB.is_null() == true) { + RCP test = Teuchos::rcp(new BlockedMultiVector(bA->getBlockedRangeMap(),rcpB)); + rcpB.swap(test); + } } // we now can guarantee that X and B are blocked multi vectors bX = Teuchos::rcp_dynamic_cast(rcpX); bB = Teuchos::rcp_dynamic_cast(rcpB); - // check the type of operator - RCP > rbA = Teuchos::rcp_dynamic_cast >(bA); + // Finally we can do a reordering of the blocked multivectors if A is a ReorderedBlockedCrsMatrix if(rbA.is_null() == false) { - // A is a ReorderedBlockedCrsMatrix + Teuchos::RCP brm = rbA->getBlockReorderManager(); // check type of X vector - if(bX->getBlockedMap()->getNumMaps() != bA->getDomainMapExtractor()->NumMaps()) { + if(bX->getBlockedMap()->getNumMaps() != bA->getDomainMapExtractor()->NumMaps() || bReorderX) { // X is a blocked multi vector but incompatible to the reordered blocked operator A - Teuchos::RCP test = - buildReorderedBlockedMultiVector(brm, bX); + Teuchos::RCP test = buildReorderedBlockedMultiVector(brm, bX); rcpX.swap(test); } - if(bB->getBlockedMap()->getNumMaps() != bA->getRangeMapExtractor()->NumMaps()) { + + if(bB->getBlockedMap()->getNumMaps() != bA->getRangeMapExtractor()->NumMaps() || bReorderB) { // B is a blocked multi vector but incompatible to the reordered blocked operator A - Teuchos::RCP test = - buildReorderedBlockedMultiVector(brm, bB); + Teuchos::RCP test = buildReorderedBlockedMultiVector(brm, bB); rcpB.swap(test); } } - - // check the type of operator - /*RCP bA = Teuchos::rcp_dynamic_cast(A_); - MUELU_TEST_FOR_EXCEPTION(bA.is_null() == true, Exceptions::RuntimeError, "MueLu::BlockedGaussSeidelSmoother::Apply(): A_ must be a BlockedCrsMatrix"); - RCP > rbA = Teuchos::rcp_dynamic_cast >(bA); - if(rbA.is_null() == false) { - // A is a ReorderedBlockedCrsMatrix - Teuchos::RCP brm = rbA->getBlockReorderManager(); - - // check type of vectors - RCP bX = Teuchos::rcp_dynamic_cast(rcpX); - RCP bB = Teuchos::rcp_dynamic_cast(rcpB); - - // check type of X vector - if(bX.is_null() == false && bX->getBlockedMap()->getNumMaps() != bA->getDomainMapExtractor()->NumMaps()) { - RCP rbX = Teuchos::rcp_dynamic_cast(bX); - if(rbX.is_null() == true) { - // X is a blocked multi vector but not reordered - // However, A is a reordered blocked operator - // We have to make sure, that A and X use compatible maps - Teuchos::RCP test = - buildReorderedBlockedMultiVector(brm, bX); - rcpX.swap(test); - } - } - if(bB.is_null() == false && bB->getBlockedMap()->getNumMaps() != bA->getRangeMapExtractor()->NumMaps()) { - RCP rbB = Teuchos::rcp_dynamic_cast(bB); - if(rbB.is_null() == true) { - // B is a blocked multi vector but not reordered - // However, A is a reordered blocked operator - // We have to make sure, that A and X use compatible maps - Teuchos::RCP test = - buildReorderedBlockedMultiVector(brm, bB); - rcpB.swap(test); - } - } - }*/ - // Throughout the rest of the algorithm rcpX and rcpB are used for solution vector and RHS RCP residual = MultiVectorFactory::Build(rcpB->getMap(), rcpB->getNumVectors()); diff --git a/packages/muelu/src/Smoothers/BlockedSmoothers/MueLu_BlockedJacobiSmoother_def.hpp b/packages/muelu/src/Smoothers/BlockedSmoothers/MueLu_BlockedJacobiSmoother_def.hpp index f4d41d8d0d1b..b2d712bda5f6 100644 --- a/packages/muelu/src/Smoothers/BlockedSmoothers/MueLu_BlockedJacobiSmoother_def.hpp +++ b/packages/muelu/src/Smoothers/BlockedSmoothers/MueLu_BlockedJacobiSmoother_def.hpp @@ -203,85 +203,73 @@ namespace MueLu { // make sure that both rcpX and rcpB are BlockedMultiVector objects bool bCopyResultX = false; + bool bReorderX = false; + bool bReorderB = false; RCP bA = Teuchos::rcp_dynamic_cast(A_); - MUELU_TEST_FOR_EXCEPTION(bA.is_null() == true, Exceptions::RuntimeError, "MueLu::BlockedJacobiSmoother::Apply(): A_ must be a BlockedCrsMatrix"); + MUELU_TEST_FOR_EXCEPTION(bA.is_null() == true, Exceptions::RuntimeError, "MueLu::BlockedGaussSeidelSmoother::Apply(): A_ must be a BlockedCrsMatrix"); RCP bX = Teuchos::rcp_dynamic_cast(rcpX); RCP bB = Teuchos::rcp_dynamic_cast(rcpB); - if(bX.is_null() == true) { - RCP test = Teuchos::rcp(new BlockedMultiVector(bA->getBlockedDomainMap(),rcpX)); - rcpX.swap(test); - bCopyResultX = true; + // check the type of operator + RCP > rbA = + Teuchos::rcp_dynamic_cast >(bA); + + if(rbA.is_null() == false) { + // A is a ReorderedBlockedCrsMatrix and thus uses nested maps, retrieve BlockedCrsMatrix and use plain blocked + // map for the construction of the blocked multivectors + + // check type of X vector + if (bX.is_null() == true) { + RCP test = Teuchos::rcp(new BlockedMultiVector(rbA->getBlockedCrsMatrix()->getBlockedDomainMap(), rcpX)); + rcpX.swap(test); + bCopyResultX = true; + bReorderX = true; + } + + // check type of B vector + if (bB.is_null() == true) { + RCP test = Teuchos::rcp(new BlockedMultiVector(rbA->getBlockedCrsMatrix()->getBlockedRangeMap(), rcpB)); + rcpB.swap(test); + bReorderB = true; + } } + else { + // A is a BlockedCrsMatrix and uses a plain blocked map + if (bX.is_null() == true) { + RCP test = Teuchos::rcp(new BlockedMultiVector(bA->getBlockedDomainMap(), rcpX)); + rcpX.swap(test); + bCopyResultX = true; + } - if(bB.is_null() == true) { - RCP test = Teuchos::rcp(new BlockedMultiVector(bA->getBlockedRangeMap(),rcpB)); - rcpB.swap(test); + if(bB.is_null() == true) { + RCP test = Teuchos::rcp(new BlockedMultiVector(bA->getBlockedRangeMap(),rcpB)); + rcpB.swap(test); + } } // we now can guarantee that X and B are blocked multi vectors bX = Teuchos::rcp_dynamic_cast(rcpX); bB = Teuchos::rcp_dynamic_cast(rcpB); - // check the type of operator - RCP > rbA = Teuchos::rcp_dynamic_cast >(bA); + // Finally we can do a reordering of the blocked multivectors if A is a ReorderedBlockedCrsMatrix if(rbA.is_null() == false) { - // A is a ReorderedBlockedCrsMatrix + Teuchos::RCP brm = rbA->getBlockReorderManager(); // check type of X vector - if(bX->getBlockedMap()->getNumMaps() != bA->getDomainMapExtractor()->NumMaps()) { + if(bX->getBlockedMap()->getNumMaps() != bA->getDomainMapExtractor()->NumMaps() || bReorderX) { // X is a blocked multi vector but incompatible to the reordered blocked operator A - Teuchos::RCP test = - buildReorderedBlockedMultiVector(brm, bX); + Teuchos::RCP test = buildReorderedBlockedMultiVector(brm, bX); rcpX.swap(test); } - if(bB->getBlockedMap()->getNumMaps() != bA->getRangeMapExtractor()->NumMaps()) { + + if(bB->getBlockedMap()->getNumMaps() != bA->getRangeMapExtractor()->NumMaps() || bReorderB) { // B is a blocked multi vector but incompatible to the reordered blocked operator A - Teuchos::RCP test = - buildReorderedBlockedMultiVector(brm, bB); + Teuchos::RCP test = buildReorderedBlockedMultiVector(brm, bB); rcpB.swap(test); } } - - // check the type of operator - /*RCP bA = Teuchos::rcp_dynamic_cast(A_); - MUELU_TEST_FOR_EXCEPTION(bA.is_null() == true, Exceptions::RuntimeError, "MueLu::BlockedJacobiSmoother::Apply(): A_ must be a BlockedCrsMatrix"); - RCP > rbA = Teuchos::rcp_dynamic_cast >(bA); - if(rbA.is_null() == false) { - // A is a ReorderedBlockedCrsMatrix - Teuchos::RCP brm = rbA->getBlockReorderManager(); - - // check type of vectors - RCP bX = Teuchos::rcp_dynamic_cast(rcpX); - RCP bB = Teuchos::rcp_dynamic_cast(rcpB); - - // check type of X vector - if(bX.is_null() == false && bX->getBlockedMap()->getNumMaps() != bA->getDomainMapExtractor()->NumMaps()) { - RCP rbX = Teuchos::rcp_dynamic_cast(bX); - if(rbX.is_null() == true) { - // X is a blocked multi vector but not reordered - // However, A is a reordered blocked operator - // We have to make sure, that A and X use compatible maps - Teuchos::RCP test = - buildReorderedBlockedMultiVector(brm, bX); - rcpX.swap(test); - } - } - if(bB.is_null() == false && bB->getBlockedMap()->getNumMaps() != bA->getRangeMapExtractor()->NumMaps()) { - RCP rbB = Teuchos::rcp_dynamic_cast(bB); - if(rbB.is_null() == true) { - // B is a blocked multi vector but not reordered - // However, A is a reordered blocked operator - // We have to make sure, that A and X use compatible maps - Teuchos::RCP test = - buildReorderedBlockedMultiVector(brm, bB); - rcpB.swap(test); - } - } - }*/ - // Throughout the rest of the algorithm rcpX and rcpB are used for solution vector and RHS RCP residual = MultiVectorFactory::Build(rcpB->getMap(), rcpB->getNumVectors()); diff --git a/packages/muelu/src/Smoothers/BlockedSmoothers/MueLu_BraessSarazinSmoother_def.hpp b/packages/muelu/src/Smoothers/BlockedSmoothers/MueLu_BraessSarazinSmoother_def.hpp index d94a11c4bc88..d33ff1a04b1a 100644 --- a/packages/muelu/src/Smoothers/BlockedSmoothers/MueLu_BraessSarazinSmoother_def.hpp +++ b/packages/muelu/src/Smoothers/BlockedSmoothers/MueLu_BraessSarazinSmoother_def.hpp @@ -177,43 +177,69 @@ namespace MueLu { // make sure that both rcpX and rcpB are BlockedMultiVector objects bool bCopyResultX = false; + bool bReorderX = false; + bool bReorderB = false; RCP bA = Teuchos::rcp_dynamic_cast(A_); MUELU_TEST_FOR_EXCEPTION(bA.is_null() == true, Exceptions::RuntimeError, "MueLu::BlockedGaussSeidelSmoother::Apply(): A_ must be a BlockedCrsMatrix"); RCP bX = Teuchos::rcp_dynamic_cast(rcpX); RCP bB = Teuchos::rcp_dynamic_cast(rcpB); - if(bX.is_null() == true) { - RCP test = Teuchos::rcp(new BlockedMultiVector(bA->getBlockedDomainMap(),rcpX)); - rcpX.swap(test); - bCopyResultX = true; + // check the type of operator + RCP > rbA = + Teuchos::rcp_dynamic_cast >(bA); + + if(rbA.is_null() == false) { + // A is a ReorderedBlockedCrsMatrix and thus uses nested maps, retrieve BlockedCrsMatrix and use plain blocked + // map for the construction of the blocked multivectors + + // check type of X vector + if (bX.is_null() == true) { + RCP test = Teuchos::rcp(new BlockedMultiVector(rbA->getBlockedCrsMatrix()->getBlockedDomainMap(), rcpX)); + rcpX.swap(test); + bCopyResultX = true; + bReorderX = true; + } + + // check type of B vector + if (bB.is_null() == true) { + RCP test = Teuchos::rcp(new BlockedMultiVector(rbA->getBlockedCrsMatrix()->getBlockedRangeMap(), rcpB)); + rcpB.swap(test); + bReorderB = true; + } } + else { + // A is a BlockedCrsMatrix and uses a plain blocked map + if (bX.is_null() == true) { + RCP test = Teuchos::rcp(new BlockedMultiVector(bA->getBlockedDomainMap(), rcpX)); + rcpX.swap(test); + bCopyResultX = true; + } - if(bB.is_null() == true) { - RCP test = Teuchos::rcp(new BlockedMultiVector(bA->getBlockedRangeMap(),rcpB)); - rcpB.swap(test); + if(bB.is_null() == true) { + RCP test = Teuchos::rcp(new BlockedMultiVector(bA->getBlockedRangeMap(),rcpB)); + rcpB.swap(test); + } } // we now can guarantee that X and B are blocked multi vectors bX = Teuchos::rcp_dynamic_cast(rcpX); bB = Teuchos::rcp_dynamic_cast(rcpB); - // check the type of operator - RCP > rbA = Teuchos::rcp_dynamic_cast >(bA); + // Finally we can do a reordering of the blocked multivectors if A is a ReorderedBlockedCrsMatrix if(rbA.is_null() == false) { - // A is a ReorderedBlockedCrsMatrix + Teuchos::RCP brm = rbA->getBlockReorderManager(); // check type of X vector - if(bX->getBlockedMap()->getNumMaps() != bA->getDomainMapExtractor()->NumMaps()) { + if(bX->getBlockedMap()->getNumMaps() != bA->getDomainMapExtractor()->NumMaps() || bReorderX) { // X is a blocked multi vector but incompatible to the reordered blocked operator A - Teuchos::RCP test = - buildReorderedBlockedMultiVector(brm, bX); + Teuchos::RCP test = buildReorderedBlockedMultiVector(brm, bX); rcpX.swap(test); } - if(bB->getBlockedMap()->getNumMaps() != bA->getRangeMapExtractor()->NumMaps()) { + + if(bB->getBlockedMap()->getNumMaps() != bA->getRangeMapExtractor()->NumMaps() || bReorderB) { // B is a blocked multi vector but incompatible to the reordered blocked operator A - Teuchos::RCP test = - buildReorderedBlockedMultiVector(brm, bB); + Teuchos::RCP test = buildReorderedBlockedMultiVector(brm, bB); rcpB.swap(test); } } diff --git a/packages/muelu/src/Smoothers/BlockedSmoothers/MueLu_IndefBlockedDiagonalSmoother_def.hpp b/packages/muelu/src/Smoothers/BlockedSmoothers/MueLu_IndefBlockedDiagonalSmoother_def.hpp index 57901b262925..59803565f338 100644 --- a/packages/muelu/src/Smoothers/BlockedSmoothers/MueLu_IndefBlockedDiagonalSmoother_def.hpp +++ b/packages/muelu/src/Smoothers/BlockedSmoothers/MueLu_IndefBlockedDiagonalSmoother_def.hpp @@ -227,43 +227,69 @@ namespace MueLu { // make sure that both rcpX and rcpB are BlockedMultiVector objects bool bCopyResultX = false; + bool bReorderX = false; + bool bReorderB = false; RCP bA = Teuchos::rcp_dynamic_cast(A_); MUELU_TEST_FOR_EXCEPTION(bA.is_null() == true, Exceptions::RuntimeError, "MueLu::BlockedGaussSeidelSmoother::Apply(): A_ must be a BlockedCrsMatrix"); RCP bX = Teuchos::rcp_dynamic_cast(rcpX); RCP bB = Teuchos::rcp_dynamic_cast(rcpB); - if(bX.is_null() == true) { - RCP test = Teuchos::rcp(new BlockedMultiVector(bA->getBlockedDomainMap(),rcpX)); - rcpX.swap(test); - bCopyResultX = true; + // check the type of operator + RCP > rbA = + Teuchos::rcp_dynamic_cast >(bA); + + if(rbA.is_null() == false) { + // A is a ReorderedBlockedCrsMatrix and thus uses nested maps, retrieve BlockedCrsMatrix and use plain blocked + // map for the construction of the blocked multivectors + + // check type of X vector + if (bX.is_null() == true) { + RCP test = Teuchos::rcp(new BlockedMultiVector(rbA->getBlockedCrsMatrix()->getBlockedDomainMap(), rcpX)); + rcpX.swap(test); + bCopyResultX = true; + bReorderX = true; + } + + // check type of B vector + if (bB.is_null() == true) { + RCP test = Teuchos::rcp(new BlockedMultiVector(rbA->getBlockedCrsMatrix()->getBlockedRangeMap(), rcpB)); + rcpB.swap(test); + bReorderB = true; + } } + else { + // A is a BlockedCrsMatrix and uses a plain blocked map + if (bX.is_null() == true) { + RCP test = Teuchos::rcp(new BlockedMultiVector(bA->getBlockedDomainMap(), rcpX)); + rcpX.swap(test); + bCopyResultX = true; + } - if(bB.is_null() == true) { - RCP test = Teuchos::rcp(new BlockedMultiVector(bA->getBlockedRangeMap(),rcpB)); - rcpB.swap(test); + if(bB.is_null() == true) { + RCP test = Teuchos::rcp(new BlockedMultiVector(bA->getBlockedRangeMap(),rcpB)); + rcpB.swap(test); + } } // we now can guarantee that X and B are blocked multi vectors bX = Teuchos::rcp_dynamic_cast(rcpX); bB = Teuchos::rcp_dynamic_cast(rcpB); - // check the type of operator - RCP > rbA = Teuchos::rcp_dynamic_cast >(bA); + // Finally we can do a reordering of the blocked multivectors if A is a ReorderedBlockedCrsMatrix if(rbA.is_null() == false) { - // A is a ReorderedBlockedCrsMatrix + Teuchos::RCP brm = rbA->getBlockReorderManager(); // check type of X vector - if(bX->getBlockedMap()->getNumMaps() != bA->getDomainMapExtractor()->NumMaps()) { + if(bX->getBlockedMap()->getNumMaps() != bA->getDomainMapExtractor()->NumMaps() || bReorderX) { // X is a blocked multi vector but incompatible to the reordered blocked operator A - Teuchos::RCP test = - buildReorderedBlockedMultiVector(brm, bX); + Teuchos::RCP test = buildReorderedBlockedMultiVector(brm, bX); rcpX.swap(test); } - if(bB->getBlockedMap()->getNumMaps() != bA->getRangeMapExtractor()->NumMaps()) { + + if(bB->getBlockedMap()->getNumMaps() != bA->getRangeMapExtractor()->NumMaps() || bReorderB) { // B is a blocked multi vector but incompatible to the reordered blocked operator A - Teuchos::RCP test = - buildReorderedBlockedMultiVector(brm, bB); + Teuchos::RCP test = buildReorderedBlockedMultiVector(brm, bB); rcpB.swap(test); } } diff --git a/packages/muelu/src/Smoothers/BlockedSmoothers/MueLu_SimpleSmoother_def.hpp b/packages/muelu/src/Smoothers/BlockedSmoothers/MueLu_SimpleSmoother_def.hpp index 974fdf3031ab..63d472853594 100644 --- a/packages/muelu/src/Smoothers/BlockedSmoothers/MueLu_SimpleSmoother_def.hpp +++ b/packages/muelu/src/Smoothers/BlockedSmoothers/MueLu_SimpleSmoother_def.hpp @@ -258,41 +258,67 @@ namespace MueLu { // make sure that both rcpX and rcpB are BlockedMultiVector objects bool bCopyResultX = false; + bool bReorderX = false; + bool bReorderB = false; RCP bA = Teuchos::rcp_dynamic_cast(A_); - MUELU_TEST_FOR_EXCEPTION(bA.is_null() == true, Exceptions::RuntimeError, - "MueLu::BlockedGaussSeidelSmoother::Apply(): A_ must be a BlockedCrsMatrix"); + MUELU_TEST_FOR_EXCEPTION(bA.is_null() == true, Exceptions::RuntimeError, "MueLu::BlockedGaussSeidelSmoother::Apply(): A_ must be a BlockedCrsMatrix"); RCP bX = Teuchos::rcp_dynamic_cast(rcpX); RCP bB = Teuchos::rcp_dynamic_cast(rcpB); - if(bX.is_null() == true) { - RCP test = Teuchos::rcp(new BlockedMultiVector(bA->getBlockedDomainMap(),rcpX)); - rcpX.swap(test); - bCopyResultX = true; + // check the type of operator + RCP > rbA = + Teuchos::rcp_dynamic_cast >(bA); + + if(rbA.is_null() == false) { + // A is a ReorderedBlockedCrsMatrix and thus uses nested maps, retrieve BlockedCrsMatrix and use plain blocked + // map for the construction of the blocked multivectors + + // check type of X vector + if (bX.is_null() == true) { + RCP test = Teuchos::rcp(new BlockedMultiVector(rbA->getBlockedCrsMatrix()->getBlockedDomainMap(), rcpX)); + rcpX.swap(test); + bCopyResultX = true; + bReorderX = true; + } + + // check type of B vector + if (bB.is_null() == true) { + RCP test = Teuchos::rcp(new BlockedMultiVector(rbA->getBlockedCrsMatrix()->getBlockedRangeMap(), rcpB)); + rcpB.swap(test); + bReorderB = true; + } } + else { + // A is a BlockedCrsMatrix and uses a plain blocked map + if (bX.is_null() == true) { + RCP test = Teuchos::rcp(new BlockedMultiVector(bA->getBlockedDomainMap(), rcpX)); + rcpX.swap(test); + bCopyResultX = true; + } - if(bB.is_null() == true) { - RCP test = Teuchos::rcp(new BlockedMultiVector(bA->getBlockedRangeMap(),rcpB)); - rcpB.swap(test); + if(bB.is_null() == true) { + RCP test = Teuchos::rcp(new BlockedMultiVector(bA->getBlockedRangeMap(),rcpB)); + rcpB.swap(test); + } } // we now can guarantee that X and B are blocked multi vectors bX = Teuchos::rcp_dynamic_cast(rcpX); bB = Teuchos::rcp_dynamic_cast(rcpB); - // check the type of operator - RCP > rbA = - Teuchos::rcp_dynamic_cast >(bA); + // Finally we can do a reordering of the blocked multivectors if A is a ReorderedBlockedCrsMatrix if(rbA.is_null() == false) { - // A is a ReorderedBlockedCrsMatrix + Teuchos::RCP brm = rbA->getBlockReorderManager(); - // check vector types - if(bX->getBlockedMap()->getNumMaps() != bA->getDomainMapExtractor()->NumMaps()) { + // check type of X vector + if(bX->getBlockedMap()->getNumMaps() != bA->getDomainMapExtractor()->NumMaps() || bReorderX) { // X is a blocked multi vector but incompatible to the reordered blocked operator A Teuchos::RCP test = buildReorderedBlockedMultiVector(brm, bX); rcpX.swap(test); } - if(bB->getBlockedMap()->getNumMaps() != bA->getRangeMapExtractor()->NumMaps()) { + + if(bB->getBlockedMap()->getNumMaps() != bA->getRangeMapExtractor()->NumMaps() || bReorderB) { // B is a blocked multi vector but incompatible to the reordered blocked operator A Teuchos::RCP test = buildReorderedBlockedMultiVector(brm, bB); rcpB.swap(test); diff --git a/packages/muelu/src/Smoothers/BlockedSmoothers/MueLu_UzawaSmoother_def.hpp b/packages/muelu/src/Smoothers/BlockedSmoothers/MueLu_UzawaSmoother_def.hpp index 8214f5c44c51..3aa6ebb05aef 100644 --- a/packages/muelu/src/Smoothers/BlockedSmoothers/MueLu_UzawaSmoother_def.hpp +++ b/packages/muelu/src/Smoothers/BlockedSmoothers/MueLu_UzawaSmoother_def.hpp @@ -191,43 +191,69 @@ namespace MueLu { // make sure that both rcpX and rcpB are BlockedMultiVector objects bool bCopyResultX = false; + bool bReorderX = false; + bool bReorderB = false; RCP bA = Teuchos::rcp_dynamic_cast(A_); MUELU_TEST_FOR_EXCEPTION(bA.is_null() == true, Exceptions::RuntimeError, "MueLu::BlockedGaussSeidelSmoother::Apply(): A_ must be a BlockedCrsMatrix"); RCP bX = Teuchos::rcp_dynamic_cast(rcpX); RCP bB = Teuchos::rcp_dynamic_cast(rcpB); - if(bX.is_null() == true) { - RCP test = Teuchos::rcp(new BlockedMultiVector(bA->getBlockedDomainMap(),rcpX)); - rcpX.swap(test); - bCopyResultX = true; + // check the type of operator + RCP > rbA = + Teuchos::rcp_dynamic_cast >(bA); + + if(rbA.is_null() == false) { + // A is a ReorderedBlockedCrsMatrix and thus uses nested maps, retrieve BlockedCrsMatrix and use plain blocked + // map for the construction of the blocked multivectors + + // check type of X vector + if (bX.is_null() == true) { + RCP test = Teuchos::rcp(new BlockedMultiVector(rbA->getBlockedCrsMatrix()->getBlockedDomainMap(), rcpX)); + rcpX.swap(test); + bCopyResultX = true; + bReorderX = true; + } + + // check type of B vector + if (bB.is_null() == true) { + RCP test = Teuchos::rcp(new BlockedMultiVector(rbA->getBlockedCrsMatrix()->getBlockedRangeMap(), rcpB)); + rcpB.swap(test); + bReorderB = true; + } } + else { + // A is a BlockedCrsMatrix and uses a plain blocked map + if (bX.is_null() == true) { + RCP test = Teuchos::rcp(new BlockedMultiVector(bA->getBlockedDomainMap(), rcpX)); + rcpX.swap(test); + bCopyResultX = true; + } - if(bB.is_null() == true) { - RCP test = Teuchos::rcp(new BlockedMultiVector(bA->getBlockedRangeMap(),rcpB)); - rcpB.swap(test); + if(bB.is_null() == true) { + RCP test = Teuchos::rcp(new BlockedMultiVector(bA->getBlockedRangeMap(),rcpB)); + rcpB.swap(test); + } } // we now can guarantee that X and B are blocked multi vectors bX = Teuchos::rcp_dynamic_cast(rcpX); bB = Teuchos::rcp_dynamic_cast(rcpB); - // check the type of operator - RCP > rbA = Teuchos::rcp_dynamic_cast >(bA); + // Finally we can do a reordering of the blocked multivectors if A is a ReorderedBlockedCrsMatrix if(rbA.is_null() == false) { - // A is a ReorderedBlockedCrsMatrix + Teuchos::RCP brm = rbA->getBlockReorderManager(); // check type of X vector - if(bX->getBlockedMap()->getNumMaps() != bA->getDomainMapExtractor()->NumMaps()) { + if(bX->getBlockedMap()->getNumMaps() != bA->getDomainMapExtractor()->NumMaps() || bReorderX) { // X is a blocked multi vector but incompatible to the reordered blocked operator A - Teuchos::RCP test = - buildReorderedBlockedMultiVector(brm, bX); + Teuchos::RCP test = buildReorderedBlockedMultiVector(brm, bX); rcpX.swap(test); } - if(bB->getBlockedMap()->getNumMaps() != bA->getRangeMapExtractor()->NumMaps()) { + + if(bB->getBlockedMap()->getNumMaps() != bA->getRangeMapExtractor()->NumMaps() || bReorderB) { // B is a blocked multi vector but incompatible to the reordered blocked operator A - Teuchos::RCP test = - buildReorderedBlockedMultiVector(brm, bB); + Teuchos::RCP test = buildReorderedBlockedMultiVector(brm, bB); rcpB.swap(test); } } diff --git a/packages/xpetra/src/BlockedCrsMatrix/Xpetra_ReorderedBlockedCrsMatrix.hpp b/packages/xpetra/src/BlockedCrsMatrix/Xpetra_ReorderedBlockedCrsMatrix.hpp index ba21356adab8..06c9eb1a6be9 100644 --- a/packages/xpetra/src/BlockedCrsMatrix/Xpetra_ReorderedBlockedCrsMatrix.hpp +++ b/packages/xpetra/src/BlockedCrsMatrix/Xpetra_ReorderedBlockedCrsMatrix.hpp @@ -264,6 +264,9 @@ namespace Xpetra { /** \brief Returns internal BlockReorderManager object */ Teuchos::RCP getBlockReorderManager() { return brm_; } + /** \brief Returns internal unmodified BlockedCrsMatrix object */ + Teuchos::RCP > getBlockedCrsMatrix() { return fullOp_; } + // @} //! @name Overridden from Teuchos::Describable From 29417e5a74c1ef64c13e61963839863b44ecef84 Mon Sep 17 00:00:00 2001 From: "Shane W. D. Hart" Date: Tue, 27 Jun 2023 08:03:08 -0400 Subject: [PATCH 03/16] Changes to Tpetra for MSVC compilation --- packages/tpetra/core/src/Tpetra_BlockView.hpp | 24 +++++++++---------- 1 file changed, 12 insertions(+), 12 deletions(-) diff --git a/packages/tpetra/core/src/Tpetra_BlockView.hpp b/packages/tpetra/core/src/Tpetra_BlockView.hpp index 2831855f92cb..cdc4addc70d4 100644 --- a/packages/tpetra/core/src/Tpetra_BlockView.hpp +++ b/packages/tpetra/core/src/Tpetra_BlockView.hpp @@ -233,7 +233,7 @@ struct SCAL { { using x_value_type = typename std::decay::type; const IndexType span = static_cast (x.span()); - x_value_type *__restrict__ x_ptr(x.data()); + x_value_type *__restrict x_ptr(x.data()); #if defined(KOKKOS_ENABLE_PRAGMA_UNROLL) #pragma unroll #endif @@ -395,8 +395,8 @@ struct AXPY { using x_value_type = typename std::decay::type; using y_value_type = typename std::decay::type; const IndexType span = static_cast (y.span()); - const x_value_type *__restrict__ x_ptr(x.data()); - y_value_type *__restrict__ y_ptr(y.data()); + const x_value_type *__restrict x_ptr(x.data()); + y_value_type *__restrict y_ptr(y.data()); #if defined(KOKKOS_ENABLE_PRAGMA_UNROLL) #pragma unroll #endif @@ -468,8 +468,8 @@ struct COPY { const IndexType span = static_cast (x.span()); using x_value_type = typename std::decay::type; using y_value_type = typename std::decay::type; - const x_value_type *__restrict__ x_ptr(x.data()); - y_value_type *__restrict__ y_ptr(y.data()); + const x_value_type *__restrict x_ptr(x.data()); + y_value_type *__restrict y_ptr(y.data()); #if defined(KOKKOS_ENABLE_PRAGMA_UNROLL) #pragma unroll @@ -530,16 +530,16 @@ struct GEMV (A.extent (0)); const IndexType numCols = static_cast (A.extent (1)); - const A_value_type *__restrict__ A_ptr(A.data()); const IndexType as1(A.stride(1)); - const x_value_type *__restrict__ x_ptr(x.data()); - y_value_type *__restrict__ y_ptr(y.data()); + const A_value_type *__restrict A_ptr(A.data()); const IndexType as1(A.stride(1)); + const x_value_type *__restrict x_ptr(x.data()); + y_value_type *__restrict y_ptr(y.data()); #if defined(KOKKOS_ENABLE_PRAGMA_UNROLL) #pragma unroll #endif for (IndexType j=0;j (A.extent (0)); const IndexType numCols = static_cast (A.extent (1)); - const A_value_type *__restrict__ A_ptr(A.data()); const IndexType as0(A.stride(0)); - const x_value_type *__restrict__ x_ptr(x.data()); - y_value_type *__restrict__ y_ptr(y.data()); + const A_value_type *__restrict A_ptr(A.data()); const IndexType as0(A.stride(0)); + const x_value_type *__restrict x_ptr(x.data()); + y_value_type *__restrict y_ptr(y.data()); #if defined(KOKKOS_ENABLE_PRAGMA_UNROLL) #pragma unroll From a42c5075e9b089a8670bf25cb020e912c2dd50a2 Mon Sep 17 00:00:00 2001 From: Chris Siefert Date: Thu, 29 Jun 2023 09:17:20 -0600 Subject: [PATCH 04/16] MueLu: Adding ML style aggregate threshold scaling --- packages/muelu/doc/UsersGuide/masterList.xml | 10 ++++++++++ .../muelu/doc/UsersGuide/options_aggregation.tex | 2 ++ packages/muelu/doc/UsersGuide/paramlist.tex | 2 ++ packages/muelu/doc/UsersGuide/paramlist_hidden.tex | 2 ++ .../MueLu_CoalesceDropFactory_def.hpp | 14 ++++++++++++-- .../MueLu_ParameterListInterpreter_def.hpp | 2 ++ packages/muelu/src/MueCentral/MueLu_MasterList.cpp | 3 +++ .../muelu/src/Operators/MueLu_Maxwell1_def.hpp | 4 ++++ 8 files changed, 37 insertions(+), 2 deletions(-) diff --git a/packages/muelu/doc/UsersGuide/masterList.xml b/packages/muelu/doc/UsersGuide/masterList.xml index 62635869a04c..0289e6bf0ae7 100644 --- a/packages/muelu/doc/UsersGuide/masterList.xml +++ b/packages/muelu/doc/UsersGuide/masterList.xml @@ -521,6 +521,16 @@ aggregation: threshold + + aggregation: use ml scaling of drop tol + bool + false + Enables ML-style scaling of drop tol, where the drop tol halves with each successive level. + parameter not existing in ML + + + + aggregation: min agg size int diff --git a/packages/muelu/doc/UsersGuide/options_aggregation.tex b/packages/muelu/doc/UsersGuide/options_aggregation.tex index 828e8c7bdafb..cd41f9039b16 100644 --- a/packages/muelu/doc/UsersGuide/options_aggregation.tex +++ b/packages/muelu/doc/UsersGuide/options_aggregation.tex @@ -13,6 +13,8 @@ \cbb{aggregation: drop tol}{double}{0.0}{Connectivity dropping threshold for a graph used in aggregation.} +\cbb{aggregation: use ml scaling of drop tol}{bool}{false}{Enables ML-style scaling of drop tol, where the drop tol halves with each successive level.} + \cbb{aggregation: min agg size}{int}{2}{Minimum size of an aggregate.} \cbb{aggregation: max agg size}{int}{-1}{Maximum size of an aggregate (-1 means unlimited).} diff --git a/packages/muelu/doc/UsersGuide/paramlist.tex b/packages/muelu/doc/UsersGuide/paramlist.tex index 05fe27de602b..c6b72a876af8 100644 --- a/packages/muelu/doc/UsersGuide/paramlist.tex +++ b/packages/muelu/doc/UsersGuide/paramlist.tex @@ -66,6 +66,8 @@ \cbb{aggregation: drop tol}{double}{0.0}{Connectivity dropping threshold for a graph used in aggregation.} +\cbb{aggregation: use ml scaling of drop tol}{bool}{false}{Enables ML-style scaling of drop tol, where the drop tol halves with each successive level.} + \cbb{aggregation: min agg size}{int}{2}{Minimum size of an aggregate.} \cbb{aggregation: max agg size}{int}{-1}{Maximum size of an aggregate (-1 means unlimited).} diff --git a/packages/muelu/doc/UsersGuide/paramlist_hidden.tex b/packages/muelu/doc/UsersGuide/paramlist_hidden.tex index 1eb82b04aa18..215e31fabcf0 100644 --- a/packages/muelu/doc/UsersGuide/paramlist_hidden.tex +++ b/packages/muelu/doc/UsersGuide/paramlist_hidden.tex @@ -95,6 +95,8 @@ \cbb{aggregation: drop tol}{double}{0.0}{Connectivity dropping threshold for a graph used in aggregation.} +\cbb{aggregation: use ml scaling of drop tol}{bool}{false}{Enables ML-style scaling of drop tol, where the drop tol halves with each successive level.} + \cbb{aggregation: min agg size}{int}{2}{Minimum size of an aggregate.} \cbb{aggregation: max agg size}{int}{-1}{Maximum size of an aggregate (-1 means unlimited).} diff --git a/packages/muelu/src/Graph/MatrixTransformation/MueLu_CoalesceDropFactory_def.hpp b/packages/muelu/src/Graph/MatrixTransformation/MueLu_CoalesceDropFactory_def.hpp index 24f454c7e99b..a26a47b7e971 100644 --- a/packages/muelu/src/Graph/MatrixTransformation/MueLu_CoalesceDropFactory_def.hpp +++ b/packages/muelu/src/Graph/MatrixTransformation/MueLu_CoalesceDropFactory_def.hpp @@ -122,6 +122,7 @@ namespace MueLu { #define SET_VALID_ENTRY(name) validParamList->setEntry(name, MasterList::getEntry(name)) SET_VALID_ENTRY("aggregation: drop tol"); + SET_VALID_ENTRY("aggregation: use ml scaling of drop tol"); SET_VALID_ENTRY("aggregation: Dirichlet threshold"); SET_VALID_ENTRY("aggregation: greedy Dirichlet"); SET_VALID_ENTRY("aggregation: row sum drop tol"); @@ -205,8 +206,10 @@ namespace MueLu { bool generateColoringGraph = false; // NOTE: If we're doing blockDiagonal, we'll not want to do rowSum twice (we'll do it - // in the block diagonalizaiton). So we'll clobber the rowSumTol with -1.0 in this case + // in the block diagonalization). So we'll clobber the rowSumTol with -1.0 in this case typename STS::magnitudeType rowSumTol = as(pL.get("aggregation: row sum drop tol")); + + RCP ghostedBlockNumber; ArrayRCP g_block_id; @@ -326,7 +329,14 @@ namespace MueLu { TEUCHOS_TEST_FOR_EXCEPTION(predrop_ != null && algo != "classical", Exceptions::RuntimeError, "Dropping function must not be provided for \"" << algo << "\" algorithm"); TEUCHOS_TEST_FOR_EXCEPTION(algo != "classical" && algo != "distance laplacian" && algo != "signed classical", Exceptions::RuntimeError, "\"algorithm\" must be one of (classical|distance laplacian|signed classical)"); - SC threshold = as(pL.get("aggregation: drop tol")); + SC threshold; + // If we're doing the ML-style halving of the drop tol at each level, we do that here. + if (pL.get("aggregation: use ml scaling of drop tol")) + threshold = pL.get("aggregation: drop tol") / pow(2.0,currentLevel.GetLevelID()); + else + threshold = as(pL.get("aggregation: drop tol")); + + std::string distanceLaplacianAlgoStr = pL.get("aggregation: distance laplacian algo"); std::string classicalAlgoStr = pL.get("aggregation: classical algo"); real_type realThreshold = STS::magnitude(threshold);// CMS: Rename this to "magnitude threshold" sometime diff --git a/packages/muelu/src/Interface/MueLu_ParameterListInterpreter_def.hpp b/packages/muelu/src/Interface/MueLu_ParameterListInterpreter_def.hpp index 002bb93108b4..ccf340ca11ed 100644 --- a/packages/muelu/src/Interface/MueLu_ParameterListInterpreter_def.hpp +++ b/packages/muelu/src/Interface/MueLu_ParameterListInterpreter_def.hpp @@ -1090,6 +1090,8 @@ namespace MueLu { MUELU_TEST_AND_SET_PARAM_2LIST(paramList, defaultList, "aggregation: row sum drop tol", double, dropParams); MUELU_TEST_AND_SET_PARAM_2LIST(paramList, defaultList, "aggregation: block diagonal: interleaved blocksize", int, dropParams); MUELU_TEST_AND_SET_PARAM_2LIST(paramList, defaultList, "aggregation: drop tol", double, dropParams); + MUELU_TEST_AND_SET_PARAM_2LIST(paramList, defaultList, "aggregation: use ml scaling of drop tol", bool, dropParams); + MUELU_TEST_AND_SET_PARAM_2LIST(paramList, defaultList, "aggregation: Dirichlet threshold", double, dropParams); MUELU_TEST_AND_SET_PARAM_2LIST(paramList, defaultList, "aggregation: greedy Dirichlet", bool, dropParams); MUELU_TEST_AND_SET_PARAM_2LIST(paramList, defaultList, "aggregation: distance laplacian algo", std::string, dropParams); diff --git a/packages/muelu/src/MueCentral/MueLu_MasterList.cpp b/packages/muelu/src/MueCentral/MueLu_MasterList.cpp index 4d1b80306bb9..262f30e274c2 100644 --- a/packages/muelu/src/MueCentral/MueLu_MasterList.cpp +++ b/packages/muelu/src/MueCentral/MueLu_MasterList.cpp @@ -221,6 +221,7 @@ namespace MueLu { "" "" "" + "" "" "" "" @@ -661,6 +662,8 @@ namespace MueLu { ("aggregation: threshold","aggregation: drop tol") + ("aggregation: use ml scaling of drop tol","aggregation: use ml scaling of drop tol") + ("aggregation: min agg size","aggregation: min agg size") ("aggregation: max agg size","aggregation: max agg size") diff --git a/packages/muelu/src/Operators/MueLu_Maxwell1_def.hpp b/packages/muelu/src/Operators/MueLu_Maxwell1_def.hpp index f20162ed5fcd..814650ac4222 100644 --- a/packages/muelu/src/Operators/MueLu_Maxwell1_def.hpp +++ b/packages/muelu/src/Operators/MueLu_Maxwell1_def.hpp @@ -112,6 +112,10 @@ namespace MueLu { newList.sublist("maxwell1: 11list").set("tentative: constant column sums", false); newList.sublist("maxwell1: 11list").set("tentative: calculate qr", false); + newList.sublist("maxwell1: 11list").set("aggregation: use ml scaling of drop tol",true); + newList.sublist("maxwell1: 22list").set("aggregation: use ml scaling of drop tol",true); + + if(list.isParameter("aggregation: damping factor") && list.get("aggregation: damping factor") == 0.0) newList.sublist("maxwell1: 11list").set("multigrid algorithm", "unsmoothed reitzinger"); else From d059523982f4bc39f0a37cefc212484681fbaef4 Mon Sep 17 00:00:00 2001 From: Christian Glusa Date: Thu, 29 Jun 2023 09:30:57 -0600 Subject: [PATCH 05/16] MueLu: match ML's drop tol in translated input decks --- .../muelu/src/Interface/MueLu_ML2MueLuParameterTranslator.cpp | 3 +++ 1 file changed, 3 insertions(+) diff --git a/packages/muelu/src/Interface/MueLu_ML2MueLuParameterTranslator.cpp b/packages/muelu/src/Interface/MueLu_ML2MueLuParameterTranslator.cpp index c87260076da6..6c84218ea930 100644 --- a/packages/muelu/src/Interface/MueLu_ML2MueLuParameterTranslator.cpp +++ b/packages/muelu/src/Interface/MueLu_ML2MueLuParameterTranslator.cpp @@ -399,6 +399,9 @@ namespace MueLu { // make sure that MueLu's phase2a matches ML's mueluss << "" << std::endl; + // make sure that MueLu's drop tol matches ML's + mueluss << "" << std::endl; + // special handling for energy minimization // TAW: this is not optimal for symmetric problems but at least works. // for symmetric problems the "energy minimization" parameter should not exist anyway... From 4a6f4be7b238e2df5f784f8f8409a5ec9dfe33a1 Mon Sep 17 00:00:00 2001 From: Brian Kelley Date: Thu, 29 Jun 2023 15:16:46 -0600 Subject: [PATCH 06/16] Tpetra: stop using deprecated Kokkos::Impl::ALL_t Use Kokkos::ALL_T instead --- .../tpetra/core/src/Tpetra_Details_WrappedDualView.hpp | 8 ++++---- packages/tpetra/core/src/Tpetra_MultiVector_def.hpp | 4 ++-- 2 files changed, 6 insertions(+), 6 deletions(-) diff --git a/packages/tpetra/core/src/Tpetra_Details_WrappedDualView.hpp b/packages/tpetra/core/src/Tpetra_Details_WrappedDualView.hpp index c2cd9932f75a..5d498707699a 100644 --- a/packages/tpetra/core/src/Tpetra_Details_WrappedDualView.hpp +++ b/packages/tpetra/core/src/Tpetra_Details_WrappedDualView.hpp @@ -202,12 +202,12 @@ class WrappedDualView { // 2D View Constructors - WrappedDualView(const WrappedDualView parent,const Kokkos::pair& rowRng, const Kokkos::Impl::ALL_t& colRng) { + WrappedDualView(const WrappedDualView parent,const Kokkos::pair& rowRng, const Kokkos::ALL_t& colRng) { originalDualView = parent.originalDualView; dualView = getSubview2D(parent.dualView,rowRng,colRng); } - WrappedDualView(const WrappedDualView parent,const Kokkos::Impl::ALL_t &rowRng, const Kokkos::pair& colRng) { + WrappedDualView(const WrappedDualView parent,const Kokkos::ALL_t &rowRng, const Kokkos::pair& colRng) { originalDualView = parent.originalDualView; dualView = getSubview2D(parent.dualView,rowRng,colRng); } @@ -589,12 +589,12 @@ class WrappedDualView { } template - ViewType getSubview2D(ViewType view, Kokkos::pair offset0, const Kokkos::Impl::ALL_t&) const { + ViewType getSubview2D(ViewType view, Kokkos::pair offset0, const Kokkos::ALL_t&) const { return Kokkos::subview(view,offset0,Kokkos::ALL()); } template - ViewType getSubview2D(ViewType view, const Kokkos::Impl::ALL_t&, Kokkos::pair offset1) const { + ViewType getSubview2D(ViewType view, const Kokkos::ALL_t&, Kokkos::pair offset1) const { return Kokkos::subview(view,Kokkos::ALL(),offset1); } diff --git a/packages/tpetra/core/src/Tpetra_MultiVector_def.hpp b/packages/tpetra/core/src/Tpetra_MultiVector_def.hpp index 84d5dbcf2c13..f9665901a0c1 100644 --- a/packages/tpetra/core/src/Tpetra_MultiVector_def.hpp +++ b/packages/tpetra/core/src/Tpetra_MultiVector_def.hpp @@ -262,7 +262,7 @@ namespace { // (anonymous) WrappedDualViewType takeSubview (const WrappedDualViewType& X, const std::pair& rowRng, - const Kokkos::Impl::ALL_t& colRng) + const Kokkos::ALL_t& colRng) { // The bug we saw below should be harder to trigger here. @@ -272,7 +272,7 @@ namespace { // (anonymous) template WrappedDualViewType takeSubview (const WrappedDualViewType& X, - const Kokkos::Impl::ALL_t& rowRng, + const Kokkos::ALL_t& rowRng, const std::pair& colRng) { using DualViewType = typename WrappedDualViewType::DVT; From 2feb5f35fb8f2cc8e45e35e162547cb490e37bdc Mon Sep 17 00:00:00 2001 From: Brian Kelley Date: Fri, 30 Jun 2023 14:56:28 -0600 Subject: [PATCH 07/16] Tpetra: use View::rank not View::Rank Kokkos::View::Rank is deprecated, use "rank" now. --- packages/tpetra/core/src/Tpetra_Details_Blas.hpp | 8 ++++---- .../core/src/Tpetra_Details_copyConvert.hpp | 16 ++++++++-------- packages/tpetra/core/src/Tpetra_Details_fill.hpp | 14 +++++++------- .../core/src/Tpetra_Details_getEntryOnHost.hpp | 4 ++-- .../Tpetra_Details_leftScaleLocalCrsMatrix.hpp | 2 +- .../tpetra/core/src/Tpetra_Details_normImpl.hpp | 4 ++-- .../Tpetra_Details_rightScaleLocalCrsMatrix.hpp | 2 +- .../tpetra/core/src/Tpetra_MultiVector_def.hpp | 2 +- .../Tpetra_computeRowAndColumnOneNorms_def.hpp | 2 +- ...ctor_Details_MultiVectorDistObjectKernels.hpp | 16 ++++++++-------- 10 files changed, 35 insertions(+), 35 deletions(-) diff --git a/packages/tpetra/core/src/Tpetra_Details_Blas.hpp b/packages/tpetra/core/src/Tpetra_Details_Blas.hpp index 2e8a1156c998..64f627a11c1d 100644 --- a/packages/tpetra/core/src/Tpetra_Details_Blas.hpp +++ b/packages/tpetra/core/src/Tpetra_Details_Blas.hpp @@ -90,7 +90,7 @@ template::value || std::is_same::value || std::is_same::value, @@ -115,7 +115,7 @@ struct GetStride1DView { static IndexType getStride (const ViewType& x) { - static_assert (ViewType::Rank == 1, "x must be a rank-1 Kokkos::View."); + static_assert (ViewType::rank == 1, "x must be a rank-1 Kokkos::View."); static_assert (std::is_same::value || std::is_same::value || std::is_same::value, @@ -137,7 +137,7 @@ struct GetStride1DView { static IndexType getStride (const ViewType&) { - static_assert (ViewType::Rank == 1, "x must be a rank-1 Kokkos::View."); + static_assert (ViewType::rank == 1, "x must be a rank-1 Kokkos::View."); static_assert (std::is_same::value, "ViewType::array_layout must be the same as array_layout."); static_assert (std::is_integral::value, @@ -153,7 +153,7 @@ struct GetStride1DView { static IndexType getStride (const ViewType&) { - static_assert (ViewType::Rank == 1, "x must be a rank-1 Kokkos::View."); + static_assert (ViewType::rank == 1, "x must be a rank-1 Kokkos::View."); static_assert (std::is_same::value, "ViewType::array_layout must be the same as array_layout."); static_assert (std::is_integral::value, diff --git a/packages/tpetra/core/src/Tpetra_Details_copyConvert.hpp b/packages/tpetra/core/src/Tpetra_Details_copyConvert.hpp index 57a15a4f0458..6a28de7f6841 100644 --- a/packages/tpetra/core/src/Tpetra_Details_copyConvert.hpp +++ b/packages/tpetra/core/src/Tpetra_Details_copyConvert.hpp @@ -130,7 +130,7 @@ namespace { // (anonymous) /// \tparam InputViewType Type of the input Kokkos::View. template (OutputViewType::Rank)> + const int rank = static_cast (OutputViewType::rank)> class CopyConvertFunctor {}; template { private: static_assert - (static_cast (OutputViewType::Rank) == 1 && - static_cast (InputViewType::Rank) == 1, + (static_cast (OutputViewType::rank) == 1 && + static_cast (InputViewType::rank) == 1, "CopyConvertFunctor (implements Tpetra::Details::copyConvert): " "OutputViewType and InputViewType must both have rank 1."); OutputViewType dst_; @@ -168,8 +168,8 @@ namespace { // (anonymous) private: static_assert - (static_cast (OutputViewType::Rank) == 2 && - static_cast (InputViewType::Rank) == 2, + (static_cast (OutputViewType::rank) == 2 && + static_cast (InputViewType::rank) == 2, "CopyConvertFunctor (implements Tpetra::Details::copyConvert): " "OutputViewType and InputViewType must both have rank 2."); OutputViewType dst_; @@ -369,8 +369,8 @@ copyConvert (const OutputViewType& dst, static_assert (std::is_same::value, "OutputViewType must be a nonconst Kokkos::View."); - static_assert (static_cast (OutputViewType::Rank) == - static_cast (InputViewType::Rank), + static_assert (static_cast (OutputViewType::rank) == + static_cast (InputViewType::rank), "src and dst must have the same rank."); if (dst.extent (0) != src.extent (0)) { @@ -381,7 +381,7 @@ copyConvert (const OutputViewType& dst, << "."; throw std::invalid_argument (os.str ()); } - if (static_cast (OutputViewType::Rank) > 1 && + if (static_cast (OutputViewType::rank) > 1 && dst.extent (1) != src.extent (1)) { std::ostringstream os; os << "Tpetra::Details::copyConvert: " diff --git a/packages/tpetra/core/src/Tpetra_Details_fill.hpp b/packages/tpetra/core/src/Tpetra_Details_fill.hpp index a8b7079d6776..4e98eb1ad4b4 100644 --- a/packages/tpetra/core/src/Tpetra_Details_fill.hpp +++ b/packages/tpetra/core/src/Tpetra_Details_fill.hpp @@ -71,7 +71,7 @@ template + const int rank = ViewType::rank> class Fill { public: static void @@ -98,7 +98,7 @@ class Fill { Fill (const ViewType& X, const ValueType& alpha) : X_ (X), alpha_ (alpha) { - static_assert (ViewType::Rank == 1, + static_assert (ViewType::rank == 1, "ViewType must be a rank-1 Kokkos::View."); static_assert (std::is_integral::value, "IndexType must be a built-in integer type."); @@ -139,7 +139,7 @@ class Fill { const IndexType numCols) : X_ (X), alpha_ (alpha), numCols_ (numCols) { - static_assert (ViewType::Rank == 2, + static_assert (ViewType::rank == 2, "ViewType must be a rank-2 Kokkos::View."); static_assert (std::is_integral::value, "IndexType must be a built-in integer type."); @@ -190,7 +190,7 @@ struct Fill::value, "IndexType must be a built-in integer type."); @@ -232,7 +232,7 @@ struct Fill::value, "IndexType must be a built-in integer type."); @@ -297,7 +297,7 @@ fill (const ExecutionSpace& execSpace, static_assert (std::is_integral::value, "IndexType must be a built-in integer type."); typedef Impl::Fill impl_type; + IndexType, ViewType::rank> impl_type; impl_type::fill (execSpace, X, alpha, numRows, numCols); } @@ -313,7 +313,7 @@ fill (const ExecutionSpace& execSpace, const IndexType numCols, const size_t whichVectors[]) { - static_assert (ViewType::Rank == 2, "ViewType must be a rank-2 " + static_assert (ViewType::rank == 2, "ViewType must be a rank-2 " "Kokkos::View in order to call the \"whichVectors\" " "specialization of fill."); static_assert (std::is_integral::value, diff --git a/packages/tpetra/core/src/Tpetra_Details_getEntryOnHost.hpp b/packages/tpetra/core/src/Tpetra_Details_getEntryOnHost.hpp index 51bbfcc0be0c..bae79227deef 100644 --- a/packages/tpetra/core/src/Tpetra_Details_getEntryOnHost.hpp +++ b/packages/tpetra/core/src/Tpetra_Details_getEntryOnHost.hpp @@ -61,7 +61,7 @@ getEntryOnHost (const ViewType& x, const IndexType ind) { using execution_space = typename ViewType::execution_space; - static_assert (ViewType::Rank == 1, "x must be a rank-1 Kokkos::View."); + static_assert (ViewType::rank == 1, "x must be a rank-1 Kokkos::View."); // Get a 0-D subview of the entry of the array, and copy to host scalar. typename ViewType::non_const_value_type val; // DEEP_COPY REVIEW - DEVICE-TO-VALUE @@ -76,7 +76,7 @@ getEntriesOnHost (const ViewType& x, const IndexType ind, const int count) { - static_assert (ViewType::Rank == 1, "x must be a rank-1 Kokkos::View."); + static_assert (ViewType::rank == 1, "x must be a rank-1 Kokkos::View."); // Get a 1-D subview of the entry of the array, and copy to host. auto subview = Kokkos::subview(x, Kokkos::make_pair(ind, ind + count)); return Kokkos::create_mirror_view_and_copy(typename ViewType::HostMirror::memory_space(), subview); diff --git a/packages/tpetra/core/src/Tpetra_Details_leftScaleLocalCrsMatrix.hpp b/packages/tpetra/core/src/Tpetra_Details_leftScaleLocalCrsMatrix.hpp index da6291febc9c..302301cd38da 100644 --- a/packages/tpetra/core/src/Tpetra_Details_leftScaleLocalCrsMatrix.hpp +++ b/packages/tpetra/core/src/Tpetra_Details_leftScaleLocalCrsMatrix.hpp @@ -73,7 +73,7 @@ class LeftScaleLocalCrsMatrix { using val_type = typename std::remove_const::type; using mag_type = typename ScalingFactorsViewType::non_const_value_type; - static_assert (ScalingFactorsViewType::Rank == 1, + static_assert (ScalingFactorsViewType::rank == 1, "scalingFactors must be a rank-1 Kokkos::View."); using device_type = typename LocalSparseMatrixType::device_type; diff --git a/packages/tpetra/core/src/Tpetra_Details_normImpl.hpp b/packages/tpetra/core/src/Tpetra_Details_normImpl.hpp index 95692326a0a3..43d90730832a 100644 --- a/packages/tpetra/core/src/Tpetra_Details_normImpl.hpp +++ b/packages/tpetra/core/src/Tpetra_Details_normImpl.hpp @@ -119,13 +119,13 @@ lclNormImpl (const RV& normsOut, using Kokkos::subview; using mag_type = typename RV::non_const_value_type; - static_assert (static_cast (RV::Rank) == 1, + static_assert (static_cast (RV::rank) == 1, "Tpetra::MultiVector::lclNormImpl: " "The first argument normsOut must have rank 1."); static_assert (Kokkos::is_view::value, "Tpetra::MultiVector::lclNormImpl: " "The second argument X is not a Kokkos::View."); - static_assert (static_cast (XMV::Rank) == 2, + static_assert (static_cast (XMV::rank) == 2, "Tpetra::MultiVector::lclNormImpl: " "The second argument X must have rank 2."); diff --git a/packages/tpetra/core/src/Tpetra_Details_rightScaleLocalCrsMatrix.hpp b/packages/tpetra/core/src/Tpetra_Details_rightScaleLocalCrsMatrix.hpp index 8d2e637ed0a7..219888164161 100644 --- a/packages/tpetra/core/src/Tpetra_Details_rightScaleLocalCrsMatrix.hpp +++ b/packages/tpetra/core/src/Tpetra_Details_rightScaleLocalCrsMatrix.hpp @@ -73,7 +73,7 @@ class RightScaleLocalCrsMatrix { using val_type = typename std::remove_const::type; using mag_type = typename ScalingFactorsViewType::non_const_value_type; - static_assert (ScalingFactorsViewType::Rank == 1, + static_assert (ScalingFactorsViewType::rank == 1, "scalingFactors must be a rank-1 Kokkos::View."); using device_type = typename LocalSparseMatrixType::device_type; diff --git a/packages/tpetra/core/src/Tpetra_MultiVector_def.hpp b/packages/tpetra/core/src/Tpetra_MultiVector_def.hpp index f9665901a0c1..72a992bba841 100644 --- a/packages/tpetra/core/src/Tpetra_MultiVector_def.hpp +++ b/packages/tpetra/core/src/Tpetra_MultiVector_def.hpp @@ -318,7 +318,7 @@ namespace { // (anonymous) // method yet, but its Views do. // NOTE: dv.stride() returns a vector of length one // more than its rank - size_t strides[WrappedOrNotDualViewType::t_dev::Rank+1]; + size_t strides[WrappedOrNotDualViewType::t_dev::rank+1]; dv.stride(strides); const size_t LDA = strides[1]; const size_t numRows = dv.extent (0); diff --git a/packages/tpetra/core/src/Tpetra_computeRowAndColumnOneNorms_def.hpp b/packages/tpetra/core/src/Tpetra_computeRowAndColumnOneNorms_def.hpp index fb68751b47d4..f7b7dfbd796e 100644 --- a/packages/tpetra/core/src/Tpetra_computeRowAndColumnOneNorms_def.hpp +++ b/packages/tpetra/core/src/Tpetra_computeRowAndColumnOneNorms_def.hpp @@ -780,7 +780,7 @@ copyMultiVectorColumnInto1DView ( template class FindZero { public: - static_assert (OneDViewType::Rank == 1, + static_assert (OneDViewType::rank == 1, "OneDViewType must be a rank-1 Kokkos::View."); static_assert (std::is_integral::value, "IndexType must be a built-in integer type."); diff --git a/packages/tpetra/core/src/kokkos_refactor/Tpetra_KokkosRefactor_Details_MultiVectorDistObjectKernels.hpp b/packages/tpetra/core/src/kokkos_refactor/Tpetra_KokkosRefactor_Details_MultiVectorDistObjectKernels.hpp index fe075ff1b9ae..e0cf496ae321 100644 --- a/packages/tpetra/core/src/kokkos_refactor/Tpetra_KokkosRefactor_Details_MultiVectorDistObjectKernels.hpp +++ b/packages/tpetra/core/src/kokkos_refactor/Tpetra_KokkosRefactor_Details_MultiVectorDistObjectKernels.hpp @@ -1643,8 +1643,8 @@ outOfBounds (const IntegerType x, const IntegerType exclusiveUpperBound) }; // clang-format on -// To do: Add enable_if<> restrictions on DstView::Rank == 1, -// SrcView::Rank == 2 +// To do: Add enable_if<> restrictions on DstView::rank == 1, +// SrcView::rank == 2 template void permute_array_multi_column(const ExecutionSpace &space, const DstView &dst, @@ -1657,8 +1657,8 @@ void permute_array_multi_column(const ExecutionSpace &space, const DstView &dst, } // clang-format off - // To do: Add enable_if<> restrictions on DstView::Rank == 1, - // SrcView::Rank == 2 + // To do: Add enable_if<> restrictions on DstView::rank == 1, + // SrcView::rank == 2 template void permute_array_multi_column(const DstView& dst, @@ -1738,8 +1738,8 @@ void permute_array_multi_column(const ExecutionSpace &space, const DstView &dst, }; // clang-format on -// To do: Add enable_if<> restrictions on DstView::Rank == 1, -// SrcView::Rank == 2 +// To do: Add enable_if<> restrictions on DstView::rank == 1, +// SrcView::rank == 2 template @@ -1756,8 +1756,8 @@ void permute_array_multi_column_variable_stride( } // clang-format off - // To do: Add enable_if<> restrictions on DstView::Rank == 1, - // SrcView::Rank == 2 + // To do: Add enable_if<> restrictions on DstView::rank == 1, + // SrcView::rank == 2 template From c5b22b457e7aacf11da3866b4d053d1089369e80 Mon Sep 17 00:00:00 2001 From: Brian Kelley Date: Fri, 30 Jun 2023 14:57:00 -0600 Subject: [PATCH 08/16] Tpetra: fix printf warning Cast values to match format specifier --- packages/tpetra/core/test/HashTable/FixedHashTableTest.cpp | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/packages/tpetra/core/test/HashTable/FixedHashTableTest.cpp b/packages/tpetra/core/test/HashTable/FixedHashTableTest.cpp index 072a9548a06c..482ec0da1233 100644 --- a/packages/tpetra/core/test/HashTable/FixedHashTableTest.cpp +++ b/packages/tpetra/core/test/HashTable/FixedHashTableTest.cpp @@ -188,12 +188,12 @@ namespace { // (anonymous) if (actualVal == Tpetra::Details::OrdinalTraits::invalid ()) { curBad = true; printf("get(key=%lld) = invalid, should have been %d\n", - key, expectedVal); + (long long) key, expectedVal); } else if (testValues && actualVal != expectedVal) { curBad = true; printf("get(key=%lld) = %d, should have been %d\n", - key, actualVal, expectedVal); + (long long) key, actualVal, expectedVal); } if (curBad) { From b9717e54818c2f09ed2aace15a23bf9bbb9170c3 Mon Sep 17 00:00:00 2001 From: Brian Kelley Date: Fri, 30 Jun 2023 14:57:18 -0600 Subject: [PATCH 09/16] Teuchos: use View::rank, not View::Rank In KokkosCompat, use View::rank because ::Rank is now deprecated. --- .../teuchos/kokkoscompat/src/KokkosCompat_View.hpp | 10 +++++----- 1 file changed, 5 insertions(+), 5 deletions(-) diff --git a/packages/teuchos/kokkoscompat/src/KokkosCompat_View.hpp b/packages/teuchos/kokkoscompat/src/KokkosCompat_View.hpp index ca7ffffd67b5..26c9166f05f2 100644 --- a/packages/teuchos/kokkoscompat/src/KokkosCompat_View.hpp +++ b/packages/teuchos/kokkoscompat/src/KokkosCompat_View.hpp @@ -129,7 +129,7 @@ namespace Kokkos { ViewType create_view (const std::string& label, size_t size) { static_assert(Kokkos::is_view::value==true,"Kokkos::Compat::create_view() called with non-view argument."); - static_assert(ViewType::Rank==1,"Kokkos::Compat::create_view() called with non-rank-1 view argument."); + static_assert(ViewType::rank==1,"Kokkos::Compat::create_view() called with non-rank-1 view argument."); return ViewType (label, size); } @@ -211,7 +211,7 @@ namespace Kokkos { const Ordinal end) { static_assert(Kokkos::is_view::value==true,"Kokkos::Compat::subview_range() called with non-view argument."); - static_assert(ViewType::Rank==1,"Kokkos::Compat::subview_range() called with non-rank-1 view argument."); + static_assert(ViewType::rank==1,"Kokkos::Compat::subview_range() called with non-rank-1 view argument."); return Kokkos::subview (view, std::make_pair (begin, end)); } @@ -222,7 +222,7 @@ namespace Kokkos { const Ordinal size) { static_assert(Kokkos::is_view::value==true,"Kokkos::Compat::subview_offset() called with non-view argument."); - static_assert(ViewType::Rank==1,"Kokkos::Compat::subview_offset() called with non-rank-1 view argument."); + static_assert(ViewType::rank==1,"Kokkos::Compat::subview_offset() called with non-rank-1 view argument."); return Kokkos::subview (view, std::make_pair (offset, offset+size)); } @@ -237,7 +237,7 @@ namespace Kokkos { { static_assert(Kokkos::is_view::value==true,"Kokkos::Compat::deep_copy_range() called with non-view argument."); static_assert(Kokkos::is_view::value==true,"Kokkos::Compat::deep_copy_range() called with non-view argument."); - static_assert(DstViewType::Rank==1 && SrcViewType::Rank==1,"Kokkos::Compat::deep_copy_range() called with non-rank-1 view argument."); + static_assert(DstViewType::rank==1 && SrcViewType::rank==1,"Kokkos::Compat::deep_copy_range() called with non-rank-1 view argument."); const Ordinal size = src_end - src_begin; const Ordinal dst_end = dst_begin + size; DstViewType dst_sub = Kokkos::subview( @@ -258,7 +258,7 @@ namespace Kokkos { { static_assert(Kokkos::is_view::value==true,"Kokkos::Compat::deep_copy_offset() called with non-view argument."); static_assert(Kokkos::is_view::value==true,"Kokkos::Compat::deep_copy_offset() called with non-view argument."); - static_assert(DstViewType::Rank==1 && SrcViewType::Rank==1,"Kokkos::Compat::deep_copy_offset() called with non-rank-1 view argument."); + static_assert(DstViewType::rank==1 && SrcViewType::rank==1,"Kokkos::Compat::deep_copy_offset() called with non-rank-1 view argument."); const Ordinal dst_end = dst_offset + size; const Ordinal src_end = src_offset + size; DstViewType dst_sub = Kokkos::subview( From 4043b28f933b3e030bb7290730613c37cfc33a04 Mon Sep 17 00:00:00 2001 From: Brian Kelley Date: Fri, 30 Jun 2023 17:31:29 -0600 Subject: [PATCH 10/16] KokkosKernels: remove non-existant include dir In CMake, don't export the binary directory kokkos-kernels/ode as an include dir. It doesn't exist since the ODE component has no auto-generated ETI files. --- packages/kokkos-kernels/ode/CMakeLists.txt | 5 ++++- 1 file changed, 4 insertions(+), 1 deletion(-) diff --git a/packages/kokkos-kernels/ode/CMakeLists.txt b/packages/kokkos-kernels/ode/CMakeLists.txt index 9d92dc07badf..b9cf089445e6 100644 --- a/packages/kokkos-kernels/ode/CMakeLists.txt +++ b/packages/kokkos-kernels/ode/CMakeLists.txt @@ -11,5 +11,8 @@ ENDIF() # Adding unit-tests -KOKKOSKERNELS_INCLUDE_DIRECTORIES(${CMAKE_CURRENT_BINARY_DIR}/ode) +# Note BMK: Since ODE has no auto-generated ETI files, this directory does not exist in a build without unit tests. +# This causes configure errors when building an app against a Trilinos install, and the unit test build dir doesn't contain any headers that need to be found. +# But uncomment the next line if ETI headers are ever added. +# KOKKOSKERNELS_INCLUDE_DIRECTORIES(${CMAKE_CURRENT_BINARY_DIR}/ode) KOKKOSKERNELS_INCLUDE_DIRECTORIES(REQUIRED_DURING_INSTALLATION_TESTING ${CMAKE_CURRENT_SOURCE_DIR}/ode) From 30f90f1f986bbee6c1e8aabaa9d2c793875fab5c Mon Sep 17 00:00:00 2001 From: maxfirmbach Date: Tue, 4 Jul 2023 20:36:58 +0200 Subject: [PATCH 11/16] Add unit testing routines --- .../unit_tests/Smoothers/BlockedSmoother.cpp | 1321 +++++++++++++++-- 1 file changed, 1168 insertions(+), 153 deletions(-) diff --git a/packages/muelu/test/unit_tests/Smoothers/BlockedSmoother.cpp b/packages/muelu/test/unit_tests/Smoothers/BlockedSmoother.cpp index bee397b6677b..b38fb4dfb530 100644 --- a/packages/muelu/test/unit_tests/Smoothers/BlockedSmoother.cpp +++ b/packages/muelu/test/unit_tests/Smoothers/BlockedSmoother.cpp @@ -63,6 +63,9 @@ #include #include +#include +#include +#include namespace MueLuTests { @@ -930,6 +933,238 @@ namespace MueLuTests { } // end UseTpetra } + TEUCHOS_UNIT_TEST_TEMPLATE_4_DECL(BlockedSmoother, NestedII30II12II_BGS_Setup_Apply2, Scalar, LocalOrdinal, GlobalOrdinal, Node) + { +# include + MUELU_TESTING_SET_OSTREAM; + MUELU_TESTING_LIMIT_SCOPE(Scalar,GlobalOrdinal,Node); + MUELU_TEST_ONLY_FOR(Xpetra::UseTpetra) { + + RCP > comm = Parameters::getDefaultComm(); + + int noBlocks = 4; + Teuchos::RCP bop = CreateBlockDiagonalExampleMatrix(noBlocks, *comm); + Teuchos::RCP Aconst = Teuchos::rcp_dynamic_cast(bop); + Teuchos::RCP< Matrix> A = Teuchos::rcp_const_cast(Aconst); + + //I don't use the testApply infrastructure because it has no provision for an initial guess. + Level level; TestHelpers::TestFactory::createSingleLevelHierarchy(level); + level.Set("A", A); + + // Test ReorderBlockAFactory + Teuchos::RCP rAFact = Teuchos::rcp(new ReorderBlockAFactory()); + rAFact->SetFactory("A",MueLu::NoFactory::getRCP()); + rAFact->SetParameter(std::string("Reorder Type"), Teuchos::ParameterEntry(std::string("[ [ 3 0 ] [1 2] ]"))); + + ////////////////////////////////////////////////////////////////////// + // global level + RCP smootherPrototype = rcp( new BlockedGaussSeidelSmoother() ); + smootherPrototype->SetParameter("Sweeps", Teuchos::ParameterEntry(1)); + smootherPrototype->SetFactory("A", rAFact); + + for(int k = 0; k < 2; k++ ) { + Teuchos::RCP sA = Teuchos::rcp(new SubBlockAFactory()); + sA->SetFactory("A",rAFact); + sA->SetParameter("block row",Teuchos::ParameterEntry(k)); // local block indices relative to size of blocked operator + sA->SetParameter("block col",Teuchos::ParameterEntry(k)); + + Teuchos::RCP sP = Teuchos::rcp( new BlockedGaussSeidelSmoother() ); + sP->SetParameter("Sweeps", Teuchos::ParameterEntry(1)); + sP->SetFactory("A", sA); + + for(int l = 0; l < 2; l++) { + std::string strInfo = std::string("{ 1 }"); + Teuchos::RCP ssA = rcp(new SubBlockAFactory()); + ssA->SetFactory("A",sA); + ssA->SetParameter("block row",Teuchos::ParameterEntry(l)); // local block indices relative to size of blocked operator + ssA->SetParameter("block col",Teuchos::ParameterEntry(l)); + ssA->SetParameter("Range map: Striding info", Teuchos::ParameterEntry(strInfo)); + ssA->SetParameter("Domain map: Striding info", Teuchos::ParameterEntry(strInfo)); + RCP ssP = rcp(new Ifpack2Smoother(std::string("RELAXATION"), Teuchos::ParameterList(), 0)); + ssP->SetFactory("A", ssA); + Teuchos::RCP ssF = Teuchos::rcp(new SmootherFactory(ssP)); + Teuchos::RCP ssM = Teuchos::rcp(new FactoryManager()); + ssM->SetFactory("A", ssA); + ssM->SetFactory("Smoother", ssF); + ssM->SetIgnoreUserData(true); + sP->AddFactoryManager(ssM,l); + } + + Teuchos::RCP sF = Teuchos::rcp(new SmootherFactory(sP)); + Teuchos::RCP sM = Teuchos::rcp(new FactoryManager()); + sM->SetFactory("A", sA); + sM->SetFactory("Smoother", sF); + sM->SetIgnoreUserData(true); + smootherPrototype->AddFactoryManager(sM,k); + } + + // create master smoother factory + RCP smootherFact = rcp( new SmootherFactory(smootherPrototype) ); + + // main factory manager + FactoryManager M; + M.SetFactory("A", rAFact); + M.SetFactory("Smoother", smootherFact); + + MueLu::SetFactoryManager SFM (Teuchos::rcpFromRef(level), Teuchos::rcpFromRef(M)); + + // request BGS smoother (and all dependencies) on level + level.Request("A", rAFact.get()); + level.Request("Smoother", smootherFact.get()); + level.Request("PreSmoother", smootherFact.get()); + level.Request("PostSmoother", smootherFact.get()); + + smootherFact->Build(level); + + RCP bgsSmoother = level.Get >("PreSmoother", smootherFact.get()); + + RCP reorderedA = level.Get >("A", rAFact.get()); + RCP reorderedbA = Teuchos::rcp_dynamic_cast(reorderedA); + + TEST_EQUALITY(reorderedbA->Rows(), 2); + TEST_EQUALITY(reorderedbA->Cols(), 2); + + RCP bX = Teuchos::rcp(new BlockedMultiVector(bop->getBlockedDomainMap(), 1, true)); + RCP bRHS = Teuchos::rcp(new BlockedMultiVector(bop->getBlockedRangeMap(), 1, true)); + + RCP X = bX->Merge(); + RCP RHS = bRHS->Merge(); + + // Random X + X->putScalar(0.0); + RHS->putScalar(1.0); + + bgsSmoother->Apply(*X, *RHS, true); //zero initial guess + + typedef typename Teuchos::ScalarTraits::magnitudeType magnitude_type; + Teuchos::Array finalNorms(1); X->norm2(finalNorms); + Teuchos::Array residualNorm1 = Utilities::ResidualNorm(*reorderedA, *X, *RHS); + out << " ||Residual_final|| = " << std::setiosflags(std::ios::fixed) << std::setprecision(20) << residualNorm1[0] << std::endl; + out << " ||X_final|| = " << std::setiosflags(std::ios::fixed) << std::setprecision(10) << finalNorms[0] << std::endl; + + magnitude_type tol = 50.*Teuchos::ScalarTraits::eps(); + + TEUCHOS_TEST_COMPARE(residualNorm1[0], <, tol, out, success); + + Teuchos::RCP doMapExtractor = reorderedbA->getDomainMapExtractor(); + + { + RCP test = Teuchos::rcp(new BlockedMultiVector(bop->getBlockedDomainMap(), X)); + bX.swap(test); + + Teuchos::RCP brm = Xpetra::blockedReorderFromString("[ [0 1] [2 3] ]"); + test = Teuchos::rcp_dynamic_cast(buildReorderedBlockedMultiVector(brm, bX)); + bX.swap(test); + } + + Teuchos::RCP v0 = doMapExtractor->ExtractVector(bX,0); + Teuchos::RCP v1 = doMapExtractor->ExtractVector(bX,1); + + Teuchos::RCP bv0 = Teuchos::rcp_dynamic_cast(v0); + Teuchos::RCP bv1 = Teuchos::rcp_dynamic_cast(v1); + TEST_EQUALITY(bv0.is_null(),false); + TEST_EQUALITY(bv1.is_null(),false); + + Teuchos::RCP bv00 = bv0->getMultiVector(0,false); + Teuchos::RCP bv01 = bv0->getMultiVector(1,false); + Teuchos::RCP bv10 = bv1->getMultiVector(0,false); + Teuchos::RCP bv11 = bv1->getMultiVector(1,false); + + TEST_EQUALITY((bv00->getData(0))[0], Teuchos::as(0.25)); + TEST_EQUALITY((bv01->getData(0))[0], Teuchos::as(0.25)); + TEST_EQUALITY((bv10->getData(0))[0], Teuchos::as(0.25)); + + Teuchos::Array n0(1); v0->norm1(n0); + Teuchos::Array n1(1); v1->norm1(n1); + + TEST_EQUALITY(n0[0], Teuchos::as(comm->getSize() * 2.5)); + TEUCHOS_TEST_COMPARE(n1[0], <, comm->getSize() * 13.34, out, success); + TEUCHOS_TEST_COMPARE(n1[0], >, comm->getSize() * 13.33, out, success); + + TEST_EQUALITY(v0->getMap()->getMinLocalIndex(), 0); + TEST_EQUALITY(v0->getMap()->getMaxLocalIndex(), 9); + TEST_EQUALITY(v0->getMap()->getMinAllGlobalIndex(), 0); + TEST_EQUALITY(v0->getMap()->getMinGlobalIndex(), comm->getRank() * 40); + TEST_EQUALITY(v0->getMap()->getMaxGlobalIndex(), comm->getRank() * 40 + 9); + TEST_EQUALITY(v0->getMap()->getMaxAllGlobalIndex(), comm->getSize() * 40 - 31); + + TEST_EQUALITY(v1->getMap()->getMinLocalIndex(), 0); + TEST_EQUALITY(v1->getMap()->getMaxLocalIndex(), 29); + TEST_EQUALITY(v1->getMap()->getMinAllGlobalIndex(), 10); + TEST_EQUALITY(v1->getMap()->getMinGlobalIndex(), comm->getRank() * 40 + 10); + TEST_EQUALITY(v1->getMap()->getMaxGlobalIndex(), comm->getRank() * 40 + 39); + TEST_EQUALITY(v1->getMap()->getMaxAllGlobalIndex(), comm->getSize() * 40 - 1); + + Teuchos::RCP b00 = Teuchos::rcp_dynamic_cast(reorderedbA->getMatrix(0,0)); + TEST_EQUALITY(b00->Rows(), 2); + TEST_EQUALITY(b00->Cols(), 2); + + Teuchos::RCP me00 = b00->getDomainMapExtractor(); + Teuchos::RCP v00 = me00->ExtractVector(v0,0); + Teuchos::RCP v01 = me00->ExtractVector(v0,1); + TEST_EQUALITY((v00->getData(0))[0], Teuchos::as(0.25)); + TEST_EQUALITY((v01->getData(0))[0], Teuchos::as(0.25)); + TEST_EQUALITY(v00->getLocalLength(), 5); + TEST_EQUALITY(v01->getLocalLength(), 5); + TEST_EQUALITY(v00->getGlobalLength(), Teuchos::as(comm->getSize() * 5)); + TEST_EQUALITY(v01->getGlobalLength(), Teuchos::as(comm->getSize() * 5)); + + TEST_EQUALITY(v00->getMap()->getMinLocalIndex(), 0); + TEST_EQUALITY(v00->getMap()->getMaxLocalIndex(), 4); + TEST_EQUALITY(v00->getMap()->getMinAllGlobalIndex(), 0); + TEST_EQUALITY(v00->getMap()->getMinGlobalIndex(), comm->getRank() * 40); + TEST_EQUALITY(v00->getMap()->getMaxGlobalIndex(), comm->getRank() * 40 + 4); + TEST_EQUALITY(v00->getMap()->getMaxAllGlobalIndex(), comm->getSize() * 40 - 36); + + TEST_EQUALITY(v01->getMap()->getMinLocalIndex(), 0); + TEST_EQUALITY(v01->getMap()->getMaxLocalIndex(), 4); + TEST_EQUALITY(v01->getMap()->getMinAllGlobalIndex(), 5); + TEST_EQUALITY(v01->getMap()->getMinGlobalIndex(), comm->getRank() * 40 + 5); + TEST_EQUALITY(v01->getMap()->getMaxGlobalIndex(), comm->getRank() * 40 + 9); + TEST_EQUALITY(v01->getMap()->getMaxAllGlobalIndex(), comm->getSize() * 40 - 31); + + + Teuchos::RCP b11 = Teuchos::rcp_dynamic_cast(reorderedbA->getMatrix(1,1)); + TEST_EQUALITY(b11->Rows(), 2); + TEST_EQUALITY(b11->Cols(), 2); + + Teuchos::RCP me11 = b11->getDomainMapExtractor(); + Teuchos::RCP v10 = me11->ExtractVector(v1,0); + Teuchos::RCP v11 = me11->ExtractVector(v1,1); + TEST_EQUALITY((v10->getData(0))[0], Teuchos::as(0.25)); + TEST_EQUALITY(v10->getLocalLength(), 10); + TEST_EQUALITY(v11->getLocalLength(), 20); + TEST_EQUALITY(v10->getGlobalLength(), Teuchos::as(comm->getSize() * 10)); + TEST_EQUALITY(v11->getGlobalLength(), Teuchos::as(comm->getSize() * 20)); + + TEST_EQUALITY(v10->getMap()->getMinLocalIndex(), 0); + TEST_EQUALITY(v10->getMap()->getMaxLocalIndex(), 9); + TEST_EQUALITY(v10->getMap()->getMinAllGlobalIndex(), 10); + TEST_EQUALITY(v10->getMap()->getMinGlobalIndex(), comm->getRank() * 40 + 10); + TEST_EQUALITY(v10->getMap()->getMaxGlobalIndex(), comm->getRank() * 40 + 19); + TEST_EQUALITY(v10->getMap()->getMaxAllGlobalIndex(), comm->getSize() * 40 - 21); + + TEST_EQUALITY(v11->getMap()->getMinLocalIndex(), 0); + TEST_EQUALITY(v11->getMap()->getMaxLocalIndex(), 19); + TEST_EQUALITY(v11->getMap()->getMinAllGlobalIndex(), 20); + TEST_EQUALITY(v11->getMap()->getMinGlobalIndex(), comm->getRank() * 40 + 20); + TEST_EQUALITY(v11->getMap()->getMaxGlobalIndex(), comm->getRank() * 40 + 39); + TEST_EQUALITY(v11->getMap()->getMaxAllGlobalIndex(), comm->getSize() * 40 - 1); + + Teuchos::Array n00(1); v00->norm1(n00); + Teuchos::Array n11(1); v01->norm1(n11); + Teuchos::Array n22(1); v10->norm1(n22); + Teuchos::Array n33(1); v11->norm1(n33); + + TEST_EQUALITY(n00[0], Teuchos::as(v00->getGlobalLength() * 0.25)); + TEST_EQUALITY(n11[0], Teuchos::as(v01->getGlobalLength() * 0.25)); + TEST_EQUALITY(n22[0], Teuchos::as(v10->getGlobalLength() * 0.25)); + TEUCHOS_TEST_COMPARE(n33[0], <, v11->getGlobalLength() * 0.54167, out, success); + TEUCHOS_TEST_COMPARE(n33[0], >, v11->getGlobalLength() * 0.54166, out, success); + + } // end UseTpetra + } + TEUCHOS_UNIT_TEST_TEMPLATE_4_DECL(BlockedSmoother, Thyra_BGS_Setup_Apply, Scalar, LocalOrdinal, GlobalOrdinal, Node) { # include @@ -1558,8 +1793,220 @@ namespace MueLuTests { RCP simpleSmoother = level.Get >("PreSmoother", smootherFact.get()); - RCP X = MultiVectorFactory::Build(A->getDomainMap(),1); - RCP RHS = MultiVectorFactory::Build(A->getRangeMap(),1); + RCP X = MultiVectorFactory::Build(A->getDomainMap(),1); + RCP RHS = MultiVectorFactory::Build(A->getRangeMap(),1); + + // Random X + X->setSeed(846930886); + X->randomize(); + + typedef typename Teuchos::ScalarTraits::magnitudeType magnitude_type; + + // Normalize X + Array norms(1); X->norm2(norms); + X->scale(1/norms[0]); + + // Compute RHS corresponding to X + A->apply(*X,*RHS, Teuchos::NO_TRANS,(SC)1.0,(SC)0.0); + + // Reset X to 0 + X->putScalar((SC) 0.0); + + RHS->norm2(norms); + out << "||RHS|| = " << std::setiosflags(std::ios::fixed) << std::setprecision(10) << norms[0] << std::endl; + + out << "solve with zero initial guess" << std::endl; + Teuchos::Array initialNorms(1); X->norm2(initialNorms); + out << " ||X_initial|| = " << std::setiosflags(std::ios::fixed) << std::setprecision(10) << initialNorms[0] << std::endl; + + simpleSmoother->Apply(*X, *RHS, true); //zero initial guess + + Teuchos::Array finalNorms(1); X->norm2(finalNorms); + Teuchos::Array residualNorm1 = Utilities::ResidualNorm(*A, *X, *RHS); + out << " ||Residual_final|| = " << std::setiosflags(std::ios::fixed) << std::setprecision(20) << residualNorm1[0] << std::endl; + out << " ||X_final|| = " << std::setiosflags(std::ios::fixed) << std::setprecision(10) << finalNorms[0] << std::endl; + + magnitude_type tol = 50.*Teuchos::ScalarTraits::eps(); + + TEUCHOS_TEST_COMPARE(residualNorm1[0], <, tol, out, success); + TEUCHOS_TEST_COMPARE(finalNorms[0] - Teuchos::ScalarTraits::magnitude(Teuchos::ScalarTraits::one()), <, tol, out, success); + + out << "solve with random initial guess" << std::endl; + X->randomize(); + X->norm2(initialNorms); + out << " ||X_initial|| = " << std::setiosflags(std::ios::fixed) << std::setprecision(10) << initialNorms[0] << std::endl; + + simpleSmoother->Apply(*X, *RHS, false); //nonzero initial guess + + X->norm2(finalNorms); + Teuchos::Array residualNorm2 = Utilities::ResidualNorm(*A, *X, *RHS); + out << " ||Residual_final|| = " << std::setiosflags(std::ios::fixed) << std::setprecision(20) << residualNorm2[0] << std::endl; + out << " ||X_final|| = " << std::setiosflags(std::ios::fixed) << std::setprecision(10) << finalNorms[0] << std::endl; + + TEUCHOS_TEST_COMPARE(residualNorm2[0], <, tol, out, success); + TEUCHOS_TEST_COMPARE(finalNorms[0] - Teuchos::ScalarTraits::magnitude(Teuchos::ScalarTraits::one()), <, tol, out, success); + + if (comm->getSize() == 1) { + TEST_EQUALITY(residualNorm1[0] != residualNorm2[0], true); + } else { + out << "Pass/Fail is only checked in serial." << std::endl; + } + } // end UseTpetra + } + + TEUCHOS_UNIT_TEST_TEMPLATE_4_DECL(BlockedSmoother, NestedII20I1I_SIMPLE_Setup_Apply, Scalar, LocalOrdinal, GlobalOrdinal, Node) + { +# include + MUELU_TESTING_SET_OSTREAM; + MUELU_TESTING_LIMIT_SCOPE(Scalar,GlobalOrdinal,Node); + // TODO test only Tpetra because of Ifpack2 smoother! + MUELU_TEST_ONLY_FOR(Xpetra::UseTpetra) { + + RCP > comm = Parameters::getDefaultComm(); + Xpetra::UnderlyingLib lib = MueLuTests::TestHelpers::Parameters::getLib(); + + Teuchos::RCP bop = TestHelpers::TestFactory::CreateBlockDiagonalExampleMatrix(lib,3, comm); + Teuchos::RCP Aconst = Teuchos::rcp_dynamic_cast(bop); + Teuchos::RCP< Matrix> A = Teuchos::rcp_const_cast(Aconst); + + //I don't use the testApply infrastructure because it has no provision for an initial guess. + Level level; TestHelpers::TestFactory::createSingleLevelHierarchy(level); + level.Set("A", A); + + // Test ReorderBlockAFactory + Teuchos::RCP rAFact = Teuchos::rcp(new ReorderBlockAFactory()); + rAFact->SetFactory("A",MueLu::NoFactory::getRCP()); + rAFact->SetParameter(std::string("Reorder Type"), Teuchos::ParameterEntry(std::string("[ [2 0] 1]"))); + + ////////////////////////////////////////////////////////////////////// + // Smoothers + RCP smootherPrototype = rcp( new SimpleSmoother() ); + smootherPrototype->SetFactory("A",rAFact); + smootherPrototype->SetParameter("Sweeps", Teuchos::ParameterEntry(Teuchos::as(1))); + smootherPrototype->SetParameter("Damping factor", Teuchos::ParameterEntry(Teuchos::as(1.0))); + smootherPrototype->SetParameter("UseSIMPLEC", Teuchos::ParameterEntry(false)); + + std::vector > sA (1, Teuchos::null); + std::vector > sF (2, Teuchos::null); + std::vector > sM (2, Teuchos::null); + + // prediction + std::string strInfo = std::string("{ 1 }"); + sA[0] = rcp(new SubBlockAFactory()); + sA[0]->SetFactory("A",rAFact); + sA[0]->SetParameter("block row",Teuchos::ParameterEntry(0)); + sA[0]->SetParameter("block col",Teuchos::ParameterEntry(0)); + sA[0]->SetParameter("Range map: Striding info", Teuchos::ParameterEntry(strInfo)); + sA[0]->SetParameter("Domain map: Striding info", Teuchos::ParameterEntry(strInfo)); + + // create a 2x2 SIMPLE for the prediction eq. + RCP smoProtoPredict = Teuchos::rcp( new SimpleSmoother() ); + smoProtoPredict->SetParameter("Sweeps", Teuchos::ParameterEntry(1)); + smoProtoPredict->SetParameter("Damping factor", Teuchos::ParameterEntry(Teuchos::as(1.0))); + smoProtoPredict->SetParameter("UseSIMPLEC", Teuchos::ParameterEntry(false)); + smoProtoPredict->SetFactory("A", sA[0]); + + for(int l = 0; l < 2; l++) { + Teuchos::RCP ssA = rcp(new SubBlockAFactory()); + ssA->SetFactory("A",sA[0]); + ssA->SetParameter("block row",Teuchos::ParameterEntry(l)); // local block indices relative to size of blocked operator + ssA->SetParameter("block col",Teuchos::ParameterEntry(l)); + ssA->SetParameter("Range map: Striding info", Teuchos::ParameterEntry(strInfo)); + ssA->SetParameter("Domain map: Striding info", Teuchos::ParameterEntry(strInfo)); + RCP ssP = rcp(new Ifpack2Smoother(std::string("RELAXATION"), Teuchos::ParameterList(), 0)); + ssP->SetFactory("A", ssA); + Teuchos::RCP ssF = Teuchos::rcp(new SmootherFactory(ssP)); + Teuchos::RCP ssM = Teuchos::rcp(new FactoryManager()); + ssM->SetFactory("A", ssA); + ssM->SetFactory("Smoother", ssF); + ssM->SetIgnoreUserData(true); + if(l == 0) smoProtoPredict->SetVelocityPredictionFactoryManager(ssM); + else smoProtoPredict->SetSchurCompFactoryManager(ssM); + } + + sF[0] = rcp( new SmootherFactory(smoProtoPredict) ); + + sM[0] = rcp(new FactoryManager()); + sM[0]->SetFactory("A", sA[0]); + sM[0]->SetFactory("Smoother", sF[0]); + sM[0]->SetIgnoreUserData(true); + + smootherPrototype->SetVelocityPredictionFactoryManager(sM[0]); + + // correction + // define SchurComplement Factory + // SchurComp gets a RCP to AFact_ which has to be the 2x2 blocked operator + // It stores the resulting SchurComplement operator as "A" generated by the SchurComplementFactory + // Instead of F^{-1} it uses the approximation \hat{F}^{-1} with \hat{F} = diag(F) + RCP AinvFact = Teuchos::rcp(new InverseApproximationFactory()); + AinvFact->SetFactory("A",rAFact); + + RCP SFact = Teuchos::rcp(new SchurComplementFactory()); + SFact->SetParameter("omega", Teuchos::ParameterEntry(Teuchos::as(1.0))); // for Simple, omega is always 1.0 in the SchurComplement + SFact->SetFactory("A",rAFact); + SFact->SetFactory("Ainv", AinvFact); + + RCP smoProtoCorrect = rcp(new Ifpack2Smoother(std::string("RELAXATION"), Teuchos::ParameterList(), 0)); + smoProtoCorrect->SetFactory("A", SFact); + sF[1] = rcp( new SmootherFactory(smoProtoCorrect) ); + + sM[1] = rcp(new FactoryManager()); + sM[1]->SetFactory("A", SFact); + sM[1]->SetFactory("Smoother", sF[1]); + sM[1]->SetIgnoreUserData(true); + + smootherPrototype->SetSchurCompFactoryManager(sM[1]); + + + RCP smootherFact = rcp( new SmootherFactory(smootherPrototype) ); + + // main factory manager + FactoryManager M; + M.SetFactory("Smoother", smootherFact); + M.SetFactory("A", rAFact); + + MueLu::SetFactoryManager SFM (Teuchos::rcpFromRef(level), Teuchos::rcpFromRef(M)); + + // request BGS smoother (and all dependencies) on level + level.Request("A", rAFact.get()); + level.Request("Smoother", smootherFact.get()); + level.Request("PreSmoother", smootherFact.get()); + level.Request("PostSmoother", smootherFact.get()); + + //smootherFact->DeclareInput(level); + smootherFact->Build(level); + + level.print(std::cout, Teuchos::VERB_EXTREME); + + RCP simpleSmoother = level.Get >("PreSmoother", smootherFact.get()); + + RCP reorderedA = level.Get >("A", rAFact.get()); + RCP reorderedbA = Teuchos::rcp_dynamic_cast(reorderedA); + + TEST_EQUALITY(reorderedbA->Rows(), 2); + TEST_EQUALITY(reorderedbA->Cols(), 2); + + RCP X = MultiVectorFactory::Build(reorderedA->getDomainMap(),1); + RCP RHS = MultiVectorFactory::Build(reorderedA->getRangeMap(),1); + + // apply simple smoother + RHS->putScalar((SC) 1.0); + X->putScalar((SC) 0.0); + + // solve system + simpleSmoother->Apply(*X, *RHS, true); //zero initial guess + + RCP bX = Teuchos::rcp_dynamic_cast(X); + TEST_EQUALITY(bX.is_null(), false); + RCP XX = bX->Merge(); + Teuchos::ArrayRCP xdata = XX->getData(0); + bool bCheck = true; + for(size_t i=0; igetLocalLength(); i++) { + if (i< 10) { if(xdata[i] != (SC) (1.0/3.0)) bCheck = false; } + if (i>=10 && i< 15) { if(xdata[i] != (SC) 1.0) bCheck = false; } + if (i>=15 && i< 20) { if(xdata[i] != (SC) 0.5) bCheck = false; } + } + TEST_EQUALITY(bCheck, true); // Random X X->setSeed(846930886); @@ -1572,7 +2019,7 @@ namespace MueLuTests { X->scale(1/norms[0]); // Compute RHS corresponding to X - A->apply(*X,*RHS, Teuchos::NO_TRANS,(SC)1.0,(SC)0.0); + reorderedA->apply(*X,*RHS, Teuchos::NO_TRANS,(SC)1.0,(SC)0.0); // Reset X to 0 X->putScalar((SC) 0.0); @@ -1587,7 +2034,7 @@ namespace MueLuTests { simpleSmoother->Apply(*X, *RHS, true); //zero initial guess Teuchos::Array finalNorms(1); X->norm2(finalNorms); - Teuchos::Array residualNorm1 = Utilities::ResidualNorm(*A, *X, *RHS); + Teuchos::Array residualNorm1 = Utilities::ResidualNorm(*reorderedA, *X, *RHS); out << " ||Residual_final|| = " << std::setiosflags(std::ios::fixed) << std::setprecision(20) << residualNorm1[0] << std::endl; out << " ||X_final|| = " << std::setiosflags(std::ios::fixed) << std::setprecision(10) << finalNorms[0] << std::endl; @@ -1595,31 +2042,10 @@ namespace MueLuTests { TEUCHOS_TEST_COMPARE(residualNorm1[0], <, tol, out, success); TEUCHOS_TEST_COMPARE(finalNorms[0] - Teuchos::ScalarTraits::magnitude(Teuchos::ScalarTraits::one()), <, tol, out, success); - - out << "solve with random initial guess" << std::endl; - X->randomize(); - X->norm2(initialNorms); - out << " ||X_initial|| = " << std::setiosflags(std::ios::fixed) << std::setprecision(10) << initialNorms[0] << std::endl; - - simpleSmoother->Apply(*X, *RHS, false); //nonzero initial guess - - X->norm2(finalNorms); - Teuchos::Array residualNorm2 = Utilities::ResidualNorm(*A, *X, *RHS); - out << " ||Residual_final|| = " << std::setiosflags(std::ios::fixed) << std::setprecision(20) << residualNorm2[0] << std::endl; - out << " ||X_final|| = " << std::setiosflags(std::ios::fixed) << std::setprecision(10) << finalNorms[0] << std::endl; - - TEUCHOS_TEST_COMPARE(residualNorm2[0], <, tol, out, success); - TEUCHOS_TEST_COMPARE(finalNorms[0] - Teuchos::ScalarTraits::magnitude(Teuchos::ScalarTraits::one()), <, tol, out, success); - - if (comm->getSize() == 1) { - TEST_EQUALITY(residualNorm1[0] != residualNorm2[0], true); - } else { - out << "Pass/Fail is only checked in serial." << std::endl; - } - } // end UseTpetra + }// end useTpetra } - TEUCHOS_UNIT_TEST_TEMPLATE_4_DECL(BlockedSmoother, NestedII20I1I_SIMPLE_Setup_Apply, Scalar, LocalOrdinal, GlobalOrdinal, Node) + TEUCHOS_UNIT_TEST_TEMPLATE_4_DECL(BlockedSmoother, NestedI2I01II_SIMPLE_Setup_Apply, Scalar, LocalOrdinal, GlobalOrdinal, Node) { # include MUELU_TESTING_SET_OSTREAM; @@ -1641,7 +2067,7 @@ namespace MueLuTests { // Test ReorderBlockAFactory Teuchos::RCP rAFact = Teuchos::rcp(new ReorderBlockAFactory()); rAFact->SetFactory("A",MueLu::NoFactory::getRCP()); - rAFact->SetParameter(std::string("Reorder Type"), Teuchos::ParameterEntry(std::string("[ [2 0] 1]"))); + rAFact->SetParameter(std::string("Reorder Type"), Teuchos::ParameterEntry(std::string("[ 2 [0 1]]"))); ////////////////////////////////////////////////////////////////////// // Smoothers @@ -1664,16 +2090,40 @@ namespace MueLuTests { sA[0]->SetParameter("Range map: Striding info", Teuchos::ParameterEntry(strInfo)); sA[0]->SetParameter("Domain map: Striding info", Teuchos::ParameterEntry(strInfo)); + RCP smoProtoCorrect = rcp(new Ifpack2Smoother(std::string("RELAXATION"), Teuchos::ParameterList(), 0)); + smoProtoCorrect->SetFactory("A", sA[0]); + sF[0] = rcp( new SmootherFactory(smoProtoCorrect) ); + + sM[0] = rcp(new FactoryManager()); + sM[0]->SetFactory("A", sA[0]); + sM[0]->SetFactory("Smoother", sF[0]); + sM[0]->SetIgnoreUserData(true); + + smootherPrototype->SetVelocityPredictionFactoryManager(sM[0]); + + // correction + // define SchurComplement Factory + // SchurComp gets a RCP to AFact_ which has to be the 2x2 blocked operator + // It stores the resulting SchurComplement operator as "A" generated by the SchurComplementFactory + // Instead of F^{-1} it uses the approximation \hat{F}^{-1} with \hat{F} = diag(F) + RCP AinvFact = Teuchos::rcp(new InverseApproximationFactory()); + AinvFact->SetFactory("A",rAFact); + + RCP SFact = Teuchos::rcp(new SchurComplementFactory()); + SFact->SetParameter("omega", Teuchos::ParameterEntry(Teuchos::as(1.0))); // for Simple, omega is always 1.0 in the SchurComplement + SFact->SetFactory("A",rAFact); + SFact->SetFactory("Ainv", AinvFact); + // create a 2x2 SIMPLE for the prediction eq. RCP smoProtoPredict = Teuchos::rcp( new SimpleSmoother() ); smoProtoPredict->SetParameter("Sweeps", Teuchos::ParameterEntry(1)); smoProtoPredict->SetParameter("Damping factor", Teuchos::ParameterEntry(Teuchos::as(1.0))); smoProtoPredict->SetParameter("UseSIMPLEC", Teuchos::ParameterEntry(false)); - smoProtoPredict->SetFactory("A", sA[0]); + smoProtoPredict->SetFactory("A", SFact); for(int l = 0; l < 2; l++) { Teuchos::RCP ssA = rcp(new SubBlockAFactory()); - ssA->SetFactory("A",sA[0]); + ssA->SetFactory("A",SFact); ssA->SetParameter("block row",Teuchos::ParameterEntry(l)); // local block indices relative to size of blocked operator ssA->SetParameter("block col",Teuchos::ParameterEntry(l)); ssA->SetParameter("Range map: Striding info", Teuchos::ParameterEntry(strInfo)); @@ -1689,31 +2139,7 @@ namespace MueLuTests { else smoProtoPredict->SetSchurCompFactoryManager(ssM); } - sF[0] = rcp( new SmootherFactory(smoProtoPredict) ); - - sM[0] = rcp(new FactoryManager()); - sM[0]->SetFactory("A", sA[0]); - sM[0]->SetFactory("Smoother", sF[0]); - sM[0]->SetIgnoreUserData(true); - - smootherPrototype->SetVelocityPredictionFactoryManager(sM[0]); - - // correction - // define SchurComplement Factory - // SchurComp gets a RCP to AFact_ which has to be the 2x2 blocked operator - // It stores the resulting SchurComplement operator as "A" generated by the SchurComplementFactory - // Instead of F^{-1} it uses the approximation \hat{F}^{-1} with \hat{F} = diag(F) - RCP AinvFact = Teuchos::rcp(new InverseApproximationFactory()); - AinvFact->SetFactory("A",rAFact); - - RCP SFact = Teuchos::rcp(new SchurComplementFactory()); - SFact->SetParameter("omega", Teuchos::ParameterEntry(Teuchos::as(1.0))); // for Simple, omega is always 1.0 in the SchurComplement - SFact->SetFactory("A",rAFact); - SFact->SetFactory("Ainv", AinvFact); - - RCP smoProtoCorrect = rcp(new Ifpack2Smoother(std::string("RELAXATION"), Teuchos::ParameterList(), 0)); - smoProtoCorrect->SetFactory("A", SFact); - sF[1] = rcp( new SmootherFactory(smoProtoCorrect) ); + sF[1] = rcp( new SmootherFactory(smoProtoPredict) ); sM[1] = rcp(new FactoryManager()); sM[1]->SetFactory("A", SFact); @@ -1722,7 +2148,6 @@ namespace MueLuTests { smootherPrototype->SetSchurCompFactoryManager(sM[1]); - RCP smootherFact = rcp( new SmootherFactory(smootherPrototype) ); // main factory manager @@ -1760,10 +2185,10 @@ namespace MueLuTests { // solve system simpleSmoother->Apply(*X, *RHS, true); //zero initial guess - RCP bX = Teuchos::rcp_dynamic_cast(X); TEST_EQUALITY(bX.is_null(), false); RCP XX = bX->Merge(); + Teuchos::ArrayRCP xdata = XX->getData(0); bool bCheck = true; for(size_t i=0; igetLocalLength(); i++) { @@ -1789,6 +2214,8 @@ namespace MueLuTests { // Reset X to 0 X->putScalar((SC) 0.0); + RHS->describe(out, Teuchos::VERB_EXTREME); + RHS->norm2(norms); out << "||RHS|| = " << std::setiosflags(std::ios::fixed) << std::setprecision(10) << norms[0] << std::endl; @@ -1810,9 +2237,9 @@ namespace MueLuTests { }// end useTpetra } - TEUCHOS_UNIT_TEST_TEMPLATE_4_DECL(BlockedSmoother, NestedI2I01II_SIMPLE_Setup_Apply, Scalar, LocalOrdinal, GlobalOrdinal, Node) + TEUCHOS_UNIT_TEST_TEMPLATE_4_DECL(BlockedSmoother, NestedI2I01II_SIMPLE_Setup_Apply2, Scalar, LocalOrdinal, GlobalOrdinal, Node) { - # include +# include MUELU_TESTING_SET_OSTREAM; MUELU_TESTING_LIMIT_SCOPE(Scalar,GlobalOrdinal,Node); // TODO test only Tpetra because of Ifpack2 smoother! @@ -1936,13 +2363,16 @@ namespace MueLuTests { RCP simpleSmoother = level.Get >("PreSmoother", smootherFact.get()); RCP reorderedA = level.Get >("A", rAFact.get()); - RCP reorderedbA = Teuchos::rcp_dynamic_cast(reorderedA); + RCP reorderedbA = Teuchos::rcp_dynamic_cast(reorderedA); TEST_EQUALITY(reorderedbA->Rows(), 2); TEST_EQUALITY(reorderedbA->Cols(), 2); - RCP X = MultiVectorFactory::Build(reorderedA->getDomainMap(),1); - RCP RHS = MultiVectorFactory::Build(reorderedA->getRangeMap(),1); + RCP bX = Teuchos::rcp(new BlockedMultiVector(bop->getBlockedDomainMap(), 1, true)); + RCP bRHS = Teuchos::rcp(new BlockedMultiVector(bop->getBlockedRangeMap(), 1, true)); + + RCP X = bX->Merge(); + RCP RHS = bRHS->Merge(); // apply simple smoother RHS->putScalar((SC) 1.0); @@ -1950,21 +2380,31 @@ namespace MueLuTests { // solve system simpleSmoother->Apply(*X, *RHS, true); //zero initial guess - RCP bX = Teuchos::rcp_dynamic_cast(X); - TEST_EQUALITY(bX.is_null(), false); - RCP XX = bX->Merge(); - Teuchos::ArrayRCP xdata = XX->getData(0); + + Teuchos::ArrayRCP xdata = X->getData(0); bool bCheck = true; - for(size_t i=0; igetLocalLength(); i++) { - if (i< 10) { if(xdata[i] != (SC) (1.0/3.0)) bCheck = false; } - if (i>=10 && i< 15) { if(xdata[i] != (SC) 1.0) bCheck = false; } - if (i>=15 && i< 20) { if(xdata[i] != (SC) 0.5) bCheck = false; } + for(size_t i=0; igetLocalLength(); i++) { + if (i < 5) { if(xdata[i] != (SC) 1.0) bCheck = false; } + if (i>=5 && i< 10) { if(xdata[i] != (SC) 0.5) bCheck = false; } + if (i>=10 && i< 20) { if(xdata[i] != (SC) (1.0/3.0)) bCheck = false; } } TEST_EQUALITY(bCheck, true); // Random X - X->setSeed(846930886); - X->randomize(); + { + Teuchos::RCP brm = Xpetra::blockedReorderFromString("[ 2 [0 1]]"); + RCP test = Teuchos::rcp_dynamic_cast(buildReorderedBlockedMultiVector(brm, bX)); + bX.swap(test); + test = Teuchos::rcp_dynamic_cast(buildReorderedBlockedMultiVector(brm, bRHS)); + bRHS.swap(test); + } + + bRHS->putScalar((SC) 1.0); + bX->setSeed(846930886); + bX->randomize(); + + RHS = bRHS->Merge(); + X = bX->Merge(); typedef typename Teuchos::ScalarTraits::magnitudeType magnitude_type; @@ -1973,7 +2413,7 @@ namespace MueLuTests { X->scale(1/norms[0]); // Compute RHS corresponding to X - reorderedA->apply(*X,*RHS, Teuchos::NO_TRANS,(SC)1.0,(SC)0.0); + bop->apply(*X,*RHS, Teuchos::NO_TRANS,(SC)1.0,(SC)0.0); // Reset X to 0 X->putScalar((SC) 0.0); @@ -1999,7 +2439,6 @@ namespace MueLuTests { }// end useTpetra } - TEUCHOS_UNIT_TEST_TEMPLATE_4_DECL(BlockedSmoother, NestedII20I1I_Thyra_SIMPLE_Setup_Apply, Scalar, LocalOrdinal, GlobalOrdinal, Node) { # include @@ -2833,18 +3272,162 @@ namespace MueLuTests { out << " ||Residual_final|| = " << std::setiosflags(std::ios::fixed) << std::setprecision(20) << residualNorm2[0] << std::endl; out << " ||X_final|| = " << std::setiosflags(std::ios::fixed) << std::setprecision(10) << finalNorms[0] << std::endl; - TEUCHOS_TEST_COMPARE(residualNorm2[0], <, tol, out, success); - TEUCHOS_TEST_COMPARE(finalNorms[0] - Teuchos::ScalarTraits::magnitude(Teuchos::ScalarTraits::one()), <, tol, out, success); + TEUCHOS_TEST_COMPARE(residualNorm2[0], <, tol, out, success); + TEUCHOS_TEST_COMPARE(finalNorms[0] - Teuchos::ScalarTraits::magnitude(Teuchos::ScalarTraits::one()), <, tol, out, success); + + if (comm->getSize() == 1) { + TEST_EQUALITY(residualNorm1[0] != residualNorm2[0], true); + } else { + out << "Pass/Fail is only checked in serial." << std::endl; + } + } // end UseTpetra + } + + TEUCHOS_UNIT_TEST_TEMPLATE_4_DECL(BlockedSmoother, NestedII20I1I_BS_Setup_Apply, Scalar, LocalOrdinal, GlobalOrdinal, Node) + { +# include + MUELU_TESTING_SET_OSTREAM; + MUELU_TESTING_LIMIT_SCOPE(Scalar,GlobalOrdinal,Node); + // TODO test only Tpetra because of Ifpack2 smoother! + MUELU_TEST_ONLY_FOR(Xpetra::UseTpetra) { + + RCP > comm = Parameters::getDefaultComm(); + Xpetra::UnderlyingLib lib = MueLuTests::TestHelpers::Parameters::getLib(); + + Teuchos::RCP bop = TestHelpers::TestFactory::CreateBlockDiagonalExampleMatrix(lib,3, comm); + Teuchos::RCP Aconst = Teuchos::rcp_dynamic_cast(bop); + Teuchos::RCP< Matrix> A = Teuchos::rcp_const_cast(Aconst); + + //I don't use the testApply infrastructure because it has no provision for an initial guess. + Level level; TestHelpers::TestFactory::createSingleLevelHierarchy(level); + level.Set("A", A); + + // Test ReorderBlockAFactory + Teuchos::RCP rAFact = Teuchos::rcp(new ReorderBlockAFactory()); + rAFact->SetFactory("A",MueLu::NoFactory::getRCP()); + rAFact->SetParameter(std::string("Reorder Type"), Teuchos::ParameterEntry(std::string("[ [2 0] 1]"))); + + ////////////////////////////////////////////////////////////////////// + // Smoothers + RCP smootherPrototype = rcp( new BraessSarazinSmoother() ); + smootherPrototype->SetFactory("A",rAFact); + smootherPrototype->SetParameter("Sweeps", Teuchos::ParameterEntry(Teuchos::as(1))); + smootherPrototype->SetParameter("Damping factor", Teuchos::ParameterEntry(Teuchos::as(1.0))); + + std::vector > sF (1, Teuchos::null); + std::vector > sM (1, Teuchos::null); + + // correction + // define SchurComplement Factory + // SchurComp gets a RCP to AFact_ which has to be the 2x2 blocked operator + // It stores the resulting SchurComplement operator as "A" generated by the SchurComplementFactory + // Instead of F^{-1} it uses the approximation \hat{F}^{-1} with \hat{F} = diag(F) + RCP AinvFact = Teuchos::rcp(new InverseApproximationFactory()); + AinvFact->SetFactory("A",rAFact); + + RCP SFact = Teuchos::rcp(new SchurComplementFactory()); + SFact->SetParameter("omega", Teuchos::ParameterEntry(Teuchos::as(1.0))); // for Simple, omega is always 1.0 in the SchurComplement + SFact->SetFactory("A",rAFact); + SFact->SetFactory("Ainv", AinvFact); + + RCP smoProtoCorrect = rcp(new Ifpack2Smoother(std::string("RELAXATION"), Teuchos::ParameterList(), 0)); + smoProtoCorrect->SetFactory("A", SFact); + sF[0] = rcp( new SmootherFactory(smoProtoCorrect) ); + + sM[0] = rcp(new FactoryManager()); + sM[0]->SetFactory("A", SFact); + sM[0]->SetFactory("Smoother", sF[0]); + sM[0]->SetIgnoreUserData(true); + + smootherPrototype->AddFactoryManager(sM[0],0); + + RCP smootherFact = rcp( new SmootherFactory(smootherPrototype) ); + + // main factory manager + FactoryManager M; + M.SetFactory("Smoother", smootherFact); + M.SetFactory("A", rAFact); + + MueLu::SetFactoryManager SFM (Teuchos::rcpFromRef(level), Teuchos::rcpFromRef(M)); + + // request BGS smoother (and all dependencies) on level + level.Request("A", rAFact.get()); + level.Request("Smoother", smootherFact.get()); + level.Request("PreSmoother", smootherFact.get()); + level.Request("PostSmoother", smootherFact.get()); + + //smootherFact->DeclareInput(level); + smootherFact->Build(level); + + level.print(std::cout, Teuchos::VERB_EXTREME); + + RCP simpleSmoother = level.Get >("PreSmoother", smootherFact.get()); + + RCP reorderedA = level.Get >("A", rAFact.get()); + RCP reorderedbA = Teuchos::rcp_dynamic_cast(reorderedA); + + TEST_EQUALITY(reorderedbA->Rows(), 2); + TEST_EQUALITY(reorderedbA->Cols(), 2); + + RCP X = MultiVectorFactory::Build(reorderedA->getDomainMap(),1); + RCP RHS = MultiVectorFactory::Build(reorderedA->getRangeMap(),1); + + // apply simple smoother + RHS->putScalar((SC) 1.0); + X->putScalar((SC) 0.0); + + // solve system + simpleSmoother->Apply(*X, *RHS, true); //zero initial guess + RCP bX = Teuchos::rcp_dynamic_cast(X); + TEST_EQUALITY(bX.is_null(), false); + RCP XX = bX->Merge(); + Teuchos::ArrayRCP xdata = XX->getData(0); + bool bCheck = true; + for(size_t i=0; igetLocalLength(); i++) { + if (i<10 ) { if(xdata[i] != (SC) (1.0/3.0)) bCheck = false; } + if (i>=10 && i<15) { if(xdata[i] != (SC) 1.0) bCheck = false; } + if (i>=15 && i<20) { if(xdata[i] != (SC) 0.5) bCheck = false; } + } + TEST_EQUALITY(bCheck, true); + + // Random X + X->setSeed(846930886); + X->randomize(); + + typedef typename Teuchos::ScalarTraits::magnitudeType magnitude_type; + + // Normalize X + Array norms(1); X->norm2(norms); + X->scale(1/norms[0]); + + // Compute RHS corresponding to X + reorderedbA->apply(*X,*RHS, Teuchos::NO_TRANS,(SC)1.0,(SC)0.0); + + // Reset X to 0 + X->putScalar((SC) 0.0); + + RHS->norm2(norms); + out << "||RHS|| = " << std::setiosflags(std::ios::fixed) << std::setprecision(10) << norms[0] << std::endl; + + out << "solve with zero initial guess" << std::endl; + Teuchos::Array initialNorms(1); X->norm2(initialNorms); + out << " ||X_initial|| = " << std::setiosflags(std::ios::fixed) << std::setprecision(10) << initialNorms[0] << std::endl; + + simpleSmoother->Apply(*X, *RHS, true); //zero initial guess + + Teuchos::Array finalNorms(1); X->norm2(finalNorms); + Teuchos::Array residualNorm1 = Utilities::ResidualNorm(*reorderedbA, *X, *RHS); + out << " ||Residual_final|| = " << std::setiosflags(std::ios::fixed) << std::setprecision(20) << residualNorm1[0] << std::endl; + out << " ||X_final|| = " << std::setiosflags(std::ios::fixed) << std::setprecision(10) << finalNorms[0] << std::endl; + + magnitude_type tol = 50.*Teuchos::ScalarTraits::eps(); - if (comm->getSize() == 1) { - TEST_EQUALITY(residualNorm1[0] != residualNorm2[0], true); - } else { - out << "Pass/Fail is only checked in serial." << std::endl; - } - } // end UseTpetra + TEUCHOS_TEST_COMPARE(residualNorm1[0], <, tol, out, success); + TEUCHOS_TEST_COMPARE(finalNorms[0] - Teuchos::ScalarTraits::magnitude(Teuchos::ScalarTraits::one()), <, tol, out, success); + }// end useTpetra } - TEUCHOS_UNIT_TEST_TEMPLATE_4_DECL(BlockedSmoother, NestedII20I1I_BS_Setup_Apply, Scalar, LocalOrdinal, GlobalOrdinal, Node) + TEUCHOS_UNIT_TEST_TEMPLATE_4_DECL(BlockedSmoother, NestedI2I10II_BS_Setup_Apply, Scalar, LocalOrdinal, GlobalOrdinal, Node) { # include MUELU_TESTING_SET_OSTREAM; @@ -2866,7 +3449,7 @@ namespace MueLuTests { // Test ReorderBlockAFactory Teuchos::RCP rAFact = Teuchos::rcp(new ReorderBlockAFactory()); rAFact->SetFactory("A",MueLu::NoFactory::getRCP()); - rAFact->SetParameter(std::string("Reorder Type"), Teuchos::ParameterEntry(std::string("[ [2 0] 1]"))); + rAFact->SetParameter(std::string("Reorder Type"), Teuchos::ParameterEntry(std::string("[ 2 [ 1 0 ]]"))); ////////////////////////////////////////////////////////////////////// // Smoothers @@ -2891,8 +3474,27 @@ namespace MueLuTests { SFact->SetFactory("A",rAFact); SFact->SetFactory("Ainv", AinvFact); - RCP smoProtoCorrect = rcp(new Ifpack2Smoother(std::string("RELAXATION"), Teuchos::ParameterList(), 0)); + RCP smoProtoCorrect = Teuchos::rcp( new BraessSarazinSmoother() ); + smoProtoCorrect->SetParameter("Sweeps", Teuchos::ParameterEntry(1)); + smoProtoCorrect->SetParameter("Damping factor", Teuchos::ParameterEntry(Teuchos::as(1.0))); smoProtoCorrect->SetFactory("A", SFact); + + std::string strInfo = std::string("{ 1 }"); + Teuchos::RCP ssA = rcp(new SubBlockAFactory()); + ssA->SetFactory("A",SFact); + ssA->SetParameter("block row",Teuchos::ParameterEntry(1)); // local block indices relative to size of blocked operator + ssA->SetParameter("block col",Teuchos::ParameterEntry(1)); + ssA->SetParameter("Range map: Striding info", Teuchos::ParameterEntry(strInfo)); + ssA->SetParameter("Domain map: Striding info", Teuchos::ParameterEntry(strInfo)); + RCP ssP = rcp(new Ifpack2Smoother(std::string("RELAXATION"), Teuchos::ParameterList(), 0)); + ssP->SetFactory("A", ssA); + Teuchos::RCP ssF = Teuchos::rcp(new SmootherFactory(ssP)); + Teuchos::RCP ssM = Teuchos::rcp(new FactoryManager()); + ssM->SetFactory("A", ssA); + ssM->SetFactory("Smoother", ssF); + ssM->SetIgnoreUserData(true); + smoProtoCorrect->AddFactoryManager(ssM,0); + sF[0] = rcp( new SmootherFactory(smoProtoCorrect) ); sM[0] = rcp(new FactoryManager()); @@ -2917,12 +3519,9 @@ namespace MueLuTests { level.Request("PreSmoother", smootherFact.get()); level.Request("PostSmoother", smootherFact.get()); - //smootherFact->DeclareInput(level); smootherFact->Build(level); - level.print(std::cout, Teuchos::VERB_EXTREME); - - RCP simpleSmoother = level.Get >("PreSmoother", smootherFact.get()); + RCP bsSmoother = level.Get >("PreSmoother", smootherFact.get()); RCP reorderedA = level.Get >("A", rAFact.get()); RCP reorderedbA = Teuchos::rcp_dynamic_cast(reorderedA); @@ -2938,7 +3537,8 @@ namespace MueLuTests { X->putScalar((SC) 0.0); // solve system - simpleSmoother->Apply(*X, *RHS, true); //zero initial guess + bsSmoother->Apply(*X, *RHS, true); //zero initial guess + RCP bX = Teuchos::rcp_dynamic_cast(X); TEST_EQUALITY(bX.is_null(), false); RCP XX = bX->Merge(); @@ -2946,8 +3546,8 @@ namespace MueLuTests { bool bCheck = true; for(size_t i=0; igetLocalLength(); i++) { if (i<10 ) { if(xdata[i] != (SC) (1.0/3.0)) bCheck = false; } - if (i>=10 && i<15) { if(xdata[i] != (SC) 1.0) bCheck = false; } - if (i>=15 && i<20) { if(xdata[i] != (SC) 0.5) bCheck = false; } + if (i>=10 && i<15) { if(xdata[i] != (SC) 0.5) bCheck = false; } + if (i>=15 && i<20) { if(xdata[i] != (SC) 1.0) bCheck = false; } } TEST_EQUALITY(bCheck, true); @@ -2974,7 +3574,7 @@ namespace MueLuTests { Teuchos::Array initialNorms(1); X->norm2(initialNorms); out << " ||X_initial|| = " << std::setiosflags(std::ios::fixed) << std::setprecision(10) << initialNorms[0] << std::endl; - simpleSmoother->Apply(*X, *RHS, true); //zero initial guess + bsSmoother->Apply(*X, *RHS, true); //zero initial guess Teuchos::Array finalNorms(1); X->norm2(finalNorms); Teuchos::Array residualNorm1 = Utilities::ResidualNorm(*reorderedbA, *X, *RHS); @@ -2988,7 +3588,7 @@ namespace MueLuTests { }// end useTpetra } - TEUCHOS_UNIT_TEST_TEMPLATE_4_DECL(BlockedSmoother, NestedI2I10II_BS_Setup_Apply, Scalar, LocalOrdinal, GlobalOrdinal, Node) + TEUCHOS_UNIT_TEST_TEMPLATE_4_DECL(BlockedSmoother, NestedI2I10II_BS_Setup_Apply2, Scalar, LocalOrdinal, GlobalOrdinal, Node) { # include MUELU_TESTING_SET_OSTREAM; @@ -3085,13 +3685,16 @@ namespace MueLuTests { RCP bsSmoother = level.Get >("PreSmoother", smootherFact.get()); RCP reorderedA = level.Get >("A", rAFact.get()); - RCP reorderedbA = Teuchos::rcp_dynamic_cast(reorderedA); + RCP reorderedbA = Teuchos::rcp_dynamic_cast(reorderedA); TEST_EQUALITY(reorderedbA->Rows(), 2); TEST_EQUALITY(reorderedbA->Cols(), 2); - RCP X = MultiVectorFactory::Build(reorderedA->getDomainMap(),1); - RCP RHS = MultiVectorFactory::Build(reorderedA->getRangeMap(),1); + RCP bX = Teuchos::rcp(new BlockedMultiVector(bop->getBlockedDomainMap(), 1, true)); + RCP bRHS = Teuchos::rcp(new BlockedMultiVector(bop->getBlockedRangeMap(), 1, true)); + + RCP X = bX->Merge(); + RCP RHS = bRHS->Merge(); // apply simple smoother RHS->putScalar((SC) 1.0); @@ -3100,21 +3703,30 @@ namespace MueLuTests { // solve system bsSmoother->Apply(*X, *RHS, true); //zero initial guess - RCP bX = Teuchos::rcp_dynamic_cast(X); - TEST_EQUALITY(bX.is_null(), false); - RCP XX = bX->Merge(); - Teuchos::ArrayRCP xdata = XX->getData(0); + Teuchos::ArrayRCP xdata = X->getData(0); bool bCheck = true; - for(size_t i=0; igetLocalLength(); i++) { - if (i<10 ) { if(xdata[i] != (SC) (1.0/3.0)) bCheck = false; } - if (i>=10 && i<15) { if(xdata[i] != (SC) 0.5) bCheck = false; } - if (i>=15 && i<20) { if(xdata[i] != (SC) 1.0) bCheck = false; } + for(size_t i=0; igetLocalLength(); i++) { + if (i < 5) { if(xdata[i] != (SC) 1.0) bCheck = false; } + if (i>=5 && i< 10) { if(xdata[i] != (SC) 0.5) bCheck = false; } + if (i>=10 && i< 20) { if(xdata[i] != (SC) (1.0/3.0)) bCheck = false; } } TEST_EQUALITY(bCheck, true); // Random X - X->setSeed(846930886); - X->randomize(); + { + Teuchos::RCP brm = Xpetra::blockedReorderFromString("[ 2 [0 1]]"); + RCP test = Teuchos::rcp_dynamic_cast(buildReorderedBlockedMultiVector(brm, bX)); + bX.swap(test); + test = Teuchos::rcp_dynamic_cast(buildReorderedBlockedMultiVector(brm, bRHS)); + bRHS.swap(test); + } + + bRHS->putScalar((SC) 1.0); + bX->setSeed(846930886); + bX->randomize(); + + RHS = bRHS->Merge(); + X = bX->Merge(); typedef typename Teuchos::ScalarTraits::magnitudeType magnitude_type; @@ -3123,7 +3735,7 @@ namespace MueLuTests { X->scale(1/norms[0]); // Compute RHS corresponding to X - reorderedbA->apply(*X,*RHS, Teuchos::NO_TRANS,(SC)1.0,(SC)0.0); + bop->apply(*X,*RHS, Teuchos::NO_TRANS,(SC)1.0,(SC)0.0); // Reset X to 0 X->putScalar((SC) 0.0); @@ -3149,7 +3761,6 @@ namespace MueLuTests { }// end useTpetra } - TEUCHOS_UNIT_TEST_TEMPLATE_4_DECL(BlockedSmoother, NestedII02I1I_Thyra_BS_Setup_Apply, Scalar, LocalOrdinal, GlobalOrdinal, Node) { # include @@ -3952,10 +4563,217 @@ namespace MueLuTests { level.print(std::cout, Teuchos::VERB_EXTREME); - RCP uzSmoother = level.Get >("PreSmoother", smootherFact.get()); + RCP uzSmoother = level.Get >("PreSmoother", smootherFact.get()); + + RCP X = MultiVectorFactory::Build(A->getDomainMap(),1); + RCP RHS = MultiVectorFactory::Build(A->getRangeMap(),1); + + // Random X + X->setSeed(846930886); + X->randomize(); + + typedef typename Teuchos::ScalarTraits::magnitudeType magnitude_type; + + // Normalize X + Array norms(1); X->norm2(norms); + X->scale(1/norms[0]); + + // Compute RHS corresponding to X + A->apply(*X,*RHS, Teuchos::NO_TRANS,(SC)1.0,(SC)0.0); + + // Reset X to 0 + X->putScalar((SC) 0.0); + + RHS->norm2(norms); + out << "||RHS|| = " << std::setiosflags(std::ios::fixed) << std::setprecision(10) << norms[0] << std::endl; + + out << "solve with zero initial guess" << std::endl; + Teuchos::Array initialNorms(1); X->norm2(initialNorms); + out << " ||X_initial|| = " << std::setiosflags(std::ios::fixed) << std::setprecision(10) << initialNorms[0] << std::endl; + + uzSmoother->Apply(*X, *RHS, true); //zero initial guess + + Teuchos::Array finalNorms(1); X->norm2(finalNorms); + Teuchos::Array residualNorm1 = Utilities::ResidualNorm(*A, *X, *RHS); + out << " ||Residual_final|| = " << std::setiosflags(std::ios::fixed) << std::setprecision(20) << residualNorm1[0] << std::endl; + out << " ||X_final|| = " << std::setiosflags(std::ios::fixed) << std::setprecision(10) << finalNorms[0] << std::endl; + + magnitude_type tol = 50.*Teuchos::ScalarTraits::eps(); + + TEUCHOS_TEST_COMPARE(residualNorm1[0], <, tol, out, success); + TEUCHOS_TEST_COMPARE(finalNorms[0] - Teuchos::ScalarTraits::magnitude(Teuchos::ScalarTraits::one()), <, tol, out, success); + + out << "solve with random initial guess" << std::endl; + X->randomize(); + X->norm2(initialNorms); + out << " ||X_initial|| = " << std::setiosflags(std::ios::fixed) << std::setprecision(10) << initialNorms[0] << std::endl; + + uzSmoother->Apply(*X, *RHS, false); //nonzero initial guess + + X->norm2(finalNorms); + Teuchos::Array residualNorm2 = Utilities::ResidualNorm(*A, *X, *RHS); + out << " ||Residual_final|| = " << std::setiosflags(std::ios::fixed) << std::setprecision(20) << residualNorm2[0] << std::endl; + out << " ||X_final|| = " << std::setiosflags(std::ios::fixed) << std::setprecision(10) << finalNorms[0] << std::endl; + + TEUCHOS_TEST_COMPARE(residualNorm2[0], <, tol, out, success); + TEUCHOS_TEST_COMPARE(finalNorms[0] - Teuchos::ScalarTraits::magnitude(Teuchos::ScalarTraits::one()), <, tol, out, success); + + if (comm->getSize() == 1) { + TEST_EQUALITY(residualNorm1[0] != residualNorm2[0], true); + } else { + out << "Pass/Fail is only checked in serial." << std::endl; + } + } // end UseTpetra + } + + TEUCHOS_UNIT_TEST_TEMPLATE_4_DECL(BlockedSmoother, NestedI2I01II_Uzawa_Setup_Apply, Scalar, LocalOrdinal, GlobalOrdinal, Node) + { +# include + MUELU_TESTING_SET_OSTREAM; + MUELU_TESTING_LIMIT_SCOPE(Scalar,GlobalOrdinal,Node); + // TODO test only Tpetra because of Ifpack2 smoother! + MUELU_TEST_ONLY_FOR(Xpetra::UseTpetra) { + + RCP > comm = Parameters::getDefaultComm(); + Xpetra::UnderlyingLib lib = MueLuTests::TestHelpers::Parameters::getLib(); + + Teuchos::RCP bop = TestHelpers::TestFactory::CreateBlockDiagonalExampleMatrix(lib,3, comm); + Teuchos::RCP Aconst = Teuchos::rcp_dynamic_cast(bop); + Teuchos::RCP< Matrix> A = Teuchos::rcp_const_cast(Aconst); + + //I don't use the testApply infrastructure because it has no provision for an initial guess. + Level level; TestHelpers::TestFactory::createSingleLevelHierarchy(level); + level.Set("A", A); + + // Test ReorderBlockAFactory + Teuchos::RCP rAFact = Teuchos::rcp(new ReorderBlockAFactory()); + rAFact->SetFactory("A",MueLu::NoFactory::getRCP()); + rAFact->SetParameter(std::string("Reorder Type"), Teuchos::ParameterEntry(std::string("[ 2 [0 1]]"))); + + ////////////////////////////////////////////////////////////////////// + // Smoothers + RCP smootherPrototype = rcp( new UzawaSmoother() ); + smootherPrototype->SetFactory("A",rAFact); + smootherPrototype->SetParameter("Sweeps", Teuchos::ParameterEntry(Teuchos::as(1))); + smootherPrototype->SetParameter("Damping factor", Teuchos::ParameterEntry(Teuchos::as(1.0))); + + std::vector > sA (1, Teuchos::null); + std::vector > sF (2, Teuchos::null); + std::vector > sM (2, Teuchos::null); + + // prediction + std::string strInfo = std::string("{ 1 }"); + sA[0] = rcp(new SubBlockAFactory()); + sA[0]->SetFactory("A",rAFact); + sA[0]->SetParameter("block row",Teuchos::ParameterEntry(0)); + sA[0]->SetParameter("block col",Teuchos::ParameterEntry(0)); + sA[0]->SetParameter("Range map: Striding info", Teuchos::ParameterEntry(strInfo)); + sA[0]->SetParameter("Domain map: Striding info", Teuchos::ParameterEntry(strInfo)); + + RCP smoProtoCorrect = rcp(new Ifpack2Smoother(std::string("RELAXATION"), Teuchos::ParameterList(), 0)); + smoProtoCorrect->SetFactory("A", sA[0]); + sF[0] = rcp( new SmootherFactory(smoProtoCorrect) ); + + sM[0] = rcp(new FactoryManager()); + sM[0]->SetFactory("A", sA[0]); + sM[0]->SetFactory("Smoother", sF[0]); + sM[0]->SetIgnoreUserData(true); + + smootherPrototype->AddFactoryManager(sM[0],0); + + // correction + // define SchurComplement Factory + // SchurComp gets a RCP to AFact_ which has to be the 2x2 blocked operator + // It stores the resulting SchurComplement operator as "A" generated by the SchurComplementFactory + // Instead of F^{-1} it uses the approximation \hat{F}^{-1} with \hat{F} = diag(F) + RCP AinvFact = Teuchos::rcp(new InverseApproximationFactory()); + AinvFact->SetFactory("A",rAFact); + + RCP SFact = Teuchos::rcp(new SchurComplementFactory()); + SFact->SetParameter("omega", Teuchos::ParameterEntry(Teuchos::as(1.0))); // for Simple, omega is always 1.0 in the SchurComplement + SFact->SetFactory("A",rAFact); + SFact->SetFactory("Ainv", AinvFact); + + // create a 2x2 SIMPLE for the prediction eq. + RCP smoProtoPredict = Teuchos::rcp( new UzawaSmoother() ); + smoProtoPredict->SetParameter("Sweeps", Teuchos::ParameterEntry(1)); + smoProtoPredict->SetParameter("Damping factor", Teuchos::ParameterEntry(Teuchos::as(1.0))); + smoProtoPredict->SetFactory("A", SFact); + + for(int l = 0; l < 2; l++) { + Teuchos::RCP ssA = rcp(new SubBlockAFactory()); + ssA->SetFactory("A",SFact); + ssA->SetParameter("block row",Teuchos::ParameterEntry(l)); // local block indices relative to size of blocked operator + ssA->SetParameter("block col",Teuchos::ParameterEntry(l)); + ssA->SetParameter("Range map: Striding info", Teuchos::ParameterEntry(strInfo)); + ssA->SetParameter("Domain map: Striding info", Teuchos::ParameterEntry(strInfo)); + RCP ssP = rcp(new Ifpack2Smoother(std::string("RELAXATION"), Teuchos::ParameterList(), 0)); + ssP->SetFactory("A", ssA); + Teuchos::RCP ssF = Teuchos::rcp(new SmootherFactory(ssP)); + Teuchos::RCP ssM = Teuchos::rcp(new FactoryManager()); + ssM->SetFactory("A", ssA); + ssM->SetFactory("Smoother", ssF); + ssM->SetIgnoreUserData(true); + smoProtoPredict->AddFactoryManager(ssM,l); + } + + sF[1] = rcp( new SmootherFactory(smoProtoPredict) ); + + sM[1] = rcp(new FactoryManager()); + sM[1]->SetFactory("A", SFact); + sM[1]->SetFactory("Smoother", sF[1]); + sM[1]->SetIgnoreUserData(true); + + smootherPrototype->AddFactoryManager(sM[1],1); + + RCP smootherFact = rcp( new SmootherFactory(smootherPrototype) ); + + // main factory manager + FactoryManager M; + M.SetFactory("Smoother", smootherFact); + M.SetFactory("A", rAFact); + + MueLu::SetFactoryManager SFM (Teuchos::rcpFromRef(level), Teuchos::rcpFromRef(M)); + + // request BGS smoother (and all dependencies) on level + level.Request("A", rAFact.get()); + level.Request("Smoother", smootherFact.get()); + level.Request("PreSmoother", smootherFact.get()); + level.Request("PostSmoother", smootherFact.get()); + + //smootherFact->DeclareInput(level); + smootherFact->Build(level); + + level.print(std::cout, Teuchos::VERB_EXTREME); + + RCP simpleSmoother = level.Get >("PreSmoother", smootherFact.get()); - RCP X = MultiVectorFactory::Build(A->getDomainMap(),1); - RCP RHS = MultiVectorFactory::Build(A->getRangeMap(),1); + RCP reorderedA = level.Get >("A", rAFact.get()); + RCP reorderedbA = Teuchos::rcp_dynamic_cast(reorderedA); + + TEST_EQUALITY(reorderedbA->Rows(), 2); + TEST_EQUALITY(reorderedbA->Cols(), 2); + + RCP X = MultiVectorFactory::Build(reorderedA->getDomainMap(),1); + RCP RHS = MultiVectorFactory::Build(reorderedA->getRangeMap(),1); + + // apply simple smoother + RHS->putScalar((SC) 1.0); + X->putScalar((SC) 0.0); + + // solve system + simpleSmoother->Apply(*X, *RHS, true); //zero initial guess + RCP bX = Teuchos::rcp_dynamic_cast(X); + TEST_EQUALITY(bX.is_null(), false); + RCP XX = bX->Merge(); + Teuchos::ArrayRCP xdata = XX->getData(0); + bool bCheck = true; + for(size_t i=0; igetLocalLength(); i++) { + if (i<10) { if(xdata[i] != (SC) (1.0/3.0)) bCheck = false; } + if (i>=10 && i< 15) { if(xdata[i] != (SC) 1.0) bCheck = false; } + if (i>=15 && i< 20) { if(xdata[i] != (SC) 0.5) bCheck = false; } + } + TEST_EQUALITY(bCheck, true); // Random X X->setSeed(846930886); @@ -3968,7 +4786,7 @@ namespace MueLuTests { X->scale(1/norms[0]); // Compute RHS corresponding to X - A->apply(*X,*RHS, Teuchos::NO_TRANS,(SC)1.0,(SC)0.0); + reorderedA->apply(*X,*RHS, Teuchos::NO_TRANS,(SC)1.0,(SC)0.0); // Reset X to 0 X->putScalar((SC) 0.0); @@ -3980,10 +4798,10 @@ namespace MueLuTests { Teuchos::Array initialNorms(1); X->norm2(initialNorms); out << " ||X_initial|| = " << std::setiosflags(std::ios::fixed) << std::setprecision(10) << initialNorms[0] << std::endl; - uzSmoother->Apply(*X, *RHS, true); //zero initial guess + simpleSmoother->Apply(*X, *RHS, true); //zero initial guess Teuchos::Array finalNorms(1); X->norm2(finalNorms); - Teuchos::Array residualNorm1 = Utilities::ResidualNorm(*A, *X, *RHS); + Teuchos::Array residualNorm1 = Utilities::ResidualNorm(*reorderedA, *X, *RHS); out << " ||Residual_final|| = " << std::setiosflags(std::ios::fixed) << std::setprecision(20) << residualNorm1[0] << std::endl; out << " ||X_final|| = " << std::setiosflags(std::ios::fixed) << std::setprecision(10) << finalNorms[0] << std::endl; @@ -3991,31 +4809,10 @@ namespace MueLuTests { TEUCHOS_TEST_COMPARE(residualNorm1[0], <, tol, out, success); TEUCHOS_TEST_COMPARE(finalNorms[0] - Teuchos::ScalarTraits::magnitude(Teuchos::ScalarTraits::one()), <, tol, out, success); - - out << "solve with random initial guess" << std::endl; - X->randomize(); - X->norm2(initialNorms); - out << " ||X_initial|| = " << std::setiosflags(std::ios::fixed) << std::setprecision(10) << initialNorms[0] << std::endl; - - uzSmoother->Apply(*X, *RHS, false); //nonzero initial guess - - X->norm2(finalNorms); - Teuchos::Array residualNorm2 = Utilities::ResidualNorm(*A, *X, *RHS); - out << " ||Residual_final|| = " << std::setiosflags(std::ios::fixed) << std::setprecision(20) << residualNorm2[0] << std::endl; - out << " ||X_final|| = " << std::setiosflags(std::ios::fixed) << std::setprecision(10) << finalNorms[0] << std::endl; - - TEUCHOS_TEST_COMPARE(residualNorm2[0], <, tol, out, success); - TEUCHOS_TEST_COMPARE(finalNorms[0] - Teuchos::ScalarTraits::magnitude(Teuchos::ScalarTraits::one()), <, tol, out, success); - - if (comm->getSize() == 1) { - TEST_EQUALITY(residualNorm1[0] != residualNorm2[0], true); - } else { - out << "Pass/Fail is only checked in serial." << std::endl; - } - } // end UseTpetra + }// end useTpetra } - TEUCHOS_UNIT_TEST_TEMPLATE_4_DECL(BlockedSmoother, NestedI2I01II_Uzawa_Setup_Apply, Scalar, LocalOrdinal, GlobalOrdinal, Node) + TEUCHOS_UNIT_TEST_TEMPLATE_4_DECL(BlockedSmoother, NestedI2I01II_Uzawa_Setup_Apply2, Scalar, LocalOrdinal, GlobalOrdinal, Node) { # include MUELU_TESTING_SET_OSTREAM; @@ -4077,7 +4874,7 @@ namespace MueLuTests { // Instead of F^{-1} it uses the approximation \hat{F}^{-1} with \hat{F} = diag(F) RCP AinvFact = Teuchos::rcp(new InverseApproximationFactory()); AinvFact->SetFactory("A",rAFact); - + RCP SFact = Teuchos::rcp(new SchurComplementFactory()); SFact->SetParameter("omega", Teuchos::ParameterEntry(Teuchos::as(1.0))); // for Simple, omega is always 1.0 in the SchurComplement SFact->SetFactory("A",rAFact); @@ -4138,13 +4935,16 @@ namespace MueLuTests { RCP simpleSmoother = level.Get >("PreSmoother", smootherFact.get()); RCP reorderedA = level.Get >("A", rAFact.get()); - RCP reorderedbA = Teuchos::rcp_dynamic_cast(reorderedA); + RCP reorderedbA = Teuchos::rcp_dynamic_cast(reorderedA); TEST_EQUALITY(reorderedbA->Rows(), 2); TEST_EQUALITY(reorderedbA->Cols(), 2); - RCP X = MultiVectorFactory::Build(reorderedA->getDomainMap(),1); - RCP RHS = MultiVectorFactory::Build(reorderedA->getRangeMap(),1); + RCP bX = Teuchos::rcp(new BlockedMultiVector(bop->getBlockedDomainMap(), 1, true)); + RCP bRHS = Teuchos::rcp(new BlockedMultiVector(bop->getBlockedRangeMap(), 1, true)); + + RCP X = bX->Merge(); + RCP RHS = bRHS->Merge(); // apply simple smoother RHS->putScalar((SC) 1.0); @@ -4152,21 +4952,32 @@ namespace MueLuTests { // solve system simpleSmoother->Apply(*X, *RHS, true); //zero initial guess - RCP bX = Teuchos::rcp_dynamic_cast(X); - TEST_EQUALITY(bX.is_null(), false); - RCP XX = bX->Merge(); - Teuchos::ArrayRCP xdata = XX->getData(0); + + Teuchos::ArrayRCP xdata = X->getData(0); bool bCheck = true; - for(size_t i=0; igetLocalLength(); i++) { - if (i<10) { if(xdata[i] != (SC) (1.0/3.0)) bCheck = false; } - if (i>=10 && i< 15) { if(xdata[i] != (SC) 1.0) bCheck = false; } - if (i>=15 && i< 20) { if(xdata[i] != (SC) 0.5) bCheck = false; } + X->describe(out, Teuchos::VERB_EXTREME); + for(size_t i=0; igetLocalLength(); i++) { + if (i < 5) { if(xdata[i] != (SC) 1.0) bCheck = false; } + if (i>=5 && i< 10) { if(xdata[i] != (SC) 0.5) bCheck = false; } + if (i>=10 && i< 20) { if(xdata[i] != (SC) (1.0/3.0)) bCheck = false; } } TEST_EQUALITY(bCheck, true); // Random X - X->setSeed(846930886); - X->randomize(); + { + Teuchos::RCP brm = Xpetra::blockedReorderFromString("[ 2 [0 1]]"); + RCP test = Teuchos::rcp_dynamic_cast(buildReorderedBlockedMultiVector(brm, bX)); + bX.swap(test); + test = Teuchos::rcp_dynamic_cast(buildReorderedBlockedMultiVector(brm, bRHS)); + bRHS.swap(test); + } + + bRHS->putScalar((SC) 1.0); + bX->setSeed(846930886); + bX->randomize(); + + RHS = bRHS->Merge(); + X = bX->Merge(); typedef typename Teuchos::ScalarTraits::magnitudeType magnitude_type; @@ -4175,7 +4986,7 @@ namespace MueLuTests { X->scale(1/norms[0]); // Compute RHS corresponding to X - reorderedA->apply(*X,*RHS, Teuchos::NO_TRANS,(SC)1.0,(SC)0.0); + bop->apply(*X,*RHS, Teuchos::NO_TRANS,(SC)1.0,(SC)0.0); // Reset X to 0 X->putScalar((SC) 0.0); @@ -5031,6 +5842,205 @@ namespace MueLuTests { }// end useTpetra } + TEUCHOS_UNIT_TEST_TEMPLATE_4_DECL(BlockedSmoother, NestedI2I01II_Indef_Setup_Apply2, Scalar, LocalOrdinal, GlobalOrdinal, Node) + { +# include + MUELU_TESTING_SET_OSTREAM; + MUELU_TESTING_LIMIT_SCOPE(Scalar,GlobalOrdinal,Node); + // TODO test only Tpetra because of Ifpack2 smoother! + MUELU_TEST_ONLY_FOR(Xpetra::UseTpetra) { + + RCP > comm = Parameters::getDefaultComm(); + Xpetra::UnderlyingLib lib = MueLuTests::TestHelpers::Parameters::getLib(); + + Teuchos::RCP bop = TestHelpers::TestFactory::CreateBlockDiagonalExampleMatrix(lib,3, comm); + Teuchos::RCP Aconst = Teuchos::rcp_dynamic_cast(bop); + Teuchos::RCP< Matrix> A = Teuchos::rcp_const_cast(Aconst); + + //I don't use the testApply infrastructure because it has no provision for an initial guess. + Level level; TestHelpers::TestFactory::createSingleLevelHierarchy(level); + level.Set("A", A); + + // Test ReorderBlockAFactory + Teuchos::RCP rAFact = Teuchos::rcp(new ReorderBlockAFactory()); + rAFact->SetFactory("A",MueLu::NoFactory::getRCP()); + rAFact->SetParameter(std::string("Reorder Type"), Teuchos::ParameterEntry(std::string("[ 2 [0 1]]"))); + + ////////////////////////////////////////////////////////////////////// + // Smoothers + RCP smootherPrototype = rcp( new IndefBlockedDiagonalSmoother() ); + smootherPrototype->SetFactory("A",rAFact); + smootherPrototype->SetParameter("Sweeps", Teuchos::ParameterEntry(Teuchos::as(1))); + smootherPrototype->SetParameter("Damping factor", Teuchos::ParameterEntry(Teuchos::as(1.0))); + + std::vector > sA (1, Teuchos::null); + std::vector > sF (2, Teuchos::null); + std::vector > sM (2, Teuchos::null); + + // prediction + std::string strInfo = std::string("{ 1 }"); + sA[0] = rcp(new SubBlockAFactory()); + sA[0]->SetFactory("A",rAFact); + sA[0]->SetParameter("block row",Teuchos::ParameterEntry(0)); + sA[0]->SetParameter("block col",Teuchos::ParameterEntry(0)); + sA[0]->SetParameter("Range map: Striding info", Teuchos::ParameterEntry(strInfo)); + sA[0]->SetParameter("Domain map: Striding info", Teuchos::ParameterEntry(strInfo)); + + RCP smoProtoCorrect = rcp(new Ifpack2Smoother(std::string("RELAXATION"), Teuchos::ParameterList(), 0)); + smoProtoCorrect->SetFactory("A", sA[0]); + sF[0] = rcp( new SmootherFactory(smoProtoCorrect) ); + + sM[0] = rcp(new FactoryManager()); + sM[0]->SetFactory("A", sA[0]); + sM[0]->SetFactory("Smoother", sF[0]); + sM[0]->SetIgnoreUserData(true); + + smootherPrototype->AddFactoryManager(sM[0],0); + + // correction + // define SchurComplement Factory + // SchurComp gets a RCP to AFact_ which has to be the 2x2 blocked operator + // It stores the resulting SchurComplement operator as "A" generated by the SchurComplementFactory + // Instead of F^{-1} it uses the approximation \hat{F}^{-1} with \hat{F} = diag(F) + RCP AinvFact = Teuchos::rcp(new InverseApproximationFactory()); + AinvFact->SetFactory("A",rAFact); + + RCP SFact = Teuchos::rcp(new SchurComplementFactory()); + SFact->SetParameter("omega", Teuchos::ParameterEntry(Teuchos::as(1.0))); // for Simple, omega is always 1.0 in the SchurComplement + SFact->SetFactory("A",rAFact); + SFact->SetFactory("Ainv", AinvFact); + + // create a 2x2 block smoother for the prediction eq. + RCP smoProtoPredict = Teuchos::rcp( new IndefBlockedDiagonalSmoother() ); + smoProtoPredict->SetParameter("Sweeps", Teuchos::ParameterEntry(1)); + smoProtoPredict->SetParameter("Damping factor", Teuchos::ParameterEntry(Teuchos::as(1.0))); + smoProtoPredict->SetFactory("A", SFact); + + for(int l = 0; l < 2; l++) { + Teuchos::RCP ssA = rcp(new SubBlockAFactory()); + ssA->SetFactory("A",SFact); + ssA->SetParameter("block row",Teuchos::ParameterEntry(l)); // local block indices relative to size of blocked operator + ssA->SetParameter("block col",Teuchos::ParameterEntry(l)); + ssA->SetParameter("Range map: Striding info", Teuchos::ParameterEntry(strInfo)); + ssA->SetParameter("Domain map: Striding info", Teuchos::ParameterEntry(strInfo)); + RCP ssP = rcp(new Ifpack2Smoother(std::string("RELAXATION"), Teuchos::ParameterList(), 0)); + ssP->SetFactory("A", ssA); + Teuchos::RCP ssF = Teuchos::rcp(new SmootherFactory(ssP)); + Teuchos::RCP ssM = Teuchos::rcp(new FactoryManager()); + ssM->SetFactory("A", ssA); + ssM->SetFactory("Smoother", ssF); + ssM->SetIgnoreUserData(true); + smoProtoPredict->AddFactoryManager(ssM,l); + } + + sF[1] = rcp( new SmootherFactory(smoProtoPredict) ); + + sM[1] = rcp(new FactoryManager()); + sM[1]->SetFactory("A", SFact); + sM[1]->SetFactory("Smoother", sF[1]); + sM[1]->SetIgnoreUserData(true); + + smootherPrototype->AddFactoryManager(sM[1],1); + + RCP smootherFact = rcp( new SmootherFactory(smootherPrototype) ); + + // main factory manager + FactoryManager M; + M.SetFactory("Smoother", smootherFact); + M.SetFactory("A", rAFact); + + MueLu::SetFactoryManager SFM (Teuchos::rcpFromRef(level), Teuchos::rcpFromRef(M)); + + // request block smoother (and all dependencies) on level + level.Request("A", rAFact.get()); + level.Request("Smoother", smootherFact.get()); + level.Request("PreSmoother", smootherFact.get()); + level.Request("PostSmoother", smootherFact.get()); + + //smootherFact->DeclareInput(level); + smootherFact->Build(level); + + level.print(std::cout, Teuchos::VERB_EXTREME); + + RCP inSmoother = level.Get >("PreSmoother", smootherFact.get()); + + RCP reorderedA = level.Get >("A", rAFact.get()); + RCP reorderedbA = Teuchos::rcp_dynamic_cast(reorderedA); + + TEST_EQUALITY(reorderedbA->Rows(), 2); + TEST_EQUALITY(reorderedbA->Cols(), 2); + + RCP bX = Teuchos::rcp(new BlockedMultiVector(bop->getBlockedDomainMap(), 1, true)); + RCP bRHS = Teuchos::rcp(new BlockedMultiVector(bop->getBlockedRangeMap(), 1, true)); + + RCP X = bX->Merge(); + RCP RHS = bRHS->Merge(); + + // apply simple smoother + RHS->putScalar((SC) 1.0); + X->putScalar((SC) 0.0); + + // solve system + inSmoother->Apply(*X, *RHS, true); //zero initial guess + + Teuchos::ArrayRCP xdata = X->getData(0); + bool bCheck = true; + for(size_t i=0; igetLocalLength(); i++) { + if (i < 5) { if(xdata[i] != (SC) 1.0) bCheck = false; } + if (i>=5 && i< 10) { if(xdata[i] != (SC) 0.5) bCheck = false; } + if (i>=10 && i< 20) { if(xdata[i] != (SC) (1.0/3.0)) bCheck = false; } + } + TEST_EQUALITY(bCheck, true); + + // Random X + { + Teuchos::RCP brm = Xpetra::blockedReorderFromString("[ 2 [0 1]]"); + RCP test = Teuchos::rcp_dynamic_cast(buildReorderedBlockedMultiVector(brm, bX)); + bX.swap(test); + test = Teuchos::rcp_dynamic_cast(buildReorderedBlockedMultiVector(brm, bRHS)); + bRHS.swap(test); + } + + bRHS->putScalar((SC) 1.0); + bX->setSeed(846930886); + bX->randomize(); + + RHS = bRHS->Merge(); + X = bX->Merge(); + + typedef typename Teuchos::ScalarTraits::magnitudeType magnitude_type; + + // Normalize X + Array norms(1); X->norm2(norms); + X->scale(1/norms[0]); + + // Compute RHS corresponding to X + bop->apply(*X,*RHS, Teuchos::NO_TRANS,(SC)1.0,(SC)0.0); + + // Reset X to 0 + X->putScalar((SC) 0.0); + + RHS->norm2(norms); + out << "||RHS|| = " << std::setiosflags(std::ios::fixed) << std::setprecision(10) << norms[0] << std::endl; + + out << "solve with zero initial guess" << std::endl; + Teuchos::Array initialNorms(1); X->norm2(initialNorms); + out << " ||X_initial|| = " << std::setiosflags(std::ios::fixed) << std::setprecision(10) << initialNorms[0] << std::endl; + + inSmoother->Apply(*X, *RHS, true); //zero initial guess + + Teuchos::Array finalNorms(1); X->norm2(finalNorms); + Teuchos::Array residualNorm1 = Utilities::ResidualNorm(*reorderedA, *X, *RHS); + out << " ||Residual_final|| = " << std::setiosflags(std::ios::fixed) << std::setprecision(20) << residualNorm1[0] << std::endl; + out << " ||X_final|| = " << std::setiosflags(std::ios::fixed) << std::setprecision(10) << finalNorms[0] << std::endl; + + magnitude_type tol = 50.*Teuchos::ScalarTraits::eps(); + + TEUCHOS_TEST_COMPARE(residualNorm1[0], <, tol, out, success); + TEUCHOS_TEST_COMPARE(finalNorms[0] - Teuchos::ScalarTraits::magnitude(Teuchos::ScalarTraits::one()), <, tol, out, success); + }// end useTpetra + } + TEUCHOS_UNIT_TEST_TEMPLATE_4_DECL(BlockedSmoother, NestedI2I01II_Thyra_Indef_Setup_Apply, Scalar, LocalOrdinal, GlobalOrdinal, Node) { # include @@ -5621,12 +6631,14 @@ namespace MueLuTests { TEUCHOS_UNIT_TEST_TEMPLATE_4_INSTANT(BlockedSmoother,BGS_Setup_Apply,SC,LO,GO,NO) \ TEUCHOS_UNIT_TEST_TEMPLATE_4_INSTANT(BlockedSmoother,Reordered_BGS_Setup_Apply,SC,LO,GO,NO) \ TEUCHOS_UNIT_TEST_TEMPLATE_4_INSTANT(BlockedSmoother,NestedII30II12II_BGS_Setup_Apply,SC,LO,GO,NO) \ + TEUCHOS_UNIT_TEST_TEMPLATE_4_INSTANT(BlockedSmoother,NestedII30II12II_BGS_Setup_Apply2,SC,LO,GO,NO) \ TEUCHOS_UNIT_TEST_TEMPLATE_4_INSTANT(BlockedSmoother,Thyra_BGS_Setup_Apply,SC,LO,GO,NO) \ TEUCHOS_UNIT_TEST_TEMPLATE_4_INSTANT(BlockedSmoother,Thyra_Nested_BGS_Setup_Apply,SC,LO,GO,NO) \ TEUCHOS_UNIT_TEST_TEMPLATE_4_INSTANT(BlockedSmoother,Thyra_Nested_BGS_Setup_Apply2,SC,LO,GO,NO) \ TEUCHOS_UNIT_TEST_TEMPLATE_4_INSTANT(BlockedSmoother,SIMPLE_Setup_Apply,SC,LO,GO,NO) \ TEUCHOS_UNIT_TEST_TEMPLATE_4_INSTANT(BlockedSmoother,NestedII20I1I_SIMPLE_Setup_Apply,SC,LO,GO,NO) \ TEUCHOS_UNIT_TEST_TEMPLATE_4_INSTANT(BlockedSmoother,NestedI2I01II_SIMPLE_Setup_Apply,SC,LO,GO,NO) \ + TEUCHOS_UNIT_TEST_TEMPLATE_4_INSTANT(BlockedSmoother,NestedI2I01II_SIMPLE_Setup_Apply2,SC,LO,GO,NO) \ TEUCHOS_UNIT_TEST_TEMPLATE_4_INSTANT(BlockedSmoother,NestedII20I1I_Thyra_SIMPLE_Setup_Apply,SC,LO,GO,NO) \ TEUCHOS_UNIT_TEST_TEMPLATE_4_INSTANT(BlockedSmoother,NestedI2I01II_Thyra_SIMPLE_Setup_Apply,SC,LO,GO,NO) \ TEUCHOS_UNIT_TEST_TEMPLATE_4_INSTANT(BlockedSmoother,NestedII01I2I_Thyra_SIMPLE_Setup_Apply2,SC,LO,GO,NO) \ @@ -5634,6 +6646,7 @@ namespace MueLuTests { TEUCHOS_UNIT_TEST_TEMPLATE_4_INSTANT(BlockedSmoother,BS_Setup_Apply,SC,LO,GO,NO) \ TEUCHOS_UNIT_TEST_TEMPLATE_4_INSTANT(BlockedSmoother,NestedII20I1I_BS_Setup_Apply,SC,LO,GO,NO) \ TEUCHOS_UNIT_TEST_TEMPLATE_4_INSTANT(BlockedSmoother,NestedI2I10II_BS_Setup_Apply,SC,LO,GO,NO) \ + TEUCHOS_UNIT_TEST_TEMPLATE_4_INSTANT(BlockedSmoother,NestedI2I10II_BS_Setup_Apply2,SC,LO,GO,NO) \ TEUCHOS_UNIT_TEST_TEMPLATE_4_INSTANT(BlockedSmoother,NestedII02I1I_Thyra_BS_Setup_Apply,SC,LO,GO,NO) \ TEUCHOS_UNIT_TEST_TEMPLATE_4_INSTANT(BlockedSmoother,NestedI2I01II_Thyra_BS_Setup_Apply,SC,LO,GO,NO) \ TEUCHOS_UNIT_TEST_TEMPLATE_4_INSTANT(BlockedSmoother,NestedI2I10II_Thyra_BS_Setup_Apply,SC,LO,GO,NO) \ @@ -5641,11 +6654,13 @@ namespace MueLuTests { TEUCHOS_UNIT_TEST_TEMPLATE_4_INSTANT(BlockedSmoother,NestedI0I21II_Thyra_BS_Setup_Apply3,SC,LO,GO,NO) \ TEUCHOS_UNIT_TEST_TEMPLATE_4_INSTANT(BlockedSmoother,Uzawa_Setup_Apply,SC,LO,GO,NO) \ TEUCHOS_UNIT_TEST_TEMPLATE_4_INSTANT(BlockedSmoother,NestedI2I01II_Uzawa_Setup_Apply,SC,LO,GO,NO) \ + TEUCHOS_UNIT_TEST_TEMPLATE_4_INSTANT(BlockedSmoother,NestedI2I01II_Uzawa_Setup_Apply2,SC,LO,GO,NO) \ TEUCHOS_UNIT_TEST_TEMPLATE_4_INSTANT(BlockedSmoother,NestedI2I01II_Thyra_Uzawa_Setup_Apply,SC,LO,GO,NO) \ TEUCHOS_UNIT_TEST_TEMPLATE_4_INSTANT(BlockedSmoother,NestedII01I2I_Thyra_Uzawa_Setup_Apply2,SC,LO,GO,NO) \ TEUCHOS_UNIT_TEST_TEMPLATE_4_INSTANT(BlockedSmoother,NestedI0I21II_Thyra_Uzawa_Setup_Apply3,SC,LO,GO,NO) \ TEUCHOS_UNIT_TEST_TEMPLATE_4_INSTANT(BlockedSmoother,Indef_Setup_Apply,SC,LO,GO,NO) \ TEUCHOS_UNIT_TEST_TEMPLATE_4_INSTANT(BlockedSmoother,NestedI2I01II_Indef_Setup_Apply,SC,LO,GO,NO) \ + TEUCHOS_UNIT_TEST_TEMPLATE_4_INSTANT(BlockedSmoother,NestedI2I01II_Indef_Setup_Apply2,SC,LO,GO,NO) \ TEUCHOS_UNIT_TEST_TEMPLATE_4_INSTANT(BlockedSmoother,NestedI2I01II_Thyra_Indef_Setup_Apply,SC,LO,GO,NO) \ TEUCHOS_UNIT_TEST_TEMPLATE_4_INSTANT(BlockedSmoother,NestedII01I2I_Thyra_Indef_Setup_Apply2,SC,LO,GO,NO) \ TEUCHOS_UNIT_TEST_TEMPLATE_4_INSTANT(BlockedSmoother,NestedI0I21II_Thyra_Indef_Setup_Apply3,SC,LO,GO,NO) \ From 35b2dc5c2198886f572dc61fc79bd8c26a5cffff Mon Sep 17 00:00:00 2001 From: maxfirmbach Date: Tue, 4 Jul 2023 20:44:42 +0200 Subject: [PATCH 12/16] Update naming scheme of temporary variables --- .../MueLu_BlockedGaussSeidelSmoother_def.hpp | 24 +++++++++---------- .../MueLu_BlockedJacobiSmoother_def.hpp | 24 +++++++++---------- .../MueLu_BraessSarazinSmoother_def.hpp | 24 +++++++++---------- ...MueLu_IndefBlockedDiagonalSmoother_def.hpp | 24 +++++++++---------- .../MueLu_SimpleSmoother_def.hpp | 24 +++++++++---------- .../MueLu_UzawaSmoother_def.hpp | 24 +++++++++---------- 6 files changed, 72 insertions(+), 72 deletions(-) diff --git a/packages/muelu/src/Smoothers/BlockedSmoothers/MueLu_BlockedGaussSeidelSmoother_def.hpp b/packages/muelu/src/Smoothers/BlockedSmoothers/MueLu_BlockedGaussSeidelSmoother_def.hpp index fea0bee9d9ff..c0b689f1d72e 100644 --- a/packages/muelu/src/Smoothers/BlockedSmoothers/MueLu_BlockedGaussSeidelSmoother_def.hpp +++ b/packages/muelu/src/Smoothers/BlockedSmoothers/MueLu_BlockedGaussSeidelSmoother_def.hpp @@ -218,30 +218,30 @@ namespace MueLu { // check type of X vector if (bX.is_null() == true) { - RCP test = Teuchos::rcp(new BlockedMultiVector(rbA->getBlockedCrsMatrix()->getBlockedDomainMap(), rcpX)); - rcpX.swap(test); + RCP vectorWithBlockedMap = Teuchos::rcp(new BlockedMultiVector(rbA->getBlockedCrsMatrix()->getBlockedDomainMap(), rcpX)); + rcpX.swap(vectorWithBlockedMap); bCopyResultX = true; bReorderX = true; } // check type of B vector if (bB.is_null() == true) { - RCP test = Teuchos::rcp(new BlockedMultiVector(rbA->getBlockedCrsMatrix()->getBlockedRangeMap(), rcpB)); - rcpB.swap(test); + RCP vectorWithBlockedMap = Teuchos::rcp(new BlockedMultiVector(rbA->getBlockedCrsMatrix()->getBlockedRangeMap(), rcpB)); + rcpB.swap(vectorWithBlockedMap); bReorderB = true; } } else { // A is a BlockedCrsMatrix and uses a plain blocked map if (bX.is_null() == true) { - RCP test = Teuchos::rcp(new BlockedMultiVector(bA->getBlockedDomainMap(), rcpX)); - rcpX.swap(test); + RCP vectorWithBlockedMap = Teuchos::rcp(new BlockedMultiVector(bA->getBlockedDomainMap(), rcpX)); + rcpX.swap(vectorWithBlockedMap); bCopyResultX = true; } if(bB.is_null() == true) { - RCP test = Teuchos::rcp(new BlockedMultiVector(bA->getBlockedRangeMap(),rcpB)); - rcpB.swap(test); + RCP vectorWithBlockedMap = Teuchos::rcp(new BlockedMultiVector(bA->getBlockedRangeMap(),rcpB)); + rcpB.swap(vectorWithBlockedMap); } } @@ -257,14 +257,14 @@ namespace MueLu { // check type of X vector if(bX->getBlockedMap()->getNumMaps() != bA->getDomainMapExtractor()->NumMaps() || bReorderX) { // X is a blocked multi vector but incompatible to the reordered blocked operator A - Teuchos::RCP test = buildReorderedBlockedMultiVector(brm, bX); - rcpX.swap(test); + Teuchos::RCP reorderedBlockedVector = buildReorderedBlockedMultiVector(brm, bX); + rcpX.swap(reorderedBlockedVector); } if(bB->getBlockedMap()->getNumMaps() != bA->getRangeMapExtractor()->NumMaps() || bReorderB) { // B is a blocked multi vector but incompatible to the reordered blocked operator A - Teuchos::RCP test = buildReorderedBlockedMultiVector(brm, bB); - rcpB.swap(test); + Teuchos::RCP reorderedBlockedVector = buildReorderedBlockedMultiVector(brm, bB); + rcpB.swap(reorderedBlockedVector); } } diff --git a/packages/muelu/src/Smoothers/BlockedSmoothers/MueLu_BlockedJacobiSmoother_def.hpp b/packages/muelu/src/Smoothers/BlockedSmoothers/MueLu_BlockedJacobiSmoother_def.hpp index b2d712bda5f6..e1e1a908987e 100644 --- a/packages/muelu/src/Smoothers/BlockedSmoothers/MueLu_BlockedJacobiSmoother_def.hpp +++ b/packages/muelu/src/Smoothers/BlockedSmoothers/MueLu_BlockedJacobiSmoother_def.hpp @@ -220,30 +220,30 @@ namespace MueLu { // check type of X vector if (bX.is_null() == true) { - RCP test = Teuchos::rcp(new BlockedMultiVector(rbA->getBlockedCrsMatrix()->getBlockedDomainMap(), rcpX)); - rcpX.swap(test); + RCP vectorWithBlockedMap = Teuchos::rcp(new BlockedMultiVector(rbA->getBlockedCrsMatrix()->getBlockedDomainMap(), rcpX)); + rcpX.swap(vectorWithBlockedMap); bCopyResultX = true; bReorderX = true; } // check type of B vector if (bB.is_null() == true) { - RCP test = Teuchos::rcp(new BlockedMultiVector(rbA->getBlockedCrsMatrix()->getBlockedRangeMap(), rcpB)); - rcpB.swap(test); + RCP vectorWithBlockedMap = Teuchos::rcp(new BlockedMultiVector(rbA->getBlockedCrsMatrix()->getBlockedRangeMap(), rcpB)); + rcpB.swap(vectorWithBlockedMap); bReorderB = true; } } else { // A is a BlockedCrsMatrix and uses a plain blocked map if (bX.is_null() == true) { - RCP test = Teuchos::rcp(new BlockedMultiVector(bA->getBlockedDomainMap(), rcpX)); - rcpX.swap(test); + RCP vectorWithBlockedMap = Teuchos::rcp(new BlockedMultiVector(bA->getBlockedDomainMap(), rcpX)); + rcpX.swap(vectorWithBlockedMap); bCopyResultX = true; } if(bB.is_null() == true) { - RCP test = Teuchos::rcp(new BlockedMultiVector(bA->getBlockedRangeMap(),rcpB)); - rcpB.swap(test); + RCP vectorWithBlockedMap = Teuchos::rcp(new BlockedMultiVector(bA->getBlockedRangeMap(),rcpB)); + rcpB.swap(vectorWithBlockedMap); } } @@ -259,14 +259,14 @@ namespace MueLu { // check type of X vector if(bX->getBlockedMap()->getNumMaps() != bA->getDomainMapExtractor()->NumMaps() || bReorderX) { // X is a blocked multi vector but incompatible to the reordered blocked operator A - Teuchos::RCP test = buildReorderedBlockedMultiVector(brm, bX); - rcpX.swap(test); + Teuchos::RCP reorderedBlockedVector = buildReorderedBlockedMultiVector(brm, bX); + rcpX.swap(reorderedBlockedVector); } if(bB->getBlockedMap()->getNumMaps() != bA->getRangeMapExtractor()->NumMaps() || bReorderB) { // B is a blocked multi vector but incompatible to the reordered blocked operator A - Teuchos::RCP test = buildReorderedBlockedMultiVector(brm, bB); - rcpB.swap(test); + Teuchos::RCP reorderedBlockedVector = buildReorderedBlockedMultiVector(brm, bB); + rcpB.swap(reorderedBlockedVector); } } diff --git a/packages/muelu/src/Smoothers/BlockedSmoothers/MueLu_BraessSarazinSmoother_def.hpp b/packages/muelu/src/Smoothers/BlockedSmoothers/MueLu_BraessSarazinSmoother_def.hpp index d33ff1a04b1a..4e8b2c8ab29d 100644 --- a/packages/muelu/src/Smoothers/BlockedSmoothers/MueLu_BraessSarazinSmoother_def.hpp +++ b/packages/muelu/src/Smoothers/BlockedSmoothers/MueLu_BraessSarazinSmoother_def.hpp @@ -194,30 +194,30 @@ namespace MueLu { // check type of X vector if (bX.is_null() == true) { - RCP test = Teuchos::rcp(new BlockedMultiVector(rbA->getBlockedCrsMatrix()->getBlockedDomainMap(), rcpX)); - rcpX.swap(test); + RCP vectorWithBlockedMap = Teuchos::rcp(new BlockedMultiVector(rbA->getBlockedCrsMatrix()->getBlockedDomainMap(), rcpX)); + rcpX.swap(vectorWithBlockedMap); bCopyResultX = true; bReorderX = true; } // check type of B vector if (bB.is_null() == true) { - RCP test = Teuchos::rcp(new BlockedMultiVector(rbA->getBlockedCrsMatrix()->getBlockedRangeMap(), rcpB)); - rcpB.swap(test); + RCP vectorWithBlockedMap = Teuchos::rcp(new BlockedMultiVector(rbA->getBlockedCrsMatrix()->getBlockedRangeMap(), rcpB)); + rcpB.swap(vectorWithBlockedMap); bReorderB = true; } } else { // A is a BlockedCrsMatrix and uses a plain blocked map if (bX.is_null() == true) { - RCP test = Teuchos::rcp(new BlockedMultiVector(bA->getBlockedDomainMap(), rcpX)); - rcpX.swap(test); + RCP vectorWithBlockedMap = Teuchos::rcp(new BlockedMultiVector(bA->getBlockedDomainMap(), rcpX)); + rcpX.swap(vectorWithBlockedMap); bCopyResultX = true; } if(bB.is_null() == true) { - RCP test = Teuchos::rcp(new BlockedMultiVector(bA->getBlockedRangeMap(),rcpB)); - rcpB.swap(test); + RCP vectorWithBlockedMap = Teuchos::rcp(new BlockedMultiVector(bA->getBlockedRangeMap(),rcpB)); + rcpB.swap(vectorWithBlockedMap); } } @@ -233,14 +233,14 @@ namespace MueLu { // check type of X vector if(bX->getBlockedMap()->getNumMaps() != bA->getDomainMapExtractor()->NumMaps() || bReorderX) { // X is a blocked multi vector but incompatible to the reordered blocked operator A - Teuchos::RCP test = buildReorderedBlockedMultiVector(brm, bX); - rcpX.swap(test); + Teuchos::RCP reorderedBlockedVector = buildReorderedBlockedMultiVector(brm, bX); + rcpX.swap(reorderedBlockedVector); } if(bB->getBlockedMap()->getNumMaps() != bA->getRangeMapExtractor()->NumMaps() || bReorderB) { // B is a blocked multi vector but incompatible to the reordered blocked operator A - Teuchos::RCP test = buildReorderedBlockedMultiVector(brm, bB); - rcpB.swap(test); + Teuchos::RCP reorderedBlockedVector = buildReorderedBlockedMultiVector(brm, bB); + rcpB.swap(reorderedBlockedVector); } } diff --git a/packages/muelu/src/Smoothers/BlockedSmoothers/MueLu_IndefBlockedDiagonalSmoother_def.hpp b/packages/muelu/src/Smoothers/BlockedSmoothers/MueLu_IndefBlockedDiagonalSmoother_def.hpp index 59803565f338..37c5eb238afa 100644 --- a/packages/muelu/src/Smoothers/BlockedSmoothers/MueLu_IndefBlockedDiagonalSmoother_def.hpp +++ b/packages/muelu/src/Smoothers/BlockedSmoothers/MueLu_IndefBlockedDiagonalSmoother_def.hpp @@ -244,30 +244,30 @@ namespace MueLu { // check type of X vector if (bX.is_null() == true) { - RCP test = Teuchos::rcp(new BlockedMultiVector(rbA->getBlockedCrsMatrix()->getBlockedDomainMap(), rcpX)); - rcpX.swap(test); + RCP vectorWithBlockedMap = Teuchos::rcp(new BlockedMultiVector(rbA->getBlockedCrsMatrix()->getBlockedDomainMap(), rcpX)); + rcpX.swap(vectorWithBlockedMap); bCopyResultX = true; bReorderX = true; } // check type of B vector if (bB.is_null() == true) { - RCP test = Teuchos::rcp(new BlockedMultiVector(rbA->getBlockedCrsMatrix()->getBlockedRangeMap(), rcpB)); - rcpB.swap(test); + RCP vectorWithBlockedMap = Teuchos::rcp(new BlockedMultiVector(rbA->getBlockedCrsMatrix()->getBlockedRangeMap(), rcpB)); + rcpB.swap(vectorWithBlockedMap); bReorderB = true; } } else { // A is a BlockedCrsMatrix and uses a plain blocked map if (bX.is_null() == true) { - RCP test = Teuchos::rcp(new BlockedMultiVector(bA->getBlockedDomainMap(), rcpX)); - rcpX.swap(test); + RCP vectorWithBlockedMap = Teuchos::rcp(new BlockedMultiVector(bA->getBlockedDomainMap(), rcpX)); + rcpX.swap(vectorWithBlockedMap); bCopyResultX = true; } if(bB.is_null() == true) { - RCP test = Teuchos::rcp(new BlockedMultiVector(bA->getBlockedRangeMap(),rcpB)); - rcpB.swap(test); + RCP vectorWithBlockedMap = Teuchos::rcp(new BlockedMultiVector(bA->getBlockedRangeMap(),rcpB)); + rcpB.swap(vectorWithBlockedMap); } } @@ -283,14 +283,14 @@ namespace MueLu { // check type of X vector if(bX->getBlockedMap()->getNumMaps() != bA->getDomainMapExtractor()->NumMaps() || bReorderX) { // X is a blocked multi vector but incompatible to the reordered blocked operator A - Teuchos::RCP test = buildReorderedBlockedMultiVector(brm, bX); - rcpX.swap(test); + Teuchos::RCP reorderedBlockedVector = buildReorderedBlockedMultiVector(brm, bX); + rcpX.swap(reorderedBlockedVector); } if(bB->getBlockedMap()->getNumMaps() != bA->getRangeMapExtractor()->NumMaps() || bReorderB) { // B is a blocked multi vector but incompatible to the reordered blocked operator A - Teuchos::RCP test = buildReorderedBlockedMultiVector(brm, bB); - rcpB.swap(test); + Teuchos::RCP reorderedBlockedVector = buildReorderedBlockedMultiVector(brm, bB); + rcpB.swap(reorderedBlockedVector); } } diff --git a/packages/muelu/src/Smoothers/BlockedSmoothers/MueLu_SimpleSmoother_def.hpp b/packages/muelu/src/Smoothers/BlockedSmoothers/MueLu_SimpleSmoother_def.hpp index 63d472853594..6e892869de0f 100644 --- a/packages/muelu/src/Smoothers/BlockedSmoothers/MueLu_SimpleSmoother_def.hpp +++ b/packages/muelu/src/Smoothers/BlockedSmoothers/MueLu_SimpleSmoother_def.hpp @@ -275,30 +275,30 @@ namespace MueLu { // check type of X vector if (bX.is_null() == true) { - RCP test = Teuchos::rcp(new BlockedMultiVector(rbA->getBlockedCrsMatrix()->getBlockedDomainMap(), rcpX)); - rcpX.swap(test); + RCP vectorWithBlockedMap = Teuchos::rcp(new BlockedMultiVector(rbA->getBlockedCrsMatrix()->getBlockedDomainMap(), rcpX)); + rcpX.swap(vectorWithBlockedMap); bCopyResultX = true; bReorderX = true; } // check type of B vector if (bB.is_null() == true) { - RCP test = Teuchos::rcp(new BlockedMultiVector(rbA->getBlockedCrsMatrix()->getBlockedRangeMap(), rcpB)); - rcpB.swap(test); + RCP vectorWithBlockedMap = Teuchos::rcp(new BlockedMultiVector(rbA->getBlockedCrsMatrix()->getBlockedRangeMap(), rcpB)); + rcpB.swap(vectorWithBlockedMap); bReorderB = true; } } else { // A is a BlockedCrsMatrix and uses a plain blocked map if (bX.is_null() == true) { - RCP test = Teuchos::rcp(new BlockedMultiVector(bA->getBlockedDomainMap(), rcpX)); - rcpX.swap(test); + RCP vectorWithBlockedMap = Teuchos::rcp(new BlockedMultiVector(bA->getBlockedDomainMap(), rcpX)); + rcpX.swap(vectorWithBlockedMap); bCopyResultX = true; } if(bB.is_null() == true) { - RCP test = Teuchos::rcp(new BlockedMultiVector(bA->getBlockedRangeMap(),rcpB)); - rcpB.swap(test); + RCP vectorWithBlockedMap = Teuchos::rcp(new BlockedMultiVector(bA->getBlockedRangeMap(),rcpB)); + rcpB.swap(vectorWithBlockedMap); } } @@ -314,14 +314,14 @@ namespace MueLu { // check type of X vector if(bX->getBlockedMap()->getNumMaps() != bA->getDomainMapExtractor()->NumMaps() || bReorderX) { // X is a blocked multi vector but incompatible to the reordered blocked operator A - Teuchos::RCP test = buildReorderedBlockedMultiVector(brm, bX); - rcpX.swap(test); + Teuchos::RCP reorderedBlockedVector = buildReorderedBlockedMultiVector(brm, bX); + rcpX.swap(reorderedBlockedVector); } if(bB->getBlockedMap()->getNumMaps() != bA->getRangeMapExtractor()->NumMaps() || bReorderB) { // B is a blocked multi vector but incompatible to the reordered blocked operator A - Teuchos::RCP test = buildReorderedBlockedMultiVector(brm, bB); - rcpB.swap(test); + Teuchos::RCP reorderedBlockedVector = buildReorderedBlockedMultiVector(brm, bB); + rcpB.swap(reorderedBlockedVector); } } diff --git a/packages/muelu/src/Smoothers/BlockedSmoothers/MueLu_UzawaSmoother_def.hpp b/packages/muelu/src/Smoothers/BlockedSmoothers/MueLu_UzawaSmoother_def.hpp index 3aa6ebb05aef..a9d1d98bca89 100644 --- a/packages/muelu/src/Smoothers/BlockedSmoothers/MueLu_UzawaSmoother_def.hpp +++ b/packages/muelu/src/Smoothers/BlockedSmoothers/MueLu_UzawaSmoother_def.hpp @@ -208,30 +208,30 @@ namespace MueLu { // check type of X vector if (bX.is_null() == true) { - RCP test = Teuchos::rcp(new BlockedMultiVector(rbA->getBlockedCrsMatrix()->getBlockedDomainMap(), rcpX)); - rcpX.swap(test); + RCP vectorWithBlockedMap = Teuchos::rcp(new BlockedMultiVector(rbA->getBlockedCrsMatrix()->getBlockedDomainMap(), rcpX)); + rcpX.swap(vectorWithBlockedMap); bCopyResultX = true; bReorderX = true; } // check type of B vector if (bB.is_null() == true) { - RCP test = Teuchos::rcp(new BlockedMultiVector(rbA->getBlockedCrsMatrix()->getBlockedRangeMap(), rcpB)); - rcpB.swap(test); + RCP vectorWithBlockedMap = Teuchos::rcp(new BlockedMultiVector(rbA->getBlockedCrsMatrix()->getBlockedRangeMap(), rcpB)); + rcpB.swap(vectorWithBlockedMap); bReorderB = true; } } else { // A is a BlockedCrsMatrix and uses a plain blocked map if (bX.is_null() == true) { - RCP test = Teuchos::rcp(new BlockedMultiVector(bA->getBlockedDomainMap(), rcpX)); - rcpX.swap(test); + RCP vectorWithBlockedMap = Teuchos::rcp(new BlockedMultiVector(bA->getBlockedDomainMap(), rcpX)); + rcpX.swap(vectorWithBlockedMap); bCopyResultX = true; } if(bB.is_null() == true) { - RCP test = Teuchos::rcp(new BlockedMultiVector(bA->getBlockedRangeMap(),rcpB)); - rcpB.swap(test); + RCP vectorWithBlockedMap = Teuchos::rcp(new BlockedMultiVector(bA->getBlockedRangeMap(),rcpB)); + rcpB.swap(vectorWithBlockedMap); } } @@ -247,14 +247,14 @@ namespace MueLu { // check type of X vector if(bX->getBlockedMap()->getNumMaps() != bA->getDomainMapExtractor()->NumMaps() || bReorderX) { // X is a blocked multi vector but incompatible to the reordered blocked operator A - Teuchos::RCP test = buildReorderedBlockedMultiVector(brm, bX); - rcpX.swap(test); + Teuchos::RCP reorderedBlockedVector = buildReorderedBlockedMultiVector(brm, bX); + rcpX.swap(reorderedBlockedVector); } if(bB->getBlockedMap()->getNumMaps() != bA->getRangeMapExtractor()->NumMaps() || bReorderB) { // B is a blocked multi vector but incompatible to the reordered blocked operator A - Teuchos::RCP test = buildReorderedBlockedMultiVector(brm, bB); - rcpB.swap(test); + Teuchos::RCP reorderedBlockedVector = buildReorderedBlockedMultiVector(brm, bB); + rcpB.swap(reorderedBlockedVector); } } From f9e3e095c0d9554714b6b5e20ab21aef043d102d Mon Sep 17 00:00:00 2001 From: Greg Sjaardema Date: Fri, 7 Jul 2023 09:02:54 -0600 Subject: [PATCH 13/16] EXODUS: Fix for mpi-io with fixed-length path limit --- packages/seacas/libraries/exodus/src/ex_create_par.c | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/packages/seacas/libraries/exodus/src/ex_create_par.c b/packages/seacas/libraries/exodus/src/ex_create_par.c index 02f568238d58..bb1f27cea32c 100644 --- a/packages/seacas/libraries/exodus/src/ex_create_par.c +++ b/packages/seacas/libraries/exodus/src/ex_create_par.c @@ -192,10 +192,10 @@ int ex_create_par_int(const char *path, int cmode, int *comp_ws, int *io_ws, MPI } #endif - /* There is an issue on some versions of mpi that limit the length of the path to <256 characters - * Check for that here and use `path` if `canon_path` is >=256 characters... + /* There is an issue on some versions of mpi that limit the length of the path to <250 characters + * Check for that here and use `path` if `canon_path` is >=250 characters... */ - if (strlen(canon_path) >= 256) { + if (strlen(canon_path) >= 250) { status = nc_create_par(path, nc_mode, comm, info, &exoid); } else { From e96c24c5825a2d0ee111b13e344f466470ce1786 Mon Sep 17 00:00:00 2001 From: Greg Sjaardema Date: Fri, 7 Jul 2023 09:09:33 -0600 Subject: [PATCH 14/16] Update ex_open_par.c --- packages/seacas/libraries/exodus/src/ex_open_par.c | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/packages/seacas/libraries/exodus/src/ex_open_par.c b/packages/seacas/libraries/exodus/src/ex_open_par.c index fa13ce31b837..aabad6c82b16 100644 --- a/packages/seacas/libraries/exodus/src/ex_open_par.c +++ b/packages/seacas/libraries/exodus/src/ex_open_par.c @@ -180,10 +180,10 @@ int ex_open_par_int(const char *path, int mode, int *comp_ws, int *io_ws, float else { nc_mode = (NC_NOWRITE | NC_SHARE | NC_MPIIO); } - /* There is an issue on some versions of mpi that limit the length of the path to <256 characters - * Check for that here and use `path` if `canon_path` is >=256 characters... + /* There is an issue on some versions of mpi that limit the length of the path to <250 characters + * Check for that here and use `path` if `canon_path` is >=250 characters... */ - if (strlen(canon_path) >= 256) { + if (strlen(canon_path) >= 250) { status = nc_open_par(path, nc_mode, comm, info, &exoid); } else { From a79aa7f7a4deff82265e24f42f7d7f2d194b5345 Mon Sep 17 00:00:00 2001 From: Chris Siefert Date: Fri, 7 Jul 2023 09:13:27 -0600 Subject: [PATCH 15/16] ML: Adding sorting support for Maxwell1 Part of the continuing ML/MueLu harmonization project --- packages/ml/src/Utils/ml_MultiLevelPreconditioner.cpp | 4 +++- 1 file changed, 3 insertions(+), 1 deletion(-) diff --git a/packages/ml/src/Utils/ml_MultiLevelPreconditioner.cpp b/packages/ml/src/Utils/ml_MultiLevelPreconditioner.cpp index b3a32d3da202..169214c05d78 100644 --- a/packages/ml/src/Utils/ml_MultiLevelPreconditioner.cpp +++ b/packages/ml/src/Utils/ml_MultiLevelPreconditioner.cpp @@ -2242,8 +2242,10 @@ ComputePreconditioner(const bool CheckPreconditioner) if(AMGSolver_ != ML_CLASSICAL_FAMILY) ML_CHK_ERR(SetupCoordinates()); - if (List_.get("RAP: sort columns",0)) // + if (List_.get("RAP: sort columns",0)) { + if (ml_nodes_) ml_nodes_->sortColumnsAfterRAP = 1; ml_->sortColumnsAfterRAP = 1; + } // ========================================================================// // Setting Repartitioning // // ========================================================================// From 9a3b26b65ac6cf3f82c72ed2a88775ebe60a1250 Mon Sep 17 00:00:00 2001 From: Chris Siefert Date: Fri, 7 Jul 2023 13:10:03 -0600 Subject: [PATCH 16/16] MueLu: Fixing minor bug --- packages/muelu/test/scaling/MatvecKernelDriver.cpp | 14 ++++++++------ 1 file changed, 8 insertions(+), 6 deletions(-) diff --git a/packages/muelu/test/scaling/MatvecKernelDriver.cpp b/packages/muelu/test/scaling/MatvecKernelDriver.cpp index 580dd47ecb65..fc499d42e8cf 100644 --- a/packages/muelu/test/scaling/MatvecKernelDriver.cpp +++ b/packages/muelu/test/scaling/MatvecKernelDriver.cpp @@ -166,12 +166,14 @@ void report_performance_models(const Teuchos::RCP & A, int nrepeat if(A->hasCrsGraph()) { auto importer = A->getCrsGraph()->getImporter(); - size_t recv_size = importer->getRemoteLIDs().size() * sizeof(SC); - size_t send_size = importer->getExportLIDs().size() * sizeof(SC); - int local_log_max = ceil(log(std::max(send_size,recv_size)) / log(2))+1; - int global_log_max=local_log_max; - Teuchos::reduceAll(*comm,Teuchos::REDUCE_MAX,1,&local_log_max,&global_log_max); - PM.halopong_make_table(nrepeat,global_log_max, importer); + if(!importer.is_null()) { + size_t recv_size = importer->getRemoteLIDs().size() * sizeof(SC); + size_t send_size = importer->getExportLIDs().size() * sizeof(SC); + int local_log_max = ceil(log(std::max(send_size,recv_size)) / log(2))+1; + int global_log_max=local_log_max; + Teuchos::reduceAll(*comm,Teuchos::REDUCE_MAX,1,&local_log_max,&global_log_max); + PM.halopong_make_table(nrepeat,global_log_max, importer); + } } if(verbose && rank == 0) {