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>
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>
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_NullspaceFactory.hpp"
69 #include "MueLu_PerfUtils.hpp"
70 #include "MueLu_Monitor.hpp"
71 #include "MueLu_Utilities.hpp"
72 
73 namespace MueLu {
74 
75  template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
77  RCP<ParameterList> validParamList = rcp(new ParameterList());
78 
79  validParamList->set< RCP<const FactoryBase> >("A", Teuchos::null, "Generating factory of the matrix A");
80  validParamList->set< RCP<const FactoryBase> >("Aggregates", Teuchos::null, "Generating factory of the aggregates");
81  validParamList->set< RCP<const FactoryBase> >("Nullspace", Teuchos::null, "Generating factory of the nullspace");
82  validParamList->set< RCP<const FactoryBase> >("UnAmalgamationInfo", Teuchos::null, "Generating factory of UnAmalgamationInfo");
83  validParamList->set< RCP<const FactoryBase> >("CoarseMap", Teuchos::null, "Generating factory of the coarse map");
84  return validParamList;
85  }
86 
87  template <class Scalar,class LocalOrdinal, class GlobalOrdinal, class Node>
89  Input(fineLevel, "A");
90  Input(fineLevel, "Aggregates");
91  Input(fineLevel, "Nullspace");
92  Input(fineLevel, "UnAmalgamationInfo");
93  Input(fineLevel, "CoarseMap");
94  }
95 
96  template <class Scalar,class LocalOrdinal, class GlobalOrdinal, class Node>
98  return BuildP(fineLevel, coarseLevel);
99  }
100 
101  template <class Scalar,class LocalOrdinal, class GlobalOrdinal, class Node>
103  FactoryMonitor m(*this, "Build", coarseLevel);
104 
105  RCP<Matrix> A = Get< RCP<Matrix> > (fineLevel, "A");
106  RCP<Aggregates> aggregates = Get< RCP<Aggregates> > (fineLevel, "Aggregates");
107  RCP<AmalgamationInfo> amalgInfo = Get< RCP<AmalgamationInfo> >(fineLevel, "UnAmalgamationInfo");
108  RCP<MultiVector> fineNullspace = Get< RCP<MultiVector> > (fineLevel, "Nullspace");
109  RCP<const Map> coarseMap = Get< RCP<const Map> > (fineLevel, "CoarseMap");
110 
111  RCP<Matrix> Ptentative;
112  RCP<MultiVector> coarseNullspace;
113  if (!aggregates->AggregatesCrossProcessors())
114  BuildPuncoupled(A, aggregates, amalgInfo, fineNullspace, coarseMap, Ptentative, coarseNullspace);
115  else
116  BuildPcoupled (A, aggregates, amalgInfo, fineNullspace, coarseMap, Ptentative, coarseNullspace);
117 
118  // If available, use striding information of fine level matrix A for range
119  // map and coarseMap as domain map; otherwise use plain range map of
120  // Ptent = plain range map of A for range map and coarseMap as domain map.
121  // NOTE:
122  // The latter is not really safe, since there is no striding information
123  // for the range map. This is not really a problem, since striding
124  // information is always available on the intermedium levels and the
125  // coarsest levels.
126  if (A->IsView("stridedMaps") == true)
127  Ptentative->CreateView("stridedMaps", A->getRowMap("stridedMaps"), coarseMap);
128  else
129  Ptentative->CreateView("stridedMaps", Ptentative->getRangeMap(), coarseMap);
130 
131  Set(coarseLevel, "Nullspace", coarseNullspace);
132  Set(coarseLevel, "P", Ptentative);
133 
134  if (IsPrint(Statistics1)) {
135  RCP<ParameterList> params = rcp(new ParameterList());
136  params->set("printLoadBalancingInfo", true);
137  GetOStream(Statistics1) << PerfUtils::PrintMatrixInfo(*Ptentative, "Ptent", params);
138  }
139  }
140 
141  template <class Scalar,class LocalOrdinal, class GlobalOrdinal, class Node>
143  BuildPuncoupled(RCP<Matrix> A, RCP<Aggregates> aggregates, RCP<AmalgamationInfo> amalgInfo, RCP<MultiVector> fineNullspace,
144  RCP<const Map> coarseMap, RCP<Matrix>& Ptentative, RCP<MultiVector>& coarseNullspace) const {
145  RCP<const Map> rowMap = A->getRowMap();
146  RCP<const Map> colMap = A->getColMap();
147 
148  const size_t numRows = rowMap->getNodeNumElements();
149 
150  typedef Teuchos::ScalarTraits<SC> STS;
151  typedef typename STS::magnitudeType Magnitude;
152  const SC zero = STS::zero();
153  const SC one = STS::one();
154  const LO INVALID = Teuchos::OrdinalTraits<LO>::invalid();
155 
156  const GO numAggs = aggregates->GetNumAggregates();
157  const size_t NSDim = fineNullspace->getNumVectors();
158 
159  // Aggregates map is based on the amalgamated column map
160  // We can skip global-to-local conversion if LIDs in row map are
161  // same as LIDs in column map
162  bool goodMap = isGoodMap(*rowMap, *colMap);
163 
164  // Create a lookup table to determine the rows (fine DOFs) that belong to a given aggregate.
165  // aggStart is a pointer into aggToRowMapLO
166  // aggStart[i]..aggStart[i+1] are indices into aggToRowMapLO
167  // aggToRowMapLO[aggStart[i]]..aggToRowMapLO[aggStart[i+1]-1] are the DOFs in aggregate i
168  ArrayRCP<LO> aggStart;
169  ArrayRCP<LO> aggToRowMapLO;
170  ArrayRCP<GO> aggToRowMapGO;
171  if (goodMap) {
172  amalgInfo->UnamalgamateAggregatesLO(*aggregates, aggStart, aggToRowMapLO);
173  GetOStream(Runtime1) << "Column map is consistent with the row map, good." << std::endl;
174 
175  } else {
176  amalgInfo->UnamalgamateAggregates(*aggregates, aggStart, aggToRowMapGO);
177  GetOStream(Warnings0) << "Column map is not consistent with the row map\n"
178  << "using GO->LO conversion with performance penalty" << std::endl;
179  }
180 
181  coarseNullspace = MultiVectorFactory::Build(coarseMap, NSDim);
182 
183  // Pull out the nullspace vectors so that we can have random access.
184  ArrayRCP<ArrayRCP<const SC> > fineNS (NSDim);
185  ArrayRCP<ArrayRCP<SC> > coarseNS(NSDim);
186  for (size_t i = 0; i < NSDim; i++) {
187  fineNS[i] = fineNullspace->getData(i);
188  if (coarseMap->getNodeNumElements() > 0)
189  coarseNS[i] = coarseNullspace->getDataNonConst(i);
190  }
191 
192  size_t nnzEstimate = numRows * NSDim;
193 
194  // Time to construct the matrix and fill in the values
195  Ptentative = rcp(new CrsMatrixWrap(rowMap, coarseMap, 0, Xpetra::StaticProfile));
196  RCP<CrsMatrix> PtentCrs = rcp_dynamic_cast<CrsMatrixWrap>(Ptentative)->getCrsMatrix();
197 
198  ArrayRCP<size_t> iaPtent;
199  ArrayRCP<LO> jaPtent;
200  ArrayRCP<SC> valPtent;
201 
202  PtentCrs->allocateAllValues(nnzEstimate, iaPtent, jaPtent, valPtent);
203 
204  ArrayView<size_t> ia = iaPtent();
205  ArrayView<LO> ja = jaPtent();
206  ArrayView<SC> val = valPtent();
207 
208  ia[0] = 0;
209  for (size_t i = 1; i <= numRows; i++)
210  ia[i] = ia[i-1] + NSDim;
211 
212  for (size_t j = 0; j < nnzEstimate; j++) {
213  ja [j] = INVALID;
214  val[j] = zero;
215  }
216 
217  for (GO agg = 0; agg < numAggs; agg++) {
218  LO aggSize = aggStart[agg+1] - aggStart[agg];
219 
220  Xpetra::global_size_t offset = agg*NSDim;
221 
222  // Extract the piece of the nullspace corresponding to the aggregate, and
223  // put it in the flat array, "localQR" (in column major format) for the
224  // QR routine.
225  Teuchos::SerialDenseMatrix<LO,SC> localQR(aggSize, NSDim);
226  if (goodMap) {
227  for (size_t j = 0; j < NSDim; j++)
228  for (LO k = 0; k < aggSize; k++)
229  localQR(k,j) = fineNS[j][aggToRowMapLO[aggStart[agg]+k]];
230  } else {
231  for (size_t j = 0; j < NSDim; j++)
232  for (LO k = 0; k < aggSize; k++)
233  localQR(k,j) = fineNS[j][rowMap->getLocalElement(aggToRowMapGO[aggStart[agg]+k])];
234  }
235 
236  // Test for zero columns
237  for (size_t j = 0; j < NSDim; j++) {
238  bool bIsZeroNSColumn = true;
239 
240  for (LO k = 0; k < aggSize; k++)
241  if (localQR(k,j) != zero)
242  bIsZeroNSColumn = false;
243 
244  TEUCHOS_TEST_FOR_EXCEPTION(bIsZeroNSColumn == true, Exceptions::RuntimeError,
245  "MueLu::TentativePFactory::MakeTentative: fine level NS part has a zero column");
246  }
247 
248  // Calculate QR decomposition (standard)
249  // NOTE: Q is stored in localQR and R is stored in coarseNS
250  if (aggSize >= Teuchos::as<LO>(NSDim)) {
251 
252  if (NSDim == 1) {
253  // Only one nullspace vector, calculate Q and R by hand
254  Magnitude norm = STS::magnitude(zero);
255  for (size_t k = 0; k < Teuchos::as<size_t>(aggSize); k++)
256  norm += STS::magnitude(localQR(k,0)*localQR(k,0));
257  norm = Teuchos::ScalarTraits<Magnitude>::squareroot(norm);
258 
259  // R = norm
260  coarseNS[0][offset] = norm;
261 
262  // Q = localQR(:,0)/norm
263  for (LO i = 0; i < aggSize; i++)
264  localQR(i,0) /= norm;
265 
266  } else {
267  Teuchos::SerialQRDenseSolver<LO,SC> qrSolver;
268  qrSolver.setMatrix(Teuchos::rcp(&localQR, false));
269  qrSolver.factor();
270 
271  // R = upper triangular part of localQR
272  for (size_t j = 0; j < NSDim; j++)
273  for (size_t k = 0; k <= j; k++)
274  coarseNS[j][offset+k] = localQR(k,j); //TODO is offset+k the correct local ID?!
275 
276  // Calculate Q, the tentative prolongator.
277  // The Lapack GEQRF call only works for myAggsize >= NSDim
278  qrSolver.formQ();
279  Teuchos::RCP<Teuchos::SerialDenseMatrix<LO,SC> > qFactor = qrSolver.getQ();
280  for (size_t j = 0; j < NSDim; j++)
281  for (size_t i = 0; i < Teuchos::as<size_t>(aggSize); i++)
282  localQR(i,j) = (*qFactor)(i,j);
283  }
284 
285  } else {
286  // Special handling for aggSize < NSDim (i.e. single node aggregates in structural mechanics)
287 
288  // The local QR decomposition is not possible in the "overconstrained"
289  // case (i.e. number of columns in localQR > number of rows), which
290  // corresponds to #DOFs in Aggregate < NSDim. For usual problems this
291  // is only possible for single node aggregates in structural mechanics.
292  // (Similar problems may arise in discontinuous Galerkin problems...)
293  // We bypass the QR decomposition and use an identity block in the
294  // tentative prolongator for the single node aggregate and transfer the
295  // corresponding fine level null space information 1-to-1 to the coarse
296  // level null space part.
297 
298  // NOTE: The resulting tentative prolongation operator has
299  // (aggSize*DofsPerNode-NSDim) zero columns leading to a singular
300  // coarse level operator A. To deal with that one has the following
301  // options:
302  // - Use the "RepairMainDiagonal" flag in the RAPFactory (default:
303  // false) to add some identity block to the diagonal of the zero rows
304  // in the coarse level operator A, such that standard level smoothers
305  // can be used again.
306  // - Use special (projection-based) level smoothers, which can deal
307  // with singular matrices (very application specific)
308  // - Adapt the code below to avoid zero columns. However, we do not
309  // support a variable number of DOFs per node in MueLu/Xpetra which
310  // makes the implementation really hard.
311 
312  // R = extended (by adding identity rows) localQR
313  for (size_t j = 0; j < NSDim; j++)
314  for (size_t k = 0; k < NSDim; k++)
315  if (k < as<size_t>(aggSize))
316  coarseNS[j][offset+k] = localQR(k,j);
317  else
318  coarseNS[j][offset+k] = (k == j ? one : zero);
319 
320  // Q = I (rectangular)
321  for (size_t i = 0; i < as<size_t>(aggSize); i++)
322  for (size_t j = 0; j < NSDim; j++)
323  localQR(i,j) = (j == i ? one : zero);
324  }
325 
326  // Process each row in the local Q factor
327  // FIXME: What happens if maps are blocked?
328  for (LO j = 0; j < aggSize; j++) {
329  LO localRow = (goodMap ? aggToRowMapLO[aggStart[agg]+j] : rowMap->getLocalElement(aggToRowMapGO[aggStart[agg]+j]));
330 
331  size_t rowStart = ia[localRow];
332  for (size_t k = 0, lnnz = 0; k < NSDim; k++) {
333  // Skip zeros (there may be plenty of them, i.e., NSDim > 1 or boundary conditions)
334  if (localQR(j,k) != zero) {
335  ja [rowStart+lnnz] = offset + k;
336  val[rowStart+lnnz] = localQR(j,k);
337  lnnz++;
338  }
339  }
340  }
341  }
342 
343  // Compress storage (remove all INVALID, which happen when we skip zeros)
344  // We do that in-place
345  size_t ia_tmp = 0, nnz = 0;
346  for (size_t i = 0; i < numRows; i++) {
347  for (size_t j = ia_tmp; j < ia[i+1]; j++)
348  if (ja[j] != INVALID) {
349  ja [nnz] = ja [j];
350  val[nnz] = val[j];
351  nnz++;
352  }
353  ia_tmp = ia[i+1];
354  ia[i+1] = nnz;
355  }
356  if (rowMap->lib() == Xpetra::UseTpetra) {
357  // - Cannot resize for Epetra, as it checks for same pointers
358  // - Need to resize for Tpetra, as it check ().size() == ia[numRows]
359  // NOTE: these invalidate ja and val views
360  jaPtent .resize(nnz);
361  valPtent.resize(nnz);
362  }
363 
364  GetOStream(Runtime1) << "TentativePFactory : aggregates do not cross process boundaries" << std::endl;
365 
366  PtentCrs->setAllValues(iaPtent, jaPtent, valPtent);
367  PtentCrs->expertStaticFillComplete(coarseMap, A->getDomainMap());
368  }
369 
370  template <class Scalar,class LocalOrdinal, class GlobalOrdinal, class Node>
372  BuildPcoupled(RCP<Matrix> A, RCP<Aggregates> aggregates, RCP<AmalgamationInfo> amalgInfo, RCP<MultiVector> fineNullspace,
373  RCP<const Map> coarseMap, RCP<Matrix>& Ptentative, RCP<MultiVector>& coarseNullspace) const {
374  typedef Teuchos::ScalarTraits<SC> STS;
375  typedef typename STS::magnitudeType Magnitude;
376  const SC zero = STS::zero();
377  const SC one = STS::one();
378 
379  // number of aggregates
380  GO numAggs = aggregates->GetNumAggregates();
381 
382  // Create a lookup table to determine the rows (fine DOFs) that belong to a given aggregate.
383  // aggStart is a pointer into aggToRowMap
384  // aggStart[i]..aggStart[i+1] are indices into aggToRowMap
385  // aggToRowMap[aggStart[i]]..aggToRowMap[aggStart[i+1]-1] are the DOFs in aggregate i
386  ArrayRCP<LO> aggStart;
387  ArrayRCP< GO > aggToRowMap;
388  amalgInfo->UnamalgamateAggregates(*aggregates, aggStart, aggToRowMap);
389 
390  // find size of the largest aggregate
391  LO maxAggSize=0;
392  for (GO i=0; i<numAggs; ++i) {
393  LO sizeOfThisAgg = aggStart[i+1] - aggStart[i];
394  if (sizeOfThisAgg > maxAggSize) maxAggSize = sizeOfThisAgg;
395  }
396 
397  // dimension of fine level nullspace
398  const size_t NSDim = fineNullspace->getNumVectors();
399 
400  // index base for coarse Dof map (usually 0)
401  GO indexBase=A->getRowMap()->getIndexBase();
402 
403  const RCP<const Map> nonUniqueMap = amalgInfo->ComputeUnamalgamatedImportDofMap(*aggregates);
404  const RCP<const Map> uniqueMap = A->getDomainMap();
405  RCP<const Import> importer = ImportFactory::Build(uniqueMap, nonUniqueMap);
406  RCP<MultiVector> fineNullspaceWithOverlap = MultiVectorFactory::Build(nonUniqueMap,NSDim);
407  fineNullspaceWithOverlap->doImport(*fineNullspace,*importer,Xpetra::INSERT);
408 
409  // Pull out the nullspace vectors so that we can have random access.
410  ArrayRCP< ArrayRCP<const SC> > fineNS(NSDim);
411  for (size_t i=0; i<NSDim; ++i)
412  fineNS[i] = fineNullspaceWithOverlap->getData(i);
413 
414  //Allocate storage for the coarse nullspace.
415  coarseNullspace = MultiVectorFactory::Build(coarseMap, NSDim);
416 
417  ArrayRCP< ArrayRCP<SC> > coarseNS(NSDim);
418  for (size_t i=0; i<NSDim; ++i)
419  if (coarseMap->getNodeNumElements() > 0) coarseNS[i] = coarseNullspace->getDataNonConst(i);
420 
421  //This makes the rowmap of Ptent the same as that of A->
422  //This requires moving some parts of some local Q's to other processors
423  //because aggregates can span processors.
424  RCP<const Map > rowMapForPtent = A->getRowMap();
425  const Map& rowMapForPtentRef = *rowMapForPtent;
426 
427  // Set up storage for the rows of the local Qs that belong to other processors.
428  // FIXME This is inefficient and could be done within the main loop below with std::vector's.
429  RCP<const Map> colMap = A->getColMap();
430 
431  RCP<const Map > ghostQMap;
432  RCP<MultiVector> ghostQvalues;
433  Array<RCP<Xpetra::Vector<GO,LO,GO,Node> > > ghostQcolumns;
434  RCP<Xpetra::Vector<GO,LO,GO,Node> > ghostQrowNums;
435  ArrayRCP< ArrayRCP<SC> > ghostQvals;
436  ArrayRCP< ArrayRCP<GO> > ghostQcols;
437  ArrayRCP< GO > ghostQrows;
438 
439  Array<GO> ghostGIDs;
440  for (LO j=0; j<numAggs; ++j) {
441  for (LO k=aggStart[j]; k<aggStart[j+1]; ++k) {
442  if (rowMapForPtentRef.isNodeGlobalElement(aggToRowMap[k]) == false) {
443  ghostGIDs.push_back(aggToRowMap[k]);
444  }
445  }
446  }
447  ghostQMap = MapFactory::Build(A->getRowMap()->lib(),
448  Teuchos::OrdinalTraits<Xpetra::global_size_t>::invalid(),
449  ghostGIDs,
450  indexBase, A->getRowMap()->getComm()); //JG:Xpetra::global_size_t>?
451  //Vector to hold bits of Q that go to other processors.
452  ghostQvalues = MultiVectorFactory::Build(ghostQMap,NSDim);
453  //Note that Epetra does not support MultiVectors templated on Scalar != double.
454  //So to work around this, we allocate an array of Vectors. This shouldn't be too
455  //expensive, as the number of Vectors is NSDim.
456  ghostQcolumns.resize(NSDim);
457  for (size_t i=0; i<NSDim; ++i)
458  ghostQcolumns[i] = Xpetra::VectorFactory<GO,LO,GO,Node>::Build(ghostQMap);
459  ghostQrowNums = Xpetra::VectorFactory<GO,LO,GO,Node>::Build(ghostQMap);
460  if (ghostQvalues->getLocalLength() > 0) {
461  ghostQvals.resize(NSDim);
462  ghostQcols.resize(NSDim);
463  for (size_t i=0; i<NSDim; ++i) {
464  ghostQvals[i] = ghostQvalues->getDataNonConst(i);
465  ghostQcols[i] = ghostQcolumns[i]->getDataNonConst(0);
466  }
467  ghostQrows = ghostQrowNums->getDataNonConst(0);
468  }
469 
470  //importer to handle moving Q
471  importer = ImportFactory::Build(ghostQMap, A->getRowMap());
472 
473  // Dense QR solver
474  Teuchos::SerialQRDenseSolver<LO,SC> qrSolver;
475 
476  //Allocate temporary storage for the tentative prolongator.
477  Array<GO> globalColPtr(maxAggSize*NSDim,0);
478  Array<LO> localColPtr(maxAggSize*NSDim,0);
479  Array<SC> valPtr(maxAggSize*NSDim,0.);
480 
481  //Create column map for Ptent, estimate local #nonzeros in Ptent, and create Ptent itself.
482  const Map& coarseMapRef = *coarseMap;
483 
484  // For the 3-arrays constructor
485  ArrayRCP<size_t> ptent_rowptr;
486  ArrayRCP<LO> ptent_colind;
487  ArrayRCP<Scalar> ptent_values;
488 
489  // Because ArrayRCPs are slow...
490  ArrayView<size_t> rowptr_v;
491  ArrayView<LO> colind_v;
492  ArrayView<Scalar> values_v;
493 
494  // For temporary usage
495  Array<size_t> rowptr_temp;
496  Array<LO> colind_temp;
497  Array<Scalar> values_temp;
498 
499  RCP<CrsMatrix> PtentCrs;
500 
501  RCP<CrsMatrixWrap> PtentCrsWrap = rcp(new CrsMatrixWrap(rowMapForPtent, NSDim, Xpetra::StaticProfile));
502  PtentCrs = PtentCrsWrap->getCrsMatrix();
503  Ptentative = PtentCrsWrap;
504 
505  //*****************************************************************
506  //Loop over all aggregates and calculate local QR decompositions.
507  //*****************************************************************
508  GO qctr=0; //for indexing into Ptent data vectors
509  const Map& nonUniqueMapRef = *nonUniqueMap;
510 
511  size_t total_nnz_count=0;
512 
513  for (GO agg=0; agg<numAggs; ++agg)
514  {
515  LO myAggSize = aggStart[agg+1]-aggStart[agg];
516  // For each aggregate, extract the corresponding piece of the nullspace and put it in the flat array,
517  // "localQR" (in column major format) for the QR routine.
518  Teuchos::SerialDenseMatrix<LO,SC> localQR(myAggSize, NSDim);
519  for (size_t j=0; j<NSDim; ++j) {
520  bool bIsZeroNSColumn = true;
521  for (LO k=0; k<myAggSize; ++k)
522  {
523  // aggToRowMap[aggPtr[i]+k] is the kth DOF in the ith aggregate
524  // fineNS[j][n] is the nth entry in the jth NS vector
525  try{
526  SC nsVal = fineNS[j][ nonUniqueMapRef.getLocalElement(aggToRowMap[aggStart[agg]+k]) ]; // extract information from fine level NS
527  localQR(k,j) = nsVal;
528  if (nsVal != zero) bIsZeroNSColumn = false;
529  }
530  catch(...) {
531  GetOStream(Runtime1,-1) << "length of fine level nsp: " << fineNullspace->getGlobalLength() << std::endl;
532  GetOStream(Runtime1,-1) << "length of fine level nsp w overlap: " << fineNullspaceWithOverlap->getGlobalLength() << std::endl;
533  GetOStream(Runtime1,-1) << "(local?) aggId=" << agg << std::endl;
534  GetOStream(Runtime1,-1) << "aggSize=" << myAggSize << std::endl;
535  GetOStream(Runtime1,-1) << "agg DOF=" << k << std::endl;
536  GetOStream(Runtime1,-1) << "NS vector j=" << j << std::endl;
537  GetOStream(Runtime1,-1) << "j*myAggSize + k = " << j*myAggSize + k << std::endl;
538  GetOStream(Runtime1,-1) << "aggToRowMap["<<agg<<"][" << k << "] = " << aggToRowMap[aggStart[agg]+k] << std::endl;
539  GetOStream(Runtime1,-1) << "id aggToRowMap[agg][k]=" << aggToRowMap[aggStart[agg]+k] << " is global element in nonUniqueMap = " <<
540  nonUniqueMapRef.isNodeGlobalElement(aggToRowMap[aggStart[agg]+k]) << std::endl;
541  GetOStream(Runtime1,-1) << "colMap local id aggToRowMap[agg][k]=" << nonUniqueMapRef.getLocalElement(aggToRowMap[aggStart[agg]+k]) << std::endl;
542  GetOStream(Runtime1,-1) << "fineNS...=" << fineNS[j][ nonUniqueMapRef.getLocalElement(aggToRowMap[aggStart[agg]+k]) ] << std::endl;
543  GetOStream(Errors,-1) << "caught an error!" << std::endl;
544  }
545  } //for (LO k=0 ...
546  TEUCHOS_TEST_FOR_EXCEPTION(bIsZeroNSColumn == true, Exceptions::RuntimeError, "MueLu::TentativePFactory::MakeTentative: fine level NS part has a zero column. Error.");
547  } //for (LO j=0 ...
548 
549  Xpetra::global_size_t offset=agg*NSDim;
550 
551  if(myAggSize >= Teuchos::as<LocalOrdinal>(NSDim)) {
552  // calculate QR decomposition (standard)
553  // R is stored in localQR (size: myAggSize x NSDim)
554 
555  // Householder multiplier
556  SC tau = localQR(0,0);
557 
558  if (NSDim == 1) {
559  // Only one nullspace vector, so normalize by hand
560  Magnitude dtemp=0;
561  for (size_t k = 0; k < Teuchos::as<size_t>(myAggSize); ++k) {
562  Magnitude tmag = STS::magnitude(localQR(k,0));
563  dtemp += tmag*tmag;
564  }
565  dtemp = Teuchos::ScalarTraits<Magnitude>::squareroot(dtemp);
566  tau = localQR(0,0);
567  localQR(0,0) = dtemp;
568  } else {
569  qrSolver.setMatrix( Teuchos::rcp(&localQR, false) );
570  qrSolver.factor();
571  }
572 
573  // Extract R, the coarse nullspace. This is stored in upper triangular part of localQR.
574  // Note: coarseNS[i][.] is the ith coarse nullspace vector, which may be counter to your intuition.
575  // This stores the (offset+k)th entry only if it is local according to the coarseMap.
576  for (size_t j=0; j<NSDim; ++j) {
577  for (size_t k=0; k<=j; ++k) {
578  try {
579  if (coarseMapRef.isNodeLocalElement(offset+k)) {
580  coarseNS[j][offset+k] = localQR(k, j); //TODO is offset+k the correct local ID?!
581  }
582  }
583  catch(...) {
584  GetOStream(Errors,-1) << "caught error in coarseNS insert, j="<<j<<", offset+k = "<<offset+k<<std::endl;
585  }
586  }
587  }
588 
589  // Calculate Q, the tentative prolongator.
590  // The Lapack GEQRF call only works for myAggsize >= NSDim
591 
592  if (NSDim == 1) {
593  // Only one nullspace vector, so calculate Q by hand
594  Magnitude dtemp = Teuchos::ScalarTraits<SC>::magnitude(localQR(0,0));
595  localQR(0,0) = tau;
596  dtemp = 1 / dtemp;
597  for (LocalOrdinal i=0; i<myAggSize; ++i) {
598  localQR(i,0) *= dtemp ;
599  }
600  } else {
601  qrSolver.formQ();
602  Teuchos::RCP<Teuchos::SerialDenseMatrix<LO,SC> > qFactor = qrSolver.getQ();
603  for (size_t j=0; j<NSDim; j++) {
604  for (size_t i = 0; i < Teuchos::as<size_t>(myAggSize); i++) {
605  localQR(i,j) = (*qFactor)(i,j);
606  }
607  }
608  }
609 
610  // end default case (myAggSize >= NSDim)
611  } else { // special handling for myAggSize < NSDim (i.e. 1pt nodes)
612  // See comments for the uncoupled case
613 
614  // R = extended (by adding identity rows) localQR
615  for (size_t j = 0; j < NSDim; j++)
616  for (size_t k = 0; k < NSDim; k++) {
617  TEUCHOS_TEST_FOR_EXCEPTION(!coarseMapRef.isNodeLocalElement(offset+k), Exceptions::RuntimeError,
618  "Caught error in coarseNS insert, j=" << j << ", offset+k = " << offset+k);
619 
620  if (k < as<size_t>(myAggSize))
621  coarseNS[j][offset+k] = localQR(k,j);
622  else
623  coarseNS[j][offset+k] = (k == j ? one : zero);
624  }
625 
626  // Q = I (rectangular)
627  for (size_t i = 0; i < as<size_t>(myAggSize); i++)
628  for (size_t j = 0; j < NSDim; j++)
629  localQR(i,j) = (j == i ? one : zero);
630  } // end else (special handling for 1pt aggregates)
631 
632  //Process each row in the local Q factor. If the row is local to the current processor
633  //according to the rowmap, insert it into Ptentative. Otherwise, save it in ghostQ
634  //to be communicated later to the owning processor.
635  //FIXME -- what happens if maps are blocked?
636  for (GO j=0; j<myAggSize; ++j) {
637  //This loop checks whether row associated with current DOF is local, according to rowMapForPtent.
638  //If it is, the row is inserted. If not, the row number, columns, and values are saved in
639  //MultiVectors that will be sent to other processors.
640  GO globalRow = aggToRowMap[aggStart[agg]+j];
641 
642  //TODO is the use of Xpetra::global_size_t below correct?
643  if (rowMapForPtentRef.isNodeGlobalElement(globalRow) == false ) {
644  ghostQrows[qctr] = globalRow;
645  for (size_t k=0; k<NSDim; ++k) {
646  ghostQcols[k][qctr] = coarseMapRef.getGlobalElement(agg*NSDim+k);
647  ghostQvals[k][qctr] = localQR(j,k);
648  }
649  ++qctr;
650  } else {
651  size_t nnz=0;
652  for (size_t k=0; k<NSDim; ++k) {
653  try{
654  if (localQR(j,k) != Teuchos::ScalarTraits<SC>::zero()) {
655  localColPtr[nnz] = agg * NSDim + k;
656  globalColPtr[nnz] = coarseMapRef.getGlobalElement(localColPtr[nnz]);
657  valPtr[nnz] = localQR(j,k);
658  ++total_nnz_count;
659  ++nnz;
660  }
661  }
662  catch(...) {
663  GetOStream(Errors,-1) << "caught error in colPtr/valPtr insert, current index="<<nnz<<std::endl;
664  }
665  } //for (size_t k=0; k<NSDim; ++k)
666 
667  try{
668  Ptentative->insertGlobalValues(globalRow,globalColPtr.view(0,nnz),valPtr.view(0,nnz));
669  }
670  catch(...) {
671  GetOStream(Errors,-1) << "pid " << A->getRowMap()->getComm()->getRank()
672  << "caught error during Ptent row insertion, global row "
673  << globalRow << std::endl;
674  }
675  }
676  } //for (GO j=0; j<myAggSize; ++j)
677 
678  } // for (LO agg=0; agg<numAggs; ++agg)
679 
680 
681  // ***********************************************************
682  // ************* end of aggregate-wise QR ********************
683  // ***********************************************************
684  GetOStream(Runtime1) << "TentativePFactory : aggregates may cross process boundaries" << std::endl;
685  // Import ghost parts of Q factors and insert into Ptentative.
686  // First import just the global row numbers.
687  RCP<Xpetra::Vector<GO,LO,GO,Node> > targetQrowNums = Xpetra::VectorFactory<GO,LO,GO,Node>::Build(rowMapForPtent);
688  targetQrowNums->putScalar(-1);
689  targetQrowNums->doImport(*ghostQrowNums,*importer,Xpetra::INSERT);
690  ArrayRCP< GO > targetQrows = targetQrowNums->getDataNonConst(0);
691 
692  // Now create map based on just the row numbers imported.
693  Array<GO> gidsToImport;
694  gidsToImport.reserve(targetQrows.size());
695  for (typename ArrayRCP<GO>::iterator r=targetQrows.begin(); r!=targetQrows.end(); ++r) {
696  if (*r > -1) {
697  gidsToImport.push_back(*r);
698  }
699  }
700  RCP<const Map > reducedMap = MapFactory::Build( A->getRowMap()->lib(),
701  Teuchos::OrdinalTraits<Xpetra::global_size_t>::invalid(),
702  gidsToImport, indexBase, A->getRowMap()->getComm() );
703 
704  // Import using the row numbers that this processor will receive.
705  importer = ImportFactory::Build(ghostQMap, reducedMap);
706 
707  Array<RCP<Xpetra::Vector<GO,LO,GO,Node> > > targetQcolumns(NSDim);
708  for (size_t i=0; i<NSDim; ++i) {
709  targetQcolumns[i] = Xpetra::VectorFactory<GO,LO,GO,Node>::Build(reducedMap);
710  targetQcolumns[i]->doImport(*(ghostQcolumns[i]),*importer,Xpetra::INSERT);
711  }
712  RCP<MultiVector> targetQvalues = MultiVectorFactory::Build(reducedMap,NSDim);
713  targetQvalues->doImport(*ghostQvalues,*importer,Xpetra::INSERT);
714 
715  ArrayRCP< ArrayRCP<SC> > targetQvals;
716  ArrayRCP<ArrayRCP<GO> > targetQcols;
717  if (targetQvalues->getLocalLength() > 0) {
718  targetQvals.resize(NSDim);
719  targetQcols.resize(NSDim);
720  for (size_t i=0; i<NSDim; ++i) {
721  targetQvals[i] = targetQvalues->getDataNonConst(i);
722  targetQcols[i] = targetQcolumns[i]->getDataNonConst(0);
723  }
724  }
725 
726  valPtr = Array<SC>(NSDim,0.);
727  globalColPtr = Array<GO>(NSDim,0);
728  for (typename Array<GO>::iterator r=gidsToImport.begin(); r!=gidsToImport.end(); ++r) {
729  if (targetQvalues->getLocalLength() > 0) {
730  for (size_t j=0; j<NSDim; ++j) {
731  valPtr[j] = targetQvals[j][reducedMap->getLocalElement(*r)];
732  globalColPtr[j] = targetQcols[j][reducedMap->getLocalElement(*r)];
733  }
734  Ptentative->insertGlobalValues(*r, globalColPtr.view(0,NSDim), valPtr.view(0,NSDim));
735  } //if (targetQvalues->getLocalLength() > 0)
736  }
737 
738  Ptentative->fillComplete(coarseMap, A->getDomainMap());
739  }
740 
741  template <class Scalar,class LocalOrdinal, class GlobalOrdinal, class Node>
742  bool TentativePFactory<Scalar, LocalOrdinal, GlobalOrdinal, Node>::isGoodMap(const Map& rowMap, const Map& colMap) const {
743  ArrayView<const GO> rowElements = rowMap.getNodeElementList();
744  ArrayView<const GO> colElements = colMap.getNodeElementList();
745 
746  const size_t numElements = rowElements.size();
747 
748  bool goodMap = true;
749  for (size_t i = 0; i < numElements; i++)
750  if (rowElements[i] != colElements[i]) {
751  goodMap = false;
752  break;
753  }
754 
755  return goodMap;
756  }
757 
758 } //namespace MueLu
759 
760 // TODO ReUse: If only P or Nullspace is missing, TentativePFactory can be smart and skip part of the computation.
761 
762 #define MUELU_TENTATIVEPFACTORY_SHORT
763 #endif // MUELU_TENTATIVEPFACTORY_DEF_HPP
Important warning messages (one line)
GlobalOrdinal GO
Timer to be used in factories. Similar to Monitor but with additional timers.
Print more statistics.
void BuildP(Level &fineLevel, Level &coarseLevel) const
Abstract Build method.
LocalOrdinal LO
RCP< const ParameterList > GetValidParameterList() const
Return a const parameter list of valid parameters that setParameterList() will accept.
Namespace for MueLu classes and methods.
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
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
static std::string PrintMatrixInfo(const Matrix &A, const std::string &msgTag, RCP< const Teuchos::ParameterList > params=Teuchos::null)
size_t global_size_t
static RCP< Vector > Build(const Teuchos::RCP< const Map > &map, bool zeroOut=true)
Scalar SC
void Build(Level &fineLevel, Level &coarseLevel) const
Build an object with this factory.
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