146 typedef Teuchos::ScalarTraits<SC> STS;
148 if (predrop_ != Teuchos::null)
149 GetOStream(
Parameters0) << predrop_->description();
151 RCP<Matrix> A = Get< RCP<Matrix> >(currentLevel,
"A");
153 const ParameterList & pL = GetParameterList();
155 LO nPDEs = A->GetFixedBlockSize();
157 RCP< MultiVector > testVecs;
158 RCP< MultiVector > nearNull;
161 testVecs = Xpetra::IO<SC,LO,GO,Node>::ReadMultiVector(
"TpetraTVecs.mm", A->getRowMap());
163 size_t numRandom= as<size_t>(pL.get<
int>(
"aggregation: number of random vectors"));
164 testVecs = MultiVectorFactory::Build(A->getRowMap(), numRandom,
true);
167 testVecs->randomize();
168 for (
size_t kk = 0; kk < testVecs->getNumVectors(); kk++ ) {
169 Teuchos::ArrayRCP< Scalar > curVec = testVecs->getDataNonConst(kk);
170 for (
size_t ii = kk; ii < as<size_t>(A->getRowMap()->getLocalNumElements()); ii++ ) curVec[ii] = Teuchos::ScalarTraits<SC>::magnitude(curVec[ii]);
172 nearNull = MultiVectorFactory::Build(A->getRowMap(), nPDEs,
true);
175 for (
size_t kk = 0; kk < nearNull->getNumVectors(); kk++ ) {
176 Teuchos::ArrayRCP< Scalar > curVec = nearNull->getDataNonConst(kk);
177 for (
size_t ii = kk; ii < as<size_t>(A->getRowMap()->getLocalNumElements()); ii += nearNull->getNumVectors() ) curVec[ii] = Teuchos::ScalarTraits<Scalar>::one();
180 RCP< MultiVector > zeroVec_TVecs;
181 RCP< MultiVector > zeroVec_Null;
183 zeroVec_TVecs = MultiVectorFactory::Build(A->getRowMap(), testVecs->getNumVectors(),
true);
184 zeroVec_Null = MultiVectorFactory::Build(A->getRowMap(), nPDEs,
true);
185 zeroVec_TVecs->putScalar(Teuchos::ScalarTraits<Scalar>::zero());
186 zeroVec_Null->putScalar( Teuchos::ScalarTraits<Scalar>::zero());
188 size_t nInvokeSmoother=as<size_t>(pL.get<
int>(
"aggregation: number of times to pre or post smooth"));
190 RCP<SmootherBase> preSmoo = currentLevel.
Get< RCP<SmootherBase> >(
"PreSmoother");
191 for (
size_t ii = 0; ii < nInvokeSmoother; ii++) preSmoo->Apply(*testVecs,*zeroVec_TVecs,
false);
192 for (
size_t ii = 0; ii < nInvokeSmoother; ii++) preSmoo->Apply(*nearNull,*zeroVec_Null,
false);
194 else if (currentLevel.
IsAvailable(
"PostSmoother")) {
195 RCP<SmootherBase> postSmoo = currentLevel.
Get< RCP<SmootherBase> >(
"PostSmoother");
196 for (
size_t ii = 0; ii < nInvokeSmoother; ii++) postSmoo->Apply(*testVecs,*zeroVec_TVecs,
false);
197 for (
size_t ii = 0; ii < nInvokeSmoother; ii++) postSmoo->Apply(*nearNull, *zeroVec_Null,
false);
202 Teuchos::ArrayRCP<Scalar> penaltyPolyCoef(5);
203 Teuchos::ArrayView<const double> inputPolyCoef;
211 if(pL.isParameter(
"aggregation: penalty parameters") && pL.get<Teuchos::Array<double> >(
"aggregation: penalty parameters").size() > 0) {
212 if (pL.get<Teuchos::Array<double> >(
"aggregation: penalty parameters").size() > penaltyPolyCoef.size())
213 TEUCHOS_TEST_FOR_EXCEPTION(
true,
Exceptions::RuntimeError,
"Number of penalty parameters must be " << penaltyPolyCoef.size() <<
" or less");
214 inputPolyCoef = pL.get<Teuchos::Array<double> >(
"aggregation: penalty parameters")();
216 for (
size_t i = 0; i < as<size_t>(inputPolyCoef.size()) ; i++) penaltyPolyCoef[i] = as<Scalar>(inputPolyCoef[i]);
217 for (
size_t i = as<size_t>(inputPolyCoef.size()); i < as<size_t>(penaltyPolyCoef.size()); i++) penaltyPolyCoef[i] = Teuchos::ScalarTraits<Scalar>::zero();
221 RCP<GraphBase> filteredGraph;
222 badGuysCoalesceDrop(*A, penaltyPolyCoef, nPDEs, *testVecs, *nearNull, filteredGraph);
228 FILE* fp = fopen(
"codeOutput",
"w");
229 fprintf(fp,
"%d %d %d\n",(
int) filteredGraph->GetNodeNumVertices(),(
int) filteredGraph->GetNodeNumVertices(),
230 (
int) filteredGraph->GetNodeNumEdges());
231 for (
size_t i = 0; i < filteredGraph->GetNodeNumVertices(); i++) {
232 ArrayView<const LO> inds = filteredGraph->getNeighborVertices(as<LO>(i));
233 for (
size_t j = 0; j < as<size_t>(inds.size()); j++) {
234 fprintf(fp,
"%d %d 1.00e+00\n",(
int) i+1,(
int) inds[j]+1);
241 Set<bool>(currentLevel,
"Filtering", (threshold != STS::zero()));
242 Set(currentLevel,
"Graph", filteredGraph);
243 Set(currentLevel,
"DofsPerNode", 1);
286 GO numMyNnz = Teuchos::as<GO>(Amat.getLocalNumEntries());
287 size_t nLoc = Amat.getRowMap()->getLocalNumElements();
289 size_t nBlks = nLoc/nPDEs;
290 if (nBlks*nPDEs != nLoc )
293 Teuchos::ArrayRCP<LO> newRowPtr(nBlks+1);
294 Teuchos::ArrayRCP<LO> newCols(numMyNnz);
296 Teuchos::ArrayRCP<LO> bcols(nBlks);
297 Teuchos::ArrayRCP<bool> keepOrNot(nBlks);
301 LO maxNzPerRow = 200;
302 Teuchos::ArrayRCP<Scalar> penalties(maxNzPerRow);
305 Teuchos::ArrayRCP<bool> keepStatus(nBlks,
true);
306 Teuchos::ArrayRCP<LO> bColList(nBlks);
316 Teuchos::ArrayRCP<bool> alreadyOnBColList(nBlks,
false);
321 Teuchos::ArrayRCP<bool> boundaryNodes(nBlks,
false);
324 for (LO i = 0; i < maxNzPerRow; i++)
331 LO nzTotal = 0, numBCols = 0, row = -1, Nbcols, bcol;
335 for (LO i = 0; i < as<LO>(nBlks); i++) {
336 newRowPtr[i+1] = newRowPtr[i];
337 for (LO j = 0; j < nPDEs; j++) {
340 Teuchos::ArrayView<const LocalOrdinal> indices;
341 Teuchos::ArrayView<const Scalar> vals;
343 Amat.getLocalRowView(row, indices, vals);
345 if (indices.size() > maxNzPerRow) {
346 LO oldSize = maxNzPerRow;
347 maxNzPerRow = indices.size() + 100;
348 penalties.resize(as<size_t>(maxNzPerRow),0.0);
349 for (LO k = oldSize; k < maxNzPerRow; k++)
356 badGuysDropfunc(row, indices, vals, testVecs, nPDEs, penalties, nearNull, bcols,keepOrNot,Nbcols,nLoc);
357 for (LO k=0; k < Nbcols; k++) {
362 if (alreadyOnBColList[bcol] ==
false) {
363 bColList[numBCols++] = bcol;
364 alreadyOnBColList[bcol] =
true;
368 if (keepOrNot[k] ==
false) keepStatus[bcol] =
false;
376 if ( numBCols < 2) boundaryNodes[i] =
true;
377 for (LO j=0; j < numBCols; j++) {
379 if (keepStatus[bcol] ==
true) {
380 newCols[nzTotal] = bColList[j];
382 nzTotal = nzTotal + 1;
384 keepStatus[bcol] =
true;
385 alreadyOnBColList[bcol] =
false;
393 Teuchos::ArrayRCP<LO> finalCols(nzTotal);
394 for (LO i = 0; i < nzTotal; i++) finalCols[i] = newCols[i];
399 RCP<const Map> rowMap = Amat.getRowMap();
401 LO nAmalgNodesOnProc = rowMap->getLocalNumElements()/nPDEs;
402 Teuchos::Array<GO> nodalGIDs(nAmalgNodesOnProc);
403 typename Teuchos::ScalarTraits<Scalar>::coordinateType temp;
404 for (
size_t i = 0; i < as<size_t>(nAmalgNodesOnProc); i++ ) {
405 GO gid = rowMap->getGlobalElement(i*nPDEs);
406 temp = ((
typename Teuchos::ScalarTraits<Scalar>::coordinateType) (gid))/((
typename Teuchos::ScalarTraits<Scalar>::coordinateType) (nPDEs));
407 nodalGIDs[i] = as<GO>(floor(temp));
409 GO nAmalgNodesGlobal = rowMap->getGlobalNumElements();
410 GO nBlkGlobal = nAmalgNodesGlobal/nPDEs;
411 if (nBlkGlobal*nPDEs != nAmalgNodesGlobal)
414 Teuchos::RCP<Map> AmalgRowMap = MapFactory::Build(rowMap->lib(), nBlkGlobal,
415 nodalGIDs(),0,rowMap->getComm());
417 filteredGraph = rcp(
new LWGraph(newRowPtr, finalCols, AmalgRowMap, AmalgRowMap,
"thresholded graph of A"));
418 filteredGraph->SetBoundaryNodeMap(boundaryNodes);
423 void SmooVecCoalesceDropFactory<Scalar, LocalOrdinal, GlobalOrdinal, Node>::badGuysDropfunc(LO row,
const Teuchos::ArrayView<const LocalOrdinal>& cols,
const Teuchos::ArrayView<const Scalar>& vals,
const MultiVector& testVecs, LO nPDEs, Teuchos::ArrayRCP<Scalar> & penalties,
const MultiVector& nearNull, Teuchos::ArrayRCP<LO>& Bcols, Teuchos::ArrayRCP<bool>& keepOrNot, LO &Nbcols, LO nLoc)
const {
424 using TST=Teuchos::ScalarTraits<Scalar>;
426 LO nLeng = cols.size();
427 typename TST::coordinateType temp;
428 temp = ((
typename TST::coordinateType) (row))/((
typename TST::coordinateType) (nPDEs));
429 LO blkRow = as<LO>(floor(temp));
430 Teuchos::ArrayRCP<Scalar> badGuy( nLeng, 0.0);
431 Teuchos::ArrayRCP<Scalar> subNull(nLeng, 0.0);
445 for (LO i = 0; i < nLeng; i++) keepOrNot[i] =
false;
449 LO rowDof = row - blkRow*nPDEs;
450 Teuchos::ArrayRCP< const Scalar > oneNull = nearNull.getData( as<size_t>(rowDof));
452 for (LO i = 0; i < nLeng; i++) {
453 if ((cols[i] < nLoc ) && (TST::magnitude(vals[i]) != 0.0)) {
454 temp = ((
typename TST::coordinateType) (cols[i]))/((
typename TST::coordinateType) (nPDEs));
455 LO colDof = cols[i] - (as<LO>(floor( temp )))*nPDEs;
456 if (colDof == rowDof) {
457 Bcols[ Nbcols] = (cols[i] - colDof)/nPDEs;
458 subNull[Nbcols] = oneNull[cols[i]];
460 if (cols[i] != row) {
461 Scalar worstRatio = -TST::one();
462 Scalar targetRatio = subNull[Nbcols]/oneNull[row];
464 for (
size_t kk = 0; kk < testVecs.getNumVectors(); kk++ ) {
465 Teuchos::ArrayRCP< const Scalar > curVec = testVecs.getData(kk);
466 actualRatio = curVec[cols[i]]/curVec[row];
467 if (TST::magnitude(actualRatio - targetRatio) > TST::magnitude(worstRatio)) {
468 badGuy[Nbcols] = actualRatio;
469 worstRatio = Teuchos::ScalarTraits<SC>::magnitude(actualRatio - targetRatio);
474 badGuy[ Nbcols] = 1.;
475 keepOrNot[Nbcols] =
true;
486 Bcols[ Nbcols] = (row - rowDof)/nPDEs;
487 subNull[ Nbcols] = 1.;
488 badGuy[ Nbcols] = 1.;
489 keepOrNot[Nbcols] =
true;
494 Scalar currentRP = oneNull[row]*oneNull[row];
495 Scalar currentRTimesBadGuy= oneNull[row]*badGuy[diagInd];
496 Scalar currentScore = penalties[0];
507 LO nKeep = 1, flag = 1, minId;
508 Scalar minFit, minFitRP = 0., minFitRTimesBadGuy = 0.;
509 Scalar newRP, newRTimesBadGuy;
518 for (LO i=0; i < Nbcols; i++) {
519 if (keepOrNot[i] ==
false) {
521 newRP = currentRP + subNull[i]*subNull[i];
522 newRTimesBadGuy= currentRTimesBadGuy + subNull[i]*badGuy[i];
523 Scalar ratio = newRTimesBadGuy/newRP;
526 for (LO k=0; k < Nbcols; k++) {
527 if (keepOrNot[k] ==
true) {
528 Scalar diff = badGuy[k] - ratio*subNull[k];
529 newFit = newFit + diff*diff;
532 if (Teuchos::ScalarTraits<SC>::magnitude(newFit) < Teuchos::ScalarTraits<SC>::magnitude(minFit)) {
536 minFitRTimesBadGuy= newRTimesBadGuy;
538 keepOrNot[i] =
false;
541 if (minId == -1) flag = 0;
543 minFit = sqrt(minFit);
544 Scalar newScore = penalties[nKeep] + minFit;
545 if (Teuchos::ScalarTraits<SC>::magnitude(newScore) < Teuchos::ScalarTraits<SC>::magnitude(currentScore)) {
547 keepOrNot[minId]=
true;
548 currentScore = newScore;
549 currentRP = minFitRP;
550 currentRTimesBadGuy= minFitRTimesBadGuy;