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>
149 using STS = Teuchos::ScalarTraits<SC>;
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()->getNodeNumElements(),(int)Acrs->getColMap()->getNodeNumElements());
267 printf(
"[%d] graph local size = %dx%d\n",rank,(
int)graph->GetDomainMap()->getNodeNumElements(),(int)graph->GetImportMap()->getNodeNumElements());
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->getNodeNumElements(),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);
375 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
377 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 {
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.getNodeNumRows();
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.getNodeNumEntries(),(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()->getNodeNumElements(),LO_INVALID);
527 ArrayRCP<LO> Pghostcol_to_Pcol;
528 if(!Pghost.is_null()) {
529 Pghostcol_to_Pcol.resize(Pghost->getColMap()->getNodeNumElements(),LO_INVALID);
530 for(LO i=0; i<(LO) Pghost->getColMap()->getNodeNumElements(); 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()->getNodeNumElements(),LO_INVALID);
538 for(LO i=0; i<(LO)Aghost->getColMap()->getNodeNumElements(); 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));
776 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
778 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 {
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.getNodeNumRows();
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.getNodeNumEntries(),(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;
886 printf(
"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());
1002 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
1004 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 {
1041 TEUCHOS_TEST_FOR_EXCEPTION(1,std::runtime_error,
"ClassicalPFactory: Ext+i not implemented");
1049 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
1055 size_t Nrows = A.getNodeNumRows();
1056 eis_rowptr.resize(Nrows+1);
1058 if(edgeIsStrong.size() == 0) {
1060 edgeIsStrong.resize(A.getNodeNumEntries(),
false);
1063 edgeIsStrong.resize(A.getNodeNumEntries(),
false);
1064 edgeIsStrong.assign(A.getNodeNumEntries(),
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;
1100 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
1102 GhostCoarseMap(
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 1157 #endif // MUELU_CLASSICALPFACTORY_DEF_HPP void GhostCoarseMap(const Matrix &A, const Import &Importer, const ArrayRCP< const LO > myPointType, const RCP< const Map > &coarseMap, RCP< const Map > &coarseColMap) const
void Build(Level &fineLevel, Level &coarseLevel) const
Build an object with this factory.
virtual size_t GetNodeNumVertices() const =0
Return number of vertices owned by the calling node.
virtual size_t GetNodeNumEdges() const =0
Return number of edges owned by the calling node.
Timer to be used in factories. Similar to Monitor but with additional timers.
Print additional debugging information.
Namespace for MueLu classes and methods.
int GetLevelID() const
Return level number.
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
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
Class that holds all level-specific information.
Timer to be used in factories. Similar to SubMonitor but adds a timer level by level.
static std::string PrintMatrixInfo(const Matrix &A, const std::string &msgTag, RCP< const Teuchos::ParameterList > params=Teuchos::null)
RCP< const ParameterList > GetValidParameterList() const
Return a const parameter list of valid parameters that setParameterList() will accept.
MueLu representation of a graph.
void DeclareInput(Level &fineLevel, Level &coarseLevel) const
Input.
#define SET_VALID_ENTRY(name)
void BuildP(Level &fineLevel, Level &coarseLevel) const
Abstract Build method.
void GenerateStrengthFlags(const Matrix &A, const GraphBase &graph, Teuchos::Array< size_t > &eis_rowptr, Teuchos::Array< bool > &edgeIsStrong) const
typename ClassicalMapFactory::point_type point_type
Exception throws to report errors in the internal logical of the program.
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
virtual Teuchos::ArrayView< const LocalOrdinal > getNeighborVertices(LocalOrdinal v) const =0
Return the list of vertices adjacent to the vertex 'v'.