47#ifndef MUELU_MATRIXANALYSISFACTORY_DEF_HPP_
48#define MUELU_MATRIXANALYSISFACTORY_DEF_HPP_
52#include <Xpetra_Map.hpp>
53#include <Xpetra_StridedMap.hpp>
54#include <Xpetra_Vector.hpp>
55#include <Xpetra_VectorFactory.hpp>
56#include <Xpetra_Matrix.hpp>
59#include "MueLu_Utilities.hpp"
63template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
67template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
70template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
72 RCP<ParameterList> validParamList = rcp(
new ParameterList());
74 validParamList->set< RCP<const FactoryBase> >(
"A", Teuchos::null,
"Generating factory of the matrix A to be permuted.");
76 return validParamList;
79template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
81 Input(fineLevel,
"A");
84template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
88 Teuchos::RCP<Matrix> A = Get< Teuchos::RCP<Matrix> > (fineLevel,
"A");
89 Teuchos::RCP<const Teuchos::Comm<int> > comm = A->getRangeMap()->getComm();
95 GetOStream(
Runtime0) <<
"~~~~~~~~ GENERAL INFORMATION ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ (Level " << fineLevel.
GetLevelID() <<
")" << std::endl;
96 GetOStream(
Runtime0) <<
"A is a " << A->getRangeMap()->getGlobalNumElements() <<
" x " << A->getDomainMap()->getGlobalNumElements() <<
" matrix." << std::endl;
98 if (A->IsView(
"stridedMaps") && Teuchos::rcp_dynamic_cast<const StridedMap>(A->getRowMap(
"stridedMaps")) != Teuchos::null) {
99 Teuchos::RCP<const StridedMap> strRowMap = Teuchos::rcp_dynamic_cast<const StridedMap>(A->getRowMap(
"stridedMaps"));
102 std::vector<size_t> stridingInfo = strRowMap->getStridingData();
103 LO dofsPerNode = Teuchos::as<LocalOrdinal>(stridingInfo[blockid]);
104 GetOStream(
Runtime0) <<
"Strided row: " << dofsPerNode <<
" dofs per node. Strided block id = " << blockid << std::endl;
106 GetOStream(
Runtime0) <<
"Row: " << strRowMap->getFixedBlockSize() <<
" dofs per node" << std::endl;
110 if (A->IsView(
"stridedMaps") && Teuchos::rcp_dynamic_cast<const StridedMap>(A->getColMap(
"stridedMaps")) != Teuchos::null) {
111 Teuchos::RCP<const StridedMap> strDomMap = Teuchos::rcp_dynamic_cast<const StridedMap>(A->getColMap(
"stridedMaps"));
114 std::vector<size_t> stridingInfo = strDomMap->getStridingData();
115 LO dofsPerNode = Teuchos::as<LocalOrdinal>(stridingInfo[blockid]);
116 GetOStream(
Runtime0) <<
"Strided column: " << dofsPerNode <<
" dofs per node. Strided block id = " << blockid << std::endl;
118 GetOStream(
Runtime0) <<
"Column: " << strDomMap->getFixedBlockSize() <<
" dofs per node" << std::endl;
123 GetOStream(
Runtime0) <<
"A is distributed over " << comm->getSize() <<
" processors" << std::endl;
126 std::vector<LO> lelePerProc(comm->getSize(),0);
127 std::vector<LO> gelePerProc(comm->getSize(),0);
128 lelePerProc[comm->getRank()] = A->getLocalNumEntries ();
129 Teuchos::reduceAll(*comm,Teuchos::REDUCE_MAX,comm->getSize(),&lelePerProc[0],&gelePerProc[0]);
130 if(comm->getRank() == 0) {
131 for(
int i=0; i<comm->getSize(); i++) {
132 if(gelePerProc[i] == 0) {
133 GetOStream(
Runtime0) <<
"Proc " << i <<
" is empty." << std::endl;
141 GetOStream(
Runtime0) <<
"~~~~~~~~ MATRIX DIAGONAL ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ (Level " << fineLevel.
GetLevelID() <<
")" << std::endl;
142 Teuchos::RCP<Vector> diagAVec = VectorFactory::Build(A->getRowMap(),
true);
143 A->getLocalDiagCopy(*diagAVec);
144 Teuchos::ArrayRCP< const Scalar > diagAVecData = diagAVec->getData(0);
152 for(
size_t i = 0; i<diagAVec->getMap()->getLocalNumElements(); ++i) {
153 if(diagAVecData[i] == Teuchos::ScalarTraits<Scalar>::zero()) {
155 }
else if(Teuchos::ScalarTraits<Scalar>::magnitude(diagAVecData[i]) < 1e-6) {
156 lnearlyzerosOnDiagonal++;
157 }
else if(Teuchos::ScalarTraits<Scalar>::isnaninf(diagAVecData[i])) {
164 MueLu_sumAll(comm, lnearlyzerosOnDiagonal, gnearlyzerosOnDiagonal);
167 if(gzerosOnDiagonal > 0) GetOStream(
Runtime0) <<
"Found " << gzerosOnDiagonal <<
" zeros on diagonal of A" << std::endl;
168 if(gnearlyzerosOnDiagonal > 0) GetOStream(
Runtime0) <<
"Found " << gnearlyzerosOnDiagonal <<
" entries with absolute value < 1.0e-6 on diagonal of A" << std::endl;
169 if(gnansOnDiagonal > 0) GetOStream(
Runtime0) <<
"Found " << gnansOnDiagonal <<
" entries with NAN or INF on diagonal of A" << std::endl;
174 GetOStream(
Runtime0) <<
"~~~~~~~~ DIAGONAL DOMINANCE ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ (Level " << fineLevel.
GetLevelID() <<
")" << std::endl;
181 double worstRatio = 99999999.9;
182 for (
size_t row = 0; row < A->getRowMap()->getLocalNumElements(); row++) {
185 size_t nnz = A->getNumEntriesInLocalRow(row);
188 Teuchos::ArrayView<const LocalOrdinal> indices;
189 Teuchos::ArrayView<const Scalar> vals;
190 A->getLocalRowView(row, indices, vals);
192 TEUCHOS_TEST_FOR_EXCEPTION(Teuchos::as<size_t>(indices.size()) != nnz,
Exceptions::RuntimeError,
"MueLu::MatrixAnalysisFactory::Build: number of nonzeros not equal to number of indices? Error.");
196 double normdiag = 0.0;
197 for (
size_t j = 0; j < Teuchos::as<size_t>(indices.size()); j++) {
198 GO gcol = A->getColMap()->getGlobalElement(indices[j]);
199 if (gcol==grow) normdiag = Teuchos::as<double>(Teuchos::ScalarTraits<Scalar>::magnitude(vals[j]));
200 else norm1 += Teuchos::as<double>(Teuchos::ScalarTraits<Scalar>::magnitude(vals[j]));
203 if (normdiag >= norm1) lnumWeakDiagDomRows++;
204 else if (normdiag > norm1) lnumStrictDiagDomRows++;
207 if (normdiag / norm1 < worstRatio) worstRatio = normdiag / norm1;
211 MueLu_sumAll(comm, lnumWeakDiagDomRows, gnumWeakDiagDomRows);
212 MueLu_sumAll(comm, lnumStrictDiagDomRows, gnumStrictDiagDomRows);
214 GetOStream(
Runtime0) <<
"A has " << gnumWeakDiagDomRows <<
"/" << A->getRangeMap()->getGlobalNumElements() <<
" weakly diagonal dominant rows. (" << Teuchos::as<Scalar>(gnumWeakDiagDomRows/A->getRangeMap()->getGlobalNumElements()*100) <<
"%)" << std::endl;
215 GetOStream(
Runtime0) <<
"A has " << gnumStrictDiagDomRows <<
"/" << A->getRangeMap()->getGlobalNumElements() <<
" strictly diagonal dominant rows. (" << Teuchos::as<Scalar>(gnumStrictDiagDomRows/A->getRangeMap()->getGlobalNumElements()*100) <<
"%)" << std::endl;
218 Teuchos::reduceAll(*comm,Teuchos::REDUCE_MIN,comm->getSize(),&worstRatio,&gworstRatio);
219 GetOStream(
Runtime0) <<
"The minimum of the ratio of diagonal element/ sum of off-diagonal elements is " << gworstRatio <<
". Values about 1.0 are ok." << std::endl;
223 SC one = Teuchos::ScalarTraits<Scalar>::one(), zero = Teuchos::ScalarTraits<Scalar>::zero();
227 GetOStream(
Runtime0) <<
"~~~~~~~~ Av with one vector ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ (Level " << fineLevel.
GetLevelID() <<
")" << std::endl;
228 Teuchos::RCP<Vector> ones = VectorFactory::Build(A->getDomainMap(), 1);
229 ones->putScalar(one);
231 Teuchos::RCP<Vector> res1 = VectorFactory::Build(A->getRangeMap(),
false);
232 A->apply(*ones, *res1, Teuchos::NO_TRANS, one, zero);
233 checkVectorEntries(res1,std::string(
"after applying the one vector to A"));
237 GetOStream(
Runtime0) <<
"~~~~~~~~ Av with random vector ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ (Level " << fineLevel.
GetLevelID() <<
")" << std::endl;
238 Teuchos::RCP<Vector> randvec = VectorFactory::Build(A->getDomainMap(), 1);
239 randvec->randomize();
241 Teuchos::RCP<Vector> resrand = VectorFactory::Build(A->getRangeMap(),
false);
242 A->apply(*randvec, *resrand, Teuchos::NO_TRANS, one, zero);
243 checkVectorEntries(resrand,std::string(
"after applying random vector to A"));
248 GetOStream(
Runtime0) <<
"~~~~~~~~ Damped Jacobi sweep (one vector) ~~~~~~~~~~~~~~~~~~~~~~~ (Level " << fineLevel.
GetLevelID() <<
")" << std::endl;
249 Teuchos::RCP<Vector> ones = VectorFactory::Build(A->getDomainMap(), 1);
250 ones->putScalar(one);
252 Teuchos::RCP<Vector> res1 = VectorFactory::Build(A->getRangeMap(),
false);
253 A->apply(*ones, *res1, Teuchos::NO_TRANS, one, zero);
254 checkVectorEntries(res1,std::string(
"after applying one vector to A"));
257 checkVectorEntries(invDiag,std::string(
"in invDiag"));
259 Teuchos::RCP<Vector> res2 = VectorFactory::Build(A->getRangeMap(),
false);
260 res2->elementWiseMultiply (0.8, *invDiag, *res1, 0.0);
261 checkVectorEntries(res2,std::string(
"after scaling Av with invDiag (with v the one vector)"));
262 res2->update(1.0, *ones, -1.0);
264 checkVectorEntries(res2,std::string(
"after applying one damped Jacobi sweep (with one vector)"));
269 GetOStream(
Runtime0) <<
"~~~~~~~~ Damped Jacobi sweep (random vector) ~~~~~~~~~~~~~~~~~~~~ (Level " << fineLevel.
GetLevelID() <<
")" << std::endl;
270 Teuchos::RCP<Vector> ones = VectorFactory::Build(A->getDomainMap(), 1);
273 Teuchos::RCP<Vector> res1 = VectorFactory::Build(A->getRangeMap(),
false);
274 A->apply(*ones, *res1, Teuchos::NO_TRANS, one, zero);
275 checkVectorEntries(res1,std::string(
"after applying a random vector to A"));
278 checkVectorEntries(invDiag,std::string(
"in invDiag"));
280 Teuchos::RCP<Vector> res2 = VectorFactory::Build(A->getRangeMap(),
false);
281 res2->elementWiseMultiply (0.8, *invDiag, *res1, 0.0);
282 checkVectorEntries(res2,std::string(
"after scaling Av with invDiag (with v a random vector)"));
283 res2->update(1.0, *ones, -1.0);
285 checkVectorEntries(res2,std::string(
"after applying one damped Jacobi sweep (with v a random vector)"));
288 GetOStream(
Runtime0) <<
"~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ (Level " << fineLevel.
GetLevelID() <<
")" << std::endl;
291template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
293 SC zero = Teuchos::ScalarTraits<Scalar>::zero();
294 Teuchos::RCP<const Teuchos::Comm<int> > comm = vec->getMap()->getComm();
295 Teuchos::ArrayRCP< const Scalar > vecData = vec->getData(0);
301 for(
size_t i = 0; i<vec->getMap()->getLocalNumElements(); ++i) {
302 if(vecData[i] == zero) {
304 }
else if(Teuchos::ScalarTraits<Scalar>::isnaninf(vecData[i])) {
313 if(gzeros > 0) GetOStream(
Runtime0) <<
"Found " << gzeros <<
" zeros " << info << std::endl;
314 if(gnans > 0) GetOStream(
Runtime0) <<
"Found " << gnans <<
" entries " << info << std::endl;
#define MueLu_sumAll(rcpComm, in, out)
MueLu::DefaultLocalOrdinal LocalOrdinal
MueLu::DefaultGlobalOrdinal GlobalOrdinal
Exception throws to report errors in the internal logical of the program.
Timer to be used in factories. Similar to Monitor but with additional timers.
Class that holds all level-specific information.
int GetLevelID() const
Return level number.
void checkVectorEntries(const Teuchos::RCP< Vector > &vec, std::string info) const
virtual ~MatrixAnalysisFactory()
Destructor.
void DeclareInput(Level &fineLevel, Level &coarseLevel) const
Specifies the data that this class needs, and the factories that generate that data.
RCP< const ParameterList > GetValidParameterList() const
Return a const parameter list of valid parameters that setParameterList() will accept.
MatrixAnalysisFactory()
Constructor.
void Build(Level &fineLevel, Level &coarseLevel) const
Build an object with this factory.
static RCP< Xpetra::Vector< Scalar, LocalOrdinal, GlobalOrdinal, Node > > GetMatrixDiagonalInverse(const Xpetra::Matrix< Scalar, LocalOrdinal, GlobalOrdinal, Node > &A, Magnitude tol=Teuchos::ScalarTraits< Scalar >::eps() *100, Scalar valReplacement=Teuchos::ScalarTraits< Scalar >::zero())
Namespace for MueLu classes and methods.
@ Runtime0
One-liner description of what is happening.