137 this->GetOStream(
Warnings0) <<
"MueLu::Amesos2Smoother::Setup(): Setup() has already been called" << std::endl;
139 RCP<Matrix> A = Factory::Get< RCP<Matrix> >(currentLevel,
"A");
142 RCP<const Map> rowMap = A->getRowMap();
144 Teuchos::ParameterList pL = this->GetParameterList();
145 if (pL.get<
bool>(
"fix nullspace")) {
146 this->GetOStream(
Runtime1) <<
"MueLu::Amesos2Smoother::Setup(): fixing nullspace" << std::endl;
148 rowMap = A->getRowMap();
149 size_t M = rowMap->getGlobalNumElements();
151 RCP<MultiVector> Nullspace = Factory::Get< RCP<MultiVector> >(currentLevel,
"Nullspace");
154 "MueLu::Amesos2Smoother::Setup Fixing nullspace for coarse matrix for Amesos2 for nullspace of dim > 1 has not been implemented yet.");
156 RCP<MultiVector> NullspaceImp;
157 RCP<const Map> colMap;
158 RCP<const Import> importer;
159 if (rowMap->getComm()->getSize() > 1) {
160 this->GetOStream(
Warnings0) <<
"MueLu::Amesos2Smoother::Setup(): Applying nullspace fix on distributed matrix. Try rebalancing to single rank!" << std::endl;
161 ArrayRCP<GO> elements_RCP;
162 elements_RCP.resize(M);
163 ArrayView<GO> elements = elements_RCP();
164 for (
size_t k = 0; k<M; k++)
165 elements[k] = Teuchos::as<GO>(k);
166 colMap = MapFactory::Build(rowMap->lib(),M*rowMap->getComm()->getSize(),elements,Teuchos::ScalarTraits<GO>::zero(),rowMap->getComm());
167 importer = ImportFactory::Build(rowMap,colMap);
168 NullspaceImp = MultiVectorFactory::Build(colMap, Nullspace->getNumVectors());
169 NullspaceImp->doImport(*Nullspace,*importer,Xpetra::INSERT);
171 NullspaceImp = Nullspace;
175 ArrayRCP<const SC> nullspaceRCP, nullspaceImpRCP;
176 RCP<CrsMatrixWrap> Acrs = rcp_dynamic_cast<CrsMatrixWrap>(A);
179 "MueLu::Amesos2Smoother::Setup Fixing nullspace for coarse matrix for Amesos2 when matrix is not a Crs matrix has not been implemented yet.");
181 ArrayRCP<const size_t> rowPointers;
182 ArrayRCP<const LO> colIndices;
183 ArrayRCP<const SC> values;
184 Acrs->getCrsMatrix()->getAllValues(rowPointers, colIndices, values);
186 ArrayRCP<size_t> newRowPointers_RCP;
187 ArrayRCP<LO> newColIndices_RCP;
188 ArrayRCP<SC> newValues_RCP;
190 size_t N = rowMap->getLocalNumElements();
191 newRowPointers_RCP.resize(N+1);
192 newColIndices_RCP.resize(N*M);
193 newValues_RCP.resize(N*M);
195 ArrayView<size_t> newRowPointers = newRowPointers_RCP();
196 ArrayView<LO> newColIndices = newColIndices_RCP();
197 ArrayView<SC> newValues = newValues_RCP();
199 SC normalization = Nullspace->getVector(0)->norm2();
200 normalization = Teuchos::ScalarTraits<Scalar>::one()/(normalization*normalization);
202 ArrayView<const SC> nullspace, nullspaceImp;
203 nullspaceRCP = Nullspace->getData(0);
204 nullspace = nullspaceRCP();
205 nullspaceImpRCP = NullspaceImp->getData(0);
206 nullspaceImp = nullspaceImpRCP();
209 for (
size_t i = 0; i < N; i++) {
210 newRowPointers[i] = i*M;
211 for (
size_t j = 0; j < M; j++) {
212 newColIndices[i*M+j] = Teuchos::as<LO>(j);
213 newValues[i*M+j] = normalization * nullspace[i]*Teuchos::ScalarTraits<Scalar>::conjugate(nullspaceImp[j]);
216 newRowPointers[N] = N*M;
219 for (
size_t i = 0; i < N; i++) {
220 for (
size_t jj = rowPointers[i]; jj < rowPointers[i+1]; jj++) {
221 LO j = colMap->getLocalElement(A->getColMap()->getGlobalElement(colIndices[jj]));
223 newValues[i*M+j] += v;
227 RCP<Matrix> newA = rcp(
new CrsMatrixWrap(rowMap, colMap, 0));
228 RCP<CrsMatrix> newAcrs = rcp_dynamic_cast<CrsMatrixWrap>(newA)->getCrsMatrix();
230 using Teuchos::arcp_const_cast;
231 newAcrs->setAllValues(arcp_const_cast<size_t>(newRowPointers_RCP), arcp_const_cast<LO>(newColIndices_RCP), arcp_const_cast<SC>(newValues_RCP));
232 newAcrs->expertStaticFillComplete(A->getDomainMap(), A->getRangeMap(),
233 importer, A->getCrsGraph()->getExporter());
242 prec_ = Amesos2::create<Tpetra_CrsMatrix,Tpetra_MultiVector>(type_, tA);
243 TEUCHOS_TEST_FOR_EXCEPTION(prec_ == Teuchos::null,
Exceptions::RuntimeError,
"Amesos2::create returns Teuchos::null");
244 RCP<Teuchos::ParameterList> amesos2_params = Teuchos::rcpFromRef(pL.sublist(
"Amesos2"));
245 amesos2_params->setName(
"Amesos2");
246 if (rowMap->getGlobalNumElements() != as<size_t>((rowMap->getMaxAllGlobalIndex() - rowMap->getMinAllGlobalIndex())+1)) {
247 if (!(amesos2_params->sublist(prec_->name()).template isType<bool>(
"IsContiguous")))
248 amesos2_params->sublist(prec_->name()).set(
"IsContiguous",
false,
"Are GIDs Contiguous");
250 prec_->setParameters(amesos2_params);
259 RCP<Tpetra_MultiVector> tX, tB;
260 if (!useTransformation_) {
265 size_t numVectors = X.getNumVectors();
266 size_t length = X.getLocalLength();
269 "MueLu::Amesos2Smoother::Apply: Fixing coarse matrix for Amesos2 for multivectors has not been implemented yet.");
270 ArrayRCP<const SC> Xdata = X. getData(0), Bdata = B. getData(0);
271 ArrayRCP<SC> X_data = X_->getDataNonConst(0), B_data = B_->getDataNonConst(0);
273 for (
size_t i = 0; i < length; i++) {
274 X_data[i] = Xdata[i];
275 B_data[i] = Bdata[i];
287 prec_->setX(Teuchos::null);
288 prec_->setB(Teuchos::null);
290 if (useTransformation_) {
292 size_t length = X.getLocalLength();
294 ArrayRCP<SC> Xdata = X. getDataNonConst(0);
295 ArrayRCP<const SC> X_data = X_->getData(0);
297 for (
size_t i = 0; i < length; i++)
298 Xdata[i] = X_data[i];