Tpetra parallel linear algebra Version of the Day
Loading...
Searching...
No Matches
Tpetra_replaceDiagonalCrsMatrix_def.hpp
Go to the documentation of this file.
1// @HEADER
2// ***********************************************************************
3//
4// Tpetra: Templated Linear Algebra Services Package
5// Copyright (2008) Sandia Corporation
6//
7// Under the terms of Contract DE-AC04-94AL85000 with Sandia Corporation,
8// the U.S. Government retains certain rights in this software.
9//
10// Redistribution and use in source and binary forms, with or without
11// modification, are permitted provided that the following conditions are
12// met:
13//
14// 1. Redistributions of source code must retain the above copyright
15// notice, this list of conditions and the following disclaimer.
16//
17// 2. Redistributions in binary form must reproduce the above copyright
18// notice, this list of conditions and the following disclaimer in the
19// documentation and/or other materials provided with the distribution.
20//
21// 3. Neither the name of the Corporation nor the names of the
22// contributors may be used to endorse or promote products derived from
23// this software without specific prior written permission.
24//
25// THIS SOFTWARE IS PROVIDED BY SANDIA CORPORATION "AS IS" AND ANY
26// EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
27// IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR
28// PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL SANDIA CORPORATION OR THE
29// CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL,
30// EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO,
31// PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR
32// PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF
33// LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING
34// NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
35// SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
36//
37// Questions? Contact Michael A. Heroux (maherou@sandia.gov)
38//
39// ************************************************************************
40// @HEADER
41
42#ifndef TPETRA_REPLACEDIAGONALCRSMATRIX_DEF_HPP
43#define TPETRA_REPLACEDIAGONALCRSMATRIX_DEF_HPP
44
47
49#include "Tpetra_CrsMatrix.hpp"
50#include "Tpetra_Vector.hpp"
51
52namespace Tpetra {
53
54template<class SC, class LO, class GO, class NT>
55LO
56replaceDiagonalCrsMatrix (CrsMatrix<SC, LO, GO, NT>& matrix,
57 const Vector<SC, LO, GO, NT>& newDiag) {
58
59 using map_type = Map<LO, GO, NT>;
60 using crs_matrix_type = CrsMatrix<SC, LO, GO, NT>;
61 // using vec_type = ::Tpetra::Vector<SC, LO, GO, NT>;
62
63 const LO oneLO = Teuchos::OrdinalTraits<LO>::one();
64
65 // Count number of successfully replaced diagonal entries
66 LO numReplacedDiagEntries = 0;
67
68 // Extract some useful data
69 auto rowMapPtr = matrix.getRowMap();
70 if (rowMapPtr.is_null() || rowMapPtr->getComm().is_null()) {
71 // Processes on which the row Map or its communicator is null
72 // don't participate. Users shouldn't even call this method on
73 // those processes.
74 return numReplacedDiagEntries;
75 }
76 auto colMapPtr = matrix.getColMap();
77
78 TEUCHOS_TEST_FOR_EXCEPTION
79 (colMapPtr.get () == nullptr, std::invalid_argument,
80 "replaceDiagonalCrsMatrix: "
81 "Input matrix must have a nonnull column Map.");
82
83 const map_type& rowMap = *rowMapPtr;
84 const map_type& colMap = *colMapPtr;
85 const LO myNumRows = static_cast<LO>(matrix.getLocalNumRows());
86 const bool isFillCompleteOnInput = matrix.isFillComplete();
87
89 TEUCHOS_TEST_FOR_EXCEPTION(!rowMap.isSameAs(*newDiag.getMap()), std::runtime_error,
90 "Row map of matrix and map of input vector do not match.");
91 }
92
93 // KJ: This fence is necessary for UVM. Views used in the row map and colmap
94 // can use UVM and they are accessed in the following routine. So, we need to
95 // make sure that the values are available for touching in host.
96 typename crs_matrix_type::execution_space().fence();
97
98 if (isFillCompleteOnInput)
99 matrix.resumeFill();
100
101 Teuchos::ArrayRCP<const SC> newDiagData = newDiag.getVector(0)->getData();
102 LO numReplacedEntriesPerRow = 0;
103
104 auto invalid = Teuchos::OrdinalTraits<LO>::invalid();
105
106 // Loop over all local rows to replace the diagonal entry row by row
107 for (LO lclRowInd = 0; lclRowInd < myNumRows; ++lclRowInd) {
108
109 // Get row and column indices of the diagonal entry
110 const GO gblInd = rowMap.getGlobalElement(lclRowInd);
111 const LO lclColInd = colMap.getLocalElement(gblInd);
112
113 // If the row map is not one-to-one, the diagonal may not be on this proc.
114 // Skip this row; some processor will have the diagonal for this row.
115 if (lclColInd == invalid) continue;
116
117 const SC vals[] = {static_cast<SC>(newDiagData[lclRowInd])};
118 const LO cols[] = {lclColInd};
119
120 // Do the actual replacement of the diagonal element, if on this proc
121 numReplacedEntriesPerRow = matrix.replaceLocalValues(lclRowInd, oneLO,
122 vals, cols);
123
124 // Check for success of replacement.
125 // numReplacedEntriesPerRow is one if the diagonal was replaced.
126 // numReplacedEntriesPerRow is zero if the diagonal is not on
127 // this processor. For example, in a 2D matrix distribution, gblInd may
128 // be in both the row and column map, but the diagonal may not be on
129 // this processor.
130 if (numReplacedEntriesPerRow == oneLO) {
131 ++numReplacedDiagEntries;
132 }
133 }
134
135 if (isFillCompleteOnInput)
136 matrix.fillComplete();
137
138 return numReplacedDiagEntries;
139}
140
141} // namespace Tpetra
142//
143// Explicit instantiation macro
144//
145// Must be expanded from within the Tpetra namespace!
146//
147
148#define TPETRA_REPLACEDIAGONALCRSMATRIX_INSTANT(SCALAR,LO,GO,NODE) \
149 \
150 template \
151 LO replaceDiagonalCrsMatrix( \
152 CrsMatrix<SCALAR, LO, GO, NODE>& matrix, \
153 const Vector<SCALAR, LO, GO, NODE>& newDiag); \
154 \
155
156#endif // #ifndef TPETRA_REPLACEDIAGONALCRSMATRIX_DEF_HPP
Declaration of Tpetra::Details::Behavior, a class that describes Tpetra's behavior.
static bool debug()
Whether Tpetra is in debug mode.
Namespace Tpetra contains the class and methods constituting the Tpetra library.
LO replaceDiagonalCrsMatrix(::Tpetra::CrsMatrix< SC, LO, GO, NT > &matrix, const ::Tpetra::Vector< SC, LO, GO, NT > &newDiag)
Replace diagonal entries of an input Tpetra::CrsMatrix matrix with values given in newDiag.