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
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
70namespace 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=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>
696 void PgPFactory<Scalar, LocalOrdinal, GlobalOrdinal, Node>::ReUseDampingParameters(bool bReuse) {
697 SetParameter("ReUseRowBasedOmegas", ParameterEntry(bReuse));
698 }
699
700} //namespace MueLu
701
702#endif /* MUELU_PGPFACTORY_DEF_HPP */
MueLu_FactoryManagerBase.hpp
MueLu_Monitor.hpp
MueLu_maxAll
#define MueLu_maxAll(rcpComm, in, out)
Definition: MueLu_PerfUtils_decl.hpp:65
MueLu_sumAll
#define MueLu_sumAll(rcpComm, in, out)
Definition: MueLu_PerfUtils_decl.hpp:61
MueLu_minAll
#define MueLu_minAll(rcpComm, in, out)
Definition: MueLu_PerfUtils_decl.hpp:63
MueLu_PgPFactory_decl.hpp
MueLu_SingleLevelFactoryBase.hpp
LocalOrdinal
MueLu::DefaultLocalOrdinal LocalOrdinal
Definition: MueLu_UseDefaultTypes.hpp:50
Scalar
MueLu::DefaultScalar Scalar
Definition: MueLu_UseDefaultTypes.hpp:49
GlobalOrdinal
MueLu::DefaultGlobalOrdinal GlobalOrdinal
Definition: MueLu_UseDefaultTypes.hpp:51
MueLu::Exceptions::RuntimeError
Exception throws to report errors in the internal logical of the program.
Definition: MueLu_Exceptions.hpp:70
MueLu::FactoryMonitor
Timer to be used in factories. Similar to Monitor but with additional timers.
Definition: MueLu_Monitor.hpp:209
MueLu::Level
Class that holds all level-specific information.
Definition: MueLu_Level.hpp:99
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
MueLu::Level::GetFactoryManager
const RCP< const FactoryManagerBase > GetFactoryManager()
returns the current factory manager
Definition: MueLu_Level.cpp:96
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)....
Definition: MueLu_Level.hpp:192
MueLu::Level::Set
void Set(const std::string &ename, const T &entry, const FactoryBase *factory=NoFactory::get())
Definition: MueLu_Level.hpp:157
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's value exist...
Definition: MueLu_Level.hpp:345
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:127
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
MueLu::PgPFactory::Build
void Build(Level &fineLevel, Level &coarseLevel) const
Build method.
Definition: MueLu_PgPFactory_def.hpp:128
MueLu::PgPFactory::ReUseDampingParameters
void ReUseDampingParameters(bool bReuse)
Definition: MueLu_PgPFactory_def.hpp:696
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:465
MueLu::PgPFactory::DeclareInput
void DeclareInput(Level &fineLevel, Level &coarseLevel) const
Input.
Definition: MueLu_PgPFactory_def.hpp:96
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::PgPFactory::SetMinimizationMode
void SetMinimizationMode(MinimizationNorm minnorm)
Set minimization mode (L2NORM for cheapest, ANORM more expensive, DINVANORM = default)
Definition: MueLu_PgPFactory_def.hpp:85
MueLu::PgPFactory::GetMinimizationMode
MinimizationNorm GetMinimizationMode()
return minimization mode
Definition: MueLu_PgPFactory_def.hpp:90
MueLu::PgPFactory::BuildP
void BuildP(Level &fineLevel, Level &coarseLevel) const
Abstract Build method.
Definition: MueLu_PgPFactory_def.hpp:691
MueLu::PgPFactory::MultiplySelfAll
void MultiplySelfAll(const RCP< Matrix > &Op, Teuchos::RCP< Xpetra::Vector< Scalar, LocalOrdinal, GlobalOrdinal, Node > > &InnerProdVec) const
Definition: MueLu_PgPFactory_def.hpp:426
MueLu::SubFactoryMonitor
Timer to be used in factories. Similar to SubMonitor but adds a timer level by level.
Definition: MueLu_Monitor.hpp:261
MueLu::Utilities::GetMatrixDiagonal
static Teuchos::ArrayRCP< Scalar > GetMatrixDiagonal(const Xpetra::Matrix< Scalar, LocalOrdinal, GlobalOrdinal, Node > &A)
Definition: MueLu_Utilities_decl.hpp:209
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:466
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
Namespace for MueLu classes and methods.
Definition: MueLu_BrickAggregationFactory_decl.hpp:76
MueLu::Statistics2
@ Statistics2
Print even more statistics.
Definition: MueLu_VerbosityLevel.hpp:72
MueLu::Statistics1
@ Statistics1
Print more statistics.
Definition: MueLu_VerbosityLevel.hpp:71
MueLu::Statistics
@ Statistics
Print all statistics.
Definition: MueLu_VerbosityLevel.hpp:92