MueLu Version of the Day
Loading...
Searching...
No Matches
MueLu_Constraint_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#ifndef MUELU_CONSTRAINT_DEF_HPP
47#define MUELU_CONSTRAINT_DEF_HPP
48
49#include <Teuchos_BLAS.hpp>
50#include <Teuchos_LAPACK.hpp>
51#include <Teuchos_SerialDenseVector.hpp>
52#include <Teuchos_SerialDenseMatrix.hpp>
53#include <Teuchos_SerialDenseHelpers.hpp>
54
55#include <Xpetra_Import_fwd.hpp>
56#include <Xpetra_ImportFactory.hpp>
57#include <Xpetra_Map.hpp>
58#include <Xpetra_Matrix.hpp>
59#include <Xpetra_MultiVectorFactory.hpp>
60#include <Xpetra_MultiVector.hpp>
61#include <Xpetra_CrsGraph.hpp>
62
63#include "MueLu_Utilities.hpp"
65
66
67namespace MueLu {
68
69 template<class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
70 void Constraint<Scalar, LocalOrdinal, GlobalOrdinal, Node>::Setup(const MultiVector& /* B */, const MultiVector& Bc, RCP<const CrsGraph> Ppattern) {
71 const size_t NSDim = Bc.getNumVectors();
72
73 Ppattern_ = Ppattern;
74
75 size_t numRows = Ppattern_->getLocalNumRows();
76 XXtInv_.resize(numRows);
77
78 RCP<const Import> importer = Ppattern_->getImporter();
79
80 X_ = MultiVectorFactory::Build(Ppattern_->getColMap(), NSDim);
81 if (!importer.is_null())
82 X_->doImport(Bc, *importer, Xpetra::INSERT);
83 else
84 *X_ = Bc;
85
86 std::vector<const SC*> Xval(NSDim);
87 for (size_t j = 0; j < NSDim; j++)
88 Xval[j] = X_->getData(j).get();
89
90 SC zero = Teuchos::ScalarTraits<SC>::zero();
91 SC one = Teuchos::ScalarTraits<SC>::one();
92
93 Teuchos::BLAS <LO,SC> blas;
94 Teuchos::LAPACK<LO,SC> lapack;
95 LO lwork = 3*NSDim;
96 ArrayRCP<LO> IPIV(NSDim);
97 ArrayRCP<SC> WORK(lwork);
98
99 for (size_t i = 0; i < numRows; i++) {
100 Teuchos::ArrayView<const LO> indices;
101 Ppattern_->getLocalRowView(i, indices);
102
103 size_t nnz = indices.size();
104
105 XXtInv_[i] = Teuchos::SerialDenseMatrix<LO,SC>(NSDim, NSDim, false/*zeroOut*/);
106 Teuchos::SerialDenseMatrix<LO,SC>& XXtInv = XXtInv_[i];
107
108 if (NSDim == 1) {
109 SC d = zero;
110 for (size_t j = 0; j < nnz; j++)
111 d += Xval[0][indices[j]] * Xval[0][indices[j]];
112 XXtInv(0,0) = one/d;
113
114 } else {
115 Teuchos::SerialDenseMatrix<LO,SC> locX(NSDim, nnz, false/*zeroOut*/);
116 for (size_t j = 0; j < nnz; j++)
117 for (size_t k = 0; k < NSDim; k++)
118 locX(k,j) = Xval[k][indices[j]];
119
120 // XXtInv_ = (locX*locX^T)^{-1}
121 blas.GEMM(Teuchos::NO_TRANS, Teuchos::CONJ_TRANS, NSDim, NSDim, nnz,
122 one, locX.values(), locX.stride(),
123 locX.values(), locX.stride(),
124 zero, XXtInv.values(), XXtInv.stride());
125
126 LO info;
127 // Compute LU factorization using partial pivoting with row exchanges
128 lapack.GETRF(NSDim, NSDim, XXtInv.values(), XXtInv.stride(), IPIV.get(), &info);
129 // Use the computed factorization to compute the inverse
130 lapack.GETRI(NSDim, XXtInv.values(), XXtInv.stride(), IPIV.get(), WORK.get(), lwork, &info);
131 }
132 }
133 }
134
136 template<class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
137 void Constraint<Scalar, LocalOrdinal, GlobalOrdinal, Node>::Apply(const Matrix& P, Matrix& Projected) const {
138 // We check only row maps. Column may be different.
139 TEUCHOS_TEST_FOR_EXCEPTION(!P.getRowMap()->isSameAs(*Projected.getRowMap()), Exceptions::Incompatible,
140 "Row maps are incompatible");
141 const size_t NSDim = X_->getNumVectors();
142 const size_t numRows = P.getLocalNumRows();
143
144 const Map& colMap = *P.getColMap();
145 const Map& PColMap = *Projected.getColMap();
146
147 Projected.resumeFill();
148
149 Teuchos::ArrayView<const LO> indices, pindices;
150 Teuchos::ArrayView<const SC> values, pvalues;
151 Teuchos::Array<SC> valuesAll(colMap.getLocalNumElements()), newValues;
152
153 LO invalid = Teuchos::OrdinalTraits<LO>::invalid();
154 LO oneLO = Teuchos::OrdinalTraits<LO>::one();
155 SC zero = Teuchos::ScalarTraits<SC> ::zero();
156 SC one = Teuchos::ScalarTraits<SC> ::one();
157
158 std::vector<const SC*> Xval(NSDim);
159 for (size_t j = 0; j < NSDim; j++)
160 Xval[j] = X_->getData(j).get();
161
162 for (size_t i = 0; i < numRows; i++) {
163 P .getLocalRowView(i, indices, values);
164 Projected.getLocalRowView(i, pindices, pvalues);
165
166 size_t nnz = indices.size(); // number of nonzeros in the supplied matrix
167 size_t pnnz = pindices.size(); // number of nonzeros in the constrained matrix
168
169 newValues.resize(pnnz);
170
171 // Step 1: fix stencil
172 // Projected *must* already have the correct stencil
173
174 // Step 2: copy correct stencil values
175 // The algorithm is very similar to the one used in the calculation of
176 // Frobenius dot product, see src/Transfers/Energy-Minimization/Solvers/MueLu_CGSolver_def.hpp
177
178 // NOTE: using extra array allows us to skip the search among indices
179 for (size_t j = 0; j < nnz; j++)
180 valuesAll[indices[j]] = values[j];
181 for (size_t j = 0; j < pnnz; j++) {
182 LO ind = colMap.getLocalElement(PColMap.getGlobalElement(pindices[j])); // FIXME: we could do that before the full loop just once
183 if (ind != invalid)
184 // index indices[j] is part of template, copy corresponding value
185 newValues[j] = valuesAll[ind];
186 else
187 newValues[j] = zero;
188 }
189 for (size_t j = 0; j < nnz; j++)
190 valuesAll[indices[j]] = zero;
191
192 // Step 3: project to the space
193 Teuchos::SerialDenseMatrix<LO,SC>& XXtInv = XXtInv_[i];
194
195 Teuchos::SerialDenseMatrix<LO,SC> locX(NSDim, pnnz, false);
196 for (size_t j = 0; j < pnnz; j++)
197 for (size_t k = 0; k < NSDim; k++)
198 locX(k,j) = Xval[k][pindices[j]];
199
200 Teuchos::SerialDenseVector<LO,SC> val(pnnz, false), val1(NSDim, false), val2(NSDim, false);
201 for (size_t j = 0; j < pnnz; j++)
202 val[j] = newValues[j];
203
204 Teuchos::BLAS<LO,SC> blas;
205 // val1 = locX * val;
206 blas.GEMV(Teuchos::NO_TRANS, NSDim, pnnz,
207 one, locX.values(), locX.stride(),
208 val.values(), oneLO,
209 zero, val1.values(), oneLO);
210 // val2 = XXtInv * val1
211 blas.GEMV(Teuchos::NO_TRANS, NSDim, NSDim,
212 one, XXtInv.values(), XXtInv.stride(),
213 val1.values(), oneLO,
214 zero, val2.values(), oneLO);
215 // val = X^T * val2
216 blas.GEMV(Teuchos::CONJ_TRANS, NSDim, pnnz,
217 one, locX.values(), locX.stride(),
218 val2.values(), oneLO,
219 zero, val.values(), oneLO);
220
221 for (size_t j = 0; j < pnnz; j++)
222 newValues[j] -= val[j];
223
224 Projected.replaceLocalValues(i, pindices, newValues);
225 }
226
227 Projected.fillComplete(Projected.getDomainMap(), Projected.getRangeMap()); //FIXME: maps needed?
228 }
229
230} // namespace MueLu
231
232#endif //ifndef MUELU_CONSTRAINT_DEF_HPP
void Setup(const MultiVector &B, const MultiVector &Bc, RCP< const CrsGraph > Ppattern)
void Apply(const Matrix &P, Matrix &Projected) const
Apply constraint.
Exception throws to report incompatible objects (like maps).
Namespace for MueLu classes and methods.