MueLu Version of the Day
Loading...
Searching...
No Matches
MueLu_UtilitiesBase_decl.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#ifndef MUELU_UTILITIESBASE_DECL_HPP
47#define MUELU_UTILITIESBASE_DECL_HPP
48
49#ifndef _WIN32
50#include <unistd.h> //necessary for "sleep" function in debugging methods (PauseForDebugging)
51#endif
52
53#include <string>
54
55#include "MueLu_ConfigDefs.hpp"
56
57#include <Teuchos_DefaultComm.hpp>
58#include <Teuchos_ScalarTraits.hpp>
59#include <Teuchos_ParameterList.hpp>
60
61#include <Xpetra_BlockedCrsMatrix_fwd.hpp>
62#include <Xpetra_CrsGraphFactory_fwd.hpp>
63#include <Xpetra_CrsGraph_fwd.hpp>
64#include <Xpetra_CrsMatrix_fwd.hpp>
65#include <Xpetra_CrsMatrixWrap_fwd.hpp>
66#include <Xpetra_Map_fwd.hpp>
67#include <Xpetra_BlockedMap_fwd.hpp>
68#include <Xpetra_MapFactory_fwd.hpp>
69#include <Xpetra_Matrix_fwd.hpp>
70#include <Xpetra_MatrixFactory_fwd.hpp>
71#include <Xpetra_MultiVector_fwd.hpp>
72#include <Xpetra_MultiVectorFactory_fwd.hpp>
73#include <Xpetra_Operator_fwd.hpp>
74#include <Xpetra_Vector_fwd.hpp>
75#include <Xpetra_BlockedMultiVector.hpp>
76#include <Xpetra_BlockedVector.hpp>
77#include <Xpetra_VectorFactory_fwd.hpp>
78#include <Xpetra_ExportFactory.hpp>
79
80#include <Xpetra_Import.hpp>
81#include <Xpetra_ImportFactory.hpp>
82#include <Xpetra_MatrixMatrix.hpp>
83#include <Xpetra_CrsGraph.hpp>
84#include <Xpetra_CrsGraphFactory.hpp>
85#include <Xpetra_CrsMatrixWrap.hpp>
86#include <Xpetra_StridedMap.hpp>
87
88#include "MueLu_Exceptions.hpp"
89
90
91namespace MueLu {
92
93// MPI helpers
94#define MueLu_sumAll(rcpComm, in, out) \
95 Teuchos::reduceAll(*rcpComm, Teuchos::REDUCE_SUM, in, Teuchos::outArg(out))
96#define MueLu_minAll(rcpComm, in, out) \
97 Teuchos::reduceAll(*rcpComm, Teuchos::REDUCE_MIN, in, Teuchos::outArg(out))
98#define MueLu_maxAll(rcpComm, in, out) \
99 Teuchos::reduceAll(*rcpComm, Teuchos::REDUCE_MAX, in, Teuchos::outArg(out))
100
108 template <class Scalar,
111 class Node = DefaultNode>
113 public:
114#undef MUELU_UTILITIESBASE_SHORT
115//#include "MueLu_UseShortNames.hpp"
116 private:
117 using CrsGraph = Xpetra::CrsGraph<LocalOrdinal,GlobalOrdinal,Node>;
118 using CrsGraphFactory = Xpetra::CrsGraphFactory<LocalOrdinal,GlobalOrdinal,Node>;
119 using CrsMatrixWrap = Xpetra::CrsMatrixWrap<Scalar,LocalOrdinal,GlobalOrdinal,Node>;
120 using CrsMatrix = Xpetra::CrsMatrix<Scalar,LocalOrdinal,GlobalOrdinal,Node>;
121 using Matrix = Xpetra::Matrix<Scalar,LocalOrdinal,GlobalOrdinal,Node>;
122 using Vector = Xpetra::Vector<Scalar,LocalOrdinal,GlobalOrdinal,Node>;
123 using BlockedVector = Xpetra::BlockedVector<Scalar,LocalOrdinal,GlobalOrdinal,Node>;
124 using MultiVector = Xpetra::MultiVector<Scalar,LocalOrdinal,GlobalOrdinal,Node>;
125 using BlockedMultiVector = Xpetra::BlockedMultiVector<Scalar,LocalOrdinal,GlobalOrdinal,Node>;
126 using BlockedMap = Xpetra::BlockedMap<LocalOrdinal,GlobalOrdinal,Node>;
127 using Map = Xpetra::Map<LocalOrdinal,GlobalOrdinal,Node>;
128 public:
129 typedef typename Teuchos::ScalarTraits<Scalar>::magnitudeType Magnitude;
130
131
132 static RCP<Matrix> Crs2Op(RCP<CrsMatrix> Op) {
133 if (Op.is_null())
134 return Teuchos::null;
135 return rcp(new CrsMatrixWrap(Op));
136 }
137
144 static RCP<CrsMatrixWrap> GetThresholdedMatrix(const RCP<Matrix>& Ain, const Scalar threshold, const bool keepDiagonal=true, const GlobalOrdinal expectedNNZperRow=-1) {
145
146 RCP<const Map> rowmap = Ain->getRowMap();
147 RCP<const Map> colmap = Ain->getColMap();
148 RCP<CrsMatrixWrap> Aout = rcp(new CrsMatrixWrap(rowmap, expectedNNZperRow <= 0 ? Ain->getGlobalMaxNumRowEntries() : expectedNNZperRow));
149 // loop over local rows
150 for(size_t row=0; row<Ain->getLocalNumRows(); row++)
151 {
152 size_t nnz = Ain->getNumEntriesInLocalRow(row);
153
154 Teuchos::ArrayView<const LocalOrdinal> indices;
155 Teuchos::ArrayView<const Scalar> vals;
156 Ain->getLocalRowView(row, indices, vals);
157
158 TEUCHOS_TEST_FOR_EXCEPTION(Teuchos::as<size_t>(indices.size()) != nnz, Exceptions::RuntimeError, "MueLu::ThresholdAFilterFactory::Build: number of nonzeros not equal to number of indices? Error.");
159
160 Teuchos::ArrayRCP<GlobalOrdinal> indout(indices.size(),Teuchos::ScalarTraits<GlobalOrdinal>::zero());
161 Teuchos::ArrayRCP<Scalar> valout(indices.size(),Teuchos::ScalarTraits<Scalar>::zero());
162 size_t nNonzeros = 0;
163 if (keepDiagonal) {
164 GlobalOrdinal glbRow = rowmap->getGlobalElement(row);
165 LocalOrdinal lclColIdx = colmap->getLocalElement(glbRow);
166 for(size_t i=0; i<(size_t)indices.size(); i++) {
167 if(Teuchos::ScalarTraits<Scalar>::magnitude(vals[i]) > Teuchos::ScalarTraits<Scalar>::magnitude(threshold) || indices[i]==lclColIdx) {
168 indout[nNonzeros] = colmap->getGlobalElement(indices[i]); // LID -> GID (column)
169 valout[nNonzeros] = vals[i];
170 nNonzeros++;
171 }
172 }
173 } else
174 for(size_t i=0; i<(size_t)indices.size(); i++) {
175 if(Teuchos::ScalarTraits<Scalar>::magnitude(vals[i]) > Teuchos::ScalarTraits<Scalar>::magnitude(threshold)) {
176 indout[nNonzeros] = colmap->getGlobalElement(indices[i]); // LID -> GID (column)
177 valout[nNonzeros] = vals[i];
178 nNonzeros++;
179 }
180 }
181
182 indout.resize(nNonzeros);
183 valout.resize(nNonzeros);
184
185 Aout->insertGlobalValues(Ain->getRowMap()->getGlobalElement(row), indout.view(0,indout.size()), valout.view(0,valout.size()));
186 }
187 Aout->fillComplete(Ain->getDomainMap(), Ain->getRangeMap());
188
189 return Aout;
190 }
191
198 static RCP<Xpetra::CrsGraph<LocalOrdinal, GlobalOrdinal, Node> > GetThresholdedGraph(const RCP<Matrix>& A, const Magnitude threshold, const GlobalOrdinal expectedNNZperRow=-1) {
199
200 using STS = Teuchos::ScalarTraits<Scalar>;
201 RCP<CrsGraph> sparsityPattern = CrsGraphFactory::Build(A->getRowMap(), expectedNNZperRow <= 0 ? A->getGlobalMaxNumRowEntries() : expectedNNZperRow);
202
203 RCP<Vector> diag = GetMatrixOverlappedDiagonal(*A);
204 ArrayRCP<const Scalar> D = diag->getData(0);
205
206 for(size_t row=0; row<A->getLocalNumRows(); row++)
207 {
208 ArrayView<const LocalOrdinal> indices;
209 ArrayView<const Scalar> vals;
210 A->getLocalRowView(row, indices, vals);
211
212 GlobalOrdinal globalRow = A->getRowMap()->getGlobalElement(row);
213 LocalOrdinal col = A->getColMap()->getLocalElement(globalRow);
214
215 const Scalar Dk = STS::magnitude(D[col]) > 0.0 ? STS::magnitude(D[col]) : 1.0;
216 Array<GlobalOrdinal> indicesNew;
217
218 for(size_t i=0; i<size_t(indices.size()); i++)
219 // keep diagonal per default
220 if(col == indices[i] || STS::magnitude(STS::squareroot(Dk)*vals[i]*STS::squareroot(Dk)) > STS::magnitude(threshold))
221 indicesNew.append(A->getColMap()->getGlobalElement(indices[i]));
222
223 sparsityPattern->insertGlobalIndices(globalRow, ArrayView<const GlobalOrdinal>(indicesNew.data(), indicesNew.length()));
224 }
225 sparsityPattern->fillComplete();
226
227 return sparsityPattern;
228 }
229
236 static Teuchos::ArrayRCP<Scalar> GetMatrixDiagonal(const Matrix& A) {
237 size_t numRows = A.getRowMap()->getLocalNumElements();
238 Teuchos::ArrayRCP<Scalar> diag(numRows);
239 Teuchos::ArrayView<const LocalOrdinal> cols;
240 Teuchos::ArrayView<const Scalar> vals;
241 for (size_t i = 0; i < numRows; ++i) {
242 A.getLocalRowView(i, cols, vals);
243 LocalOrdinal j = 0;
244 for (; j < cols.size(); ++j) {
245 if (Teuchos::as<size_t>(cols[j]) == i) {
246 diag[i] = vals[j];
247 break;
248 }
249 }
250 if (j == cols.size()) {
251 // Diagonal entry is absent
252 diag[i] = Teuchos::ScalarTraits<Scalar>::zero();
253 }
254 }
255 return diag;
256 }
257
264 static RCP<Vector> GetMatrixDiagonalInverse(const Matrix& A, Magnitude tol = Teuchos::ScalarTraits<Scalar>::eps()*100, Scalar valReplacement = Teuchos::ScalarTraits<Scalar>::zero()) {
265 Teuchos::TimeMonitor MM = *Teuchos::TimeMonitor::getNewTimer("UtilitiesBase::GetMatrixDiagonalInverse");
266
267 RCP<const Map> rowMap = A.getRowMap();
268 RCP<Vector> diag = Xpetra::VectorFactory<Scalar,LocalOrdinal,GlobalOrdinal,Node>::Build(rowMap,true);
269
270 A.getLocalDiagCopy(*diag);
271
272 RCP<Vector> inv = MueLu::UtilitiesBase<Scalar,LocalOrdinal,GlobalOrdinal,Node>::GetInverse(diag, tol, valReplacement);
273
274 return inv;
275 }
276
283 static Teuchos::RCP<Vector> GetLumpedMatrixDiagonal(Matrix const & A, const bool doReciprocal = false,
284 Magnitude tol = Teuchos::ScalarTraits<Scalar>::magnitude(Teuchos::ScalarTraits<Scalar>::zero()),
285 Scalar valReplacement = Teuchos::ScalarTraits<Scalar>::zero(),
286 const bool replaceSingleEntryRowWithZero = false,
287 const bool useAverageAbsDiagVal = false) {
288
289 typedef Teuchos::ScalarTraits<Scalar> TST;
290
291 RCP<Vector> diag = Teuchos::null;
292 const Scalar zero = TST::zero();
293 const Scalar one = TST::one();
294 const Scalar two = one + one;
295
296 Teuchos::RCP<const Matrix> rcpA = Teuchos::rcpFromRef(A);
297
298 RCP<const Xpetra::BlockedCrsMatrix<Scalar,LocalOrdinal,GlobalOrdinal,Node> > bA =
299 Teuchos::rcp_dynamic_cast<const Xpetra::BlockedCrsMatrix<Scalar,LocalOrdinal,GlobalOrdinal,Node> >(rcpA);
300 if(bA == Teuchos::null) {
301 RCP<const Map> rowMap = rcpA->getRowMap();
302 diag = Xpetra::VectorFactory<Scalar,LocalOrdinal,GlobalOrdinal,Node>::Build(rowMap,true);
303 ArrayRCP<Scalar> diagVals = diag->getDataNonConst(0);
304 Teuchos::Array<Scalar> regSum(diag->getLocalLength());
305 Teuchos::ArrayView<const LocalOrdinal> cols;
306 Teuchos::ArrayView<const Scalar> vals;
307
308 std::vector<int> nnzPerRow(rowMap->getLocalNumElements());
309
310 //FIXME 2021-10-22 JHU If this is called with doReciprocal=false, what should the correct behavior be? Currently,
311 //FIXME 2021-10-22 JHU the diagonal entry is set to be the sum of the absolute values of the row entries.
312
313 const Magnitude zeroMagn = TST::magnitude(zero);
314 Magnitude avgAbsDiagVal = TST::magnitude(zero);
315 int numDiagsEqualToOne = 0;
316 for (size_t i = 0; i < rowMap->getLocalNumElements(); ++i) {
317 nnzPerRow[i] = 0;
318 rcpA->getLocalRowView(i, cols, vals);
319 diagVals[i] = zero;
320 for (LocalOrdinal j = 0; j < cols.size(); ++j) {
321 regSum[i] += vals[j];
322 const Magnitude rowEntryMagn = TST::magnitude(vals[j]);
323 if (rowEntryMagn > zeroMagn)
324 nnzPerRow[i]++;
325 diagVals[i] += rowEntryMagn;
326 if (static_cast<size_t>(cols[j]) == i)
327 avgAbsDiagVal += rowEntryMagn;
328 }
329 if (nnzPerRow[i] == 1 && TST::magnitude(diagVals[i])==1.)
330 numDiagsEqualToOne++;
331 }
332 if (useAverageAbsDiagVal)
333 tol = TST::magnitude(100 * Teuchos::ScalarTraits<Scalar>::eps()) * (avgAbsDiagVal-numDiagsEqualToOne) / (rowMap->getLocalNumElements()-numDiagsEqualToOne);
334 if (doReciprocal) {
335 for (size_t i = 0; i < rowMap->getLocalNumElements(); ++i) {
336 if (replaceSingleEntryRowWithZero && nnzPerRow[i] <= static_cast<int>(1))
337 diagVals[i] = zero;
338 else if ((diagVals[i] != zero) && (TST::magnitude(diagVals[i]) < TST::magnitude(two*regSum[i])))
339 diagVals[i] = one / TST::magnitude((two*regSum[i]));
340 else {
341 if(TST::magnitude(diagVals[i]) > tol)
342 diagVals[i] = one / diagVals[i];
343 else {
344 diagVals[i] = valReplacement;
345 }
346 }
347 }
348 }
349 } else {
350 TEUCHOS_TEST_FOR_EXCEPTION(doReciprocal, Xpetra::Exceptions::RuntimeError,
351 "UtilitiesBase::GetLumpedMatrixDiagonal(): extracting reciprocal of diagonal of a blocked matrix is not supported");
352 diag = Xpetra::VectorFactory<Scalar,LocalOrdinal,GlobalOrdinal,Node>::Build(bA->getRangeMapExtractor()->getFullMap(),true);
353
354 for (size_t row = 0; row < bA->Rows(); ++row) {
355 for (size_t col = 0; col < bA->Cols(); ++col) {
356 if (!bA->getMatrix(row,col).is_null()) {
357 // if we are in Thyra mode, but the block (row,row) is again a blocked operator, we have to use (pseudo) Xpetra-style GIDs with offset!
358 bool bThyraMode = bA->getRangeMapExtractor()->getThyraMode() && (Teuchos::rcp_dynamic_cast<Xpetra::BlockedCrsMatrix<Scalar,LocalOrdinal,GlobalOrdinal,Node> >(bA->getMatrix(row,col)) == Teuchos::null);
359 RCP<Vector> ddtemp = bA->getRangeMapExtractor()->ExtractVector(diag,row,bThyraMode);
360 RCP<const Vector> dd = GetLumpedMatrixDiagonal(*(bA->getMatrix(row,col)));
361 ddtemp->update(Teuchos::as<Scalar>(1.0),*dd,Teuchos::as<Scalar>(1.0));
362 bA->getRangeMapExtractor()->InsertVector(ddtemp,row,diag,bThyraMode);
363 }
364 }
365 }
366
367 }
368
369 return diag;
370 }
371
378 static Teuchos::ArrayRCP<Magnitude> GetMatrixMaxMinusOffDiagonal(const Xpetra::Matrix<Scalar,LocalOrdinal,GlobalOrdinal,Node>& A) {
379 size_t numRows = A.getRowMap()->getLocalNumElements();
380 Magnitude ZERO = Teuchos::ScalarTraits<Magnitude>::zero();
381 Teuchos::ArrayRCP<Magnitude> maxvec(numRows);
382 Teuchos::ArrayView<const LocalOrdinal> cols;
383 Teuchos::ArrayView<const Scalar> vals;
384 for (size_t i = 0; i < numRows; ++i) {
385 A.getLocalRowView(i, cols, vals);
386 Magnitude mymax = ZERO;
387 for (LocalOrdinal j=0; j < cols.size(); ++j) {
388 if (Teuchos::as<size_t>(cols[j]) != i) {
389 mymax = std::max(mymax,-Teuchos::ScalarTraits<Scalar>::real(vals[j]));
390 }
391 }
392 maxvec[i] = mymax;
393 }
394 return maxvec;
395 }
396
397 static Teuchos::ArrayRCP<Magnitude> GetMatrixMaxMinusOffDiagonal(const Xpetra::Matrix<Scalar,LocalOrdinal,GlobalOrdinal,Node>& A, const Xpetra::Vector<LocalOrdinal,LocalOrdinal,GlobalOrdinal,Node> &BlockNumber) {
398 TEUCHOS_TEST_FOR_EXCEPTION(!A.getColMap()->isSameAs(*BlockNumber.getMap()),std::runtime_error,"GetMatrixMaxMinusOffDiagonal: BlockNumber must match's A's column map.");
399
400 Teuchos::ArrayRCP<const LocalOrdinal> block_id = BlockNumber.getData(0);
401
402 size_t numRows = A.getRowMap()->getLocalNumElements();
403 Magnitude ZERO = Teuchos::ScalarTraits<Magnitude>::zero();
404 Teuchos::ArrayRCP<Magnitude> maxvec(numRows);
405 Teuchos::ArrayView<const LocalOrdinal> cols;
406 Teuchos::ArrayView<const Scalar> vals;
407 for (size_t i = 0; i < numRows; ++i) {
408 A.getLocalRowView(i, cols, vals);
409 Magnitude mymax = ZERO;
410 for (LocalOrdinal j=0; j < cols.size(); ++j) {
411 if (Teuchos::as<size_t>(cols[j]) != i && block_id[i] == block_id[cols[j]]) {
412 mymax = std::max(mymax,-Teuchos::ScalarTraits<Scalar>::real(vals[j]));
413 }
414 }
415 // printf("A(%d,:) row_scale(block) = %6.4e\n",(int)i,mymax);
416
417 maxvec[i] = mymax;
418 }
419 return maxvec;
420 }
421
429 static Teuchos::RCP<Vector> GetInverse(Teuchos::RCP<const Vector> v, Magnitude tol = Teuchos::ScalarTraits<Scalar>::eps()*100, Scalar valReplacement = Teuchos::ScalarTraits<Scalar>::zero()) {
430
431 RCP<Vector> ret = Xpetra::VectorFactory<Scalar,LocalOrdinal,GlobalOrdinal,Node>::Build(v->getMap(),true);
432
433 // check whether input vector "v" is a BlockedVector
434 RCP<const BlockedVector> bv = Teuchos::rcp_dynamic_cast<const BlockedVector>(v);
435 if(bv.is_null() == false) {
436 RCP<BlockedVector> bret = Teuchos::rcp_dynamic_cast<BlockedVector>(ret);
437 TEUCHOS_TEST_FOR_EXCEPTION(bret.is_null() == true, MueLu::Exceptions::RuntimeError,"MueLu::UtilitiesBase::GetInverse: return vector should be of type BlockedVector");
438 RCP<const BlockedMap> bmap = bv->getBlockedMap();
439 for(size_t r = 0; r < bmap->getNumMaps(); ++r) {
440 RCP<const MultiVector> submvec = bv->getMultiVector(r,bmap->getThyraMode());
441 RCP<const Vector> subvec = submvec->getVector(0);
442 RCP<Vector> subvecinf = MueLu::UtilitiesBase<Scalar,LocalOrdinal,GlobalOrdinal,Node>::GetInverse(subvec,tol,valReplacement);
443 bret->setMultiVector(r, subvecinf, bmap->getThyraMode());
444 }
445 return ret;
446 }
447
448 // v is an {Epetra,Tpetra}Vector: work with the underlying raw data
449 ArrayRCP<Scalar> retVals = ret->getDataNonConst(0);
450 ArrayRCP<const Scalar> inputVals = v->getData(0);
451 for (size_t i = 0; i < v->getMap()->getLocalNumElements(); ++i) {
452 if(Teuchos::ScalarTraits<Scalar>::magnitude(inputVals[i]) > tol)
453 retVals[i] = Teuchos::ScalarTraits<Scalar>::one() / inputVals[i];
454 else
455 retVals[i] = valReplacement;
456 }
457 return ret;
458 }
459
467 static RCP<Vector> GetMatrixOverlappedDiagonal(const Matrix& A) {
468 RCP<const Map> rowMap = A.getRowMap(), colMap = A.getColMap();
469
470 // Undo block map (if we have one)
471 RCP<const BlockedMap> browMap = Teuchos::rcp_dynamic_cast<const BlockedMap>(rowMap);
472 if(!browMap.is_null()) rowMap = browMap->getMap();
473
474 RCP<Vector> localDiag = Xpetra::VectorFactory<Scalar,LocalOrdinal,GlobalOrdinal,Node>::Build(rowMap);
475 try {
476 const CrsMatrixWrap* crsOp = dynamic_cast<const CrsMatrixWrap*>(&A);
477 if (crsOp == NULL) {
478 throw Exceptions::RuntimeError("cast to CrsMatrixWrap failed");
479 }
480 Teuchos::ArrayRCP<size_t> offsets;
481 crsOp->getLocalDiagOffsets(offsets);
482 crsOp->getLocalDiagCopy(*localDiag,offsets());
483 }
484 catch (...) {
485 ArrayRCP<Scalar> localDiagVals = localDiag->getDataNonConst(0);
486 Teuchos::ArrayRCP<Scalar> diagVals = GetMatrixDiagonal(A);
487 for (LocalOrdinal i = 0; i < localDiagVals.size(); i++)
488 localDiagVals[i] = diagVals[i];
489 localDiagVals = diagVals = null;
490 }
491
492 RCP<Vector> diagonal = Xpetra::VectorFactory<Scalar,LocalOrdinal,GlobalOrdinal,Node>::Build(colMap);
493 RCP< const Xpetra::Import<LocalOrdinal,GlobalOrdinal,Node> > importer;
494 importer = A.getCrsGraph()->getImporter();
495 if (importer == Teuchos::null) {
496 importer = Xpetra::ImportFactory<LocalOrdinal,GlobalOrdinal,Node>::Build(rowMap, colMap);
497 }
498 diagonal->doImport(*localDiag, *(importer), Xpetra::INSERT);
499 return diagonal;
500 }
501
502
510 static RCP<Vector> GetMatrixOverlappedDeletedRowsum(const Matrix& A) {
511 using LO = LocalOrdinal;
512 using GO = GlobalOrdinal;
513 using SC = Scalar;
514 using STS = typename Teuchos::ScalarTraits<SC>;
515
516 // Undo block map (if we have one)
517 RCP<const Map> rowMap = A.getRowMap(), colMap = A.getColMap();
518 RCP<const BlockedMap> browMap = Teuchos::rcp_dynamic_cast<const BlockedMap>(rowMap);
519 if(!browMap.is_null()) rowMap = browMap->getMap();
520
521 RCP<Vector> local = Xpetra::VectorFactory<SC,LO,GO,Node>::Build(rowMap);
522 RCP<Vector> ghosted = Xpetra::VectorFactory<SC,LO,GO,Node>::Build(colMap,true);
523 ArrayRCP<SC> localVals = local->getDataNonConst(0);
524
525 for (LO row = 0; row < static_cast<LO>(A.getRowMap()->getLocalNumElements()); ++row) {
526 size_t nnz = A.getNumEntriesInLocalRow(row);
527 ArrayView<const LO> indices;
528 ArrayView<const SC> vals;
529 A.getLocalRowView(row, indices, vals);
530
531 SC si = STS::zero();
532
533 for (LO colID = 0; colID < static_cast<LO>(nnz); colID++) {
534 if(indices[colID] != row) {
535 si += vals[colID];
536 }
537 }
538 localVals[row] = si;
539 }
540
541 RCP< const Xpetra::Import<LO,GO,Node> > importer;
542 importer = A.getCrsGraph()->getImporter();
543 if (importer == Teuchos::null) {
544 importer = Xpetra::ImportFactory<LO,GO,Node>::Build(rowMap, colMap);
545 }
546 ghosted->doImport(*local, *(importer), Xpetra::INSERT);
547 return ghosted;
548 }
549
550
551
552 static RCP<Xpetra::Vector<Magnitude,LocalOrdinal,GlobalOrdinal,Node> >
554 RCP<const Map> rowMap = A.getRowMap(), colMap = A.getColMap();
555 using STS = typename Teuchos::ScalarTraits<Scalar>;
556 using MTS = typename Teuchos::ScalarTraits<Magnitude>;
557 using MT = Magnitude;
558 using LO = LocalOrdinal;
559 using GO = GlobalOrdinal;
560 using SC = Scalar;
561 using RealValuedVector = Xpetra::Vector<MT,LO,GO,Node>;
562
563 // Undo block map (if we have one)
564 RCP<const BlockedMap> browMap = Teuchos::rcp_dynamic_cast<const BlockedMap>(rowMap);
565 if(!browMap.is_null()) rowMap = browMap->getMap();
566
567 RCP<RealValuedVector> local = Xpetra::VectorFactory<MT,LO,GO,Node>::Build(rowMap);
568 RCP<RealValuedVector> ghosted = Xpetra::VectorFactory<MT,LO,GO,Node>::Build(colMap,true);
569 ArrayRCP<MT> localVals = local->getDataNonConst(0);
570
571 for (LO rowIdx = 0; rowIdx < static_cast<LO>(A.getRowMap()->getLocalNumElements()); ++rowIdx) {
572 size_t nnz = A.getNumEntriesInLocalRow(rowIdx);
573 ArrayView<const LO> indices;
574 ArrayView<const SC> vals;
575 A.getLocalRowView(rowIdx, indices, vals);
576
577 MT si = MTS::zero();
578
579 for (LO colID = 0; colID < static_cast<LO>(nnz); ++colID) {
580 if(indices[colID] != rowIdx) {
581 si += STS::magnitude(vals[colID]);
582 }
583 }
584 localVals[rowIdx] = si;
585 }
586
587 RCP< const Xpetra::Import<LO,GO,Node> > importer;
588 importer = A.getCrsGraph()->getImporter();
589 if (importer == Teuchos::null) {
590 importer = Xpetra::ImportFactory<LO,GO,Node>::Build(rowMap, colMap);
591 }
592 ghosted->doImport(*local, *(importer), Xpetra::INSERT);
593 return ghosted;
594 }
595
596
597
598 // TODO: should NOT return an Array. Definition must be changed to:
599 // - ArrayRCP<> ResidualNorm(Matrix const &Op, MultiVector const &X, MultiVector const &RHS)
600 // or
601 // - void ResidualNorm(Matrix const &Op, MultiVector const &X, MultiVector const &RHS, Array &)
602 static Teuchos::Array<Magnitude> ResidualNorm(const Xpetra::Operator<Scalar,LocalOrdinal,GlobalOrdinal,Node>& Op, const MultiVector& X, const MultiVector& RHS) {
603 TEUCHOS_TEST_FOR_EXCEPTION(X.getNumVectors() != RHS.getNumVectors(), Exceptions::RuntimeError, "Number of solution vectors != number of right-hand sides")
604 const size_t numVecs = X.getNumVectors();
605 RCP<MultiVector> RES = Residual(Op, X, RHS);
606 Teuchos::Array<Magnitude> norms(numVecs);
607 RES->norm2(norms);
608 return norms;
609 }
610
611 static Teuchos::Array<Magnitude> ResidualNorm(const Xpetra::Operator<Scalar,LocalOrdinal,GlobalOrdinal,Node>& Op, const MultiVector& X, const MultiVector& RHS, MultiVector & Resid) {
612 TEUCHOS_TEST_FOR_EXCEPTION(X.getNumVectors() != RHS.getNumVectors(), Exceptions::RuntimeError, "Number of solution vectors != number of right-hand sides")
613 const size_t numVecs = X.getNumVectors();
614 Residual(Op,X,RHS,Resid);
615 Teuchos::Array<Magnitude> norms(numVecs);
616 Resid.norm2(norms);
617 return norms;
618 }
619
620 static RCP<MultiVector> Residual(const Xpetra::Operator<Scalar,LocalOrdinal,GlobalOrdinal,Node>& Op, const MultiVector& X, const MultiVector& RHS) {
621 TEUCHOS_TEST_FOR_EXCEPTION(X.getNumVectors() != RHS.getNumVectors(), Exceptions::RuntimeError, "Number of solution vectors != number of right-hand sides")
622 const size_t numVecs = X.getNumVectors();
623 // TODO Op.getRangeMap should return a BlockedMap if it is a BlockedCrsOperator
624 RCP<MultiVector> RES = Xpetra::MultiVectorFactory<Scalar,LocalOrdinal,GlobalOrdinal,Node>::Build(RHS.getMap(), numVecs, false); // no need to initialize to zero
625 Op.residual(X,RHS,*RES);
626 return RES;
627 }
628
629
630 static void Residual(const Xpetra::Operator<Scalar,LocalOrdinal,GlobalOrdinal,Node>& Op, const MultiVector& X, const MultiVector& RHS, MultiVector & Resid) {
631 TEUCHOS_TEST_FOR_EXCEPTION(X.getNumVectors() != RHS.getNumVectors(), Exceptions::RuntimeError, "Number of solution vectors != number of right-hand sides");
632 TEUCHOS_TEST_FOR_EXCEPTION(Resid.getNumVectors() != RHS.getNumVectors(), Exceptions::RuntimeError, "Number of residual vectors != number of right-hand sides");
633 Op.residual(X,RHS,Resid);
634 }
635
636
648 static Scalar PowerMethod(const Matrix& A, bool scaleByDiag = true,
649 LocalOrdinal niters = 10, Magnitude tolerance = 1e-2, bool verbose = false, unsigned int seed = 123) {
650 TEUCHOS_TEST_FOR_EXCEPTION(!(A.getRangeMap()->isSameAs(*(A.getDomainMap()))), Exceptions::Incompatible,
651 "Utils::PowerMethod: operator must have domain and range maps that are equivalent.");
652
653 // power iteration
654 RCP<Vector> diagInvVec;
655 if (scaleByDiag) {
656 RCP<Vector> diagVec = Xpetra::VectorFactory<Scalar,LocalOrdinal,GlobalOrdinal,Node>::Build(A.getRowMap());
657 A.getLocalDiagCopy(*diagVec);
658 diagInvVec = Xpetra::VectorFactory<Scalar,LocalOrdinal,GlobalOrdinal,Node>::Build(A.getRowMap());
659 diagInvVec->reciprocal(*diagVec);
660 }
661
662 Scalar lambda = PowerMethod(A, diagInvVec, niters, tolerance, verbose, seed);
663 return lambda;
664 }
665
677 static Scalar PowerMethod(const Matrix& A, const RCP<Vector> &diagInvVec,
678 LocalOrdinal niters = 10, Magnitude tolerance = 1e-2, bool verbose = false, unsigned int seed = 123) {
679 TEUCHOS_TEST_FOR_EXCEPTION(!(A.getRangeMap()->isSameAs(*(A.getDomainMap()))), Exceptions::Incompatible,
680 "Utils::PowerMethod: operator must have domain and range maps that are equivalent.");
681
682 // Create three vectors, fill z with random numbers
683 RCP<Vector> q = Xpetra::VectorFactory<Scalar,LocalOrdinal,GlobalOrdinal,Node>::Build(A.getDomainMap());
684 RCP<Vector> r = Xpetra::VectorFactory<Scalar,LocalOrdinal,GlobalOrdinal,Node>::Build(A.getRangeMap());
685 RCP<Vector> z = Xpetra::VectorFactory<Scalar,LocalOrdinal,GlobalOrdinal,Node>::Build(A.getRangeMap());
686
687 z->setSeed(seed); // seed random number generator
688 z->randomize(true);// use Xpetra implementation: -> same results for Epetra and Tpetra
689
690 Teuchos::Array<Magnitude> norms(1);
691
692 typedef Teuchos::ScalarTraits<Scalar> STS;
693
694 const Scalar zero = STS::zero(), one = STS::one();
695
696 Scalar lambda = zero;
697 Magnitude residual = STS::magnitude(zero);
698
699 // power iteration
700 for (int iter = 0; iter < niters; ++iter) {
701 z->norm2(norms); // Compute 2-norm of z
702 q->update(one/norms[0], *z, zero); // Set q = z / normz
703 A.apply(*q, *z); // Compute z = A*q
704 if (diagInvVec != Teuchos::null)
705 z->elementWiseMultiply(one, *diagInvVec, *z, zero);
706 lambda = q->dot(*z); // Approximate maximum eigenvalue: lamba = dot(q,z)
707
708 if (iter % 100 == 0 || iter + 1 == niters) {
709 r->update(1.0, *z, -lambda, *q, zero); // Compute A*q - lambda*q
710 r->norm2(norms);
711 residual = STS::magnitude(norms[0] / lambda);
712 if (verbose) {
713 std::cout << "Iter = " << iter
714 << " Lambda = " << lambda
715 << " Residual of A*q - lambda*q = " << residual
716 << std::endl;
717 }
718 }
719 if (residual < tolerance)
720 break;
721 }
722 return lambda;
723 }
724
725 static RCP<Teuchos::FancyOStream> MakeFancy(std::ostream& os) {
726 RCP<Teuchos::FancyOStream> fancy = Teuchos::fancyOStream(Teuchos::rcpFromRef(os));
727 return fancy;
728 }
729
734 static typename Teuchos::ScalarTraits<Scalar>::magnitudeType Distance2(const Teuchos::Array<Teuchos::ArrayRCP<const Scalar>> &v, LocalOrdinal i0, LocalOrdinal i1) {
735 const size_t numVectors = v.size();
736
737 Scalar d = Teuchos::ScalarTraits<Scalar>::zero();
738 for (size_t j = 0; j < numVectors; j++) {
739 d += (v[j][i0] - v[j][i1])*(v[j][i0] - v[j][i1]);
740 }
741 return Teuchos::ScalarTraits<Scalar>::magnitude(d);
742 }
743
748 static typename Teuchos::ScalarTraits<Scalar>::magnitudeType Distance2(const Teuchos::ArrayView<double> & weight,const Teuchos::Array<Teuchos::ArrayRCP<const Scalar>> &v, LocalOrdinal i0, LocalOrdinal i1) {
749 const size_t numVectors = v.size();
750 using MT = typename Teuchos::ScalarTraits<Scalar>::magnitudeType;
751
752 Scalar d = Teuchos::ScalarTraits<Scalar>::zero();
753 for (size_t j = 0; j < numVectors; j++) {
754 d += Teuchos::as<MT>(weight[j])*(v[j][i0] - v[j][i1])*(v[j][i0] - v[j][i1]);
755 }
756 return Teuchos::ScalarTraits<Scalar>::magnitude(d);
757 }
758
759
772 static Teuchos::ArrayRCP<const bool> DetectDirichletRows(const Xpetra::Matrix<Scalar,LocalOrdinal,GlobalOrdinal,Node>& A, const Magnitude& tol = Teuchos::ScalarTraits<Scalar>::zero(), bool count_twos_as_dirichlet=false) {
773 LocalOrdinal numRows = A.getLocalNumRows();
774 typedef Teuchos::ScalarTraits<Scalar> STS;
775 ArrayRCP<bool> boundaryNodes(numRows, true);
776 if (count_twos_as_dirichlet) {
777 for (LocalOrdinal row = 0; row < numRows; row++) {
778 ArrayView<const LocalOrdinal> indices;
779 ArrayView<const Scalar> vals;
780 A.getLocalRowView(row, indices, vals);
781 size_t nnz = A.getNumEntriesInLocalRow(row);
782 if (nnz > 2) {
783 size_t col;
784 for (col = 0; col < nnz; col++)
785 if ( (indices[col] != row) && STS::magnitude(vals[col]) > tol) {
786 if (!boundaryNodes[row])
787 break;
788 boundaryNodes[row] = false;
789 }
790 if (col == nnz)
791 boundaryNodes[row] = true;
792 }
793 }
794 } else {
795 for (LocalOrdinal row = 0; row < numRows; row++) {
796 ArrayView<const LocalOrdinal> indices;
797 ArrayView<const Scalar> vals;
798 A.getLocalRowView(row, indices, vals);
799 size_t nnz = A.getNumEntriesInLocalRow(row);
800 if (nnz > 1)
801 for (size_t col = 0; col < nnz; col++)
802 if ( (indices[col] != row) && STS::magnitude(vals[col]) > tol) {
803 boundaryNodes[row] = false;
804 break;
805 }
806 }
807 }
808 return boundaryNodes;
809 }
810
811
824 static Teuchos::ArrayRCP<const bool> DetectDirichletRowsExt(const Xpetra::Matrix<Scalar,LocalOrdinal,GlobalOrdinal,Node>& A, bool & bHasZeroDiagonal, const Magnitude& tol = Teuchos::ScalarTraits<Scalar>::zero()) {
825
826 // assume that there is no zero diagonal in matrix
827 bHasZeroDiagonal = false;
828
829 Teuchos::RCP<Vector> diagVec = Xpetra::VectorFactory<Scalar,LocalOrdinal,GlobalOrdinal,Node>::Build(A.getRowMap());
830 A.getLocalDiagCopy(*diagVec);
831 Teuchos::ArrayRCP< const Scalar > diagVecData = diagVec->getData(0);
832
833 LocalOrdinal numRows = A.getLocalNumRows();
834 typedef Teuchos::ScalarTraits<Scalar> STS;
835 ArrayRCP<bool> boundaryNodes(numRows, false);
836 for (LocalOrdinal row = 0; row < numRows; row++) {
837 ArrayView<const LocalOrdinal> indices;
838 ArrayView<const Scalar> vals;
839 A.getLocalRowView(row, indices, vals);
840 size_t nnz = 0; // collect nonzeros in row (excluding the diagonal)
841 bool bHasDiag = false;
842 for (decltype(indices.size()) col = 0; col < indices.size(); col++) {
843 if ( indices[col] != row) {
844 if (STS::magnitude(vals[col] / STS::magnitude(sqrt(STS::magnitude(diagVecData[row]) * STS::magnitude(diagVecData[col]))) ) > tol) {
845 nnz++;
846 }
847 } else bHasDiag = true; // found a diagonal entry
848 }
849 if (bHasDiag == false) bHasZeroDiagonal = true; // we found at least one row without a diagonal
850 else if(nnz == 0) boundaryNodes[row] = true;
851 }
852 return boundaryNodes;
853 }
854
862 static void FindNonZeros(const Teuchos::ArrayRCP<const Scalar> vals,
863 Teuchos::ArrayRCP<bool> nonzeros) {
864 TEUCHOS_ASSERT(vals.size() == nonzeros.size());
865 typedef typename Teuchos::ScalarTraits<Scalar>::magnitudeType magnitudeType;
866 const magnitudeType eps = 2.0*Teuchos::ScalarTraits<magnitudeType>::eps();
867 for(size_t i=0; i<static_cast<size_t>(vals.size()); i++) {
868 nonzeros[i] = (Teuchos::ScalarTraits<Scalar>::magnitude(vals[i]) > eps);
869 }
870 }
871
880 static void DetectDirichletColsAndDomains(const Xpetra::Matrix<Scalar,LocalOrdinal,GlobalOrdinal,Node>& A,
881 const Teuchos::ArrayRCP<bool>& dirichletRows,
882 Teuchos::ArrayRCP<bool> dirichletCols,
883 Teuchos::ArrayRCP<bool> dirichletDomain) {
884 const Scalar one = Teuchos::ScalarTraits<Scalar>::one();
885 RCP<const Xpetra::Map<LocalOrdinal,GlobalOrdinal,Node> > domMap = A .getDomainMap();
886 RCP<const Xpetra::Map<LocalOrdinal,GlobalOrdinal,Node> > rowMap = A.getRowMap();
887 RCP<const Xpetra::Map<LocalOrdinal,GlobalOrdinal,Node> > colMap = A.getColMap();
888 TEUCHOS_ASSERT(static_cast<size_t>(dirichletRows.size()) == rowMap->getLocalNumElements());
889 TEUCHOS_ASSERT(static_cast<size_t>(dirichletCols.size()) == colMap->getLocalNumElements());
890 TEUCHOS_ASSERT(static_cast<size_t>(dirichletDomain.size()) == domMap->getLocalNumElements());
891 RCP<Xpetra::MultiVector<Scalar,LocalOrdinal,GlobalOrdinal,Node> > myColsToZero = Xpetra::MultiVectorFactory<Scalar,LocalOrdinal,GlobalOrdinal,Node>::Build(colMap, 1, /*zeroOut=*/true);
892 // Find all local column indices that are in Dirichlet rows, record in myColsToZero as 1.0
893 for(size_t i=0; i<(size_t) dirichletRows.size(); i++) {
894 if (dirichletRows[i]) {
895 ArrayView<const LocalOrdinal> indices;
896 ArrayView<const Scalar> values;
897 A.getLocalRowView(i,indices,values);
898 for(size_t j=0; j<static_cast<size_t>(indices.size()); j++)
899 myColsToZero->replaceLocalValue(indices[j],0,one);
900 }
901 }
902
903 RCP<Xpetra::MultiVector<Scalar,LocalOrdinal,GlobalOrdinal,Node> > globalColsToZero;
904 RCP<const Xpetra::Import<LocalOrdinal,GlobalOrdinal,Node> > importer = A.getCrsGraph()->getImporter();
905 if (!importer.is_null()) {
906 globalColsToZero = Xpetra::MultiVectorFactory<Scalar,LocalOrdinal,GlobalOrdinal,Node>::Build(domMap, 1, /*zeroOut=*/true);
907 // export to domain map
908 globalColsToZero->doExport(*myColsToZero,*importer,Xpetra::ADD);
909 // import to column map
910 myColsToZero->doImport(*globalColsToZero,*importer,Xpetra::INSERT);
911 }
912 else
913 globalColsToZero = myColsToZero;
914
915 FindNonZeros(globalColsToZero->getData(0),dirichletDomain);
916 FindNonZeros(myColsToZero->getData(0),dirichletCols);
917 }
918
919
920
932 static void ApplyRowSumCriterion(const Xpetra::Matrix<Scalar,LocalOrdinal,GlobalOrdinal,Node>& A, const Magnitude rowSumTol, Teuchos::ArrayRCP<bool>& dirichletRows) {
933 typedef Teuchos::ScalarTraits<Scalar> STS;
934 typedef typename Teuchos::ScalarTraits<Scalar>::magnitudeType MT;
935 typedef Teuchos::ScalarTraits<MT> MTS;
936 RCP<const Xpetra::Map<LocalOrdinal,GlobalOrdinal,Node>> rowmap = A.getRowMap();
937 for (LocalOrdinal row = 0; row < Teuchos::as<LocalOrdinal>(rowmap->getLocalNumElements()); ++row) {
938 size_t nnz = A.getNumEntriesInLocalRow(row);
939 ArrayView<const LocalOrdinal> indices;
940 ArrayView<const Scalar> vals;
941 A.getLocalRowView(row, indices, vals);
942
943 Scalar rowsum = STS::zero();
944 Scalar diagval = STS::zero();
945
946 for (LocalOrdinal colID = 0; colID < Teuchos::as<LocalOrdinal>(nnz); colID++) {
947 LocalOrdinal col = indices[colID];
948 if (row == col)
949 diagval = vals[colID];
950 rowsum += vals[colID];
951 }
952 // printf("A(%d,:) row_sum(point) = %6.4e\n",row,rowsum);
953 if (rowSumTol < MTS::one() && STS::magnitude(rowsum) > STS::magnitude(diagval) * rowSumTol) {
954 //printf("Row %d triggers rowsum\n",(int)row);
955 dirichletRows[row] = true;
956 }
957 }
958 }
959
960 static void ApplyRowSumCriterion(const Xpetra::Matrix<Scalar,LocalOrdinal,GlobalOrdinal,Node>& A, const Xpetra::Vector<LocalOrdinal,LocalOrdinal,GlobalOrdinal,Node> &BlockNumber, const Magnitude rowSumTol, Teuchos::ArrayRCP<bool>& dirichletRows) {
961 typedef Teuchos::ScalarTraits<Scalar> STS;
962 typedef typename Teuchos::ScalarTraits<Scalar>::magnitudeType MT;
963 typedef Teuchos::ScalarTraits<MT> MTS;
964 RCP<const Xpetra::Map<LocalOrdinal,GlobalOrdinal,Node> > rowmap = A.getRowMap();
965
966 TEUCHOS_TEST_FOR_EXCEPTION(!A.getColMap()->isSameAs(*BlockNumber.getMap()),std::runtime_error,"ApplyRowSumCriterion: BlockNumber must match's A's column map.");
967
968 Teuchos::ArrayRCP<const LocalOrdinal> block_id = BlockNumber.getData(0);
969 for (LocalOrdinal row = 0; row < Teuchos::as<LocalOrdinal>(rowmap->getLocalNumElements()); ++row) {
970 size_t nnz = A.getNumEntriesInLocalRow(row);
971 ArrayView<const LocalOrdinal> indices;
972 ArrayView<const Scalar> vals;
973 A.getLocalRowView(row, indices, vals);
974
975 Scalar rowsum = STS::zero();
976 Scalar diagval = STS::zero();
977 for (LocalOrdinal colID = 0; colID < Teuchos::as<LocalOrdinal>(nnz); colID++) {
978 LocalOrdinal col = indices[colID];
979 if (row == col)
980 diagval = vals[colID];
981 if(block_id[row] == block_id[col])
982 rowsum += vals[colID];
983 }
984
985 // printf("A(%d,:) row_sum(block) = %6.4e\n",row,rowsum);
986 if (rowSumTol < MTS::one() && STS::magnitude(rowsum) > STS::magnitude(diagval) * rowSumTol) {
987 //printf("Row %d triggers rowsum\n",(int)row);
988 dirichletRows[row] = true;
989 }
990 }
991 }
992
993
994
1005 static Teuchos::ArrayRCP<const bool> DetectDirichletCols(const Xpetra::Matrix<Scalar,LocalOrdinal,GlobalOrdinal,Node>& A,
1006 const Teuchos::ArrayRCP<const bool>& dirichletRows) {
1007 Scalar zero = Teuchos::ScalarTraits<Scalar>::zero();
1008 Scalar one = Teuchos::ScalarTraits<Scalar>::one();
1009 Teuchos::RCP<const Xpetra::Map<LocalOrdinal,GlobalOrdinal,Node> > domMap = A.getDomainMap();
1010 Teuchos::RCP<const Xpetra::Map<LocalOrdinal,GlobalOrdinal,Node> > colMap = A.getColMap();
1011 Teuchos::RCP<Xpetra::MultiVector<Scalar,LocalOrdinal,GlobalOrdinal,Node> > myColsToZero = Xpetra::MultiVectorFactory<Scalar,LocalOrdinal,GlobalOrdinal,Node>::Build(colMap,1);
1012 myColsToZero->putScalar(zero);
1013 // Find all local column indices that are in Dirichlet rows, record in myColsToZero as 1.0
1014 for(size_t i=0; i<(size_t) dirichletRows.size(); i++) {
1015 if (dirichletRows[i]) {
1016 Teuchos::ArrayView<const LocalOrdinal> indices;
1017 Teuchos::ArrayView<const Scalar> values;
1018 A.getLocalRowView(i,indices,values);
1019 for(size_t j=0; j<static_cast<size_t>(indices.size()); j++)
1020 myColsToZero->replaceLocalValue(indices[j],0,one);
1021 }
1022 }
1023
1024 Teuchos::RCP<Xpetra::MultiVector<Scalar,LocalOrdinal,GlobalOrdinal,Node> > globalColsToZero = Xpetra::MultiVectorFactory<Scalar,LocalOrdinal,GlobalOrdinal,Node>::Build(domMap,1);
1025 globalColsToZero->putScalar(zero);
1026 Teuchos::RCP<Xpetra::Export<LocalOrdinal,GlobalOrdinal,Node> > exporter = Xpetra::ExportFactory<LocalOrdinal,GlobalOrdinal,Node>::Build(colMap,domMap);
1027 // export to domain map
1028 globalColsToZero->doExport(*myColsToZero,*exporter,Xpetra::ADD);
1029 // import to column map
1030 myColsToZero->doImport(*globalColsToZero,*exporter,Xpetra::INSERT);
1031 Teuchos::ArrayRCP<const Scalar> myCols = myColsToZero->getData(0);
1032 Teuchos::ArrayRCP<bool> dirichletCols(colMap->getLocalNumElements(), true);
1033 Magnitude eps = Teuchos::ScalarTraits<Magnitude>::eps();
1034 for(size_t i=0; i<colMap->getLocalNumElements(); i++) {
1035 dirichletCols[i] = Teuchos::ScalarTraits<Scalar>::magnitude(myCols[i])>2.0*eps;
1036 }
1037 return dirichletCols;
1038 }
1039
1040
1045 static Scalar Frobenius(const Xpetra::Matrix<Scalar,LocalOrdinal,GlobalOrdinal,Node>& A, const Xpetra::Matrix<Scalar,LocalOrdinal,GlobalOrdinal,Node>& B) {
1046 // We check only row maps. Column may be different. One would hope that they are the same, as we typically
1047 // calculate frobenius norm of the specified sparsity pattern with an updated matrix from the previous step,
1048 // but matrix addition, even when one is submatrix of the other, changes column map (though change may be as
1049 // simple as couple of elements swapped)
1050 TEUCHOS_TEST_FOR_EXCEPTION(!A.getRowMap()->isSameAs(*B.getRowMap()), Exceptions::Incompatible, "MueLu::CGSolver::Frobenius: row maps are incompatible");
1051 TEUCHOS_TEST_FOR_EXCEPTION(!A.isFillComplete() || !B.isFillComplete(), Exceptions::RuntimeError, "Matrices must be fill completed");
1052
1053 const Map& AColMap = *A.getColMap();
1054 const Map& BColMap = *B.getColMap();
1055
1056 Teuchos::ArrayView<const LocalOrdinal> indA, indB;
1057 Teuchos::ArrayView<const Scalar> valA, valB;
1058 size_t nnzA = 0, nnzB = 0;
1059
1060 // We use a simple algorithm
1061 // for each row we fill valBAll array with the values in the corresponding row of B
1062 // as such, it serves as both sorted array and as storage, so we don't need to do a
1063 // tricky problem: "find a value in the row of B corresponding to the specific GID"
1064 // Once we do that, we translate LID of entries of row of A to LID of B, and multiply
1065 // corresponding entries.
1066 // The algorithm should be reasonably cheap, as it does not sort anything, provided
1067 // that getLocalElement and getGlobalElement functions are reasonably effective. It
1068 // *is* possible that the costs are hidden in those functions, but if maps are close
1069 // to linear maps, we should be fine
1070 Teuchos::Array<Scalar> valBAll(BColMap.getLocalNumElements());
1071
1072 LocalOrdinal invalid = Teuchos::OrdinalTraits<LocalOrdinal>::invalid();
1073 Scalar zero = Teuchos::ScalarTraits<Scalar> ::zero(), f = zero, gf;
1074 size_t numRows = A.getLocalNumRows();
1075 for (size_t i = 0; i < numRows; i++) {
1076 A.getLocalRowView(i, indA, valA);
1077 B.getLocalRowView(i, indB, valB);
1078 nnzA = indA.size();
1079 nnzB = indB.size();
1080
1081 // Set up array values
1082 for (size_t j = 0; j < nnzB; j++)
1083 valBAll[indB[j]] = valB[j];
1084
1085 for (size_t j = 0; j < nnzA; j++) {
1086 // The cost of the whole Frobenius dot product function depends on the
1087 // cost of the getLocalElement and getGlobalElement functions here.
1088 LocalOrdinal ind = BColMap.getLocalElement(AColMap.getGlobalElement(indA[j]));
1089 if (ind != invalid)
1090 f += valBAll[ind] * valA[j];
1091 }
1092
1093 // Clean up array values
1094 for (size_t j = 0; j < nnzB; j++)
1095 valBAll[indB[j]] = zero;
1096 }
1097
1098 MueLu_sumAll(AColMap.getComm(), f, gf);
1099
1100 return gf;
1101 }
1102
1112 static void SetRandomSeed(const Teuchos::Comm<int> &comm) {
1113 // Distribute the seeds evenly in [1,maxint-1]. This guarantees nothing
1114 // about where in random number stream we are, but avoids overflow situations
1115 // in parallel when multiplying by a PID. It would be better to use
1116 // a good parallel random number generator.
1117 double one = 1.0;
1118 int maxint = INT_MAX; //= 2^31-1 = 2147483647 for 32-bit integers
1119 int mySeed = Teuchos::as<int>((maxint-1) * (one -(comm.getRank()+1)/(comm.getSize()+one)) );
1120 if (mySeed < 1 || mySeed == maxint) {
1121 std::ostringstream errStr;
1122 errStr << "Error detected with random seed = " << mySeed << ". It should be in the interval [1,2^31-2].";
1123 throw Exceptions::RuntimeError(errStr.str());
1124 }
1125 std::srand(mySeed);
1126 // For Tpetra, we could use Kokkos' random number generator here.
1127 Teuchos::ScalarTraits<Scalar>::seedrandom(mySeed);
1128 // Epetra
1129 // MultiVector::Random() -> Epetra_Util::RandomDouble() -> Epetra_Utils::RandomInt()
1130 // Its own random number generator, based on Seed_. Seed_ is initialized in Epetra_Util constructor with std::rand()
1131 // So our setting std::srand() affects that too
1132 }
1133
1134
1135
1136 // Finds the OAZ Dirichlet rows for this matrix
1137 // so far only used in IntrepidPCoarsenFactory
1138 // TODO check whether we can use DetectDirichletRows instead
1139 static void FindDirichletRows(Teuchos::RCP<Xpetra::Matrix<Scalar,LocalOrdinal,GlobalOrdinal,Node> > &A,
1140 std::vector<LocalOrdinal>& dirichletRows, bool count_twos_as_dirichlet=false) {
1141 typedef typename Teuchos::ScalarTraits<Scalar>::magnitudeType MT;
1142 dirichletRows.resize(0);
1143 for(size_t i=0; i<A->getLocalNumRows(); i++) {
1144 Teuchos::ArrayView<const LocalOrdinal> indices;
1145 Teuchos::ArrayView<const Scalar> values;
1146 A->getLocalRowView(i,indices,values);
1147 int nnz=0;
1148 for (size_t j=0; j<(size_t)indices.size(); j++) {
1149 if (Teuchos::ScalarTraits<Scalar>::magnitude(values[j]) > Teuchos::ScalarTraits<MT>::eps()) {
1150 nnz++;
1151 }
1152 }
1153 if (nnz == 1 || (count_twos_as_dirichlet && nnz == 2)) {
1154 dirichletRows.push_back(i);
1155 }
1156 }
1157 }
1158
1159 // Applies Ones-and-Zeros to matrix rows
1160 // Takes a vector of row indices
1161 static void ApplyOAZToMatrixRows(Teuchos::RCP<Xpetra::Matrix<Scalar,LocalOrdinal,GlobalOrdinal,Node> >& A,
1162 const std::vector<LocalOrdinal>& dirichletRows) {
1163 RCP<const Map> Rmap = A->getRowMap();
1164 RCP<const Map> Cmap = A->getColMap();
1165 Scalar one =Teuchos::ScalarTraits<Scalar>::one();
1166 Scalar zero =Teuchos::ScalarTraits<Scalar>::zero();
1167
1168 for(size_t i=0; i<dirichletRows.size(); i++) {
1169 GlobalOrdinal row_gid = Rmap->getGlobalElement(dirichletRows[i]);
1170
1171 Teuchos::ArrayView<const LocalOrdinal> indices;
1172 Teuchos::ArrayView<const Scalar> values;
1173 A->getLocalRowView(dirichletRows[i],indices,values);
1174 // NOTE: This won't work with fancy node types.
1175 Scalar* valuesNC = const_cast<Scalar*>(values.getRawPtr());
1176 for(size_t j=0; j<(size_t)indices.size(); j++) {
1177 if(Cmap->getGlobalElement(indices[j])==row_gid)
1178 valuesNC[j]=one;
1179 else
1180 valuesNC[j]=zero;
1181 }
1182 }
1183 }
1184
1185 // Applies Ones-and-Zeros to matrix rows
1186 // Takes a Boolean array.
1187 static void ApplyOAZToMatrixRows(Teuchos::RCP<Xpetra::Matrix<Scalar,LocalOrdinal,GlobalOrdinal,Node> >& A,
1188 const Teuchos::ArrayRCP<const bool>& dirichletRows) {
1189 TEUCHOS_ASSERT(A->isFillComplete());
1190 RCP<const Map> domMap = A->getDomainMap();
1191 RCP<const Map> ranMap = A->getRangeMap();
1192 RCP<const Map> Rmap = A->getRowMap();
1193 RCP<const Map> Cmap = A->getColMap();
1194 TEUCHOS_ASSERT(static_cast<size_t>(dirichletRows.size()) == Rmap->getLocalNumElements());
1195 const Scalar one = Teuchos::ScalarTraits<Scalar>::one();
1196 const Scalar zero = Teuchos::ScalarTraits<Scalar>::zero();
1197 A->resumeFill();
1198 for(size_t i=0; i<(size_t) dirichletRows.size(); i++) {
1199 if (dirichletRows[i]){
1200 GlobalOrdinal row_gid = Rmap->getGlobalElement(i);
1201
1202 Teuchos::ArrayView<const LocalOrdinal> indices;
1203 Teuchos::ArrayView<const Scalar> values;
1204 A->getLocalRowView(i,indices,values);
1205
1206 Teuchos::ArrayRCP<Scalar> valuesNC(values.size());
1207 for(size_t j=0; j<(size_t)indices.size(); j++) {
1208 if(Cmap->getGlobalElement(indices[j])==row_gid)
1209 valuesNC[j]=one;
1210 else
1211 valuesNC[j]=zero;
1212 }
1213 A->replaceLocalValues(i,indices,valuesNC());
1214 }
1215 }
1216 A->fillComplete(domMap, ranMap);
1217 }
1218
1219 // Zeros out rows
1220 // Takes a vector containg Dirichlet row indices
1221 static void ZeroDirichletRows(Teuchos::RCP<Xpetra::Matrix<Scalar,LocalOrdinal,GlobalOrdinal,Node> >& A,
1222 const std::vector<LocalOrdinal>& dirichletRows,
1223 Scalar replaceWith=Teuchos::ScalarTraits<Scalar>::zero()) {
1224 for(size_t i=0; i<dirichletRows.size(); i++) {
1225 Teuchos::ArrayView<const LocalOrdinal> indices;
1226 Teuchos::ArrayView<const Scalar> values;
1227 A->getLocalRowView(dirichletRows[i],indices,values);
1228 // NOTE: This won't work with fancy node types.
1229 Scalar* valuesNC = const_cast<Scalar*>(values.getRawPtr());
1230 for(size_t j=0; j<(size_t)indices.size(); j++)
1231 valuesNC[j]=replaceWith;
1232 }
1233 }
1234
1235 // Zeros out rows
1236 // Takes a Boolean ArrayRCP
1237 static void ZeroDirichletRows(Teuchos::RCP<Xpetra::Matrix<Scalar,LocalOrdinal,GlobalOrdinal,Node> >& A,
1238 const Teuchos::ArrayRCP<const bool>& dirichletRows,
1239 Scalar replaceWith=Teuchos::ScalarTraits<Scalar>::zero()) {
1240 TEUCHOS_ASSERT(static_cast<size_t>(dirichletRows.size()) == A->getRowMap()->getLocalNumElements());
1241 for(size_t i=0; i<(size_t) dirichletRows.size(); i++) {
1242 if (dirichletRows[i]) {
1243 Teuchos::ArrayView<const LocalOrdinal> indices;
1244 Teuchos::ArrayView<const Scalar> values;
1245 A->getLocalRowView(i,indices,values);
1246 // NOTE: This won't work with fancy node types.
1247 Scalar* valuesNC = const_cast<Scalar*>(values.getRawPtr());
1248 for(size_t j=0; j<(size_t)indices.size(); j++)
1249 valuesNC[j]=replaceWith;
1250 }
1251 }
1252 }
1253
1254 // Zeros out rows
1255 // Takes a Boolean ArrayRCP
1256 static void ZeroDirichletRows(Teuchos::RCP<Xpetra::MultiVector<Scalar,LocalOrdinal,GlobalOrdinal,Node> >& X,
1257 const Teuchos::ArrayRCP<const bool>& dirichletRows,
1258 Scalar replaceWith=Teuchos::ScalarTraits<Scalar>::zero()) {
1259 TEUCHOS_ASSERT(static_cast<size_t>(dirichletRows.size()) == X->getMap()->getLocalNumElements());
1260 for(size_t i=0; i<(size_t) dirichletRows.size(); i++) {
1261 if (dirichletRows[i]) {
1262 for(size_t j=0; j<X->getNumVectors(); j++)
1263 X->replaceLocalValue(i,j,replaceWith);
1264 }
1265 }
1266 }
1267
1268 // Zeros out columns
1269 // Takes a Boolean vector
1270 static void ZeroDirichletCols(Teuchos::RCP<Matrix>& A,
1271 const Teuchos::ArrayRCP<const bool>& dirichletCols,
1272 Scalar replaceWith=Teuchos::ScalarTraits<Scalar>::zero()) {
1273 TEUCHOS_ASSERT(static_cast<size_t>(dirichletCols.size()) == A->getColMap()->getLocalNumElements());
1274 for(size_t i=0; i<A->getLocalNumRows(); i++) {
1275 Teuchos::ArrayView<const LocalOrdinal> indices;
1276 Teuchos::ArrayView<const Scalar> values;
1277 A->getLocalRowView(i,indices,values);
1278 // NOTE: This won't work with fancy node types.
1279 Scalar* valuesNC = const_cast<Scalar*>(values.getRawPtr());
1280 for(size_t j=0; j<static_cast<size_t>(indices.size()); j++)
1281 if (dirichletCols[indices[j]])
1282 valuesNC[j] = replaceWith;
1283 }
1284 }
1285
1286 // Finds the OAZ Dirichlet rows for this matrix
1287 static void FindDirichletRowsAndPropagateToCols(Teuchos::RCP<Xpetra::Matrix<Scalar,LocalOrdinal,GlobalOrdinal,Node> > &A,
1288 Teuchos::RCP<Xpetra::Vector<int,LocalOrdinal,GlobalOrdinal,Node> >& isDirichletRow,
1289 Teuchos::RCP<Xpetra::Vector<int,LocalOrdinal,GlobalOrdinal,Node> >& isDirichletCol) {
1290
1291 // Make sure A's RowMap == DomainMap
1292 if(!A->getRowMap()->isSameAs(*A->getDomainMap())) {
1293 throw std::runtime_error("UtilitiesBase::FindDirichletRowsAndPropagateToCols row and domain maps must match.");
1294 }
1295 RCP<const Xpetra::Import<LocalOrdinal,GlobalOrdinal,Node> > importer = A->getCrsGraph()->getImporter();
1296 bool has_import = !importer.is_null();
1297
1298 // Find the Dirichlet rows
1299 std::vector<LocalOrdinal> dirichletRows;
1300 FindDirichletRows(A,dirichletRows);
1301
1302#if 0
1303 printf("[%d] DirichletRow Ids = ",A->getRowMap()->getComm()->getRank());
1304 for(size_t i=0; i<(size_t) dirichletRows.size(); i++)
1305 printf("%d ",dirichletRows[i]);
1306 printf("\n");
1307 fflush(stdout);
1308#endif
1309 // Allocate all as non-Dirichlet
1310 isDirichletRow = Xpetra::VectorFactory<int,LocalOrdinal,GlobalOrdinal,Node>::Build(A->getRowMap(),true);
1311 isDirichletCol = Xpetra::VectorFactory<int,LocalOrdinal,GlobalOrdinal,Node>::Build(A->getColMap(),true);
1312
1313 {
1314 Teuchos::ArrayRCP<int> dr_rcp = isDirichletRow->getDataNonConst(0);
1315 Teuchos::ArrayView<int> dr = dr_rcp();
1316 Teuchos::ArrayRCP<int> dc_rcp = isDirichletCol->getDataNonConst(0);
1317 Teuchos::ArrayView<int> dc = dc_rcp();
1318 for(size_t i=0; i<(size_t) dirichletRows.size(); i++) {
1319 dr[dirichletRows[i]] = 1;
1320 if(!has_import) dc[dirichletRows[i]] = 1;
1321 }
1322 }
1323
1324 if(has_import)
1325 isDirichletCol->doImport(*isDirichletRow,*importer,Xpetra::CombineMode::ADD);
1326
1327 }
1328
1329 // This routine takes a BlockedMap and an Importer (assuming that the BlockedMap matches the source of the importer) and generates a BlockedMap corresponding
1330 // to the Importer's target map. We assume that the targetMap is unique (which, is not a strict requirement of an Importer, but is here and no, we don't check)
1331 // This is largely intended to be used in repartitioning of blocked matrices
1332 static RCP<const Xpetra::BlockedMap<LocalOrdinal,GlobalOrdinal,Node> > GeneratedBlockedTargetMap(const Xpetra::BlockedMap<LocalOrdinal,GlobalOrdinal,Node> & sourceBlockedMap,
1333 const Xpetra::Import<LocalOrdinal,GlobalOrdinal,Node> & Importer) {
1334 typedef Xpetra::Vector<int,LocalOrdinal,GlobalOrdinal,Node> IntVector;
1335 Xpetra::UnderlyingLib lib = sourceBlockedMap.lib();
1336
1337 // De-stride the map if we have to (might regret this later)
1338 RCP<const Map> fullMap = sourceBlockedMap.getMap();
1339 RCP<const Map> stridedMap = Teuchos::rcp_dynamic_cast<const Xpetra::StridedMap<LocalOrdinal,GlobalOrdinal,Node> >(fullMap);
1340 if(!stridedMap.is_null()) fullMap = stridedMap->getMap();
1341
1342 // Initial sanity checking for map compatibil
1343 const size_t numSubMaps = sourceBlockedMap.getNumMaps();
1344 if(!Importer.getSourceMap()->isCompatible(*fullMap))
1345 throw std::runtime_error("GenerateBlockedTargetMap(): Map compatibility error");
1346
1347 // Build an indicator vector
1348 RCP<IntVector> block_ids = Xpetra::VectorFactory<int,LocalOrdinal,GlobalOrdinal,Node>::Build(fullMap);
1349
1350 for(size_t i=0; i<numSubMaps; i++) {
1351 RCP<const Map> map = sourceBlockedMap.getMap(i);
1352
1353 for(size_t j=0; j<map->getLocalNumElements(); j++) {
1354 LocalOrdinal jj = fullMap->getLocalElement(map->getGlobalElement(j));
1355 block_ids->replaceLocalValue(jj,(int)i);
1356 }
1357 }
1358
1359 // Get the block ids for the new map
1360 RCP<const Map> targetMap = Importer.getTargetMap();
1361 RCP<IntVector> new_block_ids = Xpetra::VectorFactory<int,LocalOrdinal,GlobalOrdinal,Node>::Build(targetMap);
1362 new_block_ids->doImport(*block_ids,Importer,Xpetra::CombineMode::ADD);
1363 Teuchos::ArrayRCP<const int> dataRCP = new_block_ids->getData(0);
1364 Teuchos::ArrayView<const int> data = dataRCP();
1365
1366
1367 // Get the GIDs for each subblock
1368 Teuchos::Array<Teuchos::Array<GlobalOrdinal> > elementsInSubMap(numSubMaps);
1369 for(size_t i=0; i<targetMap->getLocalNumElements(); i++) {
1370 elementsInSubMap[data[i]].push_back(targetMap->getGlobalElement(i));
1371 }
1372
1373 // Generate the new submaps
1374 std::vector<RCP<const Map> > subMaps(numSubMaps);
1375 for(size_t i=0; i<numSubMaps; i++) {
1376 subMaps[i] = Xpetra::MapFactory<LocalOrdinal,GlobalOrdinal,Node>::Build(lib,Teuchos::OrdinalTraits<GlobalOrdinal>::invalid(),elementsInSubMap[i](),targetMap->getIndexBase(),targetMap->getComm());
1377 }
1378
1379 // Build the BlockedMap
1380 return rcp(new BlockedMap(targetMap,subMaps));
1381
1382 }
1383
1384 // Checks to see if the first chunk of the colMap is also the row map. This simiplifies a bunch of
1385 // operation in coarsening
1386 static bool MapsAreNested(const Xpetra::Map<LocalOrdinal,GlobalOrdinal,Node>& rowMap, const Xpetra::Map<LocalOrdinal,GlobalOrdinal,Node>& colMap) {
1387 ArrayView<const GlobalOrdinal> rowElements = rowMap.getLocalElementList();
1388 ArrayView<const GlobalOrdinal> colElements = colMap.getLocalElementList();
1389
1390 const size_t numElements = rowElements.size();
1391
1392 if (size_t(colElements.size()) < numElements)
1393 return false;
1394
1395 bool goodMap = true;
1396 for (size_t i = 0; i < numElements; i++)
1397 if (rowElements[i] != colElements[i]) {
1398 goodMap = false;
1399 break;
1400 }
1401
1402 return goodMap;
1403 }
1404
1405
1406
1407
1408 }; // class Utils
1409
1410
1412
1413} //namespace MueLu
1414
1415#define MUELU_UTILITIESBASE_SHORT
1416#endif // MUELU_UTILITIESBASE_DECL_HPP
#define MueLu_sumAll(rcpComm, in, out)
MueLu::DefaultLocalOrdinal LocalOrdinal
MueLu::DefaultScalar Scalar
MueLu::DefaultGlobalOrdinal GlobalOrdinal
MueLu::DefaultNode Node
Exception throws to report incompatible objects (like maps).
Exception throws to report errors in the internal logical of the program.
Xpetra::MultiVector< Scalar, LocalOrdinal, GlobalOrdinal, Node > MultiVector
static RCP< Vector > GetMatrixOverlappedDeletedRowsum(const Matrix &A)
Extract Overlapped Matrix Deleted Rowsum.
Xpetra::Map< LocalOrdinal, GlobalOrdinal, Node > Map
static void ZeroDirichletRows(Teuchos::RCP< Xpetra::MultiVector< Scalar, LocalOrdinal, GlobalOrdinal, Node > > &X, const Teuchos::ArrayRCP< const bool > &dirichletRows, Scalar replaceWith=Teuchos::ScalarTraits< Scalar >::zero())
static RCP< CrsMatrixWrap > GetThresholdedMatrix(const RCP< Matrix > &Ain, const Scalar threshold, const bool keepDiagonal=true, const GlobalOrdinal expectedNNZperRow=-1)
Threshold a matrix.
static void ZeroDirichletCols(Teuchos::RCP< Matrix > &A, const Teuchos::ArrayRCP< const bool > &dirichletCols, Scalar replaceWith=Teuchos::ScalarTraits< Scalar >::zero())
static RCP< Teuchos::FancyOStream > MakeFancy(std::ostream &os)
static void DetectDirichletColsAndDomains(const Xpetra::Matrix< Scalar, LocalOrdinal, GlobalOrdinal, Node > &A, const Teuchos::ArrayRCP< bool > &dirichletRows, Teuchos::ArrayRCP< bool > dirichletCols, Teuchos::ArrayRCP< bool > dirichletDomain)
Detects Dirichlet columns & domains from a list of Dirichlet rows.
static Teuchos::ScalarTraits< Scalar >::magnitudeType Distance2(const Teuchos::Array< Teuchos::ArrayRCP< const Scalar > > &v, LocalOrdinal i0, LocalOrdinal i1)
Squared distance between two rows in a multivector.
Teuchos::ScalarTraits< Scalar >::magnitudeType Magnitude
Xpetra::CrsMatrix< Scalar, LocalOrdinal, GlobalOrdinal, Node > CrsMatrix
static Scalar PowerMethod(const Matrix &A, const RCP< Vector > &diagInvVec, LocalOrdinal niters=10, Magnitude tolerance=1e-2, bool verbose=false, unsigned int seed=123)
Power method.
static Teuchos::ArrayRCP< Scalar > GetMatrixDiagonal(const Matrix &A)
Extract Matrix Diagonal.
static void SetRandomSeed(const Teuchos::Comm< int > &comm)
Set seed for random number generator.
static void ApplyRowSumCriterion(const Xpetra::Matrix< Scalar, LocalOrdinal, GlobalOrdinal, Node > &A, const Xpetra::Vector< LocalOrdinal, LocalOrdinal, GlobalOrdinal, Node > &BlockNumber, const Magnitude rowSumTol, Teuchos::ArrayRCP< bool > &dirichletRows)
Xpetra::Vector< Scalar, LocalOrdinal, GlobalOrdinal, Node > Vector
static Teuchos::Array< Magnitude > ResidualNorm(const Xpetra::Operator< Scalar, LocalOrdinal, GlobalOrdinal, Node > &Op, const MultiVector &X, const MultiVector &RHS)
Xpetra::BlockedVector< Scalar, LocalOrdinal, GlobalOrdinal, Node > BlockedVector
static Teuchos::ArrayRCP< const bool > DetectDirichletRows(const Xpetra::Matrix< Scalar, LocalOrdinal, GlobalOrdinal, Node > &A, const Magnitude &tol=Teuchos::ScalarTraits< Scalar >::zero(), bool count_twos_as_dirichlet=false)
Detect Dirichlet rows.
static Teuchos::ArrayRCP< Magnitude > GetMatrixMaxMinusOffDiagonal(const Xpetra::Matrix< Scalar, LocalOrdinal, GlobalOrdinal, Node > &A)
Return vector containing: max_{i\not=k}(-a_ik), for each for i in the matrix.
static RCP< Vector > GetMatrixDiagonalInverse(const Matrix &A, Magnitude tol=Teuchos::ScalarTraits< Scalar >::eps() *100, Scalar valReplacement=Teuchos::ScalarTraits< Scalar >::zero())
Extract Matrix Diagonal.
static Scalar Frobenius(const Xpetra::Matrix< Scalar, LocalOrdinal, GlobalOrdinal, Node > &A, const Xpetra::Matrix< Scalar, LocalOrdinal, GlobalOrdinal, Node > &B)
Frobenius inner product of two matrices.
static Scalar PowerMethod(const Matrix &A, bool scaleByDiag=true, LocalOrdinal niters=10, Magnitude tolerance=1e-2, bool verbose=false, unsigned int seed=123)
Power method.
static void ApplyOAZToMatrixRows(Teuchos::RCP< Xpetra::Matrix< Scalar, LocalOrdinal, GlobalOrdinal, Node > > &A, const std::vector< LocalOrdinal > &dirichletRows)
static RCP< Matrix > Crs2Op(RCP< CrsMatrix > Op)
static Teuchos::ArrayRCP< const bool > DetectDirichletRowsExt(const Xpetra::Matrix< Scalar, LocalOrdinal, GlobalOrdinal, Node > &A, bool &bHasZeroDiagonal, const Magnitude &tol=Teuchos::ScalarTraits< Scalar >::zero())
Detect Dirichlet rows (extended version)
static void FindDirichletRows(Teuchos::RCP< Xpetra::Matrix< Scalar, LocalOrdinal, GlobalOrdinal, Node > > &A, std::vector< LocalOrdinal > &dirichletRows, bool count_twos_as_dirichlet=false)
static void ApplyRowSumCriterion(const Xpetra::Matrix< Scalar, LocalOrdinal, GlobalOrdinal, Node > &A, const Magnitude rowSumTol, Teuchos::ArrayRCP< bool > &dirichletRows)
Apply Rowsum Criterion.
static RCP< const Xpetra::BlockedMap< LocalOrdinal, GlobalOrdinal, Node > > GeneratedBlockedTargetMap(const Xpetra::BlockedMap< LocalOrdinal, GlobalOrdinal, Node > &sourceBlockedMap, const Xpetra::Import< LocalOrdinal, GlobalOrdinal, Node > &Importer)
static Teuchos::RCP< Vector > GetLumpedMatrixDiagonal(Matrix const &A, const bool doReciprocal=false, Magnitude tol=Teuchos::ScalarTraits< Scalar >::magnitude(Teuchos::ScalarTraits< Scalar >::zero()), Scalar valReplacement=Teuchos::ScalarTraits< Scalar >::zero(), const bool replaceSingleEntryRowWithZero=false, const bool useAverageAbsDiagVal=false)
Extract Matrix Diagonal of lumped matrix.
static void ZeroDirichletRows(Teuchos::RCP< Xpetra::Matrix< Scalar, LocalOrdinal, GlobalOrdinal, Node > > &A, const Teuchos::ArrayRCP< const bool > &dirichletRows, Scalar replaceWith=Teuchos::ScalarTraits< Scalar >::zero())
static void FindNonZeros(const Teuchos::ArrayRCP< const Scalar > vals, Teuchos::ArrayRCP< bool > nonzeros)
Find non-zero values in an ArrayRCP Compares the value to 2 * machine epsilon.
static bool MapsAreNested(const Xpetra::Map< LocalOrdinal, GlobalOrdinal, Node > &rowMap, const Xpetra::Map< LocalOrdinal, GlobalOrdinal, Node > &colMap)
Xpetra::CrsGraphFactory< LocalOrdinal, GlobalOrdinal, Node > CrsGraphFactory
static void Residual(const Xpetra::Operator< Scalar, LocalOrdinal, GlobalOrdinal, Node > &Op, const MultiVector &X, const MultiVector &RHS, MultiVector &Resid)
static Teuchos::ScalarTraits< Scalar >::magnitudeType Distance2(const Teuchos::ArrayView< double > &weight, const Teuchos::Array< Teuchos::ArrayRCP< const Scalar > > &v, LocalOrdinal i0, LocalOrdinal i1)
Weighted squared distance between two rows in a multivector.
Xpetra::CrsMatrixWrap< Scalar, LocalOrdinal, GlobalOrdinal, Node > CrsMatrixWrap
static Teuchos::ArrayRCP< const bool > DetectDirichletCols(const Xpetra::Matrix< Scalar, LocalOrdinal, GlobalOrdinal, Node > &A, const Teuchos::ArrayRCP< const bool > &dirichletRows)
Detect Dirichlet columns based on Dirichlet rows.
static RCP< Xpetra::Vector< Magnitude, LocalOrdinal, GlobalOrdinal, Node > > GetMatrixOverlappedAbsDeletedRowsum(const Matrix &A)
static Teuchos::Array< Magnitude > ResidualNorm(const Xpetra::Operator< Scalar, LocalOrdinal, GlobalOrdinal, Node > &Op, const MultiVector &X, const MultiVector &RHS, MultiVector &Resid)
static Teuchos::ArrayRCP< Magnitude > GetMatrixMaxMinusOffDiagonal(const Xpetra::Matrix< Scalar, LocalOrdinal, GlobalOrdinal, Node > &A, const Xpetra::Vector< LocalOrdinal, LocalOrdinal, GlobalOrdinal, Node > &BlockNumber)
Xpetra::BlockedMultiVector< Scalar, LocalOrdinal, GlobalOrdinal, Node > BlockedMultiVector
Xpetra::Matrix< Scalar, LocalOrdinal, GlobalOrdinal, Node > Matrix
Xpetra::CrsGraph< LocalOrdinal, GlobalOrdinal, Node > CrsGraph
static RCP< MultiVector > Residual(const Xpetra::Operator< Scalar, LocalOrdinal, GlobalOrdinal, Node > &Op, const MultiVector &X, const MultiVector &RHS)
Xpetra::BlockedMap< LocalOrdinal, GlobalOrdinal, Node > BlockedMap
static RCP< Xpetra::CrsGraph< LocalOrdinal, GlobalOrdinal, Node > > GetThresholdedGraph(const RCP< Matrix > &A, const Magnitude threshold, const GlobalOrdinal expectedNNZperRow=-1)
Threshold a graph.
static void ApplyOAZToMatrixRows(Teuchos::RCP< Xpetra::Matrix< Scalar, LocalOrdinal, GlobalOrdinal, Node > > &A, const Teuchos::ArrayRCP< const bool > &dirichletRows)
static RCP< Vector > GetMatrixOverlappedDiagonal(const Matrix &A)
Extract Overlapped Matrix Diagonal.
static Teuchos::RCP< Vector > GetInverse(Teuchos::RCP< const Vector > v, Magnitude tol=Teuchos::ScalarTraits< Scalar >::eps() *100, Scalar valReplacement=Teuchos::ScalarTraits< Scalar >::zero())
Return vector containing inverse of input vector.
static void FindDirichletRowsAndPropagateToCols(Teuchos::RCP< Xpetra::Matrix< Scalar, LocalOrdinal, GlobalOrdinal, Node > > &A, Teuchos::RCP< Xpetra::Vector< int, LocalOrdinal, GlobalOrdinal, Node > > &isDirichletRow, Teuchos::RCP< Xpetra::Vector< int, LocalOrdinal, GlobalOrdinal, Node > > &isDirichletCol)
static void ZeroDirichletRows(Teuchos::RCP< Xpetra::Matrix< Scalar, LocalOrdinal, GlobalOrdinal, Node > > &A, const std::vector< LocalOrdinal > &dirichletRows, Scalar replaceWith=Teuchos::ScalarTraits< Scalar >::zero())
Namespace for MueLu classes and methods.
KokkosClassic::DefaultNode::DefaultNodeType DefaultNode
magnitude_type tolerance