46 #ifndef MUELU_CONSTRAINT_DEF_HPP 47 #define MUELU_CONSTRAINT_DEF_HPP 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> 63 #include "MueLu_Utilities.hpp" 69 template<
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
71 const size_t NSDim = Bc.getNumVectors();
75 size_t numRows = Ppattern_->getNodeNumRows();
76 XXtInv_.resize(numRows);
78 RCP<const Import> importer = Ppattern_->getImporter();
80 X_ = MultiVectorFactory::Build(Ppattern_->getColMap(), NSDim);
81 if (!importer.is_null())
86 std::vector<const SC*> Xval(NSDim);
87 for (
size_t j = 0; j < NSDim; j++)
88 Xval[j] = X_->getData(j).get();
90 SC zero = Teuchos::ScalarTraits<SC>::zero();
91 SC one = Teuchos::ScalarTraits<SC>::one();
93 Teuchos::BLAS <LO,SC> blas;
94 Teuchos::LAPACK<LO,SC> lapack;
96 ArrayRCP<LO> IPIV(NSDim);
97 ArrayRCP<SC> WORK(lwork);
99 for (
size_t i = 0; i < numRows; i++) {
100 Teuchos::ArrayView<const LO> indices;
101 Ppattern_->getLocalRowView(i, indices);
103 size_t nnz = indices.size();
105 XXtInv_[i] = Teuchos::SerialDenseMatrix<LO,SC>(NSDim, NSDim,
false);
106 Teuchos::SerialDenseMatrix<LO,SC>& XXtInv = XXtInv_[i];
110 for (
size_t j = 0; j < nnz; j++)
111 d += Xval[0][indices[j]] * Xval[0][indices[j]];
115 Teuchos::SerialDenseMatrix<LO,SC> locX(NSDim, nnz,
false);
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]];
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());
128 lapack.GETRF(NSDim, NSDim, XXtInv.values(), XXtInv.stride(), IPIV.get(), &info);
130 lapack.GETRI(NSDim, XXtInv.values(), XXtInv.stride(), IPIV.get(), WORK.get(), lwork, &info);
136 template<
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
140 "Row maps are incompatible");
141 const size_t NSDim = X_->getNumVectors();
142 const size_t numRows = P.getNodeNumRows();
144 const Map& colMap = *P.getColMap();
145 const Map& PColMap = *Projected.getColMap();
147 Projected.resumeFill();
149 Teuchos::ArrayView<const LO> indices, pindices;
150 Teuchos::ArrayView<const SC> values, pvalues;
151 Teuchos::Array<SC> valuesAll(colMap.getNodeNumElements()), newValues;
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();
158 std::vector<const SC*> Xval(NSDim);
159 for (
size_t j = 0; j < NSDim; j++)
160 Xval[j] = X_->getData(j).get();
162 for (
size_t i = 0; i < numRows; i++) {
163 P .getLocalRowView(i, indices, values);
164 Projected.getLocalRowView(i, pindices, pvalues);
166 size_t nnz = indices.size();
167 size_t pnnz = pindices.size();
169 newValues.resize(pnnz);
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]));
185 newValues[j] = valuesAll[ind];
189 for (
size_t j = 0; j < nnz; j++)
190 valuesAll[indices[j]] = zero;
193 Teuchos::SerialDenseMatrix<LO,SC>& XXtInv = XXtInv_[i];
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]];
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];
204 Teuchos::BLAS<LO,SC> blas;
206 blas.GEMV(Teuchos::NO_TRANS, NSDim, pnnz,
207 one, locX.values(), locX.stride(),
209 zero, val1.values(), oneLO);
211 blas.GEMV(Teuchos::NO_TRANS, NSDim, NSDim,
212 one, XXtInv.values(), XXtInv.stride(),
213 val1.values(), oneLO,
214 zero, val2.values(), oneLO);
216 blas.GEMV(Teuchos::CONJ_TRANS, NSDim, pnnz,
217 one, locX.values(), locX.stride(),
218 val2.values(), oneLO,
219 zero, val.values(), oneLO);
221 for (
size_t j = 0; j < pnnz; j++)
222 newValues[j] -= val[j];
224 Projected.replaceLocalValues(i, pindices, newValues);
227 Projected.fillComplete(Projected.getDomainMap(), Projected.getRangeMap());
232 #endif //ifndef MUELU_CONSTRAINT_DEF_HPP
void Setup(const MultiVector &B, const MultiVector &Bc, RCP< const CrsGraph > Ppattern)
Namespace for MueLu classes and methods.
Exception throws to report incompatible objects (like maps).
void Apply(const Matrix &P, Matrix &Projected) const
Apply constraint.