46 #ifndef MUELU_CONSTRAINT_DEF_HPP
47 #define MUELU_CONSTRAINT_DEF_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);
80 X_ = MultiVectorFactory::Build(Ppattern_->getColMap(), NSDim);
86 std::vector<const SC*> Xval(NSDim);
87 for (
size_t j = 0; j < NSDim; j++)
88 Xval[j] = X_->getData(j).get();
99 for (
size_t i = 0; i < numRows; i++) {
101 Ppattern_->getLocalRowView(i, indices);
103 size_t nnz = indices.
size();
110 for (
size_t j = 0; j < nnz; j++)
111 d += Xval[0][indices[j]] * Xval[0][indices[j]];
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]];
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();
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;
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]];
201 for (
size_t j = 0; j < pnnz; j++)
202 val[j] = newValues[j];
209 zero, val1.values(), oneLO);
213 val1.values(), oneLO,
214 zero, 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
ScalarType * values() const
void GEMV(ETransp trans, const OrdinalType &m, const OrdinalType &n, const alpha_type alpha, const A_type *A, const OrdinalType &lda, const x_type *x, const OrdinalType &incx, const beta_type beta, ScalarType *y, const OrdinalType &incy) const
#define TEUCHOS_TEST_FOR_EXCEPTION(throw_exception_test, Exception, msg)
void Setup(const MultiVector &B, const MultiVector &Bc, RCP< const CrsGraph > Ppattern)
void GEMM(ETransp transa, ETransp transb, const OrdinalType &m, const OrdinalType &n, const OrdinalType &k, const alpha_type alpha, const A_type *A, const OrdinalType &lda, const B_type *B, const OrdinalType &ldb, const beta_type beta, ScalarType *C, const OrdinalType &ldc) const
Exception throws to report incompatible objects (like maps).
void GETRF(const OrdinalType &m, const OrdinalType &n, ScalarType *A, const OrdinalType &lda, OrdinalType *IPIV, OrdinalType *info) const
void Apply(const Matrix &P, Matrix &Projected) const
Apply constraint.
void GETRI(const OrdinalType &n, ScalarType *A, const OrdinalType &lda, const OrdinalType *IPIV, ScalarType *WORK, const OrdinalType &lwork, OrdinalType *info) const
OrdinalType stride() const