MueLu  Version of the Day
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
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>
54 
55 #include <Xpetra_Import_fwd.hpp>
56 #include <Xpetra_ImportFactory.hpp>
57 #include <Xpetra_Map.hpp>
58 #include <Xpetra_Matrix.hpp>
60 #include <Xpetra_MultiVector.hpp>
61 #include <Xpetra_CrsGraph.hpp>
62 
63 #include "MueLu_Utilities.hpp"
65 
66 
67 namespace 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_->getNodeNumRows();
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 
92 
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++) {
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.getNodeNumRows();
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.getNodeNumElements()), newValues;
152 
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 
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
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)
size_type size() const
LocalOrdinal LO
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.
Scalar SC
void GETRI(const OrdinalType &n, ScalarType *A, const OrdinalType &lda, const OrdinalType *IPIV, ScalarType *WORK, const OrdinalType &lwork, OrdinalType *info) const
T * get() const
OrdinalType stride() const
bool is_null() const