MueLu Version of the Day
Loading...
Searching...
No Matches
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
70namespace 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 // RfromPFactory used to indicate to TogglePFactory that a factory
218 // capable or producing R can be invoked later. TogglePFactory
219 // replaces dummy value with an index into it's array of prolongators
220 // pointing to the correct prolongator factory. This is later used by
221 // RfromP_Or_TransP to invoke the prolongatorfactory in RestrictionMode
222 int dummy = 7;
223 Set(coarseLevel,"RfromPfactory", dummy);
224
225 if (IsPrint(Statistics1))
226 GetOStream(Statistics1) << PerfUtils::PrintMatrixInfo(*P_smoothed, "P", params);
227
228 // NOTE: EXPERIMENTAL
229 if (Ptent->IsView("stridedMaps"))
230 P_smoothed->CreateView("stridedMaps", Ptent);
231
232 } else {
233 // prolongation factory is in restriction mode
234 RCP<Matrix> R = Utilities::Transpose(*P_smoothed, true);
235 Set(coarseLevel, "R", R);
236
237 if (IsPrint(Statistics1))
238 GetOStream(Statistics1) << PerfUtils::PrintMatrixInfo(*R, "P", params);
239
240 // NOTE: EXPERIMENTAL
241 if (Ptent->IsView("stridedMaps"))
242 R->CreateView("stridedMaps", Ptent, true);
243 }
244
245 }
246
247 template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
248 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 {
249 FactoryMonitor m(*this, "PgPFactory::ComputeRowBasedOmega", coarseLevel);
250
251 typedef typename Teuchos::ScalarTraits<Scalar>::magnitudeType Magnitude;
252 Scalar sZero = Teuchos::ScalarTraits<Scalar>::zero();
253 Magnitude mZero = Teuchos::ScalarTraits<Scalar>::magnitude(sZero);
254
255 Teuchos::RCP<Vector > Numerator = Teuchos::null;
256 Teuchos::RCP<Vector > Denominator = Teuchos::null;
257
258 const ParameterList & pL = GetParameterList();
259 MinimizationNorm min_norm = pL.get<MinimizationNorm>("Minimization norm");
260
261 switch (min_norm)
262 {
263 case ANORM: {
264 // MUEMAT mode (=paper)
265 // Minimize with respect to the (A)' A norm.
266 // Need to be smart here to avoid the construction of A' A
267 //
268 // diag( P0' (A' A) D^{-1} A P0)
269 // omega = ------------------------------------------
270 // diag( P0' A' D^{-1}' ( A' A) D^{-1} A P0)
271 //
272 // expensive, since we have to recalculate AP0 due to the lack of an explicit scaling routine for DinvAP0
273
274 // calculate A * P0
275 bool doFillComplete=true;
276 bool optimizeStorage=false;
277 RCP<Matrix> AP0 = Xpetra::MatrixMatrix<Scalar, LocalOrdinal, GlobalOrdinal, Node>::Multiply(*A, false, *P0, false, GetOStream(Statistics2), doFillComplete, optimizeStorage);
278
279 // compute A * D^{-1} * A * P0
280 RCP<Matrix> ADinvAP0 = Xpetra::MatrixMatrix<Scalar, LocalOrdinal, GlobalOrdinal, Node>::Multiply(*A, false, *DinvAP0, false, GetOStream(Statistics2), doFillComplete, optimizeStorage);
281
282 Numerator = VectorFactory::Build(ADinvAP0->getColMap(), true);
283 Denominator = VectorFactory::Build(ADinvAP0->getColMap(), true);
284 MultiplyAll(AP0, ADinvAP0, Numerator);
285 MultiplySelfAll(ADinvAP0, Denominator);
286 }
287 break;
288 case L2NORM: {
289
290 // ML mode 1 (cheapest)
291 // Minimize with respect to L2 norm
292 // diag( P0' D^{-1} A P0)
293 // omega = -----------------------------
294 // diag( P0' A' D^{-1}' D^{-1} A P0)
295 //
296 Numerator = VectorFactory::Build(DinvAP0->getColMap(), true);
297 Denominator = VectorFactory::Build(DinvAP0->getColMap(), true);
298 MultiplyAll(P0, DinvAP0, Numerator);
299 MultiplySelfAll(DinvAP0, Denominator);
300 }
301 break;
302 case DINVANORM: {
303 // ML mode 2
304 // Minimize with respect to the (D^{-1} A)' D^{-1} A norm.
305 // Need to be smart here to avoid the construction of A' A
306 //
307 // diag( P0' ( A' D^{-1}' D^{-1} A) D^{-1} A P0)
308 // omega = --------------------------------------------------------
309 // diag( P0' A' D^{-1}' ( A' D^{-1}' D^{-1} A) D^{-1} A P0)
310 //
311
312 // compute D^{-1} * A * D^{-1} * A * P0
313 bool doFillComplete=true;
314 bool optimizeStorage=true;
315 Teuchos::ArrayRCP<Scalar> diagA = Utilities::GetMatrixDiagonal(*A);
316 RCP<Matrix> DinvADinvAP0 = Xpetra::MatrixMatrix<Scalar, LocalOrdinal, GlobalOrdinal, Node>::Multiply(*A, false, *DinvAP0, false, GetOStream(Statistics2), doFillComplete, optimizeStorage);
317 Utilities::MyOldScaleMatrix(*DinvADinvAP0, diagA, true, doFillComplete, optimizeStorage); //scale matrix with reciprocal of diag
318 diagA = Teuchos::ArrayRCP<Scalar>();
319
320 Numerator = VectorFactory::Build(DinvADinvAP0->getColMap(), true);
321 Denominator = VectorFactory::Build(DinvADinvAP0->getColMap(), true);
322 MultiplyAll(DinvAP0, DinvADinvAP0, Numerator);
323 MultiplySelfAll(DinvADinvAP0, Denominator);
324 }
325 break;
326 default:
327 TEUCHOS_TEST_FOR_EXCEPTION(true, Exceptions::RuntimeError, "MueLu::PgPFactory::Build: minimization mode not supported. error");
328 }
329
330
332 Teuchos::RCP<Vector > ColBasedOmega =
333 VectorFactory::Build(Numerator->getMap()/*DinvAP0->getColMap()*/, true);
334
335 ColBasedOmega->putScalar(-666);
336
337 Teuchos::ArrayRCP< const Scalar > Numerator_local = Numerator->getData(0);
338 Teuchos::ArrayRCP< const Scalar > Denominator_local = Denominator->getData(0);
339 Teuchos::ArrayRCP< Scalar > ColBasedOmega_local = ColBasedOmega->getDataNonConst(0);
340 GlobalOrdinal zero_local = 0; // count negative colbased omegas
341 GlobalOrdinal nan_local = 0; // count NaNs -> set them to zero
342 Magnitude min_local = 1000000.0; //Teuchos::ScalarTraits<Scalar>::one() * (Scalar) 1000000;
343 Magnitude max_local = 0.0;
344 for(LocalOrdinal i = 0; i < Teuchos::as<LocalOrdinal>(Numerator->getLocalLength()); i++) {
345 if(Teuchos::ScalarTraits<Scalar>::magnitude(Denominator_local[i]) == mZero)
346 {
347 ColBasedOmega_local[i] = 0.0; // fallback: nonsmoothed basis function since denominator == 0.0
348 nan_local++;
349 }
350 else
351 {
352 ColBasedOmega_local[i] = Numerator_local[i] / Denominator_local[i]; // default case
353 }
354
355 if(Teuchos::ScalarTraits<Scalar>::magnitude(ColBasedOmega_local[i]) < mZero) { // negative omegas are not valid. set them to zero
356 ColBasedOmega_local[i] = Teuchos::ScalarTraits<Scalar>::zero();
357 zero_local++; // count zero omegas
358 }
359
360 // handle case that Nominator == Denominator -> Dirichlet bcs in A?
361 // fallback if ColBasedOmega == 1 -> very strong smoothing may lead to zero rows in P
362 // TAW: this is somewhat nonstandard and a rough fallback strategy to avoid problems
363 // also avoid "overshooting" with omega > 0.8
364 if(Teuchos::ScalarTraits<Scalar>::magnitude(ColBasedOmega_local[i]) >= 0.8) {
365 ColBasedOmega_local[i] = 0.0;
366 }
367
368 if(Teuchos::ScalarTraits<Scalar>::magnitude(ColBasedOmega_local[i]) < min_local)
369 { min_local = Teuchos::ScalarTraits<Scalar>::magnitude(ColBasedOmega_local[i]); }
370 if(Teuchos::ScalarTraits<Scalar>::magnitude(ColBasedOmega_local[i]) > max_local)
371 { max_local = Teuchos::ScalarTraits<Scalar>::magnitude(ColBasedOmega_local[i]); }
372 }
373
374 { // be verbose
375 GlobalOrdinal zero_all;
376 GlobalOrdinal nan_all;
377 Magnitude min_all;
378 Magnitude max_all;
379 MueLu_sumAll(A->getRowMap()->getComm(), zero_local, zero_all);
380 MueLu_sumAll(A->getRowMap()->getComm(), nan_local, nan_all);
381 MueLu_minAll(A->getRowMap()->getComm(), min_local, min_all);
382 MueLu_maxAll(A->getRowMap()->getComm(), max_local, max_all);
383
384 GetOStream(MueLu::Statistics1, 0) << "PgPFactory: smoothed aggregation (scheme: ";
385 switch (min_norm) {
386 case ANORM: GetOStream(Statistics1) << "Anorm)" << std::endl; break;
387 case L2NORM: GetOStream(Statistics1) << "L2norm)" << std::endl; break;
388 case DINVANORM: GetOStream(Statistics1) << "DinvAnorm)" << std::endl; break;
389 default: GetOStream(Statistics1) << "unknown)" << std::endl; break;
390 }
391 GetOStream(Statistics1) << "Damping parameter: min = " << min_all << ", max = " << max_all << std::endl;
392 GetOStream(Statistics) << "# negative omegas: " << zero_all << " out of " << ColBasedOmega->getGlobalLength() << " column-based omegas" << std::endl;
393 GetOStream(Statistics) << "# NaNs: " << nan_all << " out of " << ColBasedOmega->getGlobalLength() << " column-based omegas" << std::endl;
394 }
395
396 if(coarseLevel.IsRequested("ColBasedOmega", this)) {
397 coarseLevel.Set("ColBasedOmega", ColBasedOmega, this);
398 }
399
401 // transform column based omegas to row based omegas
402 RowBasedOmega =
403 VectorFactory::Build(DinvAP0->getRowMap(), true);
404
405 RowBasedOmega->putScalar(-666); // TODO bad programming style
406
407 bool bAtLeastOneDefined = false;
408 Teuchos::ArrayRCP< Scalar > RowBasedOmega_local = RowBasedOmega->getDataNonConst(0);
409 for(LocalOrdinal row = 0; row<Teuchos::as<LocalOrdinal>(A->getLocalNumRows()); row++) {
410 Teuchos::ArrayView<const LocalOrdinal> lindices;
411 Teuchos::ArrayView<const Scalar> lvals;
412 DinvAP0->getLocalRowView(row, lindices, lvals);
413 bAtLeastOneDefined = false;
414 for(size_t j=0; j<Teuchos::as<size_t>(lindices.size()); j++) {
415 Scalar omega = ColBasedOmega_local[lindices[j]];
416 if (Teuchos::ScalarTraits<Scalar>::magnitude(omega) != -666) { // TODO bad programming style
417 bAtLeastOneDefined = true;
418 if(Teuchos::ScalarTraits<Scalar>::magnitude(RowBasedOmega_local[row]) == -666) RowBasedOmega_local[row] = omega;
419 else if(Teuchos::ScalarTraits<Scalar>::magnitude(omega) < Teuchos::ScalarTraits<Scalar>::magnitude(RowBasedOmega_local[row])) RowBasedOmega_local[row] = omega;
420 }
421 }
422 if(bAtLeastOneDefined == true) {
423 if(Teuchos::ScalarTraits<Scalar>::magnitude(RowBasedOmega_local[row]) < mZero)
424 RowBasedOmega_local[row] = sZero;
425 }
426 }
427
428 if(coarseLevel.IsRequested("RowBasedOmega", this)) {
429 Set(coarseLevel, "RowBasedOmega", RowBasedOmega);
430 }
431 }
432
433 template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
434 void PgPFactory<Scalar, LocalOrdinal, GlobalOrdinal, Node>::MultiplySelfAll(const RCP<Matrix>& Op, Teuchos::RCP<Vector >& InnerProdVec) const {
435
436 // note: InnerProdVec is based on column map of Op
437 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");
438
439 Teuchos::ArrayRCP< Scalar > InnerProd_local = InnerProdVec->getDataNonConst(0);
440
441 Teuchos::ArrayView<const LocalOrdinal> lindices;
442 Teuchos::ArrayView<const Scalar> lvals;
443
444 for(size_t n=0; n<Op->getLocalNumRows(); n++) {
445 Op->getLocalRowView(n, lindices, lvals);
446 for(size_t i=0; i<Teuchos::as<size_t>(lindices.size()); i++) {
447 InnerProd_local[lindices[i]] += lvals[i]*lvals[i];
448 }
449 }
450 InnerProd_local = Teuchos::ArrayRCP< Scalar >();
451
452 // exporter: overlapping map to nonoverlapping map (target map is unique)
453 Teuchos::RCP<const Export> exporter =
454 ExportFactory::Build(Op->getColMap(), Op->getDomainMap());
455
456 Teuchos::RCP<Vector > nonoverlap =
457 VectorFactory::Build(Op->getDomainMap());
458
459 nonoverlap->doExport(*InnerProdVec, *exporter, Xpetra::ADD);
460
461 // importer: nonoverlapping map to overlapping map
462
463 // importer: source -> target maps
464 Teuchos::RCP<const Import > importer =
465 ImportFactory::Build(Op->getDomainMap(), Op->getColMap());
466
467 // doImport target->doImport(*source, importer, action)
468 InnerProdVec->doImport(*nonoverlap, *importer, Xpetra::INSERT);
469
470 }
471
472 template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
473 void PgPFactory<Scalar, LocalOrdinal, GlobalOrdinal, Node>::MultiplyAll(const RCP<Matrix>& left, const RCP<Matrix>& right, Teuchos::RCP<Vector >& InnerProdVec) const {
474
475 TEUCHOS_TEST_FOR_EXCEPTION(!left->getDomainMap()->isSameAs(*right->getDomainMap()), Exceptions::RuntimeError, "MueLu::PgPFactory::MultiplyAll: domain maps of left and right do not match. Error.");
476 TEUCHOS_TEST_FOR_EXCEPTION(!left->getRowMap()->isSameAs(*right->getRowMap()), Exceptions::RuntimeError, "MueLu::PgPFactory::MultiplyAll: row maps of left and right do not match. Error.");
477#if 1 // 1=new "fast code, 0=old "slow", but safe code
478#if 0 // not necessary - remove me
479 if(InnerProdVec->getMap()->isSameAs(*left->getColMap())) {
480 // initialize NewRightLocal vector and assign all entries to
481 // left->getColMap()->getLocalNumElements() + 1
482 std::vector<LocalOrdinal> NewRightLocal(right->getColMap()->getLocalNumElements(), Teuchos::as<LocalOrdinal>(left->getColMap()->getLocalNumElements()+1));
483
484 LocalOrdinal i = 0;
485 for (size_t j=0; j < right->getColMap()->getLocalNumElements(); j++) {
486 while ( (i < Teuchos::as<LocalOrdinal>(left->getColMap()->getLocalNumElements())) &&
487 (left->getColMap()->getGlobalElement(i) < right->getColMap()->getGlobalElement(j)) ) i++;
488 if (left->getColMap()->getGlobalElement(i) == right->getColMap()->getGlobalElement(j)) {
489 NewRightLocal[j] = i;
490 }
491 }
492
493 Teuchos::ArrayRCP< Scalar > InnerProd_local = InnerProdVec->getDataNonConst(0);
494 std::vector<Scalar> temp_array(left->getColMap()->getLocalNumElements()+1, 0.0);
495
496 for(size_t n=0; n<right->getLocalNumRows(); n++) {
497 Teuchos::ArrayView<const LocalOrdinal> lindices_left;
498 Teuchos::ArrayView<const Scalar> lvals_left;
499 Teuchos::ArrayView<const LocalOrdinal> lindices_right;
500 Teuchos::ArrayView<const Scalar> lvals_right;
501
502 left->getLocalRowView (n, lindices_left, lvals_left);
503 right->getLocalRowView(n, lindices_right, lvals_right);
504
505 for(size_t j=0; j<Teuchos::as<size_t>(lindices_right.size()); j++) {
506 temp_array[NewRightLocal[lindices_right[j] ] ] = lvals_right[j];
507 }
508 for (size_t j=0; j < Teuchos::as<size_t>(lindices_left.size()); j++) {
509 InnerProd_local[lindices_left[j]] += temp_array[lindices_left[j] ]*lvals_left[j];
510 }
511 for (size_t j=0; j < Teuchos::as<size_t>(lindices_right.size()); j++) {
512 temp_array[NewRightLocal[lindices_right[j] ] ] = 0.0;
513 }
514 }
515 // exporter: overlapping map to nonoverlapping map (target map is unique)
516 Teuchos::RCP<const Export> exporter =
517 ExportFactory::Build(left->getColMap(), left->getDomainMap()); // TODO: change left to right?
518
519 Teuchos::RCP<Vector > nonoverlap =
520 VectorFactory::Build(left->getDomainMap()); // TODO: change left to right?
521
522 nonoverlap->doExport(*InnerProdVec, *exporter, Xpetra::ADD);
523
524 // importer: nonoverlapping map to overlapping map
525
526 // importer: source -> target maps
527 Teuchos::RCP<const Import > importer =
528 ImportFactory::Build(left->getDomainMap(), left->getColMap()); // TODO: change left to right?
529
530 // doImport target->doImport(*source, importer, action)
531 InnerProdVec->doImport(*nonoverlap, *importer, Xpetra::INSERT);
532
533
534 } else
535#endif // end remove me
536 if(InnerProdVec->getMap()->isSameAs(*right->getColMap())) {
537 size_t szNewLeftLocal = TEUCHOS_MAX(left->getColMap()->getLocalNumElements(), right->getColMap()->getLocalNumElements());
538 Teuchos::RCP<std::vector<LocalOrdinal> > NewLeftLocal = Teuchos::rcp(new std::vector<LocalOrdinal>(szNewLeftLocal, Teuchos::as<LocalOrdinal>(right->getColMap()->getMaxLocalIndex()+1)));
539
540 LocalOrdinal j = 0;
541 for (size_t i=0; i < left->getColMap()->getLocalNumElements(); i++) {
542 while ( (j < Teuchos::as<LocalOrdinal>(right->getColMap()->getLocalNumElements())) &&
543 (right->getColMap()->getGlobalElement(j) < left->getColMap()->getGlobalElement(i)) ) j++;
544 if (right->getColMap()->getGlobalElement(j) == left->getColMap()->getGlobalElement(i)) {
545 (*NewLeftLocal)[i] = j;
546 }
547 }
548
549 /*for (size_t i=0; i < right->getColMap()->getLocalNumElements(); i++) {
550 std::cout << "left col map: " << (*NewLeftLocal)[i] << " GID: " << left->getColMap()->getGlobalElement((*NewLeftLocal)[i]) << " GID: " << right->getColMap()->getGlobalElement(i) << " right col map: " << i << std::endl;
551 }*/
552
553 Teuchos::ArrayRCP< Scalar > InnerProd_local = InnerProdVec->getDataNonConst(0);
554 Teuchos::RCP<std::vector<Scalar> > temp_array = Teuchos::rcp(new std::vector<Scalar>(right->getColMap()->getMaxLocalIndex()+2, 0.0));
555
556 for(size_t n=0; n<left->getLocalNumRows(); n++) {
557 Teuchos::ArrayView<const LocalOrdinal> lindices_left;
558 Teuchos::ArrayView<const Scalar> lvals_left;
559 Teuchos::ArrayView<const LocalOrdinal> lindices_right;
560 Teuchos::ArrayView<const Scalar> lvals_right;
561
562 left->getLocalRowView (n, lindices_left, lvals_left);
563 right->getLocalRowView(n, lindices_right, lvals_right);
564
565 for(size_t i=0; i<Teuchos::as<size_t>(lindices_left.size()); i++) {
566 (*temp_array)[(*NewLeftLocal)[lindices_left[i] ] ] = lvals_left[i];
567 }
568 for (size_t i=0; i < Teuchos::as<size_t>(lindices_right.size()); i++) {
569 InnerProd_local[lindices_right[i]] += (*temp_array)[lindices_right[i] ] * lvals_right[i];
570 }
571 for (size_t i=0; i < Teuchos::as<size_t>(lindices_left.size()); i++) {
572 (*temp_array)[(*NewLeftLocal)[lindices_left[i] ] ] = 0.0;
573 }
574 }
575 InnerProd_local = Teuchos::ArrayRCP< Scalar >();
576 // exporter: overlapping map to nonoverlapping map (target map is unique)
577 Teuchos::RCP<const Export> exporter =
578 ExportFactory::Build(right->getColMap(), right->getDomainMap()); // TODO: change left to right?
579
580 Teuchos::RCP<Vector> nonoverlap =
581 VectorFactory::Build(right->getDomainMap()); // TODO: change left to right?
582
583 nonoverlap->doExport(*InnerProdVec, *exporter, Xpetra::ADD);
584
585 // importer: nonoverlapping map to overlapping map
586
587 // importer: source -> target maps
588 Teuchos::RCP<const Import > importer =
589 ImportFactory::Build(right->getDomainMap(), right->getColMap()); // TODO: change left to right?
590 // doImport target->doImport(*source, importer, action)
591 InnerProdVec->doImport(*nonoverlap, *importer, Xpetra::INSERT);
592 } else {
593 TEUCHOS_TEST_FOR_EXCEPTION(true, Exceptions::RuntimeError, "MueLu::PgPFactory::MultiplyAll: map of InnerProdVec must be same as column map of left operator? Error.");
594 }
595
596#else // old "safe" code
597 if(InnerProdVec->getMap()->isSameAs(*left->getColMap())) {
598
599 Teuchos::ArrayRCP< Scalar > InnerProd_local = InnerProdVec->getDataNonConst(0);
600
601 // declare variables
602 Teuchos::ArrayView<const LocalOrdinal> lindices_left;
603 Teuchos::ArrayView<const Scalar> lvals_left;
604 Teuchos::ArrayView<const LocalOrdinal> lindices_right;
605 Teuchos::ArrayView<const Scalar> lvals_right;
606
607 for(size_t n=0; n<left->getLocalNumRows(); n++)
608 {
609
610 left->getLocalRowView (n, lindices_left, lvals_left);
611 right->getLocalRowView(n, lindices_right, lvals_right);
612
613 for(size_t i=0; i<Teuchos::as<size_t>(lindices_left.size()); i++)
614 {
615 GlobalOrdinal left_gid = left->getColMap()->getGlobalElement(lindices_left[i]);
616 for(size_t j=0; j<Teuchos::as<size_t>(lindices_right.size()); j++)
617 {
618 GlobalOrdinal right_gid= right->getColMap()->getGlobalElement(lindices_right[j]);
619 if(left_gid == right_gid)
620 {
621 InnerProd_local[lindices_left[i]] += lvals_left[i]*lvals_right[j];
622 break; // skip remaining gids of right operator
623 }
624 }
625 }
626 }
627
628 // exporter: overlapping map to nonoverlapping map (target map is unique)
629 Teuchos::RCP<const Export> exporter =
630 ExportFactory::Build(left->getColMap(), left->getDomainMap()); // TODO: change left to right?
631
632 Teuchos::RCP<Vector > nonoverlap =
633 VectorFactory::Build(left->getDomainMap()); // TODO: change left to right?
634
635 nonoverlap->doExport(*InnerProdVec, *exporter, Xpetra::ADD);
636
637 // importer: nonoverlapping map to overlapping map
638
639 // importer: source -> target maps
640 Teuchos::RCP<const Import > importer =
641 ImportFactory::Build(left->getDomainMap(), left->getColMap()); // TODO: change left to right?
642
643 // doImport target->doImport(*source, importer, action)
644 InnerProdVec->doImport(*nonoverlap, *importer, Xpetra::INSERT);
645 }
646 else if(InnerProdVec->getMap()->isSameAs(*right->getColMap())) {
647 Teuchos::ArrayRCP< Scalar > InnerProd_local = InnerProdVec->getDataNonConst(0);
648
649 Teuchos::ArrayView<const LocalOrdinal> lindices_left;
650 Teuchos::ArrayView<const Scalar> lvals_left;
651 Teuchos::ArrayView<const LocalOrdinal> lindices_right;
652 Teuchos::ArrayView<const Scalar> lvals_right;
653
654 for(size_t n=0; n<left->getLocalNumRows(); n++)
655 {
656 left->getLocalRowView(n, lindices_left, lvals_left);
657 right->getLocalRowView(n, lindices_right, lvals_right);
658
659 for(size_t i=0; i<Teuchos::as<size_t>(lindices_left.size()); i++)
660 {
661 GlobalOrdinal left_gid = left->getColMap()->getGlobalElement(lindices_left[i]);
662 for(size_t j=0; j<Teuchos::as<size_t>(lindices_right.size()); j++)
663 {
664 GlobalOrdinal right_gid= right->getColMap()->getGlobalElement(lindices_right[j]);
665 if(left_gid == right_gid)
666 {
667 InnerProd_local[lindices_right[j]] += lvals_left[i]*lvals_right[j];
668 break; // skip remaining gids of right operator
669 }
670 }
671 }
672 }
673
674 // exporter: overlapping map to nonoverlapping map (target map is unique)
675 Teuchos::RCP<const Export> exporter =
676 ExportFactory::Build(right->getColMap(), right->getDomainMap()); // TODO: change left to right?
677
678 Teuchos::RCP<Vector> nonoverlap =
679 VectorFactory::Build(right->getDomainMap()); // TODO: change left to right?
680
681 nonoverlap->doExport(*InnerProdVec, *exporter, Xpetra::ADD);
682
683 // importer: nonoverlapping map to overlapping map
684
685 // importer: source -> target maps
686 Teuchos::RCP<const Import > importer =
687 ImportFactory::Build(right->getDomainMap(), right->getColMap()); // TODO: change left to right?
688
689 // doImport target->doImport(*source, importer, action)
690 InnerProdVec->doImport(*nonoverlap, *importer, Xpetra::INSERT);
691 }
692 else {
693 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.");
694 }
695#endif
696 }
697
698 template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
699 void PgPFactory<Scalar, LocalOrdinal, GlobalOrdinal, Node>::BuildP(Level &/* fineLevel */, Level &/* coarseLevel */) const {
700 std::cout << "TODO: remove me" << std::endl;
701 }
702
703 template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
705 SetParameter("ReUseRowBasedOmegas", ParameterEntry(bReuse));
706 }
707
708} //namespace MueLu
709
710#endif /* MUELU_PGPFACTORY_DEF_HPP */
#define MueLu_maxAll(rcpComm, in, out)
#define MueLu_sumAll(rcpComm, in, out)
#define MueLu_minAll(rcpComm, in, out)
MueLu::DefaultLocalOrdinal LocalOrdinal
MueLu::DefaultScalar Scalar
MueLu::DefaultGlobalOrdinal GlobalOrdinal
Exception throws to report errors in the internal logical of the program.
Timer to be used in factories. Similar to Monitor but with additional timers.
Class that holds all level-specific information.
Definition: MueLu_Level.hpp:99
void DeclareInput(const std::string &ename, const FactoryBase *factory, const FactoryBase *requestedBy=NoFactory::get())
Callback from FactoryBase::CallDeclareInput() and FactoryBase::DeclareInput()
const RCP< const FactoryManagerBase > GetFactoryManager()
returns the current factory manager
Definition: MueLu_Level.cpp:96
T & Get(const std::string &ename, const FactoryBase *factory=NoFactory::get())
Get data without decrementing associated storage counter (i.e., read-only access)....
void Set(const std::string &ename, const T &entry, const FactoryBase *factory=NoFactory::get())
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...
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.
void Build(Level &fineLevel, Level &coarseLevel) const
Build method.
void ReUseDampingParameters(bool bReuse)
void MultiplyAll(const RCP< Matrix > &left, const RCP< Matrix > &right, Teuchos::RCP< Xpetra::Vector< Scalar, LocalOrdinal, GlobalOrdinal, Node > > &InnerProdVec) const
void DeclareInput(Level &fineLevel, Level &coarseLevel) const
Input.
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
void SetMinimizationMode(MinimizationNorm minnorm)
Set minimization mode (L2NORM for cheapest, ANORM more expensive, DINVANORM = default)
MinimizationNorm GetMinimizationMode()
return minimization mode
void BuildP(Level &fineLevel, Level &coarseLevel) const
Abstract Build method.
void MultiplySelfAll(const RCP< Matrix > &Op, Teuchos::RCP< Xpetra::Vector< Scalar, LocalOrdinal, GlobalOrdinal, Node > > &InnerProdVec) const
Timer to be used in factories. Similar to SubMonitor but adds a timer level by level.
static Teuchos::ArrayRCP< Scalar > GetMatrixDiagonal(const Xpetra::Matrix< Scalar, LocalOrdinal, GlobalOrdinal, Node > &A)
static RCP< Xpetra::Matrix< Scalar, LocalOrdinal, GlobalOrdinal, Node > > Transpose(Xpetra::Matrix< Scalar, LocalOrdinal, GlobalOrdinal, Node > &Op, bool optimizeTranspose=false, const std::string &label=std::string(), const Teuchos::RCP< Teuchos::ParameterList > &params=Teuchos::null)
static void MyOldScaleMatrix(Xpetra::Matrix< Scalar, LocalOrdinal, GlobalOrdinal, Node > &Op, const Teuchos::ArrayRCP< const Scalar > &scalingVector, bool doInverse=true, bool doFillComplete=true, bool doOptimizeStorage=true)
Namespace for MueLu classes and methods.
@ Statistics2
Print even more statistics.
@ Statistics1
Print more statistics.
@ Statistics
Print all statistics.