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