71 const size_t NSDim = Bc.getNumVectors();
75 size_t numRows = Ppattern_->getLocalNumRows();
76 XXtInv_.resize(numRows);
78 RCP<const Import> importer = Ppattern_->getImporter();
80 X_ = MultiVectorFactory::Build(Ppattern_->getColMap(), NSDim);
81 if (!importer.is_null())
82 X_->doImport(Bc, *importer, Xpetra::INSERT);
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);
140 "Row maps are incompatible");
141 const size_t NSDim = X_->getNumVectors();
142 const size_t numRows = P.getLocalNumRows();
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.getLocalNumElements()), 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());