MueLu  Version of the Day
MueLu_TentativePFactory_def.hpp
Go to the documentation of this file.
1 // @HEADER
2 //
3 // ***********************************************************************
4 //
5 // MueLu: A package for multigrid based preconditioning
6 // Copyright 2012 Sandia Corporation
7 //
8 // Under the terms of Contract DE-AC04-94AL85000 with Sandia Corporation,
9 // the U.S. Government retains certain rights in this software.
10 //
11 // Redistribution and use in source and binary forms, with or without
12 // modification, are permitted provided that the following conditions are
13 // met:
14 //
15 // 1. Redistributions of source code must retain the above copyright
16 // notice, this list of conditions and the following disclaimer.
17 //
18 // 2. Redistributions in binary form must reproduce the above copyright
19 // notice, this list of conditions and the following disclaimer in the
20 // documentation and/or other materials provided with the distribution.
21 //
22 // 3. Neither the name of the Corporation nor the names of the
23 // contributors may be used to endorse or promote products derived from
24 // this software without specific prior written permission.
25 //
26 // THIS SOFTWARE IS PROVIDED BY SANDIA CORPORATION "AS IS" AND ANY
27 // EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
28 // IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR
29 // PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL SANDIA CORPORATION OR THE
30 // CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL,
31 // EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO,
32 // PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR
33 // PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF
34 // LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING
35 // NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
36 // SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
37 //
38 // Questions? Contact
39 // Jonathan Hu (jhu@sandia.gov)
40 // Andrey Prokopenko (aprokop@sandia.gov)
41 // Ray Tuminaro (rstumin@sandia.gov)
42 //
43 // ***********************************************************************
44 //
45 // @HEADER
46 #ifndef MUELU_TENTATIVEPFACTORY_DEF_HPP
47 #define MUELU_TENTATIVEPFACTORY_DEF_HPP
48 
49 #include <Xpetra_MapFactory.hpp>
50 #include <Xpetra_Map.hpp>
51 #include <Xpetra_CrsMatrix.hpp>
52 #include <Xpetra_Matrix.hpp>
53 #include <Xpetra_MultiVector.hpp>
54 #include <Xpetra_MultiVectorFactory.hpp>
55 #include <Xpetra_VectorFactory.hpp>
56 #include <Xpetra_Import.hpp>
57 #include <Xpetra_ImportFactory.hpp>
58 #include <Xpetra_CrsMatrixWrap.hpp>
59 #include <Xpetra_StridedMap.hpp>
60 #include <Xpetra_StridedMapFactory.hpp>
61 
63 
64 #include "MueLu_Aggregates.hpp"
65 #include "MueLu_AmalgamationFactory.hpp"
66 #include "MueLu_AmalgamationInfo.hpp"
67 #include "MueLu_CoarseMapFactory.hpp"
68 #include "MueLu_MasterList.hpp"
69 #include "MueLu_Monitor.hpp"
70 #include "MueLu_NullspaceFactory.hpp"
71 #include "MueLu_PerfUtils.hpp"
72 #include "MueLu_Utilities.hpp"
73 
74 namespace MueLu {
75 
76  template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
78  RCP<ParameterList> validParamList = rcp(new ParameterList());
79 
80 #define SET_VALID_ENTRY(name) validParamList->setEntry(name, MasterList::getEntry(name))
81  SET_VALID_ENTRY("tentative: calculate qr");
82  SET_VALID_ENTRY("tentative: build coarse coordinates");
83 #undef SET_VALID_ENTRY
84 
85  validParamList->set< RCP<const FactoryBase> >("A", Teuchos::null, "Generating factory of the matrix A");
86  validParamList->set< RCP<const FactoryBase> >("Aggregates", Teuchos::null, "Generating factory of the aggregates");
87  validParamList->set< RCP<const FactoryBase> >("Nullspace", Teuchos::null, "Generating factory of the nullspace");
88  validParamList->set< RCP<const FactoryBase> >("UnAmalgamationInfo", Teuchos::null, "Generating factory of UnAmalgamationInfo");
89  validParamList->set< RCP<const FactoryBase> >("CoarseMap", Teuchos::null, "Generating factory of the coarse map");
90  validParamList->set< RCP<const FactoryBase> >("Coordinates", Teuchos::null, "Generating factory of the coordinates");
91 
92  // Make sure we don't recursively validate options for the matrixmatrix kernels
93  ParameterList norecurse;
94  norecurse.disableRecursiveValidation();
95  validParamList->set<ParameterList> ("matrixmatrix: kernel params", norecurse, "MatrixMatrix kernel parameters");
96 
97  return validParamList;
98  }
99 
100  template <class Scalar,class LocalOrdinal, class GlobalOrdinal, class Node>
102 
103  const ParameterList& pL = GetParameterList();
104 
105  Input(fineLevel, "A");
106  Input(fineLevel, "Aggregates");
107  Input(fineLevel, "Nullspace");
108  Input(fineLevel, "UnAmalgamationInfo");
109  Input(fineLevel, "CoarseMap");
110  if( fineLevel.GetLevelID() == 0 &&
111  fineLevel.IsAvailable("Coordinates", NoFactory::get()) && // we have coordinates (provided by user app)
112  pL.get<bool>("tentative: build coarse coordinates") ) { // and we want coordinates on other levels
113  bTransferCoordinates_ = true; // then set the transfer coordinates flag to true
114  Input(fineLevel, "Coordinates");
115  } else if (bTransferCoordinates_) {
116  Input(fineLevel, "Coordinates");
117  }
118  }
119 
120  template <class Scalar,class LocalOrdinal, class GlobalOrdinal, class Node>
122  return BuildP(fineLevel, coarseLevel);
123  }
124 
125  template <class Scalar,class LocalOrdinal, class GlobalOrdinal, class Node>
127  FactoryMonitor m(*this, "Build", coarseLevel);
128  typedef typename Teuchos::ScalarTraits<SC>::magnitudeType real_type;
129  typedef typename Xpetra::MultiVector<real_type,LO,GO,NO> RealValuedMultiVector;
130 
131  RCP<Matrix> A = Get< RCP<Matrix> > (fineLevel, "A");
132  RCP<Aggregates> aggregates = Get< RCP<Aggregates> > (fineLevel, "Aggregates");
133  RCP<AmalgamationInfo> amalgInfo = Get< RCP<AmalgamationInfo> > (fineLevel, "UnAmalgamationInfo");
134  RCP<MultiVector> fineNullspace = Get< RCP<MultiVector> > (fineLevel, "Nullspace");
135  RCP<const Map> coarseMap = Get< RCP<const Map> > (fineLevel, "CoarseMap");
136  RCP<RealValuedMultiVector> fineCoords;
137  if(bTransferCoordinates_) {
138  fineCoords = Get< RCP<RealValuedMultiVector> >(fineLevel, "Coordinates");
139  }
140 
141  TEUCHOS_TEST_FOR_EXCEPTION(A->getRowMap()->getNodeNumElements() != fineNullspace->getMap()->getNodeNumElements(),
142  Exceptions::RuntimeError,"MueLu::TentativePFactory::MakeTentative: Size mismatch between A and Nullspace");
143 
144  RCP<Matrix> Ptentative;
145  RCP<MultiVector> coarseNullspace;
146  RCP<RealValuedMultiVector> coarseCoords;
147 
148  if(bTransferCoordinates_) {
149  //*** Create the coarse coordinates ***
150  // First create the coarse map and coarse multivector
151  ArrayView<const GO> elementAList = coarseMap->getNodeElementList();
152  LO blkSize = 1;
153  if (rcp_dynamic_cast<const StridedMap>(coarseMap) != Teuchos::null)
154  blkSize = rcp_dynamic_cast<const StridedMap>(coarseMap)->getFixedBlockSize();
155  GO indexBase = coarseMap->getIndexBase();
156  size_t numNodes = elementAList.size() / blkSize;
157  Array<GO> nodeList(numNodes);
158  const int numDimensions = fineCoords->getNumVectors();
159 
160  for (LO i = 0; i < Teuchos::as<LO>(numNodes); i++) {
161  nodeList[i] = (elementAList[i*blkSize]-indexBase)/blkSize + indexBase;
162  }
163  RCP<const Map> coarseCoordsMap = MapFactory::Build(fineCoords->getMap()->lib(),
164  Teuchos::OrdinalTraits<Xpetra::global_size_t>::invalid(),
165  nodeList,
166  indexBase,
167  fineCoords->getMap()->getComm());
168  coarseCoords = RealValuedMultiVectorFactory::Build(coarseCoordsMap, numDimensions);
169 
170  // Create overlapped fine coordinates to reduce global communication
171  RCP<RealValuedMultiVector> ghostedCoords;
172  if (aggregates->AggregatesCrossProcessors()) {
173  RCP<const Map> aggMap = aggregates->GetMap();
174  RCP<const Import> importer = ImportFactory::Build(fineCoords->getMap(), aggMap);
175 
176  ghostedCoords = RealValuedMultiVectorFactory::Build(aggMap, numDimensions);
177  ghostedCoords->doImport(*fineCoords, *importer, Xpetra::INSERT);
178  } else {
179  ghostedCoords = fineCoords;
180  }
181 
182  // Get some info about aggregates
183  int myPID = coarseCoordsMap->getComm()->getRank();
184  LO numAggs = aggregates->GetNumAggregates();
185  ArrayRCP<LO> aggSizes = aggregates->ComputeAggregateSizes();
186  const ArrayRCP<const LO> vertex2AggID = aggregates->GetVertex2AggId()->getData(0);
187  const ArrayRCP<const LO> procWinner = aggregates->GetProcWinner()->getData(0);
188 
189  // Fill in coarse coordinates
190  for (int dim = 0; dim < numDimensions; ++dim) {
191  ArrayRCP<const real_type> fineCoordsData = ghostedCoords->getData(dim);
192  ArrayRCP<real_type> coarseCoordsData = coarseCoords->getDataNonConst(dim);
193 
194  for (LO lnode = 0; lnode < Teuchos::as<LO>(numNodes); lnode++) {
195  if (procWinner[lnode] == myPID &&
196  lnode < vertex2AggID.size() &&
197  lnode < fineCoordsData.size() &&
198  vertex2AggID[lnode] < coarseCoordsData.size() &&
199  Teuchos::ScalarTraits<double>::isnaninf(fineCoordsData[lnode]) == false) {
200  coarseCoordsData[vertex2AggID[lnode]] += fineCoordsData[lnode];
201  }
202  }
203  for (LO agg = 0; agg < numAggs; agg++) {
204  coarseCoordsData[agg] /= aggSizes[agg];
205  }
206  }
207  }
208 
209  if (!aggregates->AggregatesCrossProcessors())
210  BuildPuncoupled(A, aggregates, amalgInfo, fineNullspace, coarseMap, Ptentative, coarseNullspace,coarseLevel.GetLevelID());
211  else
212  BuildPcoupled (A, aggregates, amalgInfo, fineNullspace, coarseMap, Ptentative, coarseNullspace);
213 
214 
215 
216 
217  // If available, use striding information of fine level matrix A for range
218  // map and coarseMap as domain map; otherwise use plain range map of
219  // Ptent = plain range map of A for range map and coarseMap as domain map.
220  // NOTE:
221  // The latter is not really safe, since there is no striding information
222  // for the range map. This is not really a problem, since striding
223  // information is always available on the intermedium levels and the
224  // coarsest levels.
225  if (A->IsView("stridedMaps") == true)
226  Ptentative->CreateView("stridedMaps", A->getRowMap("stridedMaps"), coarseMap);
227  else
228  Ptentative->CreateView("stridedMaps", Ptentative->getRangeMap(), coarseMap);
229 
230  if(bTransferCoordinates_) {
231  Set(coarseLevel, "Coordinates", coarseCoords);
232  }
233  Set(coarseLevel, "Nullspace", coarseNullspace);
234  Set(coarseLevel, "P", Ptentative);
235 
236  if (IsPrint(Statistics2)) {
237  RCP<ParameterList> params = rcp(new ParameterList());
238  params->set("printLoadBalancingInfo", true);
239  GetOStream(Statistics2) << PerfUtils::PrintMatrixInfo(*Ptentative, "Ptent", params);
240  }
241  }
242 
243  template <class Scalar,class LocalOrdinal, class GlobalOrdinal, class Node>
245  BuildPuncoupled(RCP<Matrix> A, RCP<Aggregates> aggregates, RCP<AmalgamationInfo> amalgInfo, RCP<MultiVector> fineNullspace,
246  RCP<const Map> coarseMap, RCP<Matrix>& Ptentative, RCP<MultiVector>& coarseNullspace, const int levelID) const {
247  RCP<const Map> rowMap = A->getRowMap();
248  RCP<const Map> colMap = A->getColMap();
249 
250  const size_t numRows = rowMap->getNodeNumElements();
251 
252  typedef Teuchos::ScalarTraits<SC> STS;
253  typedef typename STS::magnitudeType Magnitude;
254  const SC zero = STS::zero();
255  const SC one = STS::one();
256  const LO INVALID = Teuchos::OrdinalTraits<LO>::invalid();
257 
258  const GO numAggs = aggregates->GetNumAggregates();
259  const size_t NSDim = fineNullspace->getNumVectors();
260 
261  // Aggregates map is based on the amalgamated column map
262  // We can skip global-to-local conversion if LIDs in row map are
263  // same as LIDs in column map
264  bool goodMap = isGoodMap(*rowMap, *colMap);
265 
266  // Create a lookup table to determine the rows (fine DOFs) that belong to a given aggregate.
267  // aggStart is a pointer into aggToRowMapLO
268  // aggStart[i]..aggStart[i+1] are indices into aggToRowMapLO
269  // aggToRowMapLO[aggStart[i]]..aggToRowMapLO[aggStart[i+1]-1] are the DOFs in aggregate i
270  ArrayRCP<LO> aggStart;
271  ArrayRCP<LO> aggToRowMapLO;
272  ArrayRCP<GO> aggToRowMapGO;
273  if (goodMap) {
274  amalgInfo->UnamalgamateAggregatesLO(*aggregates, aggStart, aggToRowMapLO);
275  GetOStream(Runtime1) << "Column map is consistent with the row map, good." << std::endl;
276 
277  } else {
278  amalgInfo->UnamalgamateAggregates(*aggregates, aggStart, aggToRowMapGO);
279  GetOStream(Warnings0) << "Column map is not consistent with the row map\n"
280  << "using GO->LO conversion with performance penalty" << std::endl;
281  }
282 
283  coarseNullspace = MultiVectorFactory::Build(coarseMap, NSDim);
284 
285  const ParameterList& pL = GetParameterList();
286  const bool &doQRStep = pL.get<bool>("tentative: calculate qr");
287 
288  // Pull out the nullspace vectors so that we can have random access.
289  ArrayRCP<ArrayRCP<const SC> > fineNS (NSDim);
290  ArrayRCP<ArrayRCP<SC> > coarseNS(NSDim);
291  for (size_t i = 0; i < NSDim; i++) {
292  fineNS[i] = fineNullspace->getData(i);
293  if (coarseMap->getNodeNumElements() > 0)
294  coarseNS[i] = coarseNullspace->getDataNonConst(i);
295  }
296 
297  size_t nnzEstimate = numRows * NSDim;
298 
299  // Time to construct the matrix and fill in the values
300  Ptentative = rcp(new CrsMatrixWrap(rowMap, coarseMap, 0, Xpetra::StaticProfile));
301  RCP<CrsMatrix> PtentCrs = rcp_dynamic_cast<CrsMatrixWrap>(Ptentative)->getCrsMatrix();
302 
303  ArrayRCP<size_t> iaPtent;
304  ArrayRCP<LO> jaPtent;
305  ArrayRCP<SC> valPtent;
306 
307  PtentCrs->allocateAllValues(nnzEstimate, iaPtent, jaPtent, valPtent);
308 
309  ArrayView<size_t> ia = iaPtent();
310  ArrayView<LO> ja = jaPtent();
311  ArrayView<SC> val = valPtent();
312 
313  ia[0] = 0;
314  for (size_t i = 1; i <= numRows; i++)
315  ia[i] = ia[i-1] + NSDim;
316 
317  for (size_t j = 0; j < nnzEstimate; j++) {
318  ja [j] = INVALID;
319  val[j] = zero;
320  }
321 
322 
323  if (doQRStep) {
325  // Standard aggregate-wise QR //
327  for (GO agg = 0; agg < numAggs; agg++) {
328  LO aggSize = aggStart[agg+1] - aggStart[agg];
329 
330  Xpetra::global_size_t offset = agg*NSDim;
331 
332  // Extract the piece of the nullspace corresponding to the aggregate, and
333  // put it in the flat array, "localQR" (in column major format) for the
334  // QR routine.
335  Teuchos::SerialDenseMatrix<LO,SC> localQR(aggSize, NSDim);
336  if (goodMap) {
337  for (size_t j = 0; j < NSDim; j++)
338  for (LO k = 0; k < aggSize; k++)
339  localQR(k,j) = fineNS[j][aggToRowMapLO[aggStart[agg]+k]];
340  } else {
341  for (size_t j = 0; j < NSDim; j++)
342  for (LO k = 0; k < aggSize; k++)
343  localQR(k,j) = fineNS[j][rowMap->getLocalElement(aggToRowMapGO[aggStart[agg]+k])];
344  }
345 
346  // Test for zero columns
347  for (size_t j = 0; j < NSDim; j++) {
348  bool bIsZeroNSColumn = true;
349 
350  for (LO k = 0; k < aggSize; k++)
351  if (localQR(k,j) != zero)
352  bIsZeroNSColumn = false;
353 
354  TEUCHOS_TEST_FOR_EXCEPTION(bIsZeroNSColumn == true, Exceptions::RuntimeError,
355  "MueLu::TentativePFactory::MakeTentative: fine level NS part has a zero column");
356  }
357 
358  // Calculate QR decomposition (standard)
359  // NOTE: Q is stored in localQR and R is stored in coarseNS
360  if (aggSize >= Teuchos::as<LO>(NSDim)) {
361 
362  if (NSDim == 1) {
363  // Only one nullspace vector, calculate Q and R by hand
364  Magnitude norm = STS::magnitude(zero);
365  for (size_t k = 0; k < Teuchos::as<size_t>(aggSize); k++)
366  norm += STS::magnitude(localQR(k,0)*localQR(k,0));
367  norm = Teuchos::ScalarTraits<Magnitude>::squareroot(norm);
368 
369  // R = norm
370  coarseNS[0][offset] = norm;
371 
372  // Q = localQR(:,0)/norm
373  for (LO i = 0; i < aggSize; i++)
374  localQR(i,0) /= norm;
375 
376  } else {
377  Teuchos::SerialQRDenseSolver<LO,SC> qrSolver;
378  qrSolver.setMatrix(Teuchos::rcp(&localQR, false));
379  qrSolver.factor();
380 
381  // R = upper triangular part of localQR
382  for (size_t j = 0; j < NSDim; j++)
383  for (size_t k = 0; k <= j; k++)
384  coarseNS[j][offset+k] = localQR(k,j); //TODO is offset+k the correct local ID?!
385 
386  // Calculate Q, the tentative prolongator.
387  // The Lapack GEQRF call only works for myAggsize >= NSDim
388  qrSolver.formQ();
389  Teuchos::RCP<Teuchos::SerialDenseMatrix<LO,SC> > qFactor = qrSolver.getQ();
390  for (size_t j = 0; j < NSDim; j++)
391  for (size_t i = 0; i < Teuchos::as<size_t>(aggSize); i++)
392  localQR(i,j) = (*qFactor)(i,j);
393  }
394 
395  } else {
396  // Special handling for aggSize < NSDim (i.e. single node aggregates in structural mechanics)
397 
398  // The local QR decomposition is not possible in the "overconstrained"
399  // case (i.e. number of columns in localQR > number of rows), which
400  // corresponds to #DOFs in Aggregate < NSDim. For usual problems this
401  // is only possible for single node aggregates in structural mechanics.
402  // (Similar problems may arise in discontinuous Galerkin problems...)
403  // We bypass the QR decomposition and use an identity block in the
404  // tentative prolongator for the single node aggregate and transfer the
405  // corresponding fine level null space information 1-to-1 to the coarse
406  // level null space part.
407 
408  // NOTE: The resulting tentative prolongation operator has
409  // (aggSize*DofsPerNode-NSDim) zero columns leading to a singular
410  // coarse level operator A. To deal with that one has the following
411  // options:
412  // - Use the "RepairMainDiagonal" flag in the RAPFactory (default:
413  // false) to add some identity block to the diagonal of the zero rows
414  // in the coarse level operator A, such that standard level smoothers
415  // can be used again.
416  // - Use special (projection-based) level smoothers, which can deal
417  // with singular matrices (very application specific)
418  // - Adapt the code below to avoid zero columns. However, we do not
419  // support a variable number of DOFs per node in MueLu/Xpetra which
420  // makes the implementation really hard.
421 
422  // R = extended (by adding identity rows) localQR
423  for (size_t j = 0; j < NSDim; j++)
424  for (size_t k = 0; k < NSDim; k++)
425  if (k < as<size_t>(aggSize))
426  coarseNS[j][offset+k] = localQR(k,j);
427  else
428  coarseNS[j][offset+k] = (k == j ? one : zero);
429 
430  // Q = I (rectangular)
431  for (size_t i = 0; i < as<size_t>(aggSize); i++)
432  for (size_t j = 0; j < NSDim; j++)
433  localQR(i,j) = (j == i ? one : zero);
434  }
435 
436 
437  // Process each row in the local Q factor
438  // FIXME: What happens if maps are blocked?
439  for (LO j = 0; j < aggSize; j++) {
440  LO localRow = (goodMap ? aggToRowMapLO[aggStart[agg]+j] : rowMap->getLocalElement(aggToRowMapGO[aggStart[agg]+j]));
441 
442  size_t rowStart = ia[localRow];
443  for (size_t k = 0, lnnz = 0; k < NSDim; k++) {
444  // Skip zeros (there may be plenty of them, i.e., NSDim > 1 or boundary conditions)
445  if (localQR(j,k) != zero) {
446  ja [rowStart+lnnz] = offset + k;
447  val[rowStart+lnnz] = localQR(j,k);
448  lnnz++;
449  }
450  }
451  }
452  }
453 
454  } else {
455  GetOStream(Runtime1) << "TentativePFactory : bypassing local QR phase" << std::endl;
456  if (NSDim>1)
457  GetOStream(Warnings0) << "TentativePFactor : for nontrivial nullspace, this may degrade performance" << std::endl;
459  // "no-QR" option //
461  // Local Q factor is just the fine nullspace support over the current aggregate.
462  // Local R factor is the identity.
463  // TODO I have not implemented any special handling for aggregates that are too
464  // TODO small to locally support the nullspace, as is done in the standard QR
465  // TODO case above.
466  if (goodMap) {
467  for (GO agg = 0; agg < numAggs; agg++) {
468  const LO aggSize = aggStart[agg+1] - aggStart[agg];
469  Xpetra::global_size_t offset = agg*NSDim;
470 
471  // Process each row in the local Q factor
472  // FIXME: What happens if maps are blocked?
473  for (LO j = 0; j < aggSize; j++) {
474 
475  //TODO Here I do not check for a zero nullspace column on the aggregate.
476  // as is done in the standard QR case.
477 
478  const LO localRow = aggToRowMapLO[aggStart[agg]+j];
479 
480  const size_t rowStart = ia[localRow];
481 
482  for (size_t k = 0, lnnz = 0; k < NSDim; k++) {
483  // Skip zeros (there may be plenty of them, i.e., NSDim > 1 or boundary conditions)
484  const SC qr_jk = fineNS[k][aggToRowMapLO[aggStart[agg]+j]];
485  if (qr_jk != zero) {
486  ja [rowStart+lnnz] = offset + k;
487  val[rowStart+lnnz] = qr_jk;
488  lnnz++;
489  }
490  }
491  }
492  for (size_t j = 0; j < NSDim; j++)
493  coarseNS[j][offset+j] = one;
494  } //for (GO agg = 0; agg < numAggs; agg++)
495 
496  } else {
497  for (GO agg = 0; agg < numAggs; agg++) {
498  const LO aggSize = aggStart[agg+1] - aggStart[agg];
499  Xpetra::global_size_t offset = agg*NSDim;
500  for (LO j = 0; j < aggSize; j++) {
501 
502  const LO localRow = rowMap->getLocalElement(aggToRowMapGO[aggStart[agg]+j]);
503 
504  const size_t rowStart = ia[localRow];
505 
506  for (size_t k = 0, lnnz = 0; k < NSDim; k++) {
507  // Skip zeros (there may be plenty of them, i.e., NSDim > 1 or boundary conditions)
508  const SC qr_jk = fineNS[k][rowMap->getLocalElement(aggToRowMapGO[aggStart[agg]+j])];
509  if (qr_jk != zero) {
510  ja [rowStart+lnnz] = offset + k;
511  val[rowStart+lnnz] = qr_jk;
512  lnnz++;
513  }
514  }
515  }
516  for (size_t j = 0; j < NSDim; j++)
517  coarseNS[j][offset+j] = one;
518  } //for (GO agg = 0; agg < numAggs; agg++)
519 
520  } //if (goodmap) else ...
521 
522  } //if doQRStep ... else
523 
524  // Compress storage (remove all INVALID, which happen when we skip zeros)
525  // We do that in-place
526  size_t ia_tmp = 0, nnz = 0;
527  for (size_t i = 0; i < numRows; i++) {
528  for (size_t j = ia_tmp; j < ia[i+1]; j++)
529  if (ja[j] != INVALID) {
530  ja [nnz] = ja [j];
531  val[nnz] = val[j];
532  nnz++;
533  }
534  ia_tmp = ia[i+1];
535  ia[i+1] = nnz;
536  }
537  if (rowMap->lib() == Xpetra::UseTpetra) {
538  // - Cannot resize for Epetra, as it checks for same pointers
539  // - Need to resize for Tpetra, as it check ().size() == ia[numRows]
540  // NOTE: these invalidate ja and val views
541  jaPtent .resize(nnz);
542  valPtent.resize(nnz);
543  }
544 
545  GetOStream(Runtime1) << "TentativePFactory : aggregates do not cross process boundaries" << std::endl;
546 
547  PtentCrs->setAllValues(iaPtent, jaPtent, valPtent);
548 
549 
550  // Managing labels & constants for ESFC
551  RCP<ParameterList> FCparams;
552  if(pL.isSublist("matrixmatrix: kernel params"))
553  FCparams=rcp(new ParameterList(pL.sublist("matrixmatrix: kernel params")));
554  else
555  FCparams= rcp(new ParameterList);
556  // By default, we don't need global constants for TentativeP
557  FCparams->set("compute global constants",FCparams->get("compute global constants",false));
558  std::string levelIDs = toString(levelID);
559  FCparams->set("Timer Label",std::string("MueLu::TentativeP-")+levelIDs);
560  RCP<const Export> dummy_e;
561  RCP<const Import> dummy_i;
562 
563  PtentCrs->expertStaticFillComplete(coarseMap, A->getDomainMap(),dummy_i,dummy_e,FCparams);
564  }
565 
566  template <class Scalar,class LocalOrdinal, class GlobalOrdinal, class Node>
568  BuildPcoupled(RCP<Matrix> A, RCP<Aggregates> aggregates, RCP<AmalgamationInfo> amalgInfo, RCP<MultiVector> fineNullspace,
569  RCP<const Map> coarseMap, RCP<Matrix>& Ptentative, RCP<MultiVector>& coarseNullspace) const {
570  typedef Teuchos::ScalarTraits<SC> STS;
571  typedef typename STS::magnitudeType Magnitude;
572  const SC zero = STS::zero();
573  const SC one = STS::one();
574 
575  // number of aggregates
576  GO numAggs = aggregates->GetNumAggregates();
577 
578  // Create a lookup table to determine the rows (fine DOFs) that belong to a given aggregate.
579  // aggStart is a pointer into aggToRowMap
580  // aggStart[i]..aggStart[i+1] are indices into aggToRowMap
581  // aggToRowMap[aggStart[i]]..aggToRowMap[aggStart[i+1]-1] are the DOFs in aggregate i
582  ArrayRCP<LO> aggStart;
583  ArrayRCP< GO > aggToRowMap;
584  amalgInfo->UnamalgamateAggregates(*aggregates, aggStart, aggToRowMap);
585 
586  // find size of the largest aggregate
587  LO maxAggSize=0;
588  for (GO i=0; i<numAggs; ++i) {
589  LO sizeOfThisAgg = aggStart[i+1] - aggStart[i];
590  if (sizeOfThisAgg > maxAggSize) maxAggSize = sizeOfThisAgg;
591  }
592 
593  // dimension of fine level nullspace
594  const size_t NSDim = fineNullspace->getNumVectors();
595 
596  // index base for coarse Dof map (usually 0)
597  GO indexBase=A->getRowMap()->getIndexBase();
598 
599  const RCP<const Map> nonUniqueMap = amalgInfo->ComputeUnamalgamatedImportDofMap(*aggregates);
600  const RCP<const Map> uniqueMap = A->getDomainMap();
601  RCP<const Import> importer = ImportFactory::Build(uniqueMap, nonUniqueMap);
602  RCP<MultiVector> fineNullspaceWithOverlap = MultiVectorFactory::Build(nonUniqueMap,NSDim);
603  fineNullspaceWithOverlap->doImport(*fineNullspace,*importer,Xpetra::INSERT);
604 
605  // Pull out the nullspace vectors so that we can have random access.
606  ArrayRCP< ArrayRCP<const SC> > fineNS(NSDim);
607  for (size_t i=0; i<NSDim; ++i)
608  fineNS[i] = fineNullspaceWithOverlap->getData(i);
609 
610  //Allocate storage for the coarse nullspace.
611  coarseNullspace = MultiVectorFactory::Build(coarseMap, NSDim);
612 
613  ArrayRCP< ArrayRCP<SC> > coarseNS(NSDim);
614  for (size_t i=0; i<NSDim; ++i)
615  if (coarseMap->getNodeNumElements() > 0) coarseNS[i] = coarseNullspace->getDataNonConst(i);
616 
617  //This makes the rowmap of Ptent the same as that of A->
618  //This requires moving some parts of some local Q's to other processors
619  //because aggregates can span processors.
620  RCP<const Map > rowMapForPtent = A->getRowMap();
621  const Map& rowMapForPtentRef = *rowMapForPtent;
622 
623  // Set up storage for the rows of the local Qs that belong to other processors.
624  // FIXME This is inefficient and could be done within the main loop below with std::vector's.
625  RCP<const Map> colMap = A->getColMap();
626 
627  RCP<const Map > ghostQMap;
628  RCP<MultiVector> ghostQvalues;
629  Array<RCP<Xpetra::Vector<GO,LO,GO,Node> > > ghostQcolumns;
630  RCP<Xpetra::Vector<GO,LO,GO,Node> > ghostQrowNums;
631  ArrayRCP< ArrayRCP<SC> > ghostQvals;
632  ArrayRCP< ArrayRCP<GO> > ghostQcols;
633  ArrayRCP< GO > ghostQrows;
634 
635  Array<GO> ghostGIDs;
636  for (LO j=0; j<numAggs; ++j) {
637  for (LO k=aggStart[j]; k<aggStart[j+1]; ++k) {
638  if (rowMapForPtentRef.isNodeGlobalElement(aggToRowMap[k]) == false) {
639  ghostGIDs.push_back(aggToRowMap[k]);
640  }
641  }
642  }
643  ghostQMap = MapFactory::Build(A->getRowMap()->lib(),
644  Teuchos::OrdinalTraits<Xpetra::global_size_t>::invalid(),
645  ghostGIDs,
646  indexBase, A->getRowMap()->getComm()); //JG:Xpetra::global_size_t>?
647  //Vector to hold bits of Q that go to other processors.
648  ghostQvalues = MultiVectorFactory::Build(ghostQMap,NSDim);
649  //Note that Epetra does not support MultiVectors templated on Scalar != double.
650  //So to work around this, we allocate an array of Vectors. This shouldn't be too
651  //expensive, as the number of Vectors is NSDim.
652  ghostQcolumns.resize(NSDim);
653  for (size_t i=0; i<NSDim; ++i)
654  ghostQcolumns[i] = Xpetra::VectorFactory<GO,LO,GO,Node>::Build(ghostQMap);
655  ghostQrowNums = Xpetra::VectorFactory<GO,LO,GO,Node>::Build(ghostQMap);
656  if (ghostQvalues->getLocalLength() > 0) {
657  ghostQvals.resize(NSDim);
658  ghostQcols.resize(NSDim);
659  for (size_t i=0; i<NSDim; ++i) {
660  ghostQvals[i] = ghostQvalues->getDataNonConst(i);
661  ghostQcols[i] = ghostQcolumns[i]->getDataNonConst(0);
662  }
663  ghostQrows = ghostQrowNums->getDataNonConst(0);
664  }
665 
666  //importer to handle moving Q
667  importer = ImportFactory::Build(ghostQMap, A->getRowMap());
668 
669  // Dense QR solver
670  Teuchos::SerialQRDenseSolver<LO,SC> qrSolver;
671 
672  //Allocate temporary storage for the tentative prolongator.
673  Array<GO> globalColPtr(maxAggSize*NSDim,0);
674  Array<LO> localColPtr(maxAggSize*NSDim,0);
675  Array<SC> valPtr(maxAggSize*NSDim,0.);
676 
677  //Create column map for Ptent, estimate local #nonzeros in Ptent, and create Ptent itself.
678  const Map& coarseMapRef = *coarseMap;
679 
680  // For the 3-arrays constructor
681  ArrayRCP<size_t> ptent_rowptr;
682  ArrayRCP<LO> ptent_colind;
683  ArrayRCP<Scalar> ptent_values;
684 
685  // Because ArrayRCPs are slow...
686  ArrayView<size_t> rowptr_v;
687  ArrayView<LO> colind_v;
688  ArrayView<Scalar> values_v;
689 
690  // For temporary usage
691  Array<size_t> rowptr_temp;
692  Array<LO> colind_temp;
693  Array<Scalar> values_temp;
694 
695  RCP<CrsMatrix> PtentCrs;
696 
697  RCP<CrsMatrixWrap> PtentCrsWrap = rcp(new CrsMatrixWrap(rowMapForPtent, NSDim, Xpetra::StaticProfile));
698  PtentCrs = PtentCrsWrap->getCrsMatrix();
699  Ptentative = PtentCrsWrap;
700 
701  //*****************************************************************
702  //Loop over all aggregates and calculate local QR decompositions.
703  //*****************************************************************
704  GO qctr=0; //for indexing into Ptent data vectors
705  const Map& nonUniqueMapRef = *nonUniqueMap;
706 
707  size_t total_nnz_count=0;
708 
709  for (GO agg=0; agg<numAggs; ++agg)
710  {
711  LO myAggSize = aggStart[agg+1]-aggStart[agg];
712  // For each aggregate, extract the corresponding piece of the nullspace and put it in the flat array,
713  // "localQR" (in column major format) for the QR routine.
714  Teuchos::SerialDenseMatrix<LO,SC> localQR(myAggSize, NSDim);
715  for (size_t j=0; j<NSDim; ++j) {
716  bool bIsZeroNSColumn = true;
717  for (LO k=0; k<myAggSize; ++k)
718  {
719  // aggToRowMap[aggPtr[i]+k] is the kth DOF in the ith aggregate
720  // fineNS[j][n] is the nth entry in the jth NS vector
721  try{
722  SC nsVal = fineNS[j][ nonUniqueMapRef.getLocalElement(aggToRowMap[aggStart[agg]+k]) ]; // extract information from fine level NS
723  localQR(k,j) = nsVal;
724  if (nsVal != zero) bIsZeroNSColumn = false;
725  }
726  catch(...) {
727  GetOStream(Runtime1,-1) << "length of fine level nsp: " << fineNullspace->getGlobalLength() << std::endl;
728  GetOStream(Runtime1,-1) << "length of fine level nsp w overlap: " << fineNullspaceWithOverlap->getGlobalLength() << std::endl;
729  GetOStream(Runtime1,-1) << "(local?) aggId=" << agg << std::endl;
730  GetOStream(Runtime1,-1) << "aggSize=" << myAggSize << std::endl;
731  GetOStream(Runtime1,-1) << "agg DOF=" << k << std::endl;
732  GetOStream(Runtime1,-1) << "NS vector j=" << j << std::endl;
733  GetOStream(Runtime1,-1) << "j*myAggSize + k = " << j*myAggSize + k << std::endl;
734  GetOStream(Runtime1,-1) << "aggToRowMap["<<agg<<"][" << k << "] = " << aggToRowMap[aggStart[agg]+k] << std::endl;
735  GetOStream(Runtime1,-1) << "id aggToRowMap[agg][k]=" << aggToRowMap[aggStart[agg]+k] << " is global element in nonUniqueMap = " <<
736  nonUniqueMapRef.isNodeGlobalElement(aggToRowMap[aggStart[agg]+k]) << std::endl;
737  GetOStream(Runtime1,-1) << "colMap local id aggToRowMap[agg][k]=" << nonUniqueMapRef.getLocalElement(aggToRowMap[aggStart[agg]+k]) << std::endl;
738  GetOStream(Runtime1,-1) << "fineNS...=" << fineNS[j][ nonUniqueMapRef.getLocalElement(aggToRowMap[aggStart[agg]+k]) ] << std::endl;
739  GetOStream(Errors,-1) << "caught an error!" << std::endl;
740  }
741  } //for (LO k=0 ...
742  TEUCHOS_TEST_FOR_EXCEPTION(bIsZeroNSColumn == true, Exceptions::RuntimeError, "MueLu::TentativePFactory::MakeTentative: fine level NS part has a zero column. Error.");
743  } //for (LO j=0 ...
744 
745  Xpetra::global_size_t offset=agg*NSDim;
746 
747  if(myAggSize >= Teuchos::as<LocalOrdinal>(NSDim)) {
748  // calculate QR decomposition (standard)
749  // R is stored in localQR (size: myAggSize x NSDim)
750 
751  // Householder multiplier
752  SC tau = localQR(0,0);
753 
754  if (NSDim == 1) {
755  // Only one nullspace vector, so normalize by hand
756  Magnitude dtemp=0;
757  for (size_t k = 0; k < Teuchos::as<size_t>(myAggSize); ++k) {
758  Magnitude tmag = STS::magnitude(localQR(k,0));
759  dtemp += tmag*tmag;
760  }
761  dtemp = Teuchos::ScalarTraits<Magnitude>::squareroot(dtemp);
762  tau = localQR(0,0);
763  localQR(0,0) = dtemp;
764  } else {
765  qrSolver.setMatrix( Teuchos::rcp(&localQR, false) );
766  qrSolver.factor();
767  }
768 
769  // Extract R, the coarse nullspace. This is stored in upper triangular part of localQR.
770  // Note: coarseNS[i][.] is the ith coarse nullspace vector, which may be counter to your intuition.
771  // This stores the (offset+k)th entry only if it is local according to the coarseMap.
772  for (size_t j=0; j<NSDim; ++j) {
773  for (size_t k=0; k<=j; ++k) {
774  try {
775  if (coarseMapRef.isNodeLocalElement(offset+k)) {
776  coarseNS[j][offset+k] = localQR(k, j); //TODO is offset+k the correct local ID?!
777  }
778  }
779  catch(...) {
780  GetOStream(Errors,-1) << "caught error in coarseNS insert, j="<<j<<", offset+k = "<<offset+k<<std::endl;
781  }
782  }
783  }
784 
785  // Calculate Q, the tentative prolongator.
786  // The Lapack GEQRF call only works for myAggsize >= NSDim
787 
788  if (NSDim == 1) {
789  // Only one nullspace vector, so calculate Q by hand
790  Magnitude dtemp = Teuchos::ScalarTraits<SC>::magnitude(localQR(0,0));
791  localQR(0,0) = tau;
792  dtemp = 1 / dtemp;
793  for (LocalOrdinal i=0; i<myAggSize; ++i) {
794  localQR(i,0) *= dtemp ;
795  }
796  } else {
797  qrSolver.formQ();
798  Teuchos::RCP<Teuchos::SerialDenseMatrix<LO,SC> > qFactor = qrSolver.getQ();
799  for (size_t j=0; j<NSDim; j++) {
800  for (size_t i = 0; i < Teuchos::as<size_t>(myAggSize); i++) {
801  localQR(i,j) = (*qFactor)(i,j);
802  }
803  }
804  }
805 
806  // end default case (myAggSize >= NSDim)
807  } else { // special handling for myAggSize < NSDim (i.e. 1pt nodes)
808  // See comments for the uncoupled case
809 
810  // R = extended (by adding identity rows) localQR
811  for (size_t j = 0; j < NSDim; j++)
812  for (size_t k = 0; k < NSDim; k++) {
813  TEUCHOS_TEST_FOR_EXCEPTION(!coarseMapRef.isNodeLocalElement(offset+k), Exceptions::RuntimeError,
814  "Caught error in coarseNS insert, j=" << j << ", offset+k = " << offset+k);
815 
816  if (k < as<size_t>(myAggSize))
817  coarseNS[j][offset+k] = localQR(k,j);
818  else
819  coarseNS[j][offset+k] = (k == j ? one : zero);
820  }
821 
822  // Q = I (rectangular)
823  for (size_t i = 0; i < as<size_t>(myAggSize); i++)
824  for (size_t j = 0; j < NSDim; j++)
825  localQR(i,j) = (j == i ? one : zero);
826  } // end else (special handling for 1pt aggregates)
827 
828  //Process each row in the local Q factor. If the row is local to the current processor
829  //according to the rowmap, insert it into Ptentative. Otherwise, save it in ghostQ
830  //to be communicated later to the owning processor.
831  //FIXME -- what happens if maps are blocked?
832  for (GO j=0; j<myAggSize; ++j) {
833  //This loop checks whether row associated with current DOF is local, according to rowMapForPtent.
834  //If it is, the row is inserted. If not, the row number, columns, and values are saved in
835  //MultiVectors that will be sent to other processors.
836  GO globalRow = aggToRowMap[aggStart[agg]+j];
837 
838  //TODO is the use of Xpetra::global_size_t below correct?
839  if (rowMapForPtentRef.isNodeGlobalElement(globalRow) == false ) {
840  ghostQrows[qctr] = globalRow;
841  for (size_t k=0; k<NSDim; ++k) {
842  ghostQcols[k][qctr] = coarseMapRef.getGlobalElement(agg*NSDim+k);
843  ghostQvals[k][qctr] = localQR(j,k);
844  }
845  ++qctr;
846  } else {
847  size_t nnz=0;
848  for (size_t k=0; k<NSDim; ++k) {
849  try{
850  if (localQR(j,k) != Teuchos::ScalarTraits<SC>::zero()) {
851  localColPtr[nnz] = agg * NSDim + k;
852  globalColPtr[nnz] = coarseMapRef.getGlobalElement(localColPtr[nnz]);
853  valPtr[nnz] = localQR(j,k);
854  ++total_nnz_count;
855  ++nnz;
856  }
857  }
858  catch(...) {
859  GetOStream(Errors,-1) << "caught error in colPtr/valPtr insert, current index="<<nnz<<std::endl;
860  }
861  } //for (size_t k=0; k<NSDim; ++k)
862 
863  try{
864  Ptentative->insertGlobalValues(globalRow,globalColPtr.view(0,nnz),valPtr.view(0,nnz));
865  }
866  catch(...) {
867  GetOStream(Errors,-1) << "pid " << A->getRowMap()->getComm()->getRank()
868  << "caught error during Ptent row insertion, global row "
869  << globalRow << std::endl;
870  }
871  }
872  } //for (GO j=0; j<myAggSize; ++j)
873 
874  } // for (LO agg=0; agg<numAggs; ++agg)
875 
876 
877  // ***********************************************************
878  // ************* end of aggregate-wise QR ********************
879  // ***********************************************************
880  GetOStream(Runtime1) << "TentativePFactory : aggregates may cross process boundaries" << std::endl;
881  // Import ghost parts of Q factors and insert into Ptentative.
882  // First import just the global row numbers.
883  RCP<Xpetra::Vector<GO,LO,GO,Node> > targetQrowNums = Xpetra::VectorFactory<GO,LO,GO,Node>::Build(rowMapForPtent);
884  targetQrowNums->putScalar(-1);
885  targetQrowNums->doImport(*ghostQrowNums,*importer,Xpetra::INSERT);
886  ArrayRCP< GO > targetQrows = targetQrowNums->getDataNonConst(0);
887 
888  // Now create map based on just the row numbers imported.
889  Array<GO> gidsToImport;
890  gidsToImport.reserve(targetQrows.size());
891  for (typename ArrayRCP<GO>::iterator r=targetQrows.begin(); r!=targetQrows.end(); ++r) {
892  if (*r > -1) {
893  gidsToImport.push_back(*r);
894  }
895  }
896  RCP<const Map > reducedMap = MapFactory::Build( A->getRowMap()->lib(),
897  Teuchos::OrdinalTraits<Xpetra::global_size_t>::invalid(),
898  gidsToImport, indexBase, A->getRowMap()->getComm() );
899 
900  // Import using the row numbers that this processor will receive.
901  importer = ImportFactory::Build(ghostQMap, reducedMap);
902 
903  Array<RCP<Xpetra::Vector<GO,LO,GO,Node> > > targetQcolumns(NSDim);
904  for (size_t i=0; i<NSDim; ++i) {
905  targetQcolumns[i] = Xpetra::VectorFactory<GO,LO,GO,Node>::Build(reducedMap);
906  targetQcolumns[i]->doImport(*(ghostQcolumns[i]),*importer,Xpetra::INSERT);
907  }
908  RCP<MultiVector> targetQvalues = MultiVectorFactory::Build(reducedMap,NSDim);
909  targetQvalues->doImport(*ghostQvalues,*importer,Xpetra::INSERT);
910 
911  ArrayRCP< ArrayRCP<SC> > targetQvals;
912  ArrayRCP<ArrayRCP<GO> > targetQcols;
913  if (targetQvalues->getLocalLength() > 0) {
914  targetQvals.resize(NSDim);
915  targetQcols.resize(NSDim);
916  for (size_t i=0; i<NSDim; ++i) {
917  targetQvals[i] = targetQvalues->getDataNonConst(i);
918  targetQcols[i] = targetQcolumns[i]->getDataNonConst(0);
919  }
920  }
921 
922  valPtr = Array<SC>(NSDim,0.);
923  globalColPtr = Array<GO>(NSDim,0);
924  for (typename Array<GO>::iterator r=gidsToImport.begin(); r!=gidsToImport.end(); ++r) {
925  if (targetQvalues->getLocalLength() > 0) {
926  for (size_t j=0; j<NSDim; ++j) {
927  valPtr[j] = targetQvals[j][reducedMap->getLocalElement(*r)];
928  globalColPtr[j] = targetQcols[j][reducedMap->getLocalElement(*r)];
929  }
930  Ptentative->insertGlobalValues(*r, globalColPtr.view(0,NSDim), valPtr.view(0,NSDim));
931  } //if (targetQvalues->getLocalLength() > 0)
932  }
933 
934  Ptentative->fillComplete(coarseMap, A->getDomainMap());
935  }
936 
937  template <class Scalar,class LocalOrdinal, class GlobalOrdinal, class Node>
938  bool TentativePFactory<Scalar, LocalOrdinal, GlobalOrdinal, Node>::isGoodMap(const Map& rowMap, const Map& colMap) const {
939  ArrayView<const GO> rowElements = rowMap.getNodeElementList();
940  ArrayView<const GO> colElements = colMap.getNodeElementList();
941 
942  const size_t numElements = rowElements.size();
943 
944  bool goodMap = true;
945  for (size_t i = 0; i < numElements; i++)
946  if (rowElements[i] != colElements[i]) {
947  goodMap = false;
948  break;
949  }
950 
951  return goodMap;
952  }
953 
954 } //namespace MueLu
955 
956 // TODO ReUse: If only P or Nullspace is missing, TentativePFactory can be smart and skip part of the computation.
957 
958 #define MUELU_TENTATIVEPFACTORY_SHORT
959 #endif // MUELU_TENTATIVEPFACTORY_DEF_HPP
Important warning messages (one line)
std::string toString(const T &what)
Little helper function to convert non-string types to strings.
Timer to be used in factories. Similar to Monitor but with additional timers.
void BuildP(Level &fineLevel, Level &coarseLevel) const
Abstract Build method.
RCP< const ParameterList > GetValidParameterList() const
Return a const parameter list of valid parameters that setParameterList() will accept.
Namespace for MueLu classes and methods.
static const NoFactory * get()
Print even more statistics.
int GetLevelID() const
Return level number.
Definition: MueLu_Level.cpp:76
#define SET_VALID_ENTRY(name)
void BuildPcoupled(RCP< Matrix > A, RCP< Aggregates > aggregates, RCP< AmalgamationInfo > amalgInfo, RCP< MultiVector > fineNullspace, RCP< const Map > coarseMap, RCP< Matrix > &Ptentative, RCP< MultiVector > &coarseNullspace) const
Class that holds all level-specific information.
Definition: MueLu_Level.hpp:99
void BuildPuncoupled(RCP< Matrix > A, RCP< Aggregates > aggregates, RCP< AmalgamationInfo > amalgInfo, RCP< MultiVector > fineNullspace, RCP< const Map > coarseMap, RCP< Matrix > &Ptentative, RCP< MultiVector > &coarseNullspace, const int levelID) const
static std::string PrintMatrixInfo(const Matrix &A, const std::string &msgTag, RCP< const Teuchos::ParameterList > params=Teuchos::null)
void Build(Level &fineLevel, Level &coarseLevel) const
Build an object with this factory.
bool IsAvailable(const std::string &ename, const FactoryBase *factory=NoFactory::get()) const
Test whether a need&#39;s value has been saved.
Exception throws to report errors in the internal logical of the program.
Description of what is happening (more verbose)
void DeclareInput(Level &fineLevel, Level &coarseLevel) const
Input.
bool isGoodMap(const Map &rowMap, const Map &colMap) const