MueLu Version of the Day
Loading...
Searching...
No Matches
MueLu_MatrixAnalysisFactory_def.hpp
Go to the documentation of this file.
1// @HEADER
2//
3// ***********************************************************************
4//
5// MueLu: A package for multigrid based preconditioning
6// Copyright 2012 Sandia Corporation
7//
8// Under the terms of Contract DE-AC04-94AL85000 with Sandia Corporation,
9// the U.S. Government retains certain rights in this software.
10//
11// Redistribution and use in source and binary forms, with or without
12// modification, are permitted provided that the following conditions are
13// met:
14//
15// 1. Redistributions of source code must retain the above copyright
16// notice, this list of conditions and the following disclaimer.
17//
18// 2. Redistributions in binary form must reproduce the above copyright
19// notice, this list of conditions and the following disclaimer in the
20// documentation and/or other materials provided with the distribution.
21//
22// 3. Neither the name of the Corporation nor the names of the
23// contributors may be used to endorse or promote products derived from
24// this software without specific prior written permission.
25//
26// THIS SOFTWARE IS PROVIDED BY SANDIA CORPORATION "AS IS" AND ANY
27// EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
28// IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR
29// PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL SANDIA CORPORATION OR THE
30// CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL,
31// EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO,
32// PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR
33// PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF
34// LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING
35// NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
36// SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
37//
38// Questions? Contact
39// Jonathan Hu (jhu@sandia.gov)
40// Andrey Prokopenko (aprokop@sandia.gov)
41// Ray Tuminaro (rstumin@sandia.gov)
42//
43// ***********************************************************************
44//
45// @HEADER
46
47#ifndef MUELU_MATRIXANALYSISFACTORY_DEF_HPP_
48#define MUELU_MATRIXANALYSISFACTORY_DEF_HPP_
49
51
52#include <Xpetra_Map.hpp>
53#include <Xpetra_StridedMap.hpp> // for nDofsPerNode...
54#include <Xpetra_Vector.hpp>
55#include <Xpetra_VectorFactory.hpp>
56#include <Xpetra_Matrix.hpp>
57
58#include "MueLu_Level.hpp"
59#include "MueLu_Utilities.hpp"
60#include "MueLu_Monitor.hpp"
61
62namespace MueLu {
63template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
66
67template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
69
70template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
72 RCP<ParameterList> validParamList = rcp(new ParameterList());
73
74 validParamList->set< RCP<const FactoryBase> >("A", Teuchos::null, "Generating factory of the matrix A to be permuted.");
75
76 return validParamList;
77}
78
79template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
81 Input(fineLevel, "A");
82}
83
84template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
86 FactoryMonitor m(*this, "MatrixAnalysis Factory ", fineLevel);
87
88 Teuchos::RCP<Matrix> A = Get< Teuchos::RCP<Matrix> > (fineLevel, "A");
89 Teuchos::RCP<const Teuchos::Comm<int> > comm = A->getRangeMap()->getComm();
90
91 //const ParameterList & pL = GetParameterList();
92
93 // General information
94 {
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;
97
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"));
100 LocalOrdinal blockid = strRowMap->getStridedBlockId();
101 if (blockid > -1) {
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;
105 } else {
106 GetOStream(Runtime0) << "Row: " << strRowMap->getFixedBlockSize() << " dofs per node" << std::endl;
107 }
108 }
109
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"));
112 LocalOrdinal blockid = strDomMap->getStridedBlockId();
113 if (blockid > -1) {
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;
117 } else {
118 GetOStream(Runtime0) << "Column: " << strDomMap->getFixedBlockSize() << " dofs per node" << std::endl;
119 }
120 }
121
122
123 GetOStream(Runtime0) << "A is distributed over " << comm->getSize() << " processors" << std::endl;
124
125 // empty processors
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;
134 }
135 }
136 }
137 }
138
139 // Matrix diagonal
140 {
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);
145 GlobalOrdinal lzerosOnDiagonal = 0;
146 GlobalOrdinal gzerosOnDiagonal = 0;
147 GlobalOrdinal lnearlyzerosOnDiagonal = 0;
148 GlobalOrdinal gnearlyzerosOnDiagonal = 0;
149 GlobalOrdinal lnansOnDiagonal = 0;
150 GlobalOrdinal gnansOnDiagonal = 0;
151
152 for(size_t i = 0; i<diagAVec->getMap()->getLocalNumElements(); ++i) {
153 if(diagAVecData[i] == Teuchos::ScalarTraits<Scalar>::zero()) {
154 lzerosOnDiagonal++;
155 } else if(Teuchos::ScalarTraits<Scalar>::magnitude(diagAVecData[i]) < 1e-6) {
156 lnearlyzerosOnDiagonal++;
157 } else if(Teuchos::ScalarTraits<Scalar>::isnaninf(diagAVecData[i])) {
158 lnansOnDiagonal++;
159 }
160 }
161
162 // sum up all entries in multipleColRequests over all processors
163 MueLu_sumAll(comm, lzerosOnDiagonal, gzerosOnDiagonal);
164 MueLu_sumAll(comm, lnearlyzerosOnDiagonal, gnearlyzerosOnDiagonal);
165 MueLu_sumAll(comm, lnansOnDiagonal, gnansOnDiagonal);
166
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;
170 }
171
172 // Diagonal dominance?
173 {
174 GetOStream(Runtime0) << "~~~~~~~~ DIAGONAL DOMINANCE ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ (Level " << fineLevel.GetLevelID() << ")" << std::endl;
175 // loop over all local rows in matrix A and keep diagonal entries if corresponding
176 // matrix rows are not contained in permRowMap
177 GlobalOrdinal lnumWeakDiagDomRows = 0;
178 GlobalOrdinal gnumWeakDiagDomRows = 0;
179 GlobalOrdinal lnumStrictDiagDomRows = 0;
180 GlobalOrdinal gnumStrictDiagDomRows = 0;
181 double worstRatio = 99999999.9;
182 for (size_t row = 0; row < A->getRowMap()->getLocalNumElements(); row++) {
183 GlobalOrdinal grow = A->getRowMap()->getGlobalElement(row);
184
185 size_t nnz = A->getNumEntriesInLocalRow(row);
186
187 // extract local row information from matrix
188 Teuchos::ArrayView<const LocalOrdinal> indices;
189 Teuchos::ArrayView<const Scalar> vals;
190 A->getLocalRowView(row, indices, vals);
191
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.");
193
194 // find column entry with max absolute value
195 double norm1 = 0.0;
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]));
201 }
202
203 if (normdiag >= norm1) lnumWeakDiagDomRows++;
204 else if (normdiag > norm1) lnumStrictDiagDomRows++;
205
206 if (norm1 != 0.0) {
207 if (normdiag / norm1 < worstRatio) worstRatio = normdiag / norm1;
208 }
209 }
210
211 MueLu_sumAll(comm, lnumWeakDiagDomRows, gnumWeakDiagDomRows);
212 MueLu_sumAll(comm, lnumStrictDiagDomRows, gnumStrictDiagDomRows);
213
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;
216
217 double gworstRatio;
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;
220 }
221
222
223 SC one = Teuchos::ScalarTraits<Scalar>::one(), zero = Teuchos::ScalarTraits<Scalar>::zero();
224
225 // multiply with one vector
226 {
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);
230
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"));
234 }
235
236 {
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();
240
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"));
244 }
245
246 // apply Jacobi sweep
247 {
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);
251
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"));
255
256 Teuchos::RCP<Vector> invDiag = Utilities::GetMatrixDiagonalInverse(*A);
257 checkVectorEntries(invDiag,std::string("in invDiag"));
258
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);
263
264 checkVectorEntries(res2,std::string("after applying one damped Jacobi sweep (with one vector)"));
265 }
266
267 // apply Jacobi sweep
268 {
269 GetOStream(Runtime0) << "~~~~~~~~ Damped Jacobi sweep (random vector) ~~~~~~~~~~~~~~~~~~~~ (Level " << fineLevel.GetLevelID() << ")" << std::endl;
270 Teuchos::RCP<Vector> ones = VectorFactory::Build(A->getDomainMap(), 1);
271 ones->randomize();
272
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"));
276
277 Teuchos::RCP<Vector> invDiag = Utilities::GetMatrixDiagonalInverse(*A);
278 checkVectorEntries(invDiag,std::string("in invDiag"));
279
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);
284
285 checkVectorEntries(res2,std::string("after applying one damped Jacobi sweep (with v a random vector)"));
286 }
287
288 GetOStream(Runtime0) << "~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ (Level " << fineLevel.GetLevelID() << ")" << std::endl;
289}
290
291template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
292void MatrixAnalysisFactory<Scalar, LocalOrdinal, GlobalOrdinal, Node>::checkVectorEntries(const Teuchos::RCP<Vector>& vec, std::string info) const {
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);
296 GO lzeros = 0;
297 GO gzeros = 0;
298 GO lnans = 0;
299 GO gnans = 0;
300
301 for(size_t i = 0; i<vec->getMap()->getLocalNumElements(); ++i) {
302 if(vecData[i] == zero) {
303 lzeros++;
304 } else if(Teuchos::ScalarTraits<Scalar>::isnaninf(vecData[i])) {
305 lnans++;
306 }
307 }
308
309 // sum up all entries in multipleColRequests over all processors
310 MueLu_sumAll(comm, lzeros, gzeros);
311 MueLu_sumAll(comm, lnans, gnans);
312
313 if(gzeros > 0) GetOStream(Runtime0) << "Found " << gzeros << " zeros " << info << std::endl;
314 if(gnans > 0) GetOStream(Runtime0) << "Found " << gnans << " entries " << info << std::endl;
315}
316
317} // namespace MueLu
318
319
320#endif /* MUELU_MATRIXANALYSISFACTORY_DEF_HPP_ */
#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
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.
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.