46 #ifndef MUELU_FILTEREDAFACTORY_DEF_HPP 47 #define MUELU_FILTEREDAFACTORY_DEF_HPP 49 #include <Xpetra_Matrix.hpp> 50 #include <Xpetra_MatrixFactory.hpp> 51 #include <Xpetra_IO.hpp> 55 #include "MueLu_FactoryManager.hpp" 59 #include "MueLu_Aggregates.hpp" 60 #include "MueLu_AmalgamationInfo.hpp" 61 #include "MueLu_Utilities.hpp" 64 #define MUELU_FILTEREDAFACTORY_LOTS_OF_PRINTING 0 71 std::sort(array.begin(),array.end());
72 std::unique(array.begin(),array.end());
77 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
79 RCP<ParameterList> validParamList = rcp(
new ParameterList());
81 #define SET_VALID_ENTRY(name) validParamList->setEntry(name, MasterList::getEntry(name)) 87 SET_VALID_ENTRY(
"filtered matrix: spread lumping diag dom growth factor");
90 #undef SET_VALID_ENTRY 92 validParamList->set< RCP<const FactoryBase> >(
"A", Teuchos::null,
"Generating factory of the matrix A used for filtering");
93 validParamList->set< RCP<const FactoryBase> >(
"Graph", Teuchos::null,
"Generating factory for coalesced filtered graph");
94 validParamList->set< RCP<const FactoryBase> >(
"Filtering", Teuchos::null,
"Generating factory for filtering boolean");
98 validParamList->set< RCP<const FactoryBase> >(
"Aggregates", Teuchos::null,
"Generating factory of the aggregates");
99 validParamList->set< RCP<const FactoryBase> >(
"UnAmalgamationInfo", Teuchos::null,
"Generating factory of UnAmalgamationInfo");
100 return validParamList;
103 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
105 Input(currentLevel,
"A");
106 Input(currentLevel,
"Filtering");
107 Input(currentLevel,
"Graph");
108 const ParameterList& pL = GetParameterList();
109 if(pL.isParameter(
"filtered matrix: use root stencil") && pL.get<
bool>(
"filtered matrix: use root stencil") ==
true){
110 Input(currentLevel,
"Aggregates");
111 Input(currentLevel,
"UnAmalgamationInfo");
115 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
119 RCP<Matrix> A = Get< RCP<Matrix> >(currentLevel,
"A");
120 if (Get<bool>(currentLevel,
"Filtering") ==
false) {
121 GetOStream(
Runtime0) <<
"Filtered matrix is not being constructed as no filtering is being done" << std::endl;
122 Set(currentLevel,
"A", A);
126 const ParameterList& pL = GetParameterList();
127 bool lumping = pL.get<
bool>(
"filtered matrix: use lumping");
129 GetOStream(
Runtime0) <<
"Lumping dropped entries" << std::endl;
131 bool use_spread_lumping = pL.get<
bool>(
"filtered matrix: use spread lumping");
132 if (use_spread_lumping && (!lumping) )
133 throw std::runtime_error(
"Must also request 'filtered matrix: use lumping' in order to use spread lumping");
135 if (use_spread_lumping) {
136 GetOStream(
Runtime0) <<
"using spread lumping " << std::endl;
139 double DdomAllowGrowthRate = 1.1;
140 double DdomCap = 2.0;
141 if (use_spread_lumping) {
142 DdomAllowGrowthRate = pL.get<
double>(
"filtered matrix: spread lumping diag dom growth factor");
143 DdomCap = pL.get<
double>(
"filtered matrix: spread lumping diag dom cap");
145 bool use_root_stencil = lumping && pL.get<
bool>(
"filtered matrix: use root stencil");
146 if (use_root_stencil)
147 GetOStream(
Runtime0) <<
"Using root stencil for dropping" << std::endl;
148 double dirichlet_threshold = pL.get<
double>(
"filtered matrix: Dirichlet threshold");
149 if(dirichlet_threshold >= 0.0)
150 GetOStream(
Runtime0) <<
"Filtering Dirichlet threshold of "<<dirichlet_threshold << std::endl;
152 if(use_root_stencil || pL.get<
bool>(
"filtered matrix: reuse graph"))
153 GetOStream(
Runtime0) <<
"Reusing graph"<<std::endl;
155 GetOStream(
Runtime0) <<
"Generating new graph"<<std::endl;
158 RCP<GraphBase> G = Get< RCP<GraphBase> >(currentLevel,
"Graph");
161 FILE * f = fopen(
"graph.dat",
"w");
162 size_t numGRows = G->GetNodeNumVertices();
163 for (
size_t i = 0; i < numGRows; i++) {
165 ArrayView<const LO> indsG = G->getNeighborVertices(i);
166 for(
size_t j=0; j<(size_t)indsG.size(); j++) {
167 fprintf(f,
"%d %d 1.0\n",(
int)i,(
int)indsG[j]);
173 RCP<ParameterList> fillCompleteParams(
new ParameterList);
174 fillCompleteParams->set(
"No Nonlocal Changes",
true);
176 RCP<Matrix> filteredA;
177 if(use_root_stencil) {
178 filteredA = MatrixFactory::Build(A->getCrsGraph());
179 filteredA->fillComplete(fillCompleteParams);
180 filteredA->resumeFill();
181 BuildNewUsingRootStencil(*A, *G, dirichlet_threshold, currentLevel,*filteredA, use_spread_lumping,DdomAllowGrowthRate, DdomCap);
182 filteredA->fillComplete(fillCompleteParams);
185 else if (pL.get<
bool>(
"filtered matrix: reuse graph")) {
186 filteredA = MatrixFactory::Build(A->getCrsGraph());
187 filteredA->resumeFill();
188 BuildReuse(*A, *G, (lumping != use_spread_lumping), dirichlet_threshold,*filteredA);
192 if (use_spread_lumping) ExperimentalLumping(*A, *filteredA, DdomAllowGrowthRate, DdomCap);
193 filteredA->fillComplete(fillCompleteParams);
197 filteredA = MatrixFactory::Build(A->getRowMap(), A->getColMap(), A->getNodeMaxNumRowEntries());
198 BuildNew(*A, *G, (lumping != use_spread_lumping), dirichlet_threshold,*filteredA);
201 if (use_spread_lumping) ExperimentalLumping(*A, *filteredA, DdomAllowGrowthRate, DdomCap);
202 filteredA->fillComplete(A->getDomainMap(), A->getRangeMap(), fillCompleteParams);
209 Xpetra::IO<SC,LO,GO,NO>::Write(
"filteredA.dat", *filteredA);
212 Xpetra::IO<SC,LO,GO,NO>::Write(
"A.dat", *A);
213 RCP<Matrix> origFilteredA = MatrixFactory::Build(A->getRowMap(), A->getColMap(), A->getNodeMaxNumRowEntries());
214 BuildNew(*A, *G, lumping, dirichlet_threshold,*origFilteredA);
215 if (use_spread_lumping) ExperimentalLumping(*A, *origFilteredA, DdomAllowGrowthRate, DdomCap);
216 origFilteredA->fillComplete(A->getDomainMap(), A->getRangeMap(), fillCompleteParams);
217 Xpetra::IO<SC,LO,GO,NO>::Write(
"origFilteredA.dat", *origFilteredA);
221 filteredA->SetFixedBlockSize(A->GetFixedBlockSize());
223 if (pL.get<
bool>(
"filtered matrix: reuse eigenvalue")) {
228 filteredA->SetMaxEigenvalueEstimate(A->GetMaxEigenvalueEstimate());
231 Set(currentLevel,
"A", filteredA);
240 #define ASSUME_DIRECT_ACCESS_TO_ROW 249 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
251 BuildReuse(
const Matrix& A,
const GraphBase& G,
const bool lumping,
double dirichletThresh, Matrix& filteredA)
const {
252 using TST =
typename Teuchos::ScalarTraits<SC>;
253 SC zero = TST::zero();
256 size_t blkSize = A.GetFixedBlockSize();
258 ArrayView<const LO> inds;
259 ArrayView<const SC> valsA;
260 #ifdef ASSUME_DIRECT_ACCESS_TO_ROW 266 Array<char> filter( std::max(blkSize*G.
GetImportMap()->getNodeNumElements(),
267 A.getColMap()->getNodeNumElements()),
271 for (
size_t i = 0; i < numGRows; i++) {
274 for (
size_t j = 0; j < as<size_t>(indsG.size()); j++)
275 for (
size_t k = 0; k < blkSize; k++)
276 filter[indsG[j]*blkSize+k] = 1;
278 for (
size_t k = 0; k < blkSize; k++) {
279 LO row = i*blkSize + k;
281 A.getLocalRowView(row, inds, valsA);
283 size_t nnz = inds.size();
287 #ifdef ASSUME_DIRECT_ACCESS_TO_ROW 289 ArrayView<const SC> vals1;
290 filteredA.getLocalRowView(row, inds, vals1);
291 vals = ArrayView<SC>(
const_cast<SC*
>(vals1.getRawPtr()), nnz);
293 memcpy(vals.getRawPtr(), valsA.getRawPtr(), nnz*
sizeof(SC));
295 vals = Array<SC>(valsA);
298 SC ZERO = Teuchos::ScalarTraits<SC>::zero();
300 SC A_rowsum = ZERO, F_rowsum = ZERO;
301 for(LO l = 0; l < (LO)inds.size(); l++)
302 A_rowsum += valsA[l];
304 if (lumping ==
false) {
305 for (
size_t j = 0; j < nnz; j++)
306 if (!filter[inds[j]])
313 for (
size_t j = 0; j < nnz; j++) {
314 if (filter[inds[j]]) {
315 if (inds[j] == row) {
322 diagExtra += vals[j];
332 if (diagIndex != -1) {
334 vals[diagIndex] += diagExtra;
335 if(dirichletThresh >= 0.0 && TST::real(vals[diagIndex]) <= dirichletThresh) {
338 for(LO l = 0; l < (LO)nnz; l++)
341 vals[diagIndex] = TST::one();
347 #ifndef ASSUME_DIRECT_ACCESS_TO_ROW 350 filteredA.replaceLocalValues(row, inds, vals);
355 for (
size_t j = 0; j < as<size_t> (indsG.size()); j++)
356 for (
size_t k = 0; k < blkSize; k++)
357 filter[indsG[j]*blkSize+k] = 0;
361 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
363 BuildNew(
const Matrix& A,
const GraphBase& G,
const bool lumping,
double dirichletThresh, Matrix& filteredA)
const {
364 using TST =
typename Teuchos::ScalarTraits<SC>;
365 SC zero = Teuchos::ScalarTraits<SC>::zero();
367 size_t blkSize = A.GetFixedBlockSize();
369 ArrayView<const LO> indsA;
370 ArrayView<const SC> valsA;
374 Array<char> filter(blkSize * G.
GetImportMap()->getNodeNumElements(), 0);
377 for (
size_t i = 0; i < numGRows; i++) {
380 for (
size_t j = 0; j < as<size_t>(indsG.size()); j++)
381 for (
size_t k = 0; k < blkSize; k++)
382 filter[indsG[j]*blkSize+k] = 1;
384 for (
size_t k = 0; k < blkSize; k++) {
385 LO row = i*blkSize + k;
387 A.getLocalRowView(row, indsA, valsA);
389 size_t nnz = indsA.size();
393 inds.resize(indsA.size());
394 vals.resize(valsA.size());
397 if (lumping ==
false) {
398 for (
size_t j = 0; j < nnz; j++)
399 if (filter[indsA[j]]) {
400 inds[numInds] = indsA[j];
401 vals[numInds] = valsA[j];
409 for (
size_t j = 0; j < nnz; j++) {
410 if (filter[indsA[j]]) {
411 inds[numInds] = indsA[j];
412 vals[numInds] = valsA[j];
415 if (inds[numInds] == row)
421 diagExtra += valsA[j];
429 if (diagIndex != -1) {
430 vals[diagIndex] += diagExtra;
431 if(dirichletThresh >= 0.0 && TST::real(vals[diagIndex]) <= dirichletThresh) {
437 vals[diagIndex] = TST::one();
442 inds.resize(numInds);
443 vals.resize(numInds);
449 filteredA.insertLocalValues(row, inds, vals);
453 for (
size_t j = 0; j < as<size_t> (indsG.size()); j++)
454 for (
size_t k = 0; k < blkSize; k++)
455 filter[indsG[j]*blkSize+k] = 0;
459 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
461 BuildNewUsingRootStencil(
const Matrix& A,
const GraphBase& G,
double dirichletThresh,
Level& currentLevel, Matrix& filteredA,
bool use_spread_lumping,
double DdomAllowGrowthRate,
double DdomCap)
const {
462 using TST =
typename Teuchos::ScalarTraits<SC>;
463 using Teuchos::arcp_const_cast;
464 SC ZERO = Teuchos::ScalarTraits<SC>::zero();
465 SC ONE = Teuchos::ScalarTraits<SC>::one();
466 LO INVALID = Teuchos::OrdinalTraits<LO>::invalid();
469 size_t blkSize = A.GetFixedBlockSize();
470 size_t numRows = A.getMap()->getNodeNumElements();
471 ArrayView<const LO> indsA;
472 ArrayView<const SC> valsA;
473 ArrayRCP<const size_t> rowptr;
474 ArrayRCP<const LO> inds;
475 ArrayRCP<const SC> vals_const;
481 RCP<CrsMatrix> filteredAcrs =
dynamic_cast<const CrsMatrixWrap*
>(&filteredA)->getCrsMatrix();
482 filteredAcrs->getAllValues(rowptr,inds,vals_const);
483 vals = arcp_const_cast<SC>(vals_const);
484 Array<bool> vals_dropped_indicator(vals.size(),
false);
488 RCP<Aggregates> aggregates = Get< RCP<Aggregates> > (currentLevel,
"Aggregates");
489 RCP<AmalgamationInfo> amalgInfo = Get< RCP<AmalgamationInfo> > (currentLevel,
"UnAmalgamationInfo");
490 LO numAggs = aggregates->GetNumAggregates();
491 Teuchos::ArrayRCP<const LO> vertex2AggId = aggregates->GetVertex2AggId()->getData(0);
494 RCP<const Map> rowMap = A.getRowMap();
495 RCP<const Map> colMap = A.getColMap();
500 Array<LO> diagIndex(numRows,INVALID);
501 Array<SC> diagExtra(numRows,ZERO);
505 Array<LO> ptr,nodes, unaggregated;
507 aggregates->ComputeNodesInAggregate(nodesInAgg.ptr, nodesInAgg.nodes, nodesInAgg.unaggregated);
508 LO graphNumCols = G.
GetImportMap()->getNodeNumElements();
509 Array<bool> filter(graphNumCols,
false);
512 for(LO i=0; i<nodesInAgg.unaggregated.size(); i++) {
513 for (LO m = 0; m < (LO)blkSize; m++) {
514 LO row = amalgInfo->ComputeLocalDOF(nodesInAgg.unaggregated[i],m);
515 if (row >= (LO)numRows)
continue;
516 size_t index_start = rowptr[row];
517 A.getLocalRowView(row, indsA, valsA);
518 for(LO k=0; k<(LO)indsA.size(); k++) {
519 if(row == indsA[k]) {
520 vals[index_start+k] = ONE;
524 vals[index_start+k] = ZERO;
530 std::vector<LO> badCount(numAggs,0);
534 for(LO i=0; i<numAggs; i++)
535 maxAggSize = std::max(maxAggSize,nodesInAgg.ptr[i+1] - nodesInAgg.ptr[i]);
542 size_t numNewDrops=0;
543 size_t numOldDrops=0;
544 size_t numFixedDiags=0;
545 size_t numSymDrops = 0;
547 for(LO i=0; i<numAggs; i++) {
548 LO numNodesInAggregate = nodesInAgg.ptr[i+1] - nodesInAgg.ptr[i];
549 if(numNodesInAggregate == 0)
continue;
552 LO root_node = INVALID;
553 for (LO k=nodesInAgg.ptr[i]; k < nodesInAgg.ptr[i+1]; k++) {
554 if(aggregates->IsRoot(nodesInAgg.nodes[k])) {
555 root_node = nodesInAgg.nodes[k];
break;
559 TEUCHOS_TEST_FOR_EXCEPTION(root_node == INVALID,
566 goodAggNeighbors.resize(0);
567 for(LO k=0; k<(LO) goodNodeNeighbors.size(); k++) {
568 goodAggNeighbors.push_back(vertex2AggId[goodNodeNeighbors[k]]);
575 badAggNeighbors.resize(0);
576 for(LO j = 0; j < (LO)blkSize; j++) {
577 LO row = amalgInfo->ComputeLocalDOF(root_node,j);
578 if (row >= (LO)numRows)
continue;
579 A.getLocalRowView(row, indsA, valsA);
580 for(LO k=0; k<(LO)indsA.size(); k++) {
581 if ( (indsA[k] < (LO)numRows) && (TST::magnitude(valsA[k]) != TST::magnitude(ZERO))) {
582 LO node = amalgInfo->ComputeLocalNode(indsA[k]);
583 LO agg = vertex2AggId[node];
584 if(!std::binary_search(goodAggNeighbors.begin(),goodAggNeighbors.end(),agg))
585 badAggNeighbors.push_back(agg);
594 for (LO k=nodesInAgg.ptr[i]; k < nodesInAgg.ptr[i+1]; k++) {
596 for (LO kk=0; kk < nodeNeighbors.size(); kk++) {
597 if ( (vertex2AggId[nodeNeighbors[kk]] >= 0) && (vertex2AggId[nodeNeighbors[kk]] < numAggs))
598 (badCount[vertex2AggId[nodeNeighbors[kk]]])++;
602 reallyBadAggNeighbors.resize(0);
603 for (LO k=0; k < (LO) badAggNeighbors.size(); k++) {
604 if (badCount[badAggNeighbors[k]] <= 1 ) reallyBadAggNeighbors.push_back(badAggNeighbors[k]);
606 for (LO k=nodesInAgg.ptr[i]; k < nodesInAgg.ptr[i+1]; k++) {
608 for (LO kk=0; kk < nodeNeighbors.size(); kk++) {
609 if ( (vertex2AggId[nodeNeighbors[kk]] >= 0) && (vertex2AggId[nodeNeighbors[kk]] < numAggs))
610 badCount[vertex2AggId[nodeNeighbors[kk]]] = 0;
616 for(LO b=0; b<(LO)reallyBadAggNeighbors.size(); b++) {
617 LO bad_agg = reallyBadAggNeighbors[b];
618 for (LO k=nodesInAgg.ptr[bad_agg]; k < nodesInAgg.ptr[bad_agg+1]; k++) {
619 LO bad_node = nodesInAgg.nodes[k];
620 for(LO j = 0; j < (LO)blkSize; j++) {
621 LO bad_row = amalgInfo->ComputeLocalDOF(bad_node,j);
622 if (bad_row >= (LO)numRows)
continue;
623 size_t index_start = rowptr[bad_row];
624 A.getLocalRowView(bad_row, indsA, valsA);
625 for(LO l = 0; l < (LO)indsA.size(); l++) {
626 if(indsA[l] < (LO)numRows && vertex2AggId[amalgInfo->ComputeLocalNode(indsA[l])] == i && vals_dropped_indicator[index_start+l] ==
false) {
627 vals_dropped_indicator[index_start + l] =
true;
628 vals[index_start + l] = ZERO;
629 diagExtra[bad_row] += valsA[l];
640 for(LO k=nodesInAgg.ptr[i]; k<nodesInAgg.ptr[i+1]; k++) {
641 LO row_node = nodesInAgg.nodes[k];
645 for (
size_t j = 0; j < as<size_t>(indsG.size()); j++)
646 filter[indsG[j]]=
true;
648 for (LO m = 0; m < (LO)blkSize; m++) {
649 LO row = amalgInfo->ComputeLocalDOF(row_node,m);
650 if (row >= (LO)numRows)
continue;
651 size_t index_start = rowptr[row];
652 A.getLocalRowView(row, indsA, valsA);
654 for(LO l = 0; l < (LO)indsA.size(); l++) {
655 int col_node = amalgInfo->ComputeLocalNode(indsA[l]);
656 bool is_good = filter[col_node];
657 if (indsA[l] == row) {
659 vals[index_start + l] = valsA[l];
664 if(vals_dropped_indicator[index_start +l] ==
true) {
665 if(is_good) numOldDrops++;
673 if(is_good && indsA[l] < (LO)numRows) {
674 int agg = vertex2AggId[col_node];
675 if(std::binary_search(reallyBadAggNeighbors.begin(),reallyBadAggNeighbors.end(),agg))
680 vals[index_start+l] = valsA[l];
683 if(!filter[col_node]) numOldDrops++;
685 diagExtra[row] += valsA[l];
686 vals[index_start+l]=ZERO;
687 vals_dropped_indicator[index_start+l]=
true;
694 for (
size_t j = 0; j < as<size_t>(indsG.size()); j++)
695 filter[indsG[j]]=
false;
700 if (!use_spread_lumping) {
702 for(LO row=0; row <(LO)numRows; row++) {
703 if (diagIndex[row] != INVALID) {
704 size_t index_start = rowptr[row];
705 size_t diagIndexInMatrix = index_start + diagIndex[row];
707 vals[diagIndexInMatrix] += diagExtra[row];
708 SC A_rowsum=ZERO, A_absrowsum = ZERO, F_rowsum = ZERO;
711 if( (dirichletThresh >= 0.0 && TST::real(vals[diagIndexInMatrix]) <= dirichletThresh) || TST::real(vals[diagIndexInMatrix]) == ZERO) {
714 A.getLocalRowView(row, indsA, valsA);
718 for(LO l = 0; l < (LO)indsA.size(); l++) {
719 A_rowsum += valsA[l];
720 A_absrowsum+=std::abs(valsA[l]);
722 for(LO l = 0; l < (LO)indsA.size(); l++)
723 F_rowsum += vals[index_start+l];
737 for(
size_t l=rowptr[row]; l<rowptr[row+1]; l++) {
740 vals[diagIndexInMatrix] = TST::one();
745 GetOStream(
Runtime0)<<
"WARNING: Row "<<row<<
" has no diagonal "<<std::endl;
751 for(LO row=0; row<(LO)numRows; row++) {
752 filteredA.replaceLocalValues(row, inds(rowptr[row],rowptr[row+1]-rowptr[row]), vals(rowptr[row],rowptr[row+1]-rowptr[row]));
754 if (use_spread_lumping) ExperimentalLumping(A, filteredA, DdomAllowGrowthRate, DdomCap);
756 size_t g_newDrops = 0, g_oldDrops = 0, g_fixedDiags=0;
758 MueLu_sumAll(A.getRowMap()->getComm(), numNewDrops, g_newDrops);
759 MueLu_sumAll(A.getRowMap()->getComm(), numOldDrops, g_oldDrops);
760 MueLu_sumAll(A.getRowMap()->getComm(), numFixedDiags, g_fixedDiags);
761 GetOStream(
Runtime0)<<
"Filtering out "<<g_newDrops<<
" edges, in addition to the "<<g_oldDrops<<
" edges dropped earlier"<<std::endl;
762 GetOStream(
Runtime0)<<
"Fixing "<< g_fixedDiags<<
" zero diagonal values" <<std::endl;
776 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
779 using TST =
typename Teuchos::ScalarTraits<SC>;
780 SC zero = TST::zero();
783 ArrayView<const LO> inds;
784 ArrayView<const SC> vals;
785 ArrayView<const LO> finds;
788 SC PosOffSum, NegOffSum, PosOffDropSum, NegOffDropSum;
789 SC diag, gamma, alpha;
790 LO NumPosKept, NumNegKept;
794 SC PosFilteredSum, NegFilteredSum;
797 SC rho = as<Scalar>(irho);
798 SC rho2 = as<Scalar>(irho2);
800 for (LO row = 0; row < (LO) A.getRowMap()->getNodeNumElements(); row++) {
801 noLumpDdom = as<Scalar>(10000.0);
814 ArrayView<const SC> tvals;
815 A.getLocalRowView(row, inds, vals);
816 size_t nnz = inds.size();
817 if (nnz == 0)
continue;
818 filteredA.getLocalRowView(row, finds, tvals);
820 fvals = ArrayView<SC>(
const_cast<SC*
>(tvals.getRawPtr()), nnz);
822 LO diagIndex = -1, fdiagIndex = -1;
824 PosOffSum=zero; NegOffSum=zero; PosOffDropSum=zero; NegOffDropSum=zero;
825 diag=zero; NumPosKept=0; NumNegKept=0;
828 for (
size_t j = 0; j < nnz; j++) {
829 if (inds[j] == row) {
834 if (TST::real(vals[j]) > TST::real(zero) ) PosOffSum += vals[j];
835 else NegOffSum += vals[j];
838 PosOffDropSum = PosOffSum;
839 NegOffDropSum = NegOffSum;
843 for (
size_t jj = 0; jj < (size_t) finds.size(); jj++) {
844 while( inds[j] != finds[jj] ) j++;
846 if (finds[jj] == row) fdiagIndex = jj;
848 if (TST::real(vals[j]) > TST::real(zero) ) {
849 PosOffDropSum -= fvals[jj];
850 if (TST::real(fvals[jj]) != TST::real(zero) ) NumPosKept++;
853 NegOffDropSum -= fvals[jj];
854 if (TST::real(fvals[jj]) != TST::real(zero) ) NumNegKept++;
860 if (TST::magnitude(diag) != TST::magnitude(zero) )
861 noLumpDdom = (PosOffSum - NegOffSum)/diag;
866 Target = rho*noLumpDdom;
867 if (TST::magnitude(Target) <= TST::magnitude(rho)) Target = rho2;
869 PosFilteredSum = PosOffSum - PosOffDropSum;
870 NegFilteredSum = NegOffSum - NegOffDropSum;
880 diag += PosOffDropSum;
883 gamma = -NegOffDropSum - PosFilteredSum;
885 if (TST::real(gamma) < TST::real(zero) ) {
893 if (fdiagIndex != -1) fvals[fdiagIndex] = diag;
895 for(LO jj = 0; jj < (LO)finds.size(); jj++) {
896 while( inds[j] != finds[jj] ) j++;
898 if ((j != diagIndex)&&(TST::real(vals[j]) > TST::real(zero) ) && (TST::magnitude(fvals[jj]) != TST::magnitude(zero)))
899 fvals[jj] = -gamma*(vals[j]/PosFilteredSum);
913 bool flipPosOffDiagsToNeg =
false;
925 if (( TST::real(diag) > TST::real(gamma)) &&
926 ( TST::real((-NegFilteredSum)/(diag - gamma)) <= TST::real(Target))) {
933 if (fdiagIndex != -1) fvals[fdiagIndex] = diag - gamma;
935 else if (NumNegKept > 0) {
941 numer = -NegFilteredSum - Target*(diag-gamma);
942 denom = gamma*(Target - TST::one());
964 if ( TST::magnitude(denom) < TST::magnitude(numer) ) alpha = TST::one();
965 else alpha = numer/denom;
966 if ( TST::real(alpha) < TST::real(zero)) alpha = zero;
967 if ( TST::real(diag) < TST::real((one-alpha)*gamma) ) alpha = TST::one();
971 if (fdiagIndex != -1) fvals[fdiagIndex] = diag - (one-alpha)*gamma;
981 SC temp = (NegFilteredSum+alpha*gamma)/NegFilteredSum;
983 for(LO jj = 0; jj < (LO)finds.size(); jj++) {
984 while( inds[j] != finds[jj] ) j++;
986 if ( (jj != fdiagIndex)&&(TST::magnitude(fvals[jj]) != TST::magnitude(zero) ) &&
987 ( TST::real(vals[j]) < TST::real(zero) ) )
988 fvals[jj] = temp*vals[j];
994 if (NumPosKept > 0) {
997 flipPosOffDiagsToNeg =
true;
1000 for(LO jj = 0; jj < (LO)finds.size(); jj++) {
1001 while( inds[j] != finds[jj] ) j++;
1003 if ( (j != diagIndex)&&(TST::magnitude(fvals[jj]) != TST::magnitude(zero) ) &&
1004 (TST::real(vals[j]) > TST::real(zero) ))
1005 fvals[jj] = -gamma/( (SC) NumPosKept);
1011 if (!flipPosOffDiagsToNeg) {
1016 for(LO jj = 0; jj < (LO)finds.size(); jj++) {
1017 while( inds[j] != finds[jj] ) j++;
1019 if ((jj != fdiagIndex)&& (TST::real(vals[j]) > TST::real(zero))) fvals[jj] = zero;
1030 #endif // MUELU_FILTEREDAFACTORY_DEF_HPP void sort_and_unique(T &array)
#define MueLu_sumAll(rcpComm, in, out)
void DeclareInput(Level ¤tLevel) const
Input.
virtual size_t getNodeMaxNumRowEntries() const =0
virtual size_t GetNodeNumVertices() const =0
Return number of vertices owned by the calling node.
Timer to be used in factories. Similar to Monitor but with additional timers.
One-liner description of what is happening.
Namespace for MueLu classes and methods.
void ExperimentalLumping(const Matrix &A, Matrix &filteredA, double rho, double rho2) const
static bool MapsAreNested(const Xpetra::Map< LocalOrdinal, GlobalOrdinal, Node > &rowMap, const Xpetra::Map< LocalOrdinal, GlobalOrdinal, Node > &colMap)
#define MUELU_FILTEREDAFACTORY_LOTS_OF_PRINTING
Class that holds all level-specific information.
void Build(Level ¤tLevel) const
Build method.
void BuildNew(const Matrix &A, const GraphBase &G, const bool lumping, double dirichletThresh, Matrix &filteredA) const
void BuildNewUsingRootStencil(const Matrix &A, const GraphBase &G, double dirichletThresh, Level ¤tLevel, Matrix &filteredA, bool use_spread_lumping, double DdomAllowGrowthRate, double DdomCap) const
RCP< const ParameterList > GetValidParameterList() const
Return a const parameter list of valid parameters that setParameterList() will accept.
MueLu representation of a graph.
void BuildReuse(const Matrix &A, const GraphBase &G, const bool lumping, double dirichletThresh, Matrix &filteredA) const
#define SET_VALID_ENTRY(name)
virtual const RCP< const Map > GetImportMap() const =0
Exception throws to report errors in the internal logical of the program.
virtual Teuchos::ArrayView< const LocalOrdinal > getNeighborVertices(LocalOrdinal v) const =0
Return the list of vertices adjacent to the vertex 'v'.