MueLu  Version of the Day
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
MueLu_PgPFactory_def.hpp
Go to the documentation of this file.
1 // @HEADER
2 // *****************************************************************************
3 // MueLu: A package for multigrid based preconditioning
4 //
5 // Copyright 2012 NTESS and the MueLu contributors.
6 // SPDX-License-Identifier: BSD-3-Clause
7 // *****************************************************************************
8 // @HEADER
9 
10 #ifndef MUELU_PGPFACTORY_DEF_HPP
11 #define MUELU_PGPFACTORY_DEF_HPP
12 
13 #include <vector>
14 
15 #include <Xpetra_Vector.hpp>
16 #include <Xpetra_MultiVectorFactory.hpp>
17 #include <Xpetra_VectorFactory.hpp>
18 #include <Xpetra_Import.hpp>
19 #include <Xpetra_ImportFactory.hpp>
20 #include <Xpetra_Export.hpp>
21 #include <Xpetra_ExportFactory.hpp>
22 #include <Xpetra_Matrix.hpp>
23 #include <Xpetra_MatrixMatrix.hpp>
24 
26 #include "MueLu_Monitor.hpp"
27 #include "MueLu_PerfUtils.hpp"
29 #include "MueLu_Utilities.hpp"
30 
31 namespace MueLu {
32 
33 template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
35  RCP<ParameterList> validParamList = rcp(new ParameterList());
36 
37  validParamList->set<RCP<const FactoryBase> >("A", Teuchos::null, "Generating factory of the matrix A used during the prolongator smoothing process");
38  validParamList->set<RCP<const FactoryBase> >("P", Teuchos::null, "Tentative prolongator factory");
39  validParamList->set<MinimizationNorm>("Minimization norm", DINVANORM, "Norm to be minimized");
40  validParamList->set<bool>("ReUseRowBasedOmegas", false, "Reuse omegas for prolongator for restrictor");
41 
42  return validParamList;
43 }
44 
45 template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
47  SetParameter("Minimization norm", ParameterEntry(minnorm)); // revalidate
48 }
49 
50 template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
52  const ParameterList& pL = GetParameterList();
53  return pL.get<MueLu::MinimizationNorm>("Minimization norm");
54 }
55 
56 template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
58  Input(fineLevel, "A");
59 
60  // Get default tentative prolongator factory
61  // Getting it that way ensure that the same factory instance will be used for both SaPFactory and NullspaceFactory.
62  // -- Warning: Do not use directly initialPFact_. Use initialPFact instead everywhere!
63  RCP<const FactoryBase> initialPFact = GetFactory("P");
64  if (initialPFact == Teuchos::null) {
65  initialPFact = coarseLevel.GetFactoryManager()->GetFactory("Ptent");
66  }
67  coarseLevel.DeclareInput("P", initialPFact.get(), this);
68 
69  /* If PgPFactory is reusing the row based damping parameters omega for
70  * restriction, it has to request the data here.
71  * we have the following scenarios:
72  * 1) Reuse omegas:
73  * PgPFactory.DeclareInput for prolongation mode requests A and P0
74  * PgPFactory.DeclareInput for restriction mode requests A, P0 and RowBasedOmega (call triggered by GenericRFactory)
75  * PgPFactory.Build for prolongation mode calculates RowBasedOmega and stores it (as requested)
76  * PgPFactory.Build for restriction mode reuses RowBasedOmega (and Releases the data with the Get call)
77  * 2) do not reuse omegas
78  * PgPFactory.DeclareInput for prolongation mode requests A and P0
79  * PgPFactory.DeclareInput for restriction mode requests A and P0
80  * PgPFactory.Build for prolongation mode calculates RowBasedOmega for prolongation operator
81  * PgPFactory.Build for restriction mode calculates RowBasedOmega for restriction operator
82  */
83  const ParameterList& pL = GetParameterList();
84  bool bReUseRowBasedOmegas = pL.get<bool>("ReUseRowBasedOmegas");
85  if (bReUseRowBasedOmegas == true && restrictionMode_ == true) {
86  coarseLevel.DeclareInput("RowBasedOmega", this, this); // RowBasedOmega is calculated by this PgPFactory and requested by this PgPFactory
87  }
88 }
89 
90 template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
92  FactoryMonitor m(*this, "Prolongator smoothing (PG-AMG)", coarseLevel);
93 
94  // Level Get
95  RCP<Matrix> A = Get<RCP<Matrix> >(fineLevel, "A");
96 
97  // Get default tentative prolongator factory
98  // Getting it that way ensure that the same factory instance will be used for both SaPFactory and NullspaceFactory.
99  // -- Warning: Do not use directly initialPFact_. Use initialPFact instead everywhere!
100  RCP<const FactoryBase> initialPFact = GetFactory("P");
101  if (initialPFact == Teuchos::null) {
102  initialPFact = coarseLevel.GetFactoryManager()->GetFactory("Ptent");
103  }
104  RCP<Matrix> Ptent = coarseLevel.Get<RCP<Matrix> >("P", initialPFact.get());
105 
107  if (restrictionMode_) {
108  SubFactoryMonitor m2(*this, "Transpose A", coarseLevel);
109  A = Utilities::Transpose(*A, true); // build transpose of A explicitely
110  }
111 
113  bool doFillComplete = true;
114  bool optimizeStorage = true;
115  RCP<Matrix> DinvAP0 = Xpetra::MatrixMatrix<Scalar, LocalOrdinal, GlobalOrdinal, Node>::Multiply(*A, false, *Ptent, false, GetOStream(Statistics2), doFillComplete, optimizeStorage);
116 
117  const auto rowMap = A->getRowMap();
119  A->getLocalDiagCopy(*diag);
120  diag->reciprocal(*diag);
121  DinvAP0->leftScale(*diag);
122 
124 
125  Teuchos::RCP<Vector> RowBasedOmega = Teuchos::null;
126 
127  const ParameterList& pL = GetParameterList();
128  bool bReUseRowBasedOmegas = pL.get<bool>("ReUseRowBasedOmegas");
129  if (restrictionMode_ == false || bReUseRowBasedOmegas == false) {
130  // if in prolongation mode: calculate row based omegas
131  // if in restriction mode: calculate omegas only if row based omegas are not used from prolongation mode
132  ComputeRowBasedOmega(fineLevel, coarseLevel, A, Ptent, DinvAP0, RowBasedOmega);
133  } // if(bReUseRowBasedOmegas == false)
134  else {
135  // reuse row based omegas, calculated by this factory in the run before (with restrictionMode_ == false)
136  RowBasedOmega = coarseLevel.Get<Teuchos::RCP<Vector> >("RowBasedOmega", this);
137 
138  // RowBasedOmega is now based on row map of A (not transposed)
139  // for restriction we use A^T instead of A
140  // -> recommunicate row based omega
141 
142  // exporter: overlapping row map to nonoverlapping domain map (target map is unique)
143  // since A is already transposed we use the RangeMap of A
144  Teuchos::RCP<const Export> exporter =
145  ExportFactory::Build(RowBasedOmega->getMap(), A->getRangeMap());
146 
147  Teuchos::RCP<Vector> noRowBasedOmega =
148  VectorFactory::Build(A->getRangeMap());
149 
150  noRowBasedOmega->doExport(*RowBasedOmega, *exporter, Xpetra::INSERT);
151 
152  // importer: nonoverlapping map to overlapping map
153 
154  // importer: source -> target maps
155  Teuchos::RCP<const Import> importer =
156  ImportFactory::Build(A->getRangeMap(), A->getRowMap());
157 
158  // doImport target->doImport(*source, importer, action)
159  RowBasedOmega->doImport(*noRowBasedOmega, *importer, Xpetra::INSERT);
160  }
161 
162  DinvAP0->leftScale(*RowBasedOmega);
163 
164  RCP<Matrix> P_smoothed = Teuchos::null;
165 
167  *DinvAP0, false, -Teuchos::ScalarTraits<Scalar>::one(),
168  P_smoothed, GetOStream(Statistics2));
169  P_smoothed->fillComplete(Ptent->getDomainMap(), Ptent->getRangeMap());
170 
172 
173  RCP<ParameterList> params = rcp(new ParameterList());
174  params->set("printLoadBalancingInfo", true);
175 
176  // Level Set
177  if (!restrictionMode_) {
178  // prolongation factory is in prolongation mode
179  Set(coarseLevel, "P", P_smoothed);
180 
181  // RfromPFactory used to indicate to TogglePFactory that a factory
182  // capable or producing R can be invoked later. TogglePFactory
183  // replaces dummy value with an index into it's array of prolongators
184  // pointing to the correct prolongator factory. This is later used by
185  // RfromP_Or_TransP to invoke the prolongatorfactory in RestrictionMode
186  int dummy = 7;
187  Set(coarseLevel, "RfromPfactory", dummy);
188 
189  if (IsPrint(Statistics1))
190  GetOStream(Statistics1) << PerfUtils::PrintMatrixInfo(*P_smoothed, "P", params);
191 
192  // NOTE: EXPERIMENTAL
193  if (Ptent->IsView("stridedMaps"))
194  P_smoothed->CreateView("stridedMaps", Ptent);
195 
196  } else {
197  // prolongation factory is in restriction mode
198  RCP<Matrix> R = Utilities::Transpose(*P_smoothed, true);
199  Set(coarseLevel, "R", R);
200 
201  if (IsPrint(Statistics1))
202  GetOStream(Statistics1) << PerfUtils::PrintMatrixInfo(*R, "P", params);
203 
204  // NOTE: EXPERIMENTAL
205  if (Ptent->IsView("stridedMaps"))
206  R->CreateView("stridedMaps", Ptent, true);
207  }
208 }
209 
210 template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
211 void PgPFactory<Scalar, LocalOrdinal, GlobalOrdinal, Node>::ComputeRowBasedOmega(Level& /* fineLevel */, Level& coarseLevel, const RCP<Matrix>& A, const RCP<Matrix>& P0, const RCP<Matrix>& DinvAP0, RCP<Vector>& RowBasedOmega) const {
212  FactoryMonitor m(*this, "PgPFactory::ComputeRowBasedOmega", coarseLevel);
213 
214  typedef typename Teuchos::ScalarTraits<Scalar>::magnitudeType Magnitude;
216  Magnitude mZero = Teuchos::ScalarTraits<Scalar>::magnitude(sZero);
217 
218  Teuchos::RCP<Vector> Numerator = Teuchos::null;
219  Teuchos::RCP<Vector> Denominator = Teuchos::null;
220 
221  const ParameterList& pL = GetParameterList();
222  MinimizationNorm min_norm = pL.get<MinimizationNorm>("Minimization norm");
223 
224  switch (min_norm) {
225  case ANORM: {
226  // MUEMAT mode (=paper)
227  // Minimize with respect to the (A)' A norm.
228  // Need to be smart here to avoid the construction of A' A
229  //
230  // diag( P0' (A' A) D^{-1} A P0)
231  // omega = ------------------------------------------
232  // diag( P0' A' D^{-1}' ( A' A) D^{-1} A P0)
233  //
234  // expensive, since we have to recalculate AP0 due to the lack of an explicit scaling routine for DinvAP0
235 
236  // calculate A * P0
237  bool doFillComplete = true;
238  bool optimizeStorage = false;
239  RCP<Matrix> AP0 = Xpetra::MatrixMatrix<Scalar, LocalOrdinal, GlobalOrdinal, Node>::Multiply(*A, false, *P0, false, GetOStream(Statistics2), doFillComplete, optimizeStorage);
240 
241  // compute A * D^{-1} * A * P0
242  RCP<Matrix> ADinvAP0 = Xpetra::MatrixMatrix<Scalar, LocalOrdinal, GlobalOrdinal, Node>::Multiply(*A, false, *DinvAP0, false, GetOStream(Statistics2), doFillComplete, optimizeStorage);
243 
244  Numerator = VectorFactory::Build(ADinvAP0->getColMap(), true);
245  Denominator = VectorFactory::Build(ADinvAP0->getColMap(), true);
246  MultiplyAll(AP0, ADinvAP0, Numerator);
247  MultiplySelfAll(ADinvAP0, Denominator);
248  } break;
249  case L2NORM: {
250  // ML mode 1 (cheapest)
251  // Minimize with respect to L2 norm
252  // diag( P0' D^{-1} A P0)
253  // omega = -----------------------------
254  // diag( P0' A' D^{-1}' D^{-1} A P0)
255  //
256  Numerator = VectorFactory::Build(DinvAP0->getColMap(), true);
257  Denominator = VectorFactory::Build(DinvAP0->getColMap(), true);
258  MultiplyAll(P0, DinvAP0, Numerator);
259  MultiplySelfAll(DinvAP0, Denominator);
260  } break;
261  case DINVANORM: {
262  // ML mode 2
263  // Minimize with respect to the (D^{-1} A)' D^{-1} A norm.
264  // Need to be smart here to avoid the construction of A' A
265  //
266  // diag( P0' ( A' D^{-1}' D^{-1} A) D^{-1} A P0)
267  // omega = --------------------------------------------------------
268  // diag( P0' A' D^{-1}' ( A' D^{-1}' D^{-1} A) D^{-1} A P0)
269  //
270 
271  // compute D^{-1} * A * D^{-1} * A * P0
272  bool doFillComplete = true;
273  bool optimizeStorage = true;
275  RCP<Matrix> DinvADinvAP0 = Xpetra::MatrixMatrix<Scalar, LocalOrdinal, GlobalOrdinal, Node>::Multiply(*A, false, *DinvAP0, false, GetOStream(Statistics2), doFillComplete, optimizeStorage);
276 
277  const auto rowMap = A->getRowMap();
279  A->getLocalDiagCopy(*diag);
280  diag->reciprocal(*diag);
281  DinvADinvAP0->leftScale(*diag);
282 
283  Numerator = VectorFactory::Build(DinvADinvAP0->getColMap(), true);
284  Denominator = VectorFactory::Build(DinvADinvAP0->getColMap(), true);
285  MultiplyAll(DinvAP0, DinvADinvAP0, Numerator);
286  MultiplySelfAll(DinvADinvAP0, Denominator);
287  } break;
288  default:
289  TEUCHOS_TEST_FOR_EXCEPTION(true, Exceptions::RuntimeError, "MueLu::PgPFactory::Build: minimization mode not supported. error");
290  }
291 
293  Teuchos::RCP<Vector> ColBasedOmega =
294  VectorFactory::Build(Numerator->getMap() /*DinvAP0->getColMap()*/, true);
295 
296  ColBasedOmega->putScalar(-666 );
297 
298  Teuchos::ArrayRCP<const Scalar> Numerator_local = Numerator->getData(0);
299  Teuchos::ArrayRCP<const Scalar> Denominator_local = Denominator->getData(0);
300  Teuchos::ArrayRCP<Scalar> ColBasedOmega_local = ColBasedOmega->getDataNonConst(0);
301  GlobalOrdinal zero_local = 0; // count negative colbased omegas
302  GlobalOrdinal nan_local = 0; // count NaNs -> set them to zero
303  Magnitude min_local = 1000000.0; // Teuchos::ScalarTraits<Scalar>::one() * (Scalar) 1000000;
304  Magnitude max_local = 0.0;
305  for (LocalOrdinal i = 0; i < Teuchos::as<LocalOrdinal>(Numerator->getLocalLength()); i++) {
306  if (Teuchos::ScalarTraits<Scalar>::magnitude(Denominator_local[i]) == mZero) {
307  ColBasedOmega_local[i] = 0.0; // fallback: nonsmoothed basis function since denominator == 0.0
308  nan_local++;
309  } else {
310  ColBasedOmega_local[i] = Numerator_local[i] / Denominator_local[i]; // default case
311  }
312 
313  if (Teuchos::ScalarTraits<Scalar>::magnitude(ColBasedOmega_local[i]) < mZero) { // negative omegas are not valid. set them to zero
314  ColBasedOmega_local[i] = Teuchos::ScalarTraits<Scalar>::zero();
315  zero_local++; // count zero omegas
316  }
317 
318  // handle case that Nominator == Denominator -> Dirichlet bcs in A?
319  // fallback if ColBasedOmega == 1 -> very strong smoothing may lead to zero rows in P
320  // TAW: this is somewhat nonstandard and a rough fallback strategy to avoid problems
321  // also avoid "overshooting" with omega > 0.8
322  if (Teuchos::ScalarTraits<Scalar>::magnitude(ColBasedOmega_local[i]) >= 0.8) {
323  ColBasedOmega_local[i] = 0.0;
324  }
325 
326  if (Teuchos::ScalarTraits<Scalar>::magnitude(ColBasedOmega_local[i]) < min_local) {
327  min_local = Teuchos::ScalarTraits<Scalar>::magnitude(ColBasedOmega_local[i]);
328  }
329  if (Teuchos::ScalarTraits<Scalar>::magnitude(ColBasedOmega_local[i]) > max_local) {
330  max_local = Teuchos::ScalarTraits<Scalar>::magnitude(ColBasedOmega_local[i]);
331  }
332  }
333 
334  { // be verbose
335  GlobalOrdinal zero_all;
336  GlobalOrdinal nan_all;
337  Magnitude min_all;
338  Magnitude max_all;
339  MueLu_sumAll(A->getRowMap()->getComm(), zero_local, zero_all);
340  MueLu_sumAll(A->getRowMap()->getComm(), nan_local, nan_all);
341  MueLu_minAll(A->getRowMap()->getComm(), min_local, min_all);
342  MueLu_maxAll(A->getRowMap()->getComm(), max_local, max_all);
343 
344  GetOStream(MueLu::Statistics1, 0) << "PgPFactory: smoothed aggregation (scheme: ";
345  switch (min_norm) {
346  case ANORM: GetOStream(Statistics1) << "Anorm)" << std::endl; break;
347  case L2NORM: GetOStream(Statistics1) << "L2norm)" << std::endl; break;
348  case DINVANORM: GetOStream(Statistics1) << "DinvAnorm)" << std::endl; break;
349  default: GetOStream(Statistics1) << "unknown)" << std::endl; break;
350  }
351  GetOStream(Statistics1) << "Damping parameter: min = " << min_all << ", max = " << max_all << std::endl;
352  GetOStream(Statistics) << "# negative omegas: " << zero_all << " out of " << ColBasedOmega->getGlobalLength() << " column-based omegas" << std::endl;
353  GetOStream(Statistics) << "# NaNs: " << nan_all << " out of " << ColBasedOmega->getGlobalLength() << " column-based omegas" << std::endl;
354  }
355 
356  if (coarseLevel.IsRequested("ColBasedOmega", this)) {
357  coarseLevel.Set("ColBasedOmega", ColBasedOmega, this);
358  }
359 
361  // transform column based omegas to row based omegas
362  RowBasedOmega =
363  VectorFactory::Build(DinvAP0->getRowMap(), true);
364 
365  RowBasedOmega->putScalar(-666); // TODO bad programming style
366 
367  bool bAtLeastOneDefined = false;
368  Teuchos::ArrayRCP<Scalar> RowBasedOmega_local = RowBasedOmega->getDataNonConst(0);
369  for (LocalOrdinal row = 0; row < Teuchos::as<LocalOrdinal>(A->getLocalNumRows()); row++) {
372  DinvAP0->getLocalRowView(row, lindices, lvals);
373  bAtLeastOneDefined = false;
374  for (size_t j = 0; j < Teuchos::as<size_t>(lindices.size()); j++) {
375  Scalar omega = ColBasedOmega_local[lindices[j]];
376  if (Teuchos::ScalarTraits<Scalar>::magnitude(omega) != -666) { // TODO bad programming style
377  bAtLeastOneDefined = true;
378  if (Teuchos::ScalarTraits<Scalar>::magnitude(RowBasedOmega_local[row]) == -666)
379  RowBasedOmega_local[row] = omega;
381  RowBasedOmega_local[row] = omega;
382  }
383  }
384  if (bAtLeastOneDefined == true) {
385  if (Teuchos::ScalarTraits<Scalar>::magnitude(RowBasedOmega_local[row]) < mZero)
386  RowBasedOmega_local[row] = sZero;
387  }
388  }
389 
390  if (coarseLevel.IsRequested("RowBasedOmega", this)) {
391  Set(coarseLevel, "RowBasedOmega", RowBasedOmega);
392  }
393 }
394 
395 template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
397  // note: InnerProdVec is based on column map of Op
398  TEUCHOS_TEST_FOR_EXCEPTION(!InnerProdVec->getMap()->isSameAs(*Op->getColMap()), Exceptions::RuntimeError, "MueLu::PgPFactory::MultiplySelfAll: map of InnerProdVec must be same as column map of operator. error");
399 
400  Teuchos::ArrayRCP<Scalar> InnerProd_local = InnerProdVec->getDataNonConst(0);
401 
404 
405  for (size_t n = 0; n < Op->getLocalNumRows(); n++) {
406  Op->getLocalRowView(n, lindices, lvals);
407  for (size_t i = 0; i < Teuchos::as<size_t>(lindices.size()); i++) {
408  InnerProd_local[lindices[i]] += lvals[i] * lvals[i];
409  }
410  }
411  InnerProd_local = Teuchos::ArrayRCP<Scalar>();
412 
413  // exporter: overlapping map to nonoverlapping map (target map is unique)
414  Teuchos::RCP<const Export> exporter =
415  ExportFactory::Build(Op->getColMap(), Op->getDomainMap());
416 
417  Teuchos::RCP<Vector> nonoverlap =
418  VectorFactory::Build(Op->getDomainMap());
419 
420  nonoverlap->doExport(*InnerProdVec, *exporter, Xpetra::ADD);
421 
422  // importer: nonoverlapping map to overlapping map
423 
424  // importer: source -> target maps
425  Teuchos::RCP<const Import> importer =
426  ImportFactory::Build(Op->getDomainMap(), Op->getColMap());
427 
428  // doImport target->doImport(*source, importer, action)
429  InnerProdVec->doImport(*nonoverlap, *importer, Xpetra::INSERT);
430 }
431 
432 template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
434  TEUCHOS_TEST_FOR_EXCEPTION(!left->getDomainMap()->isSameAs(*right->getDomainMap()), Exceptions::RuntimeError, "MueLu::PgPFactory::MultiplyAll: domain maps of left and right do not match. Error.");
435  TEUCHOS_TEST_FOR_EXCEPTION(!left->getRowMap()->isSameAs(*right->getRowMap()), Exceptions::RuntimeError, "MueLu::PgPFactory::MultiplyAll: row maps of left and right do not match. Error.");
436 #if 1 // 1=new "fast code, 0=old "slow", but safe code
437 #if 0 // not necessary - remove me
438  if(InnerProdVec->getMap()->isSameAs(*left->getColMap())) {
439  // initialize NewRightLocal vector and assign all entries to
440  // left->getColMap()->getLocalNumElements() + 1
441  std::vector<LocalOrdinal> NewRightLocal(right->getColMap()->getLocalNumElements(), Teuchos::as<LocalOrdinal>(left->getColMap()->getLocalNumElements()+1));
442 
443  LocalOrdinal i = 0;
444  for (size_t j=0; j < right->getColMap()->getLocalNumElements(); j++) {
445  while ( (i < Teuchos::as<LocalOrdinal>(left->getColMap()->getLocalNumElements())) &&
446  (left->getColMap()->getGlobalElement(i) < right->getColMap()->getGlobalElement(j)) ) i++;
447  if (left->getColMap()->getGlobalElement(i) == right->getColMap()->getGlobalElement(j)) {
448  NewRightLocal[j] = i;
449  }
450  }
451 
452  Teuchos::ArrayRCP< Scalar > InnerProd_local = InnerProdVec->getDataNonConst(0);
453  std::vector<Scalar> temp_array(left->getColMap()->getLocalNumElements()+1, 0.0);
454 
455  for(size_t n=0; n<right->getLocalNumRows(); n++) {
460 
461  left->getLocalRowView (n, lindices_left, lvals_left);
462  right->getLocalRowView(n, lindices_right, lvals_right);
463 
464  for(size_t j=0; j<Teuchos::as<size_t>(lindices_right.size()); j++) {
465  temp_array[NewRightLocal[lindices_right[j] ] ] = lvals_right[j];
466  }
467  for (size_t j=0; j < Teuchos::as<size_t>(lindices_left.size()); j++) {
468  InnerProd_local[lindices_left[j]] += temp_array[lindices_left[j] ]*lvals_left[j];
469  }
470  for (size_t j=0; j < Teuchos::as<size_t>(lindices_right.size()); j++) {
471  temp_array[NewRightLocal[lindices_right[j] ] ] = 0.0;
472  }
473  }
474  // exporter: overlapping map to nonoverlapping map (target map is unique)
475  Teuchos::RCP<const Export> exporter =
476  ExportFactory::Build(left->getColMap(), left->getDomainMap()); // TODO: change left to right?
477 
478  Teuchos::RCP<Vector > nonoverlap =
479  VectorFactory::Build(left->getDomainMap()); // TODO: change left to right?
480 
481  nonoverlap->doExport(*InnerProdVec, *exporter, Xpetra::ADD);
482 
483  // importer: nonoverlapping map to overlapping map
484 
485  // importer: source -> target maps
486  Teuchos::RCP<const Import > importer =
487  ImportFactory::Build(left->getDomainMap(), left->getColMap()); // TODO: change left to right?
488 
489  // doImport target->doImport(*source, importer, action)
490  InnerProdVec->doImport(*nonoverlap, *importer, Xpetra::INSERT);
491 
492 
493  } else
494 #endif // end remove me
495  if (InnerProdVec->getMap()->isSameAs(*right->getColMap())) {
496  size_t szNewLeftLocal = TEUCHOS_MAX(left->getColMap()->getLocalNumElements(), right->getColMap()->getLocalNumElements());
497  Teuchos::RCP<std::vector<LocalOrdinal> > NewLeftLocal = Teuchos::rcp(new std::vector<LocalOrdinal>(szNewLeftLocal, Teuchos::as<LocalOrdinal>(right->getColMap()->getMaxLocalIndex() + 1)));
498 
499  LocalOrdinal j = 0;
500  for (size_t i = 0; i < left->getColMap()->getLocalNumElements(); i++) {
501  while ((j < Teuchos::as<LocalOrdinal>(right->getColMap()->getLocalNumElements())) &&
502  (right->getColMap()->getGlobalElement(j) < left->getColMap()->getGlobalElement(i))) j++;
503  if (right->getColMap()->getGlobalElement(j) == left->getColMap()->getGlobalElement(i)) {
504  (*NewLeftLocal)[i] = j;
505  }
506  }
507 
508  /*for (size_t i=0; i < right->getColMap()->getLocalNumElements(); i++) {
509  std::cout << "left col map: " << (*NewLeftLocal)[i] << " GID: " << left->getColMap()->getGlobalElement((*NewLeftLocal)[i]) << " GID: " << right->getColMap()->getGlobalElement(i) << " right col map: " << i << std::endl;
510  }*/
511 
512  Teuchos::ArrayRCP<Scalar> InnerProd_local = InnerProdVec->getDataNonConst(0);
513  Teuchos::RCP<std::vector<Scalar> > temp_array = Teuchos::rcp(new std::vector<Scalar>(right->getColMap()->getMaxLocalIndex() + 2, 0.0));
514 
515  for (size_t n = 0; n < left->getLocalNumRows(); n++) {
520 
521  left->getLocalRowView(n, lindices_left, lvals_left);
522  right->getLocalRowView(n, lindices_right, lvals_right);
523 
524  for (size_t i = 0; i < Teuchos::as<size_t>(lindices_left.size()); i++) {
525  (*temp_array)[(*NewLeftLocal)[lindices_left[i]]] = lvals_left[i];
526  }
527  for (size_t i = 0; i < Teuchos::as<size_t>(lindices_right.size()); i++) {
528  InnerProd_local[lindices_right[i]] += (*temp_array)[lindices_right[i]] * lvals_right[i];
529  }
530  for (size_t i = 0; i < Teuchos::as<size_t>(lindices_left.size()); i++) {
531  (*temp_array)[(*NewLeftLocal)[lindices_left[i]]] = 0.0;
532  }
533  }
534  InnerProd_local = Teuchos::ArrayRCP<Scalar>();
535  // exporter: overlapping map to nonoverlapping map (target map is unique)
536  Teuchos::RCP<const Export> exporter =
537  ExportFactory::Build(right->getColMap(), right->getDomainMap()); // TODO: change left to right?
538 
539  Teuchos::RCP<Vector> nonoverlap =
540  VectorFactory::Build(right->getDomainMap()); // TODO: change left to right?
541 
542  nonoverlap->doExport(*InnerProdVec, *exporter, Xpetra::ADD);
543 
544  // importer: nonoverlapping map to overlapping map
545 
546  // importer: source -> target maps
547  Teuchos::RCP<const Import> importer =
548  ImportFactory::Build(right->getDomainMap(), right->getColMap()); // TODO: change left to right?
549  // doImport target->doImport(*source, importer, action)
550  InnerProdVec->doImport(*nonoverlap, *importer, Xpetra::INSERT);
551  } else {
552  TEUCHOS_TEST_FOR_EXCEPTION(true, Exceptions::RuntimeError, "MueLu::PgPFactory::MultiplyAll: map of InnerProdVec must be same as column map of left operator? Error.");
553  }
554 
555 #else // old "safe" code
556  if (InnerProdVec->getMap()->isSameAs(*left->getColMap())) {
557  Teuchos::ArrayRCP<Scalar> InnerProd_local = InnerProdVec->getDataNonConst(0);
558 
559  // declare variables
564 
565  for (size_t n = 0; n < left->getLocalNumRows(); n++) {
566  left->getLocalRowView(n, lindices_left, lvals_left);
567  right->getLocalRowView(n, lindices_right, lvals_right);
568 
569  for (size_t i = 0; i < Teuchos::as<size_t>(lindices_left.size()); i++) {
570  GlobalOrdinal left_gid = left->getColMap()->getGlobalElement(lindices_left[i]);
571  for (size_t j = 0; j < Teuchos::as<size_t>(lindices_right.size()); j++) {
572  GlobalOrdinal right_gid = right->getColMap()->getGlobalElement(lindices_right[j]);
573  if (left_gid == right_gid) {
574  InnerProd_local[lindices_left[i]] += lvals_left[i] * lvals_right[j];
575  break; // skip remaining gids of right operator
576  }
577  }
578  }
579  }
580 
581  // exporter: overlapping map to nonoverlapping map (target map is unique)
582  Teuchos::RCP<const Export> exporter =
583  ExportFactory::Build(left->getColMap(), left->getDomainMap()); // TODO: change left to right?
584 
585  Teuchos::RCP<Vector> nonoverlap =
586  VectorFactory::Build(left->getDomainMap()); // TODO: change left to right?
587 
588  nonoverlap->doExport(*InnerProdVec, *exporter, Xpetra::ADD);
589 
590  // importer: nonoverlapping map to overlapping map
591 
592  // importer: source -> target maps
593  Teuchos::RCP<const Import> importer =
594  ImportFactory::Build(left->getDomainMap(), left->getColMap()); // TODO: change left to right?
595 
596  // doImport target->doImport(*source, importer, action)
597  InnerProdVec->doImport(*nonoverlap, *importer, Xpetra::INSERT);
598  } else if (InnerProdVec->getMap()->isSameAs(*right->getColMap())) {
599  Teuchos::ArrayRCP<Scalar> InnerProd_local = InnerProdVec->getDataNonConst(0);
600 
605 
606  for (size_t n = 0; n < left->getLocalNumRows(); n++) {
607  left->getLocalRowView(n, lindices_left, lvals_left);
608  right->getLocalRowView(n, lindices_right, lvals_right);
609 
610  for (size_t i = 0; i < Teuchos::as<size_t>(lindices_left.size()); i++) {
611  GlobalOrdinal left_gid = left->getColMap()->getGlobalElement(lindices_left[i]);
612  for (size_t j = 0; j < Teuchos::as<size_t>(lindices_right.size()); j++) {
613  GlobalOrdinal right_gid = right->getColMap()->getGlobalElement(lindices_right[j]);
614  if (left_gid == right_gid) {
615  InnerProd_local[lindices_right[j]] += lvals_left[i] * lvals_right[j];
616  break; // skip remaining gids of right operator
617  }
618  }
619  }
620  }
621 
622  // exporter: overlapping map to nonoverlapping map (target map is unique)
623  Teuchos::RCP<const Export> exporter =
624  ExportFactory::Build(right->getColMap(), right->getDomainMap()); // TODO: change left to right?
625 
626  Teuchos::RCP<Vector> nonoverlap =
627  VectorFactory::Build(right->getDomainMap()); // TODO: change left to right?
628 
629  nonoverlap->doExport(*InnerProdVec, *exporter, Xpetra::ADD);
630 
631  // importer: nonoverlapping map to overlapping map
632 
633  // importer: source -> target maps
634  Teuchos::RCP<const Import> importer =
635  ImportFactory::Build(right->getDomainMap(), right->getColMap()); // TODO: change left to right?
636 
637  // doImport target->doImport(*source, importer, action)
638  InnerProdVec->doImport(*nonoverlap, *importer, Xpetra::INSERT);
639  } else {
640  TEUCHOS_TEST_FOR_EXCEPTION(true, Exceptions::RuntimeError, "MueLu::PgPFactory::MultiplyAll: map of InnerProdVec must be same as column map of left or right operator? Error.");
641  }
642 #endif
643 }
644 
645 template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
646 void PgPFactory<Scalar, LocalOrdinal, GlobalOrdinal, Node>::BuildP(Level& /* fineLevel */, Level& /* coarseLevel */) const {
647  std::cout << "TODO: remove me" << std::endl;
648 }
649 
650 template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
652  SetParameter("ReUseRowBasedOmegas", ParameterEntry(bReuse));
653 }
654 
655 } // namespace MueLu
656 
657 #endif /* MUELU_PGPFACTORY_DEF_HPP */
#define MueLu_sumAll(rcpComm, in, out)
MueLu::DefaultLocalOrdinal LocalOrdinal
T & Get(const std::string &ename, const FactoryBase *factory=NoFactory::get())
Get data without decrementing associated storage counter (i.e., read-only access). Usage: Level-&gt;Get&lt; RCP&lt;Matrix&gt; &gt;(&quot;A&quot;, factory) if factory == NULL =&gt; use default factory.
#define MueLu_maxAll(rcpComm, in, out)
void DeclareInput(Level &fineLevel, Level &coarseLevel) const
Input.
T & get(const std::string &name, T def_value)
void BuildP(Level &fineLevel, Level &coarseLevel) const
Abstract Build method.
void ComputeRowBasedOmega(Level &fineLevel, Level &coarseLevel, const RCP< Matrix > &A, const RCP< Matrix > &P0, const RCP< Matrix > &DinvAP0, RCP< Xpetra::Vector< Scalar, LocalOrdinal, GlobalOrdinal, Node > > &RowBasedOmega) const
Timer to be used in factories. Similar to Monitor but with additional timers.
#define TEUCHOS_TEST_FOR_EXCEPTION(throw_exception_test, Exception, msg)
Print more statistics.
size_type size() const
T * get() const
ParameterList & set(std::string const &name, T &&value, std::string const &docString="", RCP< const ParameterEntryValidator > const &validator=null)
#define MueLu_minAll(rcpComm, in, out)
Print all statistics.
Print even more statistics.
void SetMinimizationMode(MinimizationNorm minnorm)
Set minimization mode (L2NORM for cheapest, ANORM more expensive, DINVANORM = default) ...
TEUCHOS_DEPRECATED RCP< T > rcp(T *p, Dealloc_T dealloc, bool owns_mem)
MueLu::DefaultScalar Scalar
MueLu::DefaultGlobalOrdinal GlobalOrdinal
void Build(Level &fineLevel, Level &coarseLevel) const
Build method.
Class that holds all level-specific information.
Definition: MueLu_Level.hpp:63
Timer to be used in factories. Similar to SubMonitor but adds a timer level by level.
void MultiplySelfAll(const RCP< Matrix > &Op, Teuchos::RCP< Xpetra::Vector< Scalar, LocalOrdinal, GlobalOrdinal, Node > > &InnerProdVec) const
static Teuchos::ArrayRCP< Scalar > GetMatrixDiagonal_arcp(const Matrix &A)
Extract Matrix Diagonal.
static RCP< Xpetra::Matrix< Scalar, LocalOrdinal, GlobalOrdinal, Node > > Transpose(Xpetra::Matrix< Scalar, LocalOrdinal, GlobalOrdinal, Node > &Op, bool optimizeTranspose=false, const std::string &label=std::string(), const Teuchos::RCP< Teuchos::ParameterList > &params=Teuchos::null)
static std::string PrintMatrixInfo(const Matrix &A, const std::string &msgTag, RCP< const Teuchos::ParameterList > params=Teuchos::null)
MinimizationNorm GetMinimizationMode()
return minimization mode
void Set(const std::string &ename, const T &entry, const FactoryBase *factory=NoFactory::get())
static RCP< Vector > Build(const Teuchos::RCP< const Map > &map, bool zeroOut=true)
static void TwoMatrixAdd(const Matrix &A, bool transposeA, SC alpha, Matrix &B, SC beta)
static magnitudeType magnitude(T a)
void ReUseDampingParameters(bool bReuse)
RCP< const ParameterList > GetValidParameterList() const
Return a const parameter list of valid parameters that setParameterList() will accept.
static void Multiply(const Matrix &A, bool transposeA, const Matrix &B, bool transposeB, Matrix &C, bool call_FillComplete_on_result=true, bool doOptimizeStorage=true, const std::string &label=std::string(), const RCP< ParameterList > &params=null)
const RCP< const FactoryManagerBase > GetFactoryManager()
returns the current factory manager
Definition: MueLu_Level.cpp:73
Exception throws to report errors in the internal logical of the program.
bool IsRequested(const std::string &ename, const FactoryBase *factory=NoFactory::get()) const
Test whether a need has been requested. Note: this tells nothing about whether the need&#39;s value exist...
void DeclareInput(const std::string &ename, const FactoryBase *factory, const FactoryBase *requestedBy=NoFactory::get())
Callback from FactoryBase::CallDeclareInput() and FactoryBase::DeclareInput()
void MultiplyAll(const RCP< Matrix > &left, const RCP< Matrix > &right, Teuchos::RCP< Xpetra::Vector< Scalar, LocalOrdinal, GlobalOrdinal, Node > > &InnerProdVec) const