46 #ifndef MUELU_CLASSICALPFACTORY_DEF_HPP
47 #define MUELU_CLASSICALPFACTORY_DEF_HPP
49 #include <Xpetra_MultiVectorFactory.hpp>
53 #include <Xpetra_Map.hpp>
54 #include <Xpetra_MapFactory.hpp>
57 #include <Xpetra_ImportUtils.hpp>
59 #include <Xpetra_StridedMapFactory.hpp>
65 #include "MueLu_PerfUtils.hpp"
67 #include "MueLu_ClassicalMapFactory.hpp"
68 #include "MueLu_AmalgamationInfo.hpp"
69 #include "MueLu_LWGraph.hpp"
80 if (STS::real(val) > MT_ZERO)
82 else if (STS::real(val) < MT_ZERO)
92 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
95 #define SET_VALID_ENTRY(name) validParamList->setEntry(name, MasterList::getEntry(name))
106 #undef SET_VALID_ENTRY
110 validParamList->
set<
RCP<const FactoryBase> >(
"DofsPerNode", null,
"Generating factory for variable \'DofsPerNode\', usually the same as for \'Graph\'");
116 return validParamList;
119 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
121 Input(fineLevel,
"A");
122 Input(fineLevel,
"Graph");
123 Input(fineLevel,
"DofsPerNode");
125 Input(fineLevel,
"CoarseMap");
126 Input(fineLevel,
"FC Splitting");
129 std::string drop_algo = pL.
get<std::string>(
"aggregation: drop scheme");
130 if (drop_algo.find(
"block diagonal") != std::string::npos || drop_algo ==
"signed classical") {
131 Input(fineLevel,
"BlockNumber");
135 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
137 return BuildP(fineLevel, coarseLevel);
140 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
151 RCP<const Map> ownedCoarseMap = Get<RCP<const Map> >(fineLevel,
"CoarseMap");
174 std::string scheme = pL.
get<std::string>(
"aggregation: classical scheme");
175 bool need_ghost_rows =
false;
176 if (scheme ==
"ext+i")
177 need_ghost_rows =
true;
178 else if (scheme ==
"direct")
179 need_ghost_rows =
false;
180 else if (scheme ==
"classical modified")
181 need_ghost_rows =
true;
185 GetOStream(
Statistics1) <<
"ClassicalPFactory: scheme = " << scheme << std::endl;
191 if (Importer.is_null()) {
192 fc_splitting = owned_fc_splitting;
194 RCP<LocalOrdinalVector> fc_splitting_nonconst = LocalOrdinalVectorFactory::Build(A->getCrsGraph()->getColMap());
195 fc_splitting_nonconst->doImport(*owned_fc_splitting, *Importer,
Xpetra::INSERT);
196 fc_splitting = fc_splitting_nonconst;
198 myPointType = fc_splitting->getData(0);
205 if (need_ghost_rows && !Importer.is_null()) {
207 size_t numRemote = Importer->getNumRemoteIDs();
209 for (
size_t i = 0; i < numRemote; i++)
210 remoteRows[i] = Importer->getTargetMap()->getGlobalElement(remoteLIDs[i]);
213 A->getDomainMap()->getIndexBase(), A->getDomainMap()->getComm());
215 remoteOnlyImporter = Importer->createRemoteOnlyImport(remoteRowMap);
217 RCP<CrsMatrix> Aghost_crs = CrsMatrixFactory::Build(Acrs, *remoteOnlyImporter, A->getDomainMap(), remoteOnlyImporter->getTargetMap());
218 Aghost =
rcp(
new CrsMatrixWrap(Aghost_crs));
232 if (Importer.is_null())
233 coarseMap = ownedCoarseMap;
236 GhostCoarseMap(*A, *Importer, myPointType, ownedCoarseMap, coarseMap);
241 std::string drop_algo = pL.
get<std::string>(
"aggregation: drop scheme");
242 if (drop_algo.find(
"block diagonal") != std::string::npos || drop_algo ==
"signed classical") {
244 OwnedBlockNumber = Get<RCP<LocalOrdinalVector> >(fineLevel,
"BlockNumber");
245 if (Importer.is_null())
246 BlockNumber = OwnedBlockNumber;
248 BlockNumber = LocalOrdinalVectorFactory::Build(A->getColMap());
249 BlockNumber->doImport(*OwnedBlockNumber, *Importer,
Xpetra::INSERT);
253 #if defined(CMS_DEBUG) || defined(CMS_DUMP)
256 int rank = A->getRowMap()->getComm()->getRank();
257 printf(
"[%d] A local size = %dx%d\n", rank, (
int)Acrs->getRowMap()->getLocalNumElements(), (int)Acrs->getColMap()->getLocalNumElements());
259 printf(
"[%d] graph local size = %dx%d\n", rank, (
int)graph->GetDomainMap()->getLocalNumElements(), (int)graph->GetImportMap()->getLocalNumElements());
261 std::ofstream ofs(std::string(
"dropped_graph_") + std::to_string(fineLevel.GetLevelID()) + std::string(
"_") + std::to_string(rank) + std::string(
".dat"), std::ofstream::out);
263 graph->print(*fancy,
Debug);
264 std::string out_fc = std::string(
"fc_splitting_") + std::to_string(fineLevel.GetLevelID()) + std::string(
"_") + std::to_string(rank) + std::string(
".dat");
265 std::string out_block = std::string(
"block_number_") + std::to_string(fineLevel.GetLevelID()) + std::string(
"_") + std::to_string(rank) + std::string(
".dat");
277 for (
LO i = 0; i < (
LO)fc_data.
size(); i++)
278 mv_data[i] = Teuchos::as<real_type>(fc_data[i]);
286 for (
LO i = 0; i < (
LO)b_data.
size(); i++) {
287 mv_data2[i] = Teuchos::as<real_type>(b_data[i]);
302 Array<LO> pcol2cpoint(coarseMap->getLocalNumElements(), LO_INVALID);
305 for (
LO i = 0; i < (
LO)myPointType.
size(); i++) {
306 if (myPointType[i] == C_PT) {
307 cpoint2pcol[i] = num_c_points;
309 }
else if (myPointType[i] == F_PT)
312 for (
LO i = 0; i < (
LO)cpoint2pcol.size(); i++) {
313 if (cpoint2pcol[i] != LO_INVALID)
314 pcol2cpoint[cpoint2pcol[i]] = i;
323 GenerateStrengthFlags(*A, *graph, eis_rowptr, edgeIsStrong);
329 if (scheme ==
"ext+i") {
331 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);
332 }
else if (scheme ==
"direct") {
334 Coarsen_Direct(*A, Aghost, *graph, coarseColMap, coarseDomainMap, num_c_points, num_f_points, myPointType(), myPointType_ghost(), cpoint2pcol, pcol2cpoint, eis_rowptr, edgeIsStrong, BlockNumber, P);
335 }
else if (scheme ==
"classical modified") {
337 Coarsen_ClassicalModified(*A, Aghost, *graph, coarseColMap, coarseDomainMap, num_c_points, num_f_points, myPointType(), myPointType_ghost(), cpoint2pcol, pcol2cpoint, eis_rowptr, edgeIsStrong, BlockNumber, remoteOnlyImporter, P);
347 Set(coarseLevel,
"P", P);
349 Set(coarseLevel,
"P Graph", pg);
357 params->
set(
"printLoadBalancingInfo",
true);
362 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
364 Coarsen_ClassicalModified(
const Matrix& A,
const RCP<const Matrix>& Aghost,
const LWGraph& 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 {
405 SC SC_ZERO = STS::zero();
407 int rank = A.getRowMap()->getComm()->getRank();
413 block_id = BlockNumber->getData(0);
418 size_t Nrows = A.getLocalNumRows();
419 double c_point_density = (double)num_c_points / (num_c_points + num_f_points);
422 double nnz_per_row_est = c_point_density * mean_strong_neighbors_per_row;
424 size_t nnz_est = std::max(Nrows, std::min((
size_t)A.getLocalNumEntries(), (size_t)(nnz_per_row_est * Nrows)));
432 for (
LO row = 0; row < (
LO)Nrows; row++) {
433 size_t row_start = eis_rowptr[row];
437 if (myPointType[row] == DIRICHLET_PT) {
439 }
else if (myPointType[row] == C_PT) {
441 C_hat.insert(cpoint2pcol[row]);
446 A.getLocalRowView(row, indices, vals);
447 for (
LO j = 0; j < indices.
size(); j++)
448 if (myPointType[indices[j]] == C_PT && edgeIsStrong[row_start + j])
449 C_hat.insert(cpoint2pcol[indices[j]]);
453 if (ct + (
size_t)C_hat.size() > (size_t)tmp_colind.
size()) {
454 tmp_colind.
resize(std::max(ct + (
size_t)C_hat.size(), (size_t)2 * tmp_colind.
size()));
458 std::copy(C_hat.begin(), C_hat.end(), tmp_colind.
begin() + ct);
460 tmp_rowptr[row + 1] = tmp_rowptr[row] + C_hat.
size();
463 tmp_colind.
resize(tmp_rowptr[Nrows]);
465 RCP<CrsMatrix> Pcrs = CrsMatrixFactory::Build(A.getRowMap(), coarseColMap, 0);
471 P_rowptr = Teuchos::arcpFromArray(tmp_rowptr);
472 P_colind = Teuchos::arcpFromArray(tmp_colind);
473 P_values.
resize(P_rowptr[Nrows]);
478 Pcrs->allocateAllValues(tmp_rowptr[Nrows], P_rowptr, P_colind, P_values);
479 std::copy(tmp_rowptr.
begin(), tmp_rowptr.
end(), P_rowptr.begin());
480 std::copy(tmp_colind.
begin(), tmp_colind.
end(), P_colind.
begin());
481 Pcrs->setAllValues(P_rowptr, P_colind, P_values);
482 Pcrs->expertStaticFillComplete( coarseDomainMap, A.getDomainMap());
492 if (!remoteOnlyImporter.
is_null()) {
494 RCP<CrsGraph> tempPgraph = CrsGraphFactory::Build(A.getRowMap(), coarseColMap, P_rowptr, P_colind);
495 tempPgraph->fillComplete(coarseDomainMap, A.getDomainMap());
499 Pgraph = Pcrs->getCrsGraph();
502 Pghost = CrsGraphFactory::Build(Pgraph, *remoteOnlyImporter, Pgraph->getDomainMap(), remoteOnlyImporter->getTargetMap());
503 Pghost->getAllIndices(Pghost_rowptr, Pghost_colind);
507 ArrayRCP<LO> Acol_to_Pcol(A.getColMap()->getLocalNumElements(), LO_INVALID);
512 Pghostcol_to_Pcol.
resize(Pghost->getColMap()->getLocalNumElements(), LO_INVALID);
513 for (
LO i = 0; i < (
LO)Pghost->getColMap()->getLocalNumElements(); i++)
514 Pghostcol_to_Pcol[i] = Pgraph->getColMap()->getLocalElement(Pghost->getColMap()->getGlobalElement(i));
520 Aghostcol_to_Acol.
resize(Aghost->getColMap()->getLocalNumElements(), LO_INVALID);
521 for (
LO i = 0; i < (
LO)Aghost->getColMap()->getLocalNumElements(); i++)
522 Aghostcol_to_Acol[i] = A.getColMap()->getLocalElement(Aghost->getColMap()->getGlobalElement(i));
526 for (
LO i = 0; i < (
LO)Nrows; i++) {
527 if (myPointType[i] == DIRICHLET_PT) {
531 printf(
"[%d] ** A(%d,:) is a Dirichlet-Point.\n", rank, i);
533 }
else if (myPointType[i] == C_PT) {
538 printf(
"[%d] ** A(%d,:) is a C-Point.\n", rank, i);
544 printf(
"[%d] ** A(%d,:) is a F-Point.\n", rank, i);
550 A.getLocalRowView(i, A_indices_i, A_vals_i);
551 size_t row_start = eis_rowptr[i];
554 ArrayView<SC> P_vals_i = P_values.
view(P_rowptr[i], P_rowptr[i + 1] - P_rowptr[i]);
557 for (
LO j = 0; j < (
LO)P_vals_i.
size(); j++)
558 P_vals_i[j] = SC_ZERO;
562 for (
LO j = 0; j < (
LO)P_indices_i.
size(); j++) {
563 Acol_to_Pcol[pcol2cpoint[P_indices_i[j]]] = P_rowptr[i] + j;
567 SC first_denominator = SC_ZERO;
571 for (
LO k0 = 0; k0 < (
LO)A_indices_i.
size(); k0++) {
572 LO k = A_indices_i[k0];
573 SC a_ik = A_vals_i[k0];
574 LO pcol_k = Acol_to_Pcol[k];
579 first_denominator += a_ik;
582 printf(
"- A(%d,%d) is the diagonal\n", i, k);
585 }
else if (myPointType[k] == DIRICHLET_PT) {
588 printf(
"- A(%d,%d) is a Dirichlet point\n", i, k);
591 }
else if (pcol_k != LO_INVALID && pcol_k >= (
LO)P_rowptr[i]) {
593 P_values[pcol_k] += a_ik;
595 printf(
"- A(%d,%d) is a strong C-Point\n", i, k);
597 }
else if (!edgeIsStrong[row_start + k0]) {
599 if (block_id.
size() == 0 || block_id[i] == block_id[k]) {
600 first_denominator += a_ik;
602 printf(
"- A(%d,%d) is weak adding to diagonal(%d,%d) (%6.4e)\n", i, k, block_id[i], block_id[k], a_ik);
604 printf(
"- A(%d,%d) is weak but does not match blocks (%d,%d), discarding\n", i, k, block_id[i], block_id[k]);
612 printf(
"- A(%d,%d) is a strong F-Point\n", i, k);
616 SC second_denominator = SC_ZERO;
621 A.getLocalRowView(k, A_indices_k, A_vals_k);
622 for (
LO m0 = 0; m0 < (
LO)A_indices_k.
size(); m0++) {
623 LO m = A_indices_k[m0];
631 sign_akk = Sign(a_kk);
632 for (
LO m0 = 0; m0 < (
LO)A_indices_k.
size(); m0++) {
633 LO m = A_indices_k[m0];
634 if (m != k && Acol_to_Pcol[A_indices_k[m0]] >= (
LO)P_rowptr[i]) {
635 SC a_km = A_vals_k[m0];
636 second_denominator += (Sign(a_km) == sign_akk ? SC_ZERO : a_km);
642 if (second_denominator != SC_ZERO) {
643 for (
LO j0 = 0; j0 < (
LO)A_indices_k.
size(); j0++) {
644 LO j = A_indices_k[j0];
647 if (Acol_to_Pcol[j] >= (
LO)P_rowptr[i]) {
648 SC a_kj = A_vals_k[j0];
649 SC sign_akj_val = sign_akk == Sign(a_kj) ? SC_ZERO : a_kj;
650 P_values[Acol_to_Pcol[j]] += a_ik * sign_akj_val / second_denominator;
652 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]]);
659 printf(
"- - A(%d,%d) second denominator is zero\n", i, k);
661 if (block_id.
size() == 0 || block_id[i] == block_id[k])
662 first_denominator += a_ik;
667 LO kless = k - Nrows;
671 Aghost->getLocalRowView(kless, A_indices_k, A_vals_k);
672 GO k_g = Aghost->getRowMap()->getGlobalElement(kless);
673 for (
LO m0 = 0; m0 < (
LO)A_indices_k.
size(); m0++) {
674 GO m_g = Aghost->getColMap()->getGlobalElement(A_indices_k[m0]);
682 sign_akk = Sign(a_kk);
683 for (
LO m0 = 0; m0 < (
LO)A_indices_k.
size(); m0++) {
684 GO m_g = Aghost->getColMap()->getGlobalElement(A_indices_k[m0]);
685 LO mghost = A_indices_k[m0];
686 LO m = Aghostcol_to_Acol[mghost];
687 if (m_g != k_g && m != LO_INVALID && Acol_to_Pcol[m] >= (
LO)P_rowptr[i]) {
688 SC a_km = A_vals_k[m0];
689 second_denominator += (Sign(a_km) == sign_akk ? SC_ZERO : a_km);
695 if (second_denominator != SC_ZERO) {
696 for (
LO j0 = 0; j0 < (
LO)A_indices_k.
size(); j0++) {
697 LO jghost = A_indices_k[j0];
698 LO j = Aghostcol_to_Acol[jghost];
700 if ((j != LO_INVALID) && (Acol_to_Pcol[j] >= (
LO)P_rowptr[i])) {
701 SC a_kj = A_vals_k[j0];
702 SC sign_akj_val = sign_akk == Sign(a_kj) ? SC_ZERO : a_kj;
703 P_values[Acol_to_Pcol[j]] += a_ik * sign_akj_val / second_denominator;
705 printf(
"- - Unscaled P(%d,A-%d) += %6.4e\n", i, j, a_ik * sign_akj_val / second_denominator);
713 printf(
"- - A(%d,%d) second denominator is zero\n", i, k);
715 if (block_id.
size() == 0 || block_id[i] == block_id[k])
716 first_denominator += a_ik;
724 if (first_denominator != SC_ZERO) {
725 for (
LO j0 = 0; j0 < (
LO)P_indices_i.
size(); j0++) {
727 SC old_pij = P_vals_i[j0];
728 P_vals_i[j0] /= -first_denominator;
729 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);
731 P_vals_i[j0] /= -first_denominator;
741 Pcrs->setAllValues(P_rowptr, P_colind, P_values);
742 Pcrs->expertStaticFillComplete( coarseDomainMap, A.getDomainMap());
744 P =
rcp(
new CrsMatrixWrap(Pcrs));
749 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
751 Coarsen_Direct(
const Matrix& A,
const RCP<const Matrix>& Aghost,
const LWGraph& 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 {
797 using MT =
typename STS::magnitudeType;
799 size_t Nrows = A.getLocalNumRows();
800 double c_point_density = (double)num_c_points / (num_c_points + num_f_points);
803 double nnz_per_row_est = c_point_density * mean_strong_neighbors_per_row;
805 size_t nnz_est = std::max(Nrows, std::min((
size_t)A.getLocalNumEntries(), (size_t)(nnz_per_row_est * Nrows)));
806 SC SC_ZERO = STS::zero();
807 MT MT_ZERO = MTS::zero();
815 for (
LO row = 0; row < (
LO)Nrows; row++) {
816 size_t row_start = eis_rowptr[row];
820 if (myPointType[row] == DIRICHLET_PT) {
822 }
else if (myPointType[row] == C_PT) {
824 C_hat.insert(cpoint2pcol[row]);
829 A.getLocalRowView(row, indices, vals);
830 for (
LO j = 0; j < indices.
size(); j++)
831 if (myPointType[indices[j]] == C_PT && edgeIsStrong[row_start + j])
832 C_hat.insert(cpoint2pcol[indices[j]]);
836 if (ct + (
size_t)C_hat.size() > (size_t)tmp_colind.
size()) {
837 tmp_colind.
resize(std::max(ct + (
size_t)C_hat.size(), (size_t)2 * tmp_colind.
size()));
841 std::copy(C_hat.begin(), C_hat.end(), tmp_colind.
begin() + ct);
843 tmp_rowptr[row + 1] = tmp_rowptr[row] + C_hat.
size();
846 tmp_colind.
resize(tmp_rowptr[Nrows]);
849 P =
rcp(
new CrsMatrixWrap(A.getRowMap(), coarseColMap, 0));
850 RCP<CrsMatrix> PCrs = rcp_dynamic_cast<CrsMatrixWrap>(P)->getCrsMatrix();
856 printf(
"CMS: Allocating P w/ %d nonzeros\n", (
int)tmp_rowptr[Nrows]);
858 PCrs->allocateAllValues(tmp_rowptr[Nrows], P_rowptr, P_colind, P_values);
862 for (
LO i = 0; i < (
LO)Nrows + 1; i++)
863 P_rowptr[i] = tmp_rowptr[i];
864 for (
LO i = 0; i < (
LO)tmp_rowptr[Nrows]; i++)
865 P_colind[i] = tmp_colind[i];
868 for (
LO i = 0; i < (
LO)Nrows; i++) {
869 if (myPointType[i] == DIRICHLET_PT) {
873 printf(
"** A(%d,:) is a Dirichlet-Point.\n", i);
875 }
else if (myPointType[i] == C_PT) {
880 printf(
"** A(%d,:) is a C-Point.\n", i);
890 A.getLocalRowView(i, A_indices_i, A_vals_i);
891 size_t row_start = eis_rowptr[i];
893 ArrayView<LO> P_indices_i = P_colind.
view(P_rowptr[i], P_rowptr[i + 1] - P_rowptr[i]);
894 ArrayView<SC> P_vals_i = P_values.
view(P_rowptr[i], P_rowptr[i + 1] - P_rowptr[i]);
899 char mylabel[5] =
"FUCD";
901 printf(
"** A(%d,:) = ", i);
902 for (
LO j = 0; j < (
LO)A_indices_i.
size(); j++) {
903 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]]);
910 SC pos_numerator = SC_ZERO, neg_numerator = SC_ZERO;
911 SC pos_denominator = SC_ZERO, neg_denominator = SC_ZERO;
913 for (
LO j = 0; j < (
LO)A_indices_i.
size(); j++) {
914 SC a_ik = A_vals_i[j];
915 LO k = A_indices_i[j];
922 if (myPointType[k] == C_PT && edgeIsStrong[row_start + j]) {
923 if (STS::real(a_ik) > MT_ZERO)
924 pos_denominator += a_ik;
926 neg_denominator += a_ik;
932 if (STS::real(a_ik) > MT_ZERO)
933 pos_numerator += a_ik;
935 neg_numerator += a_ik;
938 SC alpha = (neg_denominator == MT_ZERO) ? SC_ZERO : (neg_numerator / neg_denominator);
939 SC beta = (pos_denominator == MT_ZERO) ? SC_ZERO : (pos_numerator / pos_denominator);
944 for (
LO p_j = 0; p_j < (
LO)P_indices_i.
size(); p_j++) {
945 LO P_col = pcol2cpoint[P_indices_i[p_j]];
950 for (
LO a_j = 0; a_j < (
LO)A_indices_i.
size(); a_j++) {
951 if (A_indices_i[a_j] == P_col) {
952 a_ij = A_vals_i[a_j];
956 SC w_ij = (STS::real(a_ij) < 0) ? (alpha * a_ij) : (beta * a_ij);
958 SC alpha_or_beta = (STS::real(a_ij) < 0) ? alpha : beta;
959 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);
961 P_vals_i[p_j] = w_ij;
967 PCrs->setAllValues(P_rowptr, P_colind, P_values);
968 PCrs->expertStaticFillComplete( coarseDomainMap, A.getDomainMap());
972 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
974 Coarsen_Ext_Plus_I(
const Matrix& A,
const RCP<const Matrix>& Aghost,
const LWGraph& 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 {
1013 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
1019 size_t Nrows = A.getLocalNumRows();
1020 eis_rowptr.
resize(Nrows + 1);
1022 if (edgeIsStrong.
size() == 0) {
1024 edgeIsStrong.
resize(A.getLocalNumEntries(),
false);
1026 edgeIsStrong.
resize(A.getLocalNumEntries(),
false);
1027 edgeIsStrong.
assign(A.getLocalNumEntries(),
false);
1031 for (
LO i = 0; i < (
LO)Nrows; i++) {
1032 LO rowstart = eis_rowptr[i];
1035 A.getLocalRowView(i, A_indices, A_values);
1039 LO G_size = (
LO)G_indices.length;
1043 for (
LO j = 0; j < A_size - 1; j++)
1044 if (A_indices[j] >= A_indices[j + 1]) {
1048 for (
LO j = 0; j < G_size - 1; j++)
1049 if (G_indices(j) >= G_indices(j + 1)) {
1056 for (
LO g_idx = 0, a_idx = 0; g_idx < G_size; g_idx++) {
1057 LO col = G_indices(g_idx);
1058 while (A_indices[a_idx] != col && a_idx < A_size) a_idx++;
1059 if (a_idx == A_size) {
1063 edgeIsStrong[rowstart + a_idx] =
true;
1066 eis_rowptr[i + 1] = eis_rowptr[i] + A_size;
1071 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
1080 for (
LO i = 0; i < (
LO)d_data.size(); i++) {
1081 if (myPointType[i] == C_PT) {
1082 d_data[i] = coarseMap->getGlobalElement(ct);
1085 d_data[i] = GO_INVALID;
1100 for (
LO i = 0; i < (
LO)c_data.
size(); i++) {
1101 if (c_data[i] != GO_INVALID) {
1102 c_gids[count] = c_data[i];
1107 std::vector<size_t> stridingInfo_(1);
1108 stridingInfo_[0] = 1;
1109 GO domainGIDOffset = 0;
1111 coarseColMap = StridedMapFactory::Build(coarseMap->lib(),
1113 c_gids.view(0, count),
1114 coarseMap->getIndexBase(),
1116 coarseMap->getComm(),
1122 #define MUELU_CLASSICALPFACTORY_SHORT
1123 #endif // MUELU_CLASSICALPFACTORY_DEF_HPP
void GenerateStrengthFlags(const Matrix &A, const LWGraph &graph, Teuchos::Array< size_t > &eis_rowptr, Teuchos::Array< bool > &edgeIsStrong) const
void setValidator(RCP< const ParameterEntryValidator > const &validator)
T & get(const std::string &name, T def_value)
ParameterList & set(std::string const &name, T const &value, std::string const &docString="", RCP< const ParameterEntryValidator > const &validator=null)
Timer to be used in factories. Similar to Monitor but with additional timers.
#define TEUCHOS_TEST_FOR_EXCEPTION(throw_exception_test, Exception, msg)
void BuildP(Level &fineLevel, Level &coarseLevel) const
Abstract Build method.
static void Write(const std::string &fileName, const Xpetra::Map< LocalOrdinal, GlobalOrdinal, Node > &M)
void GhostCoarseMap(const Matrix &A, const Import &Importer, const ArrayRCP< const LO > myPointType, const RCP< const Map > &coarseMap, RCP< const Map > &coarseColMap) const
Print additional debugging information.
void Build(Level &fineLevel, Level &coarseLevel) const
Build an object with this factory.
void Coarsen_ClassicalModified(const Matrix &A, const RCP< const Matrix > &Aghost, const LWGraph &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
RCP< const ParameterList > GetValidParameterList() const
Return a const parameter list of valid parameters that setParameterList() will accept.
KOKKOS_INLINE_FUNCTION size_type GetNodeNumVertices() const
Return number of graph vertices.
KOKKOS_INLINE_FUNCTION size_type GetNodeNumEdges() const
Return number of graph edges.
void DeclareInput(Level &fineLevel, Level &coarseLevel) const
Input.
void resize(const size_type n, const T &val=T())
TEUCHOS_DEPRECATED RCP< T > rcp(T *p, Dealloc_T dealloc, bool owns_mem)
Class that holds all level-specific information.
Timer to be used in factories. Similar to SubMonitor but adds a timer level by level.
void Coarsen_Ext_Plus_I(const Matrix &A, const RCP< const Matrix > &Aghost, const LWGraph &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 resize(size_type new_size, const value_type &x=value_type())
static std::string PrintMatrixInfo(const Matrix &A, const std::string &msgTag, RCP< const Teuchos::ParameterList > params=Teuchos::null)
KOKKOS_INLINE_FUNCTION neighbor_vertices_type getNeighborVertices(LO i) const
Return the list of vertices adjacent to the vertex 'v'.
#define SET_VALID_ENTRY(name)
Lightweight MueLu representation of a compressed row storage graph.
void Coarsen_Direct(const Matrix &A, const RCP< const Matrix > &Aghost, const LWGraph &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
int GetLevelID() const
Return level number.
typename ClassicalMapFactory::point_type point_type
Exception throws to report errors in the internal logical of the program.
#define TEUCHOS_ASSERT(assertion_test)
void assign(size_type n, const value_type &val)
ParameterEntry & getEntry(const std::string &name)
ArrayView< T > view(size_type lowerOffset, size_type size) const