MueLu  Version of the Day
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 
62 #include "MueLu_Monitor.hpp"
63 #include "MueLu_PerfUtils.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>
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>
86  SetParameter("Minimization norm", ParameterEntry(minnorm)); // revalidate
87  }
88 
89  template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
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>
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>
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=true;
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  diagA = Teuchos::ArrayRCP<Scalar>();
311 
312  Numerator = VectorFactory::Build(DinvADinvAP0->getColMap(), true);
313  Denominator = VectorFactory::Build(DinvADinvAP0->getColMap(), true);
314  MultiplyAll(DinvAP0, DinvADinvAP0, Numerator);
315  MultiplySelfAll(DinvADinvAP0, Denominator);
316  }
317  break;
318  default:
319  TEUCHOS_TEST_FOR_EXCEPTION(true, Exceptions::RuntimeError, "MueLu::PgPFactory::Build: minimization mode not supported. error");
320  }
321 
322 
324  Teuchos::RCP<Vector > ColBasedOmega =
325  VectorFactory::Build(Numerator->getMap()/*DinvAP0->getColMap()*/, true);
326 
327  ColBasedOmega->putScalar(-666);
328 
329  Teuchos::ArrayRCP< const Scalar > Numerator_local = Numerator->getData(0);
330  Teuchos::ArrayRCP< const Scalar > Denominator_local = Denominator->getData(0);
331  Teuchos::ArrayRCP< Scalar > ColBasedOmega_local = ColBasedOmega->getDataNonConst(0);
332  GlobalOrdinal zero_local = 0; // count negative colbased omegas
333  GlobalOrdinal nan_local = 0; // count NaNs -> set them to zero
334  Magnitude min_local = 1000000.0; //Teuchos::ScalarTraits<Scalar>::one() * (Scalar) 1000000;
335  Magnitude max_local = 0.0;
336  for(LocalOrdinal i = 0; i < Teuchos::as<LocalOrdinal>(Numerator->getLocalLength()); i++) {
337  if(Teuchos::ScalarTraits<Scalar>::magnitude(Denominator_local[i]) == mZero)
338  {
339  ColBasedOmega_local[i] = 0.0; // fallback: nonsmoothed basis function since denominator == 0.0
340  nan_local++;
341  }
342  else
343  {
344  ColBasedOmega_local[i] = Numerator_local[i] / Denominator_local[i]; // default case
345  }
346 
347  if(Teuchos::ScalarTraits<Scalar>::magnitude(ColBasedOmega_local[i]) < mZero) { // negative omegas are not valid. set them to zero
348  ColBasedOmega_local[i] = Teuchos::ScalarTraits<Scalar>::zero();
349  zero_local++; // count zero omegas
350  }
351 
352  // handle case that Nominator == Denominator -> Dirichlet bcs in A?
353  // fallback if ColBasedOmega == 1 -> very strong smoothing may lead to zero rows in P
354  // TAW: this is somewhat nonstandard and a rough fallback strategy to avoid problems
355  // also avoid "overshooting" with omega > 0.8
356  if(Teuchos::ScalarTraits<Scalar>::magnitude(ColBasedOmega_local[i]) >= 0.8) {
357  ColBasedOmega_local[i] = 0.0;
358  }
359 
360  if(Teuchos::ScalarTraits<Scalar>::magnitude(ColBasedOmega_local[i]) < min_local)
361  { min_local = Teuchos::ScalarTraits<Scalar>::magnitude(ColBasedOmega_local[i]); }
362  if(Teuchos::ScalarTraits<Scalar>::magnitude(ColBasedOmega_local[i]) > max_local)
363  { max_local = Teuchos::ScalarTraits<Scalar>::magnitude(ColBasedOmega_local[i]); }
364  }
365 
366  { // be verbose
367  GlobalOrdinal zero_all;
368  GlobalOrdinal nan_all;
369  Magnitude min_all;
370  Magnitude max_all;
371  MueLu_sumAll(A->getRowMap()->getComm(), zero_local, zero_all);
372  MueLu_sumAll(A->getRowMap()->getComm(), nan_local, nan_all);
373  MueLu_minAll(A->getRowMap()->getComm(), min_local, min_all);
374  MueLu_maxAll(A->getRowMap()->getComm(), max_local, max_all);
375 
376  GetOStream(MueLu::Statistics1, 0) << "PgPFactory: smoothed aggregation (scheme: ";
377  switch (min_norm) {
378  case ANORM: GetOStream(Statistics1) << "Anorm)" << std::endl; break;
379  case L2NORM: GetOStream(Statistics1) << "L2norm)" << std::endl; break;
380  case DINVANORM: GetOStream(Statistics1) << "DinvAnorm)" << std::endl; break;
381  default: GetOStream(Statistics1) << "unknown)" << std::endl; break;
382  }
383  GetOStream(Statistics1) << "Damping parameter: min = " << min_all << ", max = " << max_all << std::endl;
384  GetOStream(Statistics) << "# negative omegas: " << zero_all << " out of " << ColBasedOmega->getGlobalLength() << " column-based omegas" << std::endl;
385  GetOStream(Statistics) << "# NaNs: " << nan_all << " out of " << ColBasedOmega->getGlobalLength() << " column-based omegas" << std::endl;
386  }
387 
388  if(coarseLevel.IsRequested("ColBasedOmega", this)) {
389  coarseLevel.Set("ColBasedOmega", ColBasedOmega, this);
390  }
391 
393  // transform column based omegas to row based omegas
394  RowBasedOmega =
395  VectorFactory::Build(DinvAP0->getRowMap(), true);
396 
397  RowBasedOmega->putScalar(-666); // TODO bad programming style
398 
399  bool bAtLeastOneDefined = false;
400  Teuchos::ArrayRCP< Scalar > RowBasedOmega_local = RowBasedOmega->getDataNonConst(0);
401  for(LocalOrdinal row = 0; row<Teuchos::as<LocalOrdinal>(A->getNodeNumRows()); row++) {
402  Teuchos::ArrayView<const LocalOrdinal> lindices;
403  Teuchos::ArrayView<const Scalar> lvals;
404  DinvAP0->getLocalRowView(row, lindices, lvals);
405  bAtLeastOneDefined = false;
406  for(size_t j=0; j<Teuchos::as<size_t>(lindices.size()); j++) {
407  Scalar omega = ColBasedOmega_local[lindices[j]];
408  if (Teuchos::ScalarTraits<Scalar>::magnitude(omega) != -666) { // TODO bad programming style
409  bAtLeastOneDefined = true;
410  if(Teuchos::ScalarTraits<Scalar>::magnitude(RowBasedOmega_local[row]) == -666) RowBasedOmega_local[row] = omega;
411  else if(Teuchos::ScalarTraits<Scalar>::magnitude(omega) < Teuchos::ScalarTraits<Scalar>::magnitude(RowBasedOmega_local[row])) RowBasedOmega_local[row] = omega;
412  }
413  }
414  if(bAtLeastOneDefined == true) {
415  if(Teuchos::ScalarTraits<Scalar>::magnitude(RowBasedOmega_local[row]) < mZero)
416  RowBasedOmega_local[row] = sZero;
417  }
418  }
419 
420  if(coarseLevel.IsRequested("RowBasedOmega", this)) {
421  Set(coarseLevel, "RowBasedOmega", RowBasedOmega);
422  }
423  }
424 
425  template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
426  void PgPFactory<Scalar, LocalOrdinal, GlobalOrdinal, Node>::MultiplySelfAll(const RCP<Matrix>& Op, Teuchos::RCP<Vector >& InnerProdVec) const {
427 
428  // note: InnerProdVec is based on column map of Op
429  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");
430 
431  Teuchos::ArrayRCP< Scalar > InnerProd_local = InnerProdVec->getDataNonConst(0);
432 
433  Teuchos::ArrayView<const LocalOrdinal> lindices;
434  Teuchos::ArrayView<const Scalar> lvals;
435 
436  for(size_t n=0; n<Op->getNodeNumRows(); n++) {
437  Op->getLocalRowView(n, lindices, lvals);
438  for(size_t i=0; i<Teuchos::as<size_t>(lindices.size()); i++) {
439  InnerProd_local[lindices[i]] += lvals[i]*lvals[i];
440  }
441  }
442  InnerProd_local = Teuchos::ArrayRCP< Scalar >();
443 
444  // exporter: overlapping map to nonoverlapping map (target map is unique)
445  Teuchos::RCP<const Export> exporter =
446  ExportFactory::Build(Op->getColMap(), Op->getDomainMap());
447 
448  Teuchos::RCP<Vector > nonoverlap =
449  VectorFactory::Build(Op->getDomainMap());
450 
451  nonoverlap->doExport(*InnerProdVec, *exporter, Xpetra::ADD);
452 
453  // importer: nonoverlapping map to overlapping map
454 
455  // importer: source -> target maps
456  Teuchos::RCP<const Import > importer =
457  ImportFactory::Build(Op->getDomainMap(), Op->getColMap());
458 
459  // doImport target->doImport(*source, importer, action)
460  InnerProdVec->doImport(*nonoverlap, *importer, Xpetra::INSERT);
461 
462  }
463 
464  template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
465  void PgPFactory<Scalar, LocalOrdinal, GlobalOrdinal, Node>::MultiplyAll(const RCP<Matrix>& left, const RCP<Matrix>& right, Teuchos::RCP<Vector >& InnerProdVec) const {
466 
467  TEUCHOS_TEST_FOR_EXCEPTION(!left->getDomainMap()->isSameAs(*right->getDomainMap()), Exceptions::RuntimeError, "MueLu::PgPFactory::MultiplyAll: domain maps of left and right do not match. Error.");
468  TEUCHOS_TEST_FOR_EXCEPTION(!left->getRowMap()->isSameAs(*right->getRowMap()), Exceptions::RuntimeError, "MueLu::PgPFactory::MultiplyAll: row maps of left and right do not match. Error.");
469 #if 1 // 1=new "fast code, 0=old "slow", but safe code
470 #if 0 // not necessary - remove me
471  if(InnerProdVec->getMap()->isSameAs(*left->getColMap())) {
472  // initialize NewRightLocal vector and assign all entries to
473  // left->getColMap()->getNodeNumElements() + 1
474  std::vector<LocalOrdinal> NewRightLocal(right->getColMap()->getNodeNumElements(), Teuchos::as<LocalOrdinal>(left->getColMap()->getNodeNumElements()+1));
475 
476  LocalOrdinal i = 0;
477  for (size_t j=0; j < right->getColMap()->getNodeNumElements(); j++) {
478  while ( (i < Teuchos::as<LocalOrdinal>(left->getColMap()->getNodeNumElements())) &&
479  (left->getColMap()->getGlobalElement(i) < right->getColMap()->getGlobalElement(j)) ) i++;
480  if (left->getColMap()->getGlobalElement(i) == right->getColMap()->getGlobalElement(j)) {
481  NewRightLocal[j] = i;
482  }
483  }
484 
485  Teuchos::ArrayRCP< Scalar > InnerProd_local = InnerProdVec->getDataNonConst(0);
486  std::vector<Scalar> temp_array(left->getColMap()->getNodeNumElements()+1, 0.0);
487 
488  for(size_t n=0; n<right->getNodeNumRows(); n++) {
489  Teuchos::ArrayView<const LocalOrdinal> lindices_left;
490  Teuchos::ArrayView<const Scalar> lvals_left;
491  Teuchos::ArrayView<const LocalOrdinal> lindices_right;
492  Teuchos::ArrayView<const Scalar> lvals_right;
493 
494  left->getLocalRowView (n, lindices_left, lvals_left);
495  right->getLocalRowView(n, lindices_right, lvals_right);
496 
497  for(size_t j=0; j<Teuchos::as<size_t>(lindices_right.size()); j++) {
498  temp_array[NewRightLocal[lindices_right[j] ] ] = lvals_right[j];
499  }
500  for (size_t j=0; j < Teuchos::as<size_t>(lindices_left.size()); j++) {
501  InnerProd_local[lindices_left[j]] += temp_array[lindices_left[j] ]*lvals_left[j];
502  }
503  for (size_t j=0; j < Teuchos::as<size_t>(lindices_right.size()); j++) {
504  temp_array[NewRightLocal[lindices_right[j] ] ] = 0.0;
505  }
506  }
507  // exporter: overlapping map to nonoverlapping map (target map is unique)
508  Teuchos::RCP<const Export> exporter =
509  ExportFactory::Build(left->getColMap(), left->getDomainMap()); // TODO: change left to right?
510 
511  Teuchos::RCP<Vector > nonoverlap =
512  VectorFactory::Build(left->getDomainMap()); // TODO: change left to right?
513 
514  nonoverlap->doExport(*InnerProdVec, *exporter, Xpetra::ADD);
515 
516  // importer: nonoverlapping map to overlapping map
517 
518  // importer: source -> target maps
519  Teuchos::RCP<const Import > importer =
520  ImportFactory::Build(left->getDomainMap(), left->getColMap()); // TODO: change left to right?
521 
522  // doImport target->doImport(*source, importer, action)
523  InnerProdVec->doImport(*nonoverlap, *importer, Xpetra::INSERT);
524 
525 
526  } else
527 #endif // end remove me
528  if(InnerProdVec->getMap()->isSameAs(*right->getColMap())) {
529  size_t szNewLeftLocal = TEUCHOS_MAX(left->getColMap()->getNodeNumElements(), right->getColMap()->getNodeNumElements());
530  Teuchos::RCP<std::vector<LocalOrdinal> > NewLeftLocal = Teuchos::rcp(new std::vector<LocalOrdinal>(szNewLeftLocal, Teuchos::as<LocalOrdinal>(right->getColMap()->getMaxLocalIndex()+1)));
531 
532  LocalOrdinal j = 0;
533  for (size_t i=0; i < left->getColMap()->getNodeNumElements(); i++) {
534  while ( (j < Teuchos::as<LocalOrdinal>(right->getColMap()->getNodeNumElements())) &&
535  (right->getColMap()->getGlobalElement(j) < left->getColMap()->getGlobalElement(i)) ) j++;
536  if (right->getColMap()->getGlobalElement(j) == left->getColMap()->getGlobalElement(i)) {
537  (*NewLeftLocal)[i] = j;
538  }
539  }
540 
541  /*for (size_t i=0; i < right->getColMap()->getNodeNumElements(); i++) {
542  std::cout << "left col map: " << (*NewLeftLocal)[i] << " GID: " << left->getColMap()->getGlobalElement((*NewLeftLocal)[i]) << " GID: " << right->getColMap()->getGlobalElement(i) << " right col map: " << i << std::endl;
543  }*/
544 
545  Teuchos::ArrayRCP< Scalar > InnerProd_local = InnerProdVec->getDataNonConst(0);
546  Teuchos::RCP<std::vector<Scalar> > temp_array = Teuchos::rcp(new std::vector<Scalar>(right->getColMap()->getMaxLocalIndex()+2, 0.0));
547 
548  for(size_t n=0; n<left->getNodeNumRows(); n++) {
549  Teuchos::ArrayView<const LocalOrdinal> lindices_left;
550  Teuchos::ArrayView<const Scalar> lvals_left;
551  Teuchos::ArrayView<const LocalOrdinal> lindices_right;
552  Teuchos::ArrayView<const Scalar> lvals_right;
553 
554  left->getLocalRowView (n, lindices_left, lvals_left);
555  right->getLocalRowView(n, lindices_right, lvals_right);
556 
557  for(size_t i=0; i<Teuchos::as<size_t>(lindices_left.size()); i++) {
558  (*temp_array)[(*NewLeftLocal)[lindices_left[i] ] ] = lvals_left[i];
559  }
560  for (size_t i=0; i < Teuchos::as<size_t>(lindices_right.size()); i++) {
561  InnerProd_local[lindices_right[i]] += (*temp_array)[lindices_right[i] ] * lvals_right[i];
562  }
563  for (size_t i=0; i < Teuchos::as<size_t>(lindices_left.size()); i++) {
564  (*temp_array)[(*NewLeftLocal)[lindices_left[i] ] ] = 0.0;
565  }
566  }
567  InnerProd_local = Teuchos::ArrayRCP< Scalar >();
568  // exporter: overlapping map to nonoverlapping map (target map is unique)
569  Teuchos::RCP<const Export> exporter =
570  ExportFactory::Build(right->getColMap(), right->getDomainMap()); // TODO: change left to right?
571 
572  Teuchos::RCP<Vector> nonoverlap =
573  VectorFactory::Build(right->getDomainMap()); // TODO: change left to right?
574 
575  nonoverlap->doExport(*InnerProdVec, *exporter, Xpetra::ADD);
576 
577  // importer: nonoverlapping map to overlapping map
578 
579  // importer: source -> target maps
580  Teuchos::RCP<const Import > importer =
581  ImportFactory::Build(right->getDomainMap(), right->getColMap()); // TODO: change left to right?
582  // doImport target->doImport(*source, importer, action)
583  InnerProdVec->doImport(*nonoverlap, *importer, Xpetra::INSERT);
584  } else {
585  TEUCHOS_TEST_FOR_EXCEPTION(true, Exceptions::RuntimeError, "MueLu::PgPFactory::MultiplyAll: map of InnerProdVec must be same as column map of left operator? Error.");
586  }
587 
588 #else // old "safe" code
589  if(InnerProdVec->getMap()->isSameAs(*left->getColMap())) {
590 
591  Teuchos::ArrayRCP< Scalar > InnerProd_local = InnerProdVec->getDataNonConst(0);
592 
593  // declare variables
594  Teuchos::ArrayView<const LocalOrdinal> lindices_left;
595  Teuchos::ArrayView<const Scalar> lvals_left;
596  Teuchos::ArrayView<const LocalOrdinal> lindices_right;
597  Teuchos::ArrayView<const Scalar> lvals_right;
598 
599  for(size_t n=0; n<left->getNodeNumRows(); n++)
600  {
601 
602  left->getLocalRowView (n, lindices_left, lvals_left);
603  right->getLocalRowView(n, lindices_right, lvals_right);
604 
605  for(size_t i=0; i<Teuchos::as<size_t>(lindices_left.size()); i++)
606  {
607  GlobalOrdinal left_gid = left->getColMap()->getGlobalElement(lindices_left[i]);
608  for(size_t j=0; j<Teuchos::as<size_t>(lindices_right.size()); j++)
609  {
610  GlobalOrdinal right_gid= right->getColMap()->getGlobalElement(lindices_right[j]);
611  if(left_gid == right_gid)
612  {
613  InnerProd_local[lindices_left[i]] += lvals_left[i]*lvals_right[j];
614  break; // skip remaining gids of right operator
615  }
616  }
617  }
618  }
619 
620  // exporter: overlapping map to nonoverlapping map (target map is unique)
621  Teuchos::RCP<const Export> exporter =
622  ExportFactory::Build(left->getColMap(), left->getDomainMap()); // TODO: change left to right?
623 
624  Teuchos::RCP<Vector > nonoverlap =
625  VectorFactory::Build(left->getDomainMap()); // TODO: change left to right?
626 
627  nonoverlap->doExport(*InnerProdVec, *exporter, Xpetra::ADD);
628 
629  // importer: nonoverlapping map to overlapping map
630 
631  // importer: source -> target maps
632  Teuchos::RCP<const Import > importer =
633  ImportFactory::Build(left->getDomainMap(), left->getColMap()); // TODO: change left to right?
634 
635  // doImport target->doImport(*source, importer, action)
636  InnerProdVec->doImport(*nonoverlap, *importer, Xpetra::INSERT);
637  }
638  else if(InnerProdVec->getMap()->isSameAs(*right->getColMap())) {
639  Teuchos::ArrayRCP< Scalar > InnerProd_local = InnerProdVec->getDataNonConst(0);
640 
641  Teuchos::ArrayView<const LocalOrdinal> lindices_left;
642  Teuchos::ArrayView<const Scalar> lvals_left;
643  Teuchos::ArrayView<const LocalOrdinal> lindices_right;
644  Teuchos::ArrayView<const Scalar> lvals_right;
645 
646  for(size_t n=0; n<left->getNodeNumRows(); n++)
647  {
648  left->getLocalRowView(n, lindices_left, lvals_left);
649  right->getLocalRowView(n, lindices_right, lvals_right);
650 
651  for(size_t i=0; i<Teuchos::as<size_t>(lindices_left.size()); i++)
652  {
653  GlobalOrdinal left_gid = left->getColMap()->getGlobalElement(lindices_left[i]);
654  for(size_t j=0; j<Teuchos::as<size_t>(lindices_right.size()); j++)
655  {
656  GlobalOrdinal right_gid= right->getColMap()->getGlobalElement(lindices_right[j]);
657  if(left_gid == right_gid)
658  {
659  InnerProd_local[lindices_right[j]] += lvals_left[i]*lvals_right[j];
660  break; // skip remaining gids of right operator
661  }
662  }
663  }
664  }
665 
666  // exporter: overlapping map to nonoverlapping map (target map is unique)
667  Teuchos::RCP<const Export> exporter =
668  ExportFactory::Build(right->getColMap(), right->getDomainMap()); // TODO: change left to right?
669 
670  Teuchos::RCP<Vector> nonoverlap =
671  VectorFactory::Build(right->getDomainMap()); // TODO: change left to right?
672 
673  nonoverlap->doExport(*InnerProdVec, *exporter, Xpetra::ADD);
674 
675  // importer: nonoverlapping map to overlapping map
676 
677  // importer: source -> target maps
678  Teuchos::RCP<const Import > importer =
679  ImportFactory::Build(right->getDomainMap(), right->getColMap()); // TODO: change left to right?
680 
681  // doImport target->doImport(*source, importer, action)
682  InnerProdVec->doImport(*nonoverlap, *importer, Xpetra::INSERT);
683  }
684  else {
685  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.");
686  }
687 #endif
688  }
689 
690  template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
691  void PgPFactory<Scalar, LocalOrdinal, GlobalOrdinal, Node>::BuildP(Level &/* fineLevel */, Level &/* coarseLevel */) const {
692  std::cout << "TODO: remove me" << std::endl;
693  }
694 
695  template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
697  SetParameter("ReUseRowBasedOmegas", ParameterEntry(bReuse));
698  }
699 
700 } //namespace MueLu
701 
702 #endif /* MUELU_PGPFACTORY_DEF_HPP */
static Teuchos::ArrayRCP< Scalar > GetMatrixDiagonal(const Xpetra::Matrix< Scalar, LocalOrdinal, GlobalOrdinal, Node > &A)
void DeclareInput(Level &fineLevel, Level &coarseLevel) const
Input.
#define MueLu_sumAll(rcpComm, in, out)
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
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->Get< RCP<Matrix> >("A", factory) if factory == NULL => use default factory.
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)
#define MueLu_maxAll(rcpComm, in, out)
Timer to be used in factories. Similar to Monitor but with additional timers.
Print more statistics.
Namespace for MueLu classes and methods.
#define MueLu_minAll(rcpComm, in, out)
Print all statistics.
Print even more statistics.
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 SetMinimizationMode(MinimizationNorm minnorm)
Set minimization mode (L2NORM for cheapest, ANORM more expensive, DINVANORM = default) ...
void Build(Level &fineLevel, Level &coarseLevel) const
Build method.
MueLu::DefaultScalar Scalar
MueLu::DefaultGlobalOrdinal GlobalOrdinal
Class that holds all level-specific information.
Definition: MueLu_Level.hpp:99
Timer to be used in factories. Similar to SubMonitor but adds a timer level by level.
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)
void MultiplySelfAll(const RCP< Matrix > &Op, Teuchos::RCP< Xpetra::Vector< Scalar, LocalOrdinal, GlobalOrdinal, Node > > &InnerProdVec) const
static std::string PrintMatrixInfo(const Matrix &A, const std::string &msgTag, RCP< const Teuchos::ParameterList > params=Teuchos::null)
RCP< const ParameterList > GetValidParameterList() const
Return a const parameter list of valid parameters that setParameterList() will accept.
MinimizationNorm GetMinimizationMode()
return minimization mode
void Set(const std::string &ename, const T &entry, const FactoryBase *factory=NoFactory::get())
void ReUseDampingParameters(bool bReuse)
void BuildP(Level &fineLevel, Level &coarseLevel) const
Abstract Build method.
const RCP< const FactoryManagerBase > GetFactoryManager()
returns the current factory manager
Definition: MueLu_Level.cpp:96
Exception throws to report errors in the internal logical of the program.
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