46#ifndef MUELU_CLASSICALPFACTORY_DEF_HPP
47#define MUELU_CLASSICALPFACTORY_DEF_HPP
49#include <Xpetra_MultiVectorFactory.hpp>
50#include <Xpetra_VectorFactory.hpp>
51#include <Xpetra_CrsGraphFactory.hpp>
52#include <Xpetra_Matrix.hpp>
53#include <Xpetra_Map.hpp>
54#include <Xpetra_Map.hpp>
55#include <Xpetra_MapFactory.hpp>
56#include <Xpetra_Vector.hpp>
57#include <Xpetra_Import.hpp>
58#include <Xpetra_ImportUtils.hpp>
59#include <Xpetra_IO.hpp>
60#include <Xpetra_StridedMapFactory.hpp>
62#include <Teuchos_OrdinalTraits.hpp>
66#include "MueLu_PerfUtils.hpp"
68#include "MueLu_ClassicalMapFactory.hpp"
69#include "MueLu_Utilities.hpp"
70#include "MueLu_AmalgamationInfo.hpp"
81 using STS =
typename Teuchos::ScalarTraits<SC>;
82 typename STS::magnitudeType MT_ZERO = Teuchos::ScalarTraits<typename STS::magnitudeType>::zero();
83 if(STS::real(val) > MT_ZERO)
return 1;
84 else if(STS::real(val) < MT_ZERO)
return -1;
95 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
97 RCP<ParameterList> validParamList = rcp(
new ParameterList());
98#define SET_VALID_ENTRY(name) validParamList->setEntry(name, MasterList::getEntry(name))
106 typedef Teuchos::StringToIntegralParameterEntryValidator<int> validatorType;
107 validParamList->getEntry(
"aggregation: classical scheme").setValidator(rcp(
new validatorType(Teuchos::tuple<std::string>(
"direct",
"ext+i",
"classical modified"),
"aggregation: classical scheme")));
111#undef SET_VALID_ENTRY
112 validParamList->set< RCP<const FactoryBase> >(
"A", Teuchos::null,
"Generating factory of the matrix A");
114 validParamList->set< RCP<const FactoryBase> >(
"Graph", null,
"Generating factory of the graph");
115 validParamList->set< RCP<const FactoryBase> >(
"DofsPerNode", null,
"Generating factory for variable \'DofsPerNode\', usually the same as for \'Graph\'");
116 validParamList->set< RCP<const FactoryBase> >(
"CoarseMap", Teuchos::null,
"Generating factory of the CoarseMap");
117 validParamList->set< RCP<const FactoryBase> >(
"FC Splitting", Teuchos::null,
"Generating factory of the FC Splitting");
118 validParamList->set< RCP<const FactoryBase> >(
"BlockNumber", Teuchos::null,
"Generating factory for Block Number");
121 return validParamList;
124 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
126 Input(fineLevel,
"A");
127 Input(fineLevel,
"Graph");
128 Input(fineLevel,
"DofsPerNode");
130 Input(fineLevel,
"CoarseMap");
131 Input(fineLevel,
"FC Splitting");
133 const ParameterList& pL = GetParameterList();
134 std::string drop_algo = pL.get<std::string>(
"aggregation: drop scheme");
135 if (drop_algo.find(
"block diagonal") != std::string::npos || drop_algo ==
"signed classical") {
136 Input(fineLevel,
"BlockNumber");
141 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
143 return BuildP(fineLevel, coarseLevel);
146 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
156 RCP<const Matrix> A = Get< RCP<Matrix> >(fineLevel,
"A");
157 RCP<const Map> ownedCoarseMap = Get<RCP<const Map> >(fineLevel,
"CoarseMap");
158 RCP<const LocalOrdinalVector> owned_fc_splitting = Get<RCP<LocalOrdinalVector> >(fineLevel,
"FC Splitting");
159 RCP<const GraphBase> graph = Get< RCP<GraphBase> >(fineLevel,
"Graph");
162 RCP<const Import> Importer = A->getCrsGraph()->getImporter();
163 Xpetra::UnderlyingLib lib = ownedCoarseMap->lib();
168 LO LO_INVALID = Teuchos::OrdinalTraits<LO>::invalid();
171 const ParameterList& pL = GetParameterList();
174 TEUCHOS_TEST_FOR_EXCEPTION(A->GetFixedBlockSize() != 1,
Exceptions::RuntimeError,
"ClassicalPFactory: Multiple PDEs per node not supported yet");
177 TEUCHOS_TEST_FOR_EXCEPTION(!A->getRowMap()->isSameAs(*A->getDomainMap()),
Exceptions::RuntimeError,
"ClassicalPFactory: MPI Ranks > 1 not supported yet");
180 std::string scheme = pL.get<std::string>(
"aggregation: classical scheme");
181 bool need_ghost_rows =
false;
182 if(scheme ==
"ext+i")
183 need_ghost_rows=
true;
184 else if(scheme ==
"direct")
185 need_ghost_rows=
false;
186 else if(scheme ==
"classical modified")
187 need_ghost_rows=
true;
191 GetOStream(
Statistics1) <<
"ClassicalPFactory: scheme = "<<scheme<<std::endl;
195 RCP<const LocalOrdinalVector> fc_splitting;
196 ArrayRCP<const LO> myPointType;
197 if(Importer.is_null()) {
198 fc_splitting = owned_fc_splitting;
201 RCP<LocalOrdinalVector> fc_splitting_nonconst = LocalOrdinalVectorFactory::Build(A->getCrsGraph()->getColMap());
202 fc_splitting_nonconst->doImport(*owned_fc_splitting,*Importer,Xpetra::INSERT);
203 fc_splitting = fc_splitting_nonconst;
205 myPointType = fc_splitting->getData(0);
209 RCP<const Matrix> Aghost;
210 RCP<const LocalOrdinalVector> fc_splitting_ghost;
211 ArrayRCP<const LO> myPointType_ghost;
212 RCP<const Import> remoteOnlyImporter;
213 if(need_ghost_rows && !Importer.is_null()){
214 ArrayView<const LO> remoteLIDs = Importer->getRemoteLIDs();
215 size_t numRemote = Importer->getNumRemoteIDs();
216 Array<GO> remoteRows(numRemote);
217 for (
size_t i = 0; i < numRemote; i++)
218 remoteRows[i] = Importer->getTargetMap()->getGlobalElement(remoteLIDs[i]);
220 RCP<const Map> remoteRowMap = MapFactory::Build(lib,Teuchos::OrdinalTraits<Xpetra::global_size_t>::invalid(), remoteRows(),
221 A->getDomainMap()->getIndexBase(), A->getDomainMap()->getComm());
223 remoteOnlyImporter = Importer->createRemoteOnlyImport(remoteRowMap);
224 RCP<const CrsMatrix> Acrs = rcp_dynamic_cast<const CrsMatrixWrap>(A)->getCrsMatrix();
225 RCP<CrsMatrix> Aghost_crs = CrsMatrixFactory::Build(Acrs,*remoteOnlyImporter,A->getDomainMap(),remoteOnlyImporter->getTargetMap());
226 Aghost = rcp(
new CrsMatrixWrap(Aghost_crs));
239 RCP<const Map> coarseMap;
240 if(Importer.is_null())
241 coarseMap = ownedCoarseMap;
244 GhostCoarseMap(*A,*Importer,myPointType,ownedCoarseMap,coarseMap);
248 RCP<LocalOrdinalVector> BlockNumber;
249 std::string drop_algo = pL.get<std::string>(
"aggregation: drop scheme");
250 if (drop_algo.find(
"block diagonal") != std::string::npos || drop_algo ==
"signed classical") {
251 RCP<LocalOrdinalVector> OwnedBlockNumber;
252 OwnedBlockNumber = Get<RCP<LocalOrdinalVector> >(fineLevel,
"BlockNumber");
253 if(Importer.is_null())
254 BlockNumber = OwnedBlockNumber;
256 BlockNumber = LocalOrdinalVectorFactory::Build(A->getColMap());
257 BlockNumber->doImport(*OwnedBlockNumber,*Importer,Xpetra::INSERT);
261#if defined(CMS_DEBUG) || defined(CMS_DUMP)
263 RCP<const CrsMatrix> Acrs = rcp_dynamic_cast<const CrsMatrixWrap>(A)->getCrsMatrix();
264 int rank = A->getRowMap()->getComm()->getRank();
265 printf(
"[%d] A local size = %dx%d\n",rank,(
int)Acrs->getRowMap()->getLocalNumElements(),(
int)Acrs->getColMap()->getLocalNumElements());
267 printf(
"[%d] graph local size = %dx%d\n",rank,(
int)graph->GetDomainMap()->getLocalNumElements(),(
int)graph->GetImportMap()->getLocalNumElements());
269 std::ofstream ofs(std::string(
"dropped_graph_") + std::to_string(fineLevel.
GetLevelID())+std::string(
"_") + std::to_string(rank) + std::string(
".dat"),std::ofstream::out);
270 RCP<Teuchos::FancyOStream> fancy = Teuchos::fancyOStream(Teuchos::rcpFromRef(ofs));
271 graph->print(*fancy,
Debug);
272 std::string out_fc = std::string(
"fc_splitting_") + std::to_string(fineLevel.
GetLevelID()) +std::string(
"_") + std::to_string(rank)+ std::string(
".dat");
273 std::string out_block = std::string(
"block_number_") + std::to_string(fineLevel.
GetLevelID()) +std::string(
"_") + std::to_string(rank)+ std::string(
".dat");
276 using real_type =
typename Teuchos::ScalarTraits<SC>::magnitudeType;
277 using RealValuedMultiVector =
typename Xpetra::MultiVector<real_type,LO,GO,NO>;
278 typedef Xpetra::MultiVectorFactory<real_type,LO,GO,NO> RealValuedMultiVectorFactory;
280 RCP<RealValuedMultiVector> mv = RealValuedMultiVectorFactory::Build(fc_splitting->getMap(),1);
281 ArrayRCP<real_type> mv_data= mv->getDataNonConst(0);
284 ArrayRCP<const LO> fc_data= fc_splitting->getData(0);
285 for(LO i=0; i<(LO)fc_data.size(); i++)
286 mv_data[i] = Teuchos::as<real_type>(fc_data[i]);
287 Xpetra::IO<real_type,LO,GO,NO>::Write(out_fc,*mv);
290 if(!BlockNumber.is_null()) {
291 RCP<RealValuedMultiVector> mv2 = RealValuedMultiVectorFactory::Build(BlockNumber->getMap(),1);
292 ArrayRCP<real_type> mv_data2= mv2->getDataNonConst(0);
293 ArrayRCP<const LO> b_data= BlockNumber->getData(0);
294 for(LO i=0; i<(LO)b_data.size(); i++) {
295 mv_data2[i] = Teuchos::as<real_type>(b_data[i]);
297 Xpetra::IO<real_type,LO,GO,NO>::Write(out_block,*mv2);
311 Array<LO> cpoint2pcol(myPointType.size(),LO_INVALID);
312 Array<LO> pcol2cpoint(coarseMap->getLocalNumElements(),LO_INVALID);
315 for(LO i=0; i<(LO) myPointType.size(); i++) {
316 if(myPointType[i] == C_PT) {
317 cpoint2pcol[i] = num_c_points;
320 else if (myPointType[i] == F_PT)
323 for(LO i=0; i<(LO)cpoint2pcol.size(); i++) {
324 if(cpoint2pcol[i] != LO_INVALID)
325 pcol2cpoint[cpoint2pcol[i]] =i;
330 Teuchos::Array<size_t> eis_rowptr;
331 Teuchos::Array<bool> edgeIsStrong;
334 GenerateStrengthFlags(*A,*graph,eis_rowptr,edgeIsStrong);
338 RCP<const Map> coarseColMap = coarseMap;
339 RCP<const Map> coarseDomainMap = ownedCoarseMap;
340 if(scheme ==
"ext+i") {
342 Coarsen_Ext_Plus_I(*A,Aghost,*graph,coarseColMap,coarseDomainMap,num_c_points,num_f_points,myPointType(),myPointType_ghost(),cpoint2pcol,pcol2cpoint,eis_rowptr,edgeIsStrong,BlockNumber,P);
344 else if(scheme ==
"direct") {
346 Coarsen_Direct(*A,Aghost,*graph,coarseColMap,coarseDomainMap,num_c_points,num_f_points,myPointType(),myPointType_ghost(),cpoint2pcol,pcol2cpoint,eis_rowptr,edgeIsStrong,BlockNumber,P);
348 else if(scheme ==
"classical modified") {
350 Coarsen_ClassicalModified(*A,Aghost,*graph,coarseColMap,coarseDomainMap,num_c_points,num_f_points,myPointType(),myPointType_ghost(),cpoint2pcol,pcol2cpoint,eis_rowptr,edgeIsStrong,BlockNumber,remoteOnlyImporter,P);
355 Xpetra::IO<SC,LO,GO,NO>::Write(
"classical_p.mat."+std::to_string(coarseLevel.
GetLevelID()), *P);
360 Set(coarseLevel,
"P",P);
361 RCP<const CrsGraph> pg = P->getCrsGraph();
362 Set(coarseLevel,
"P Graph",pg);
369 RCP<ParameterList> params = rcp(
new ParameterList());
370 params->set(
"printLoadBalancingInfo",
true);
375template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
377Coarsen_ClassicalModified(
const Matrix & A,
const RCP<const Matrix> & Aghost,
const GraphBase & graph, RCP<const Map> & coarseColMap, RCP<const Map> & coarseDomainMap, LO num_c_points, LO num_f_points,
const Teuchos::ArrayView<const LO> & myPointType,
const Teuchos::ArrayView<const LO> & myPointType_ghost,
const Teuchos::Array<LO> & cpoint2pcol,
const Teuchos::Array<LO> & pcol2cpoint, Teuchos::Array<size_t> & eis_rowptr, Teuchos::Array<bool> & edgeIsStrong, RCP<LocalOrdinalVector> & BlockNumber, RCP<const Import> remoteOnlyImporter,RCP<Matrix> & P)
const {
417 using STS =
typename Teuchos::ScalarTraits<SC>;
418 LO LO_INVALID = Teuchos::OrdinalTraits<LO>::invalid();
420 SC SC_ZERO = STS::zero();
422 int rank = A.getRowMap()->getComm()->getRank();
426 ArrayRCP<const LO> block_id;
427 if(!BlockNumber.is_null())
428 block_id = BlockNumber->getData(0);
433 size_t Nrows = A.getLocalNumRows();
434 double c_point_density = (double)num_c_points / (num_c_points+num_f_points);
437 double nnz_per_row_est = c_point_density*mean_strong_neighbors_per_row;
439 size_t nnz_est = std::max(Nrows,std::min((
size_t)A.getLocalNumEntries(),(
size_t)(nnz_per_row_est*Nrows)));
440 Array<size_t> tmp_rowptr(Nrows+1);
441 Array<LO> tmp_colind(nnz_est);
447 for(LO row=0; row < (LO) Nrows; row++) {
448 size_t row_start = eis_rowptr[row];
449 ArrayView<const LO> indices;
450 ArrayView<const SC> vals;
452 if(myPointType[row] == DIRICHLET_PT) {
455 else if(myPointType[row] == C_PT) {
457 C_hat.insert(cpoint2pcol[row]);
463 A.getLocalRowView(row, indices, vals);
464 for(LO j=0; j<indices.size(); j++)
465 if(myPointType[indices[j]] == C_PT && edgeIsStrong[row_start + j])
466 C_hat.insert(cpoint2pcol[indices[j]]);
470 if(ct + (
size_t)C_hat.size() > (size_t)tmp_colind.size()) {
471 tmp_colind.resize(std::max(ct+(
size_t)C_hat.size(),(size_t)2*tmp_colind.size()));
475 std::copy(C_hat.begin(), C_hat.end(),tmp_colind.begin()+ct);
477 tmp_rowptr[row+1] = tmp_rowptr[row] + C_hat.size();
480 tmp_colind.resize(tmp_rowptr[Nrows]);
482 RCP<CrsMatrix> Pcrs = CrsMatrixFactory::Build(A.getRowMap(),coarseColMap,0);
483 ArrayRCP<size_t> P_rowptr;
484 ArrayRCP<LO> P_colind;
485 ArrayRCP<SC> P_values;
487 if (A.getRowMap()->lib() == Xpetra::UseTpetra) {
488 P_rowptr = Teuchos::arcpFromArray(tmp_rowptr);
489 P_colind = Teuchos::arcpFromArray(tmp_colind);
490 P_values.resize(P_rowptr[Nrows]);
495 Pcrs->allocateAllValues(tmp_rowptr[Nrows],P_rowptr,P_colind,P_values);
496 std::copy(tmp_rowptr.begin(),tmp_rowptr.end(), P_rowptr.begin());
497 std::copy(tmp_colind.begin(),tmp_colind.end(), P_colind.begin());
498 Pcrs->setAllValues(P_rowptr,P_colind,P_values);
499 Pcrs->expertStaticFillComplete(coarseDomainMap, A.getDomainMap());
503 RCP<const CrsGraph> Pgraph;
504 RCP<const CrsGraph> Pghost;
507 ArrayRCP<const size_t> Pghost_rowptr;
508 ArrayRCP<const LO> Pghost_colind;
509 if(!remoteOnlyImporter.is_null()) {
510 if (A.getRowMap()->lib() == Xpetra::UseTpetra) {
511 RCP<CrsGraph> tempPgraph = CrsGraphFactory::Build(A.getRowMap(),coarseColMap,P_rowptr,P_colind);
512 tempPgraph->fillComplete(coarseDomainMap,A.getDomainMap());
516 Pgraph = Pcrs->getCrsGraph();
518 TEUCHOS_ASSERT(!Pgraph.is_null());
519 Pghost = CrsGraphFactory::Build(Pgraph,*remoteOnlyImporter,Pgraph->getDomainMap(),remoteOnlyImporter->getTargetMap());
520 Pghost->getAllIndices(Pghost_rowptr,Pghost_colind);
524 ArrayRCP<LO> Acol_to_Pcol(A.getColMap()->getLocalNumElements(),LO_INVALID);
527 ArrayRCP<LO> Pghostcol_to_Pcol;
528 if(!Pghost.is_null()) {
529 Pghostcol_to_Pcol.resize(Pghost->getColMap()->getLocalNumElements(),LO_INVALID);
530 for(LO i=0; i<(LO) Pghost->getColMap()->getLocalNumElements(); i++)
531 Pghostcol_to_Pcol[i] = Pgraph->getColMap()->getLocalElement(Pghost->getColMap()->getGlobalElement(i));
535 ArrayRCP<LO> Aghostcol_to_Acol;
536 if(!Aghost.is_null()) {
537 Aghostcol_to_Acol.resize(Aghost->getColMap()->getLocalNumElements(),LO_INVALID);
538 for(LO i=0; i<(LO)Aghost->getColMap()->getLocalNumElements(); i++)
539 Aghostcol_to_Acol[i] = A.getColMap()->getLocalElement(Aghost->getColMap()->getGlobalElement(i));
544 for(LO i=0; i < (LO)Nrows; i++) {
545 if(myPointType[i] == DIRICHLET_PT) {
549 printf(
"[%d] ** A(%d,:) is a Dirichlet-Point.\n",rank,i);
552 else if (myPointType[i] == C_PT) {
554 P_values[P_rowptr[i]] = Teuchos::ScalarTraits<SC>::one();
557 printf(
"[%d] ** A(%d,:) is a C-Point.\n",rank,i);
564 printf(
"[%d] ** A(%d,:) is a F-Point.\n",rank,i);
568 ArrayView<const LO> A_indices_i, A_indices_k;
569 ArrayView<const SC> A_vals_i, A_vals_k;
570 A.getLocalRowView(i, A_indices_i, A_vals_i);
571 size_t row_start = eis_rowptr[i];
573 ArrayView<const LO> P_indices_i = P_colind.view(P_rowptr[i],P_rowptr[i+1] - P_rowptr[i]);
574 ArrayView<SC> P_vals_i = P_values.view(P_rowptr[i],P_rowptr[i+1] - P_rowptr[i]);
577 for(LO j=0; j<(LO)P_vals_i.size(); j++)
578 P_vals_i[j] = SC_ZERO;
582 for(LO j=0; j<(LO)P_indices_i.size(); j++) {
583 Acol_to_Pcol[pcol2cpoint[P_indices_i[j]]] = P_rowptr[i] + j;
587 SC first_denominator = SC_ZERO;
591 for(LO k0=0; k0<(LO)A_indices_i.size(); k0++) {
592 LO k = A_indices_i[k0];
593 SC a_ik = A_vals_i[k0];
594 LO pcol_k = Acol_to_Pcol[k];
599 first_denominator += a_ik;
602 printf(
"- A(%d,%d) is the diagonal\n",i,k);
606 else if(myPointType[k] == DIRICHLET_PT) {
609 printf(
"- A(%d,%d) is a Dirichlet point\n",i,k);
613 else if(pcol_k != LO_INVALID && pcol_k >= (LO)P_rowptr[i]) {
615 P_values[pcol_k] += a_ik;
617 printf(
"- A(%d,%d) is a strong C-Point\n",i,k);
620 else if (!edgeIsStrong[row_start + k0]) {
622 if(block_id.size() == 0 || block_id[i] == block_id[k]) {
623 first_denominator += a_ik;
625 printf(
"- A(%d,%d) is weak adding to diagonal(%d,%d) (%6.4e)\n",i,k,block_id[i],block_id[k],a_ik);
628 printf(
"- A(%d,%d) is weak but does not match blocks (%d,%d), discarding\n",i,k,block_id[i],block_id[k]);
637 printf(
"- A(%d,%d) is a strong F-Point\n",i,k);
641 SC second_denominator = SC_ZERO;
646 A.getLocalRowView(k, A_indices_k, A_vals_k);
647 for(LO m0=0; m0<(LO)A_indices_k.size(); m0++){
648 LO m = A_indices_k[m0];
656 sign_akk = Sign(a_kk);
657 for(LO m0=0; m0<(LO)A_indices_k.size(); m0++){
658 LO m = A_indices_k[m0];
659 if(m != k && Acol_to_Pcol[A_indices_k[m0]] >= (LO)P_rowptr[i]) {
660 SC a_km = A_vals_k[m0];
661 second_denominator+= (Sign(a_km) == sign_akk ? SC_ZERO : a_km);
667 if(second_denominator != SC_ZERO) {
668 for(LO j0=0; j0<(LO)A_indices_k.size(); j0++) {
669 LO j = A_indices_k[j0];
672 if(Acol_to_Pcol[j] >= (LO)P_rowptr[i]) {
673 SC a_kj = A_vals_k[j0];
674 SC sign_akj_val = sign_akk == Sign(a_kj) ? SC_ZERO : a_kj;
675 P_values[Acol_to_Pcol[j]] += a_ik * sign_akj_val / second_denominator;
677 printf(
"- - Unscaled P(%d,A-%d) += %6.4e = %5.4e\n",i,j,a_ik * sign_akj_val / second_denominator,P_values[Acol_to_Pcol[j]]);
684 printf(
"- - A(%d,%d) second denominator is zero\n",i,k);
686 if(block_id.size() == 0 || block_id[i] == block_id[k])
687 first_denominator += a_ik;
696 Aghost->getLocalRowView(kless, A_indices_k, A_vals_k);
697 GO k_g = Aghost->getRowMap()->getGlobalElement(kless);
698 for(LO m0=0; m0<(LO)A_indices_k.size(); m0++){
699 GO m_g = Aghost->getColMap()->getGlobalElement(A_indices_k[m0]);
707 sign_akk = Sign(a_kk);
708 for(LO m0=0; m0<(LO)A_indices_k.size(); m0++){
709 GO m_g = Aghost->getColMap()->getGlobalElement(A_indices_k[m0]);
710 LO mghost = A_indices_k[m0];
711 LO m = Aghostcol_to_Acol[mghost];
712 if(m_g != k_g && m != LO_INVALID && Acol_to_Pcol[m] >= (LO)P_rowptr[i]) {
713 SC a_km = A_vals_k[m0];
714 second_denominator+= (Sign(a_km) == sign_akk ? SC_ZERO : a_km);
721 if(second_denominator != SC_ZERO) {
722 for(LO j0=0; j0<(LO)A_indices_k.size(); j0++){
723 LO jghost = A_indices_k[j0];
724 LO j = Aghostcol_to_Acol[jghost];
726 if((j != LO_INVALID) && (Acol_to_Pcol[j] >= (LO)P_rowptr[i])) {
727 SC a_kj = A_vals_k[j0];
728 SC sign_akj_val = sign_akk == Sign(a_kj) ? SC_ZERO : a_kj;
729 P_values[Acol_to_Pcol[j]] += a_ik * sign_akj_val / second_denominator;
731 printf(
"- - Unscaled P(%d,A-%d) += %6.4e\n",i,j,a_ik * sign_akj_val / second_denominator);
739 printf(
"- - A(%d,%d) second denominator is zero\n",i,k);
741 if(block_id.size() == 0 || block_id[i] == block_id[k])
742 first_denominator += a_ik;
750 if(first_denominator != SC_ZERO) {
751 for(LO j0=0; j0<(LO)P_indices_i.size(); j0++) {
753 SC old_pij = P_vals_i[j0];
754 P_vals_i[j0] /= -first_denominator;
755 printf(
"P(%d,%d) = %6.4e = %6.4e / (%6.4e + %6.4e)\n",i,P_indices_i[j0],P_vals_i[j0],old_pij,a_ii,first_denominator - a_ii);
757 P_vals_i[j0] /= -first_denominator;
767 Pcrs->setAllValues(P_rowptr,P_colind,P_values);
768 Pcrs->expertStaticFillComplete(coarseDomainMap, A.getDomainMap());
770 P = rcp(
new CrsMatrixWrap(Pcrs));
776template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
778Coarsen_Direct(
const Matrix & A,
const RCP<const Matrix> & Aghost,
const GraphBase & graph, RCP<const Map> & coarseColMap, RCP<const Map> & coarseDomainMap, LO num_c_points, LO num_f_points,
const Teuchos::ArrayView<const LO> & myPointType,
const Teuchos::ArrayView<const LO> & myPointType_ghost,
const Teuchos::Array<LO> & cpoint2pcol,
const Teuchos::Array<LO> & pcol2cpoint, Teuchos::Array<size_t> & eis_rowptr, Teuchos::Array<bool> & edgeIsStrong, RCP<LocalOrdinalVector> & BlockNumber, RCP<Matrix> & P)
const {
824 using STS =
typename Teuchos::ScalarTraits<SC>;
825 using MT =
typename STS::magnitudeType;
826 using MTS =
typename Teuchos::ScalarTraits<MT>;
827 size_t Nrows = A.getLocalNumRows();
828 double c_point_density = (double)num_c_points / (num_c_points+num_f_points);
831 double nnz_per_row_est = c_point_density*mean_strong_neighbors_per_row;
833 size_t nnz_est = std::max(Nrows,std::min((
size_t)A.getLocalNumEntries(),(
size_t)(nnz_per_row_est*Nrows)));
834 SC SC_ZERO = STS::zero();
835 MT MT_ZERO = MTS::zero();
836 Array<size_t> tmp_rowptr(Nrows+1);
837 Array<LO> tmp_colind(nnz_est);
843 for(LO row=0; row < (LO) Nrows; row++) {
844 size_t row_start = eis_rowptr[row];
845 ArrayView<const LO> indices;
846 ArrayView<const SC> vals;
848 if(myPointType[row] == DIRICHLET_PT) {
851 else if(myPointType[row] == C_PT) {
853 C_hat.insert(cpoint2pcol[row]);
859 A.getLocalRowView(row, indices, vals);
860 for(LO j=0; j<indices.size(); j++)
861 if(myPointType[indices[j]] == C_PT && edgeIsStrong[row_start + j])
862 C_hat.insert(cpoint2pcol[indices[j]]);
866 if(ct + (
size_t)C_hat.size() > (size_t)tmp_colind.size()) {
867 tmp_colind.resize(std::max(ct+(
size_t)C_hat.size(),(size_t)2*tmp_colind.size()));
871 std::copy(C_hat.begin(), C_hat.end(),tmp_colind.begin()+ct);
873 tmp_rowptr[row+1] = tmp_rowptr[row] + C_hat.size();
876 tmp_colind.resize(tmp_rowptr[Nrows]);
879 P = rcp(
new CrsMatrixWrap(A.getRowMap(), coarseColMap, 0));
880 RCP<CrsMatrix> PCrs = rcp_dynamic_cast<CrsMatrixWrap>(P)->getCrsMatrix();
881 ArrayRCP<size_t> P_rowptr;
882 ArrayRCP<LO> P_colind;
883 ArrayRCP<SC> P_values;
886printf(
"CMS: Allocating P w/ %d nonzeros\n",(
int)tmp_rowptr[Nrows]);
888 PCrs->allocateAllValues(tmp_rowptr[Nrows], P_rowptr, P_colind, P_values);
889 TEUCHOS_TEST_FOR_EXCEPTION(tmp_rowptr.size() !=P_rowptr.size(),
Exceptions::RuntimeError,
"ClassicalPFactory: Allocation size error (rowptr)");
890 TEUCHOS_TEST_FOR_EXCEPTION(tmp_colind.size() !=P_colind.size(),
Exceptions::RuntimeError,
"ClassicalPFactory: Allocation size error (colind)");
892 for(LO i=0; i<(LO)Nrows+1; i++)
893 P_rowptr[i] = tmp_rowptr[i];
894 for(LO i=0; i<(LO)tmp_rowptr[Nrows]; i++)
895 P_colind[i] = tmp_colind[i];
899 for(LO i=0; i < (LO)Nrows; i++) {
900 if(myPointType[i] == DIRICHLET_PT) {
904 printf(
"** A(%d,:) is a Dirichlet-Point.\n",i);
907 else if (myPointType[i] == C_PT) {
909 P_values[P_rowptr[i]] = Teuchos::ScalarTraits<SC>::one();
912 printf(
"** A(%d,:) is a C-Point.\n",i);
921 ArrayView<const LO> A_indices_i, A_incides_k;
922 ArrayView<const SC> A_vals_i, A_indices_k;
923 A.getLocalRowView(i, A_indices_i, A_vals_i);
924 size_t row_start = eis_rowptr[i];
926 ArrayView<LO> P_indices_i = P_colind.view(P_rowptr[i],P_rowptr[i+1] - P_rowptr[i]);
927 ArrayView<SC> P_vals_i = P_values.view(P_rowptr[i],P_rowptr[i+1] - P_rowptr[i]);
932 char mylabel[5]=
"FUCD";
934 printf(
"** A(%d,:) = ",i);
935 for(LO j=0; j<(LO)A_indices_i.size(); j++){
936 printf(
"%6.4e(%d-%c%c) ",A_vals_i[j],A_indices_i[j],mylabel[1+myPointType[A_indices_i[j]]],sw[(
int)edgeIsStrong[row_start+j]]);
943 SC pos_numerator = SC_ZERO, neg_numerator = SC_ZERO;
944 SC pos_denominator = SC_ZERO, neg_denominator = SC_ZERO;
946 for(LO j=0; j<(LO)A_indices_i.size(); j++) {
947 SC a_ik = A_vals_i[j];
948 LO k = A_indices_i[j];
955 if(myPointType[k] == C_PT && edgeIsStrong[row_start + j]) {
956 if(STS::real(a_ik) > MT_ZERO) pos_denominator += a_ik;
957 else neg_denominator += a_ik;
963 if(STS::real(a_ik) > MT_ZERO) pos_numerator += a_ik;
964 else neg_numerator += a_ik;
967 SC alpha = (neg_denominator == MT_ZERO) ? SC_ZERO : (neg_numerator / neg_denominator);
968 SC beta = (pos_denominator == MT_ZERO) ? SC_ZERO : (pos_numerator / pos_denominator);
973 for(LO p_j=0; p_j<(LO)P_indices_i.size(); p_j++){
974 LO P_col = pcol2cpoint[P_indices_i[p_j]];
979 for(LO a_j =0; a_j<(LO)A_indices_i.size(); a_j++) {
980 if(A_indices_i[a_j] == P_col) {
981 a_ij = A_vals_i[a_j];
985 SC w_ij = (STS::real(a_ij) < 0 ) ? (alpha * a_ij) : (beta * a_ij);
987 SC alpha_or_beta = (STS::real(a_ij) < 0 ) ? alpha : beta;
988 printf(
"P(%d,%d/%d) = - %6.4e * %6.4e = %6.4e\n",i,P_indices_i[p_j],pcol2cpoint[P_indices_i[p_j]],alpha_or_beta,a_ij,w_ij);
990 P_vals_i[p_j] = w_ij;
996 PCrs->setAllValues(P_rowptr, P_colind, P_values);
997 PCrs->expertStaticFillComplete(coarseDomainMap, A.getDomainMap());
1002template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
1004Coarsen_Ext_Plus_I(
const Matrix & A,
const RCP<const Matrix> & Aghost,
const GraphBase & graph, RCP<const Map> & coarseColMap, RCP<const Map> & coarseDomainMap, LO num_c_points, LO num_f_points,
const Teuchos::ArrayView<const LO> & myPointType,
const Teuchos::ArrayView<const LO> & myPointType_ghost,
const Teuchos::Array<LO> & cpoint2pcol,
const Teuchos::Array<LO> & pcol2cpoint, Teuchos::Array<size_t> & eis_rowptr, Teuchos::Array<bool> & edgeIsStrong, RCP<LocalOrdinalVector> & BlockNumber, RCP<Matrix> & P)
const {
1041 TEUCHOS_TEST_FOR_EXCEPTION(1,std::runtime_error,
"ClassicalPFactory: Ext+i not implemented");
1049template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
1055 size_t Nrows = A.getLocalNumRows();
1056 eis_rowptr.resize(Nrows+1);
1058 if(edgeIsStrong.size() == 0) {
1060 edgeIsStrong.resize(A.getLocalNumEntries(),
false);
1063 edgeIsStrong.resize(A.getLocalNumEntries(),
false);
1064 edgeIsStrong.assign(A.getLocalNumEntries(),
false);
1068 for (LO i=0; i<(LO)Nrows; i++) {
1069 LO rowstart = eis_rowptr[i];
1070 ArrayView<const LO> A_indices;
1071 ArrayView<const SC> A_values;
1072 A.getLocalRowView(i, A_indices, A_values);
1073 LO A_size = (LO) A_indices.size();
1076 LO G_size = (LO) G_indices.size();
1080 for(LO j=0; j<A_size-1; j++)
1081 if (A_indices[j] >= A_indices[j+1]) { is_ok=
false;
break;}
1082 for(LO j=0; j<G_size-1; j++)
1083 if (G_indices[j] >= G_indices[j+1]) { is_ok=
false;
break;}
1084 TEUCHOS_TEST_FOR_EXCEPTION(!is_ok,
Exceptions::RuntimeError,
"ClassicalPFactory: Exected A and Graph to be sorted");
1087 for(LO g_idx=0, a_idx=0; g_idx < G_size; g_idx++) {
1088 LO col = G_indices[g_idx];
1089 while (A_indices[a_idx] != col && a_idx < A_size) a_idx++;
1090 if(a_idx == A_size) {is_ok=
false;
break;}
1091 edgeIsStrong[rowstart+a_idx] =
true;
1094 eis_rowptr[i+1] = eis_rowptr[i] + A_size;
1100template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
1102GhostCoarseMap(
const Matrix &A,
const Import & Importer,
const ArrayRCP<const LO> myPointType,
const RCP<const Map> & coarseMap, RCP<const Map> & coarseColMap)
const {
1104 const GO GO_INVALID = Teuchos::OrdinalTraits<GO>::invalid();
1105 RCP<GlobalOrdinalVector> d_coarseIds = GlobalOrdinalVectorFactory::Build(A.getRowMap());
1106 ArrayRCP<GO> d_data = d_coarseIds->getDataNonConst(0);
1109 for(LO i=0; i<(LO)d_data.size(); i++) {
1110 if(myPointType[i] == C_PT) {
1111 d_data[i] = coarseMap->getGlobalElement(ct);
1115 d_data[i] = GO_INVALID;
1119 RCP<GlobalOrdinalVector> c_coarseIds = GlobalOrdinalVectorFactory::Build(A.getColMap());
1120 c_coarseIds->doImport(*d_coarseIds,Importer,Xpetra::INSERT);
1125 ArrayRCP<GO> c_data = c_coarseIds->getDataNonConst(0);
1127 Array<GO> c_gids(c_data.size());
1130 for(LO i=0; i<(LO)c_data.size(); i++) {
1131 if(c_data[i] != GO_INVALID) {
1132 c_gids[count] = c_data[i];
1137 std::vector<size_t> stridingInfo_(1);
1139 GO domainGIDOffset = 0;
1141 coarseColMap = StridedMapFactory::Build(coarseMap->lib(),
1142 Teuchos::OrdinalTraits<Xpetra::global_size_t>::invalid(),
1143 c_gids.view(0,count),
1144 coarseMap->getIndexBase(),
1146 coarseMap->getComm(),
1156#define MUELU_CLASSICALPFACTORY_SHORT
#define SET_VALID_ENTRY(name)
void BuildP(Level &fineLevel, Level &coarseLevel) const
Abstract Build method.
void Coarsen_Direct(const Matrix &A, const RCP< const Matrix > &Aghost, const GraphBase &graph, RCP< const Map > &coarseColMap, RCP< const Map > &coarseDomainMap, LO num_c_points, LO num_f_points, const Teuchos::ArrayView< const LO > &myPointType, const Teuchos::ArrayView< const LO > &myPointType_ghost, const Teuchos::Array< LO > &cpoint2pcol, const Teuchos::Array< LO > &pcol2cpoint, Teuchos::Array< size_t > &eis_rowptr, Teuchos::Array< bool > &edgeIsStrong, RCP< LocalOrdinalVector > &BlockNumber, RCP< Matrix > &P) const
void GenerateStrengthFlags(const Matrix &A, const GraphBase &graph, Teuchos::Array< size_t > &eis_rowptr, Teuchos::Array< bool > &edgeIsStrong) const
void Build(Level &fineLevel, Level &coarseLevel) const
Build an object with this factory.
void DeclareInput(Level &fineLevel, Level &coarseLevel) const
Input.
void Coarsen_Ext_Plus_I(const Matrix &A, const RCP< const Matrix > &Aghost, const GraphBase &graph, RCP< const Map > &coarseColMap, RCP< const Map > &coarseDomainMap, LO num_c_points, LO num_f_points, const Teuchos::ArrayView< const LO > &myPointType, const Teuchos::ArrayView< const LO > &myPointType_ghost, const Teuchos::Array< LO > &cpoint2pcol, const Teuchos::Array< LO > &pcol2cpoint, Teuchos::Array< size_t > &eis_rowptr, Teuchos::Array< bool > &edgeIsStrong, RCP< LocalOrdinalVector > &BlockNumber, RCP< Matrix > &P) const
RCP< const ParameterList > GetValidParameterList() const
Return a const parameter list of valid parameters that setParameterList() will accept.
typename ClassicalMapFactory::point_type point_type
void Coarsen_ClassicalModified(const Matrix &A, const RCP< const Matrix > &Aghost, const GraphBase &graph, RCP< const Map > &coarseColMap, RCP< const Map > &coarseDomainMap, LO num_c_points, LO num_f_points, const Teuchos::ArrayView< const LO > &myPointType, const Teuchos::ArrayView< const LO > &myPointType_ghost, const Teuchos::Array< LO > &cpoint2pcol, const Teuchos::Array< LO > &pcol2cpoint, Teuchos::Array< size_t > &eis_rowptr, Teuchos::Array< bool > &edgeIsStrong, RCP< LocalOrdinalVector > &BlockNumber, RCP< const Import > remoteOnlyImporter, RCP< Matrix > &P) const
void GhostCoarseMap(const Matrix &A, const Import &Importer, const ArrayRCP< const LO > myPointType, const RCP< const Map > &coarseMap, RCP< const Map > &coarseColMap) const
Exception throws to report errors in the internal logical of the program.
Timer to be used in factories. Similar to Monitor but with additional timers.
MueLu representation of a graph.
virtual size_t GetNodeNumEdges() const =0
Return number of edges owned by the calling node.
virtual Teuchos::ArrayView< const LocalOrdinal > getNeighborVertices(LocalOrdinal v) const =0
Return the list of vertices adjacent to the vertex 'v'.
virtual size_t GetNodeNumVertices() const =0
Return number of vertices owned by the calling node.
Class that holds all level-specific information.
int GetLevelID() const
Return level number.
static std::string PrintMatrixInfo(const Matrix &A, const std::string &msgTag, RCP< const Teuchos::ParameterList > params=Teuchos::null)
Timer to be used in factories. Similar to SubMonitor but adds a timer level by level.
Namespace for MueLu classes and methods.
@ Debug
Print additional debugging information.
@ Statistics1
Print more statistics.