MueLu  Version of the Day
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
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 
62 namespace MueLu {
63 template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
65  { }
66 
67 template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
69 
70 template <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 
79 template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
81  Input(fineLevel, "A");
82 }
83 
84 template <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->getNodeNumEntries ();
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()->getNodeNumElements(); ++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()->getNodeNumElements(); row++) {
183  GlobalOrdinal grow = A->getRowMap()->getGlobalElement(row);
184 
185  size_t nnz = A->getNumEntriesInLocalRow(row);
186 
187  // extract local row information from matrix
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 
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 
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 
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 
291 template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
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()->getNodeNumElements(); ++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)
void Build(Level &fineLevel, Level &coarseLevel) const
Build an object with this factory.
GlobalOrdinal GO
ParameterList & set(std::string const &name, T const &value, std::string const &docString="", RCP< const ParameterEntryValidator > const &validator=null)
Timer to be used in factories. Similar to Monitor but with additional timers.
#define TEUCHOS_TEST_FOR_EXCEPTION(throw_exception_test, Exception, msg)
size_type size() const
LocalOrdinal LO
One-liner description of what is happening.
static RCP< Xpetra::Vector< Scalar, LocalOrdinal, GlobalOrdinal, Node > > GetMatrixDiagonalInverse(const Xpetra::Matrix< Scalar, LocalOrdinal, GlobalOrdinal, Node > &A, Magnitude tol=Teuchos::ScalarTraits< Scalar >::eps()*100)
void checkVectorEntries(const Teuchos::RCP< Vector > &vec, std::string info) const
TEUCHOS_DEPRECATED RCP< T > rcp(T *p, Dealloc_T dealloc, bool owns_mem)
Class that holds all level-specific information.
Definition: MueLu_Level.hpp:99
RCP< const ParameterList > GetValidParameterList() const
Return a const parameter list of valid parameters that setParameterList() will accept.
TEUCHOS_DEPRECATED void reduceAll(const Comm< Ordinal > &comm, const EReductionType reductType, const Packet &send, Packet *globalReduct)
static magnitudeType magnitude(T a)
Scalar SC
void DeclareInput(Level &fineLevel, Level &coarseLevel) const
Specifies the data that this class needs, and the factories that generate that data.
Exception throws to report errors in the internal logical of the program.