MueLu Version of the Day
Loading...
Searching...
No Matches
MueLu_ReitzingerPFactory_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_REITZINGERPFACTORY_DEF_HPP
47#define MUELU_REITZINGERPFACTORY_DEF_HPP
48
49#include <Xpetra_MapFactory.hpp>
50#include <Xpetra_Map.hpp>
51#include <Xpetra_CrsMatrix.hpp>
52#include <Xpetra_Matrix.hpp>
53#include <Xpetra_MatrixMatrix.hpp>
54#include <Xpetra_MultiVector.hpp>
55#include <Xpetra_MultiVectorFactory.hpp>
56#include <Xpetra_VectorFactory.hpp>
57#include <Xpetra_Import.hpp>
58#include <Xpetra_ImportUtils.hpp>
59#include <Xpetra_ImportFactory.hpp>
60#include <Xpetra_CrsMatrixWrap.hpp>
61#include <Xpetra_StridedMap.hpp>
62#include <Xpetra_StridedMapFactory.hpp>
63#include <Xpetra_IO.hpp>
64
66
67#include "MueLu_Aggregates.hpp"
68#include "MueLu_AmalgamationFactory.hpp"
69#include "MueLu_AmalgamationInfo.hpp"
70#include "MueLu_CoarseMapFactory.hpp"
71#include "MueLu_MasterList.hpp"
72#include "MueLu_Monitor.hpp"
73#include "MueLu_NullspaceFactory.hpp"
74#include "MueLu_PerfUtils.hpp"
75#include "MueLu_Utilities.hpp"
76
77namespace MueLu {
78
79 template <class Scalar,class LocalOrdinal, class GlobalOrdinal, class Node>
81 RCP<ParameterList> validParamList = rcp(new ParameterList());
82
83#define SET_VALID_ENTRY(name) validParamList->setEntry(name, MasterList::getEntry(name))
84 SET_VALID_ENTRY("repartition: enable");
85 SET_VALID_ENTRY("repartition: use subcommunicators");
86 SET_VALID_ENTRY("tentative: calculate qr");
87 SET_VALID_ENTRY("tentative: constant column sums");
88#undef SET_VALID_ENTRY
89
90 validParamList->set< RCP<const FactoryBase> >("A", Teuchos::null, "Generating factory of the matrix A");
91 validParamList->set< RCP<const FactoryBase> >("D0", Teuchos::null, "Generating factory of the matrix D0");
92 validParamList->set< RCP<const FactoryBase> >("NodeAggMatrix", Teuchos::null, "Generating factory of the matrix NodeAggMatrix");
93 validParamList->set< RCP<const FactoryBase> >("Pnodal", Teuchos::null, "Generating factory of the matrix P");
94 validParamList->set< RCP<const FactoryBase> >("NodeImporter", Teuchos::null, "Generating factory of the matrix NodeImporter");
95
96 // Make sure we don't recursively validate options for the matrixmatrix kernels
97 ParameterList norecurse;
98 norecurse.disableRecursiveValidation();
99 validParamList->set<ParameterList> ("matrixmatrix: kernel params", norecurse, "MatrixMatrix kernel parameters");
100
101 return validParamList;
102 }
103
104 template <class Scalar,class LocalOrdinal, class GlobalOrdinal, class Node>
106 Input(fineLevel, "A");
107 Input(fineLevel, "D0");
108 Input(fineLevel, "NodeAggMatrix");
109 Input(coarseLevel, "NodeAggMatrix");
110 Input(coarseLevel, "Pnodal");
111 // Input(coarseLevel, "NodeImporter");
112
113 }
114
115 template <class Scalar,class LocalOrdinal, class GlobalOrdinal, class Node>
117 return BuildP(fineLevel, coarseLevel);
118 }
119
120 template <class Scalar,class LocalOrdinal, class GlobalOrdinal, class Node>
122 FactoryMonitor m(*this, "Build", coarseLevel);
123 using Teuchos::arcp_const_cast;
124 using MT = typename Teuchos::ScalarTraits<SC>::magnitudeType;
125 using XMM = Xpetra::MatrixMatrix<SC,LO,GO,NO>;
126 Teuchos::FancyOStream& out0=GetBlackHole();
127 const ParameterList& pL = GetParameterList();
128
129 bool update_communicators = pL.get<bool>("repartition: enable") && pL.get<bool>("repartition: use subcommunicators");
130
131 // If these are set correctly we assume that the nodal P contains only ones
132 bool nodal_p_is_all_ones = !pL.get<bool>("tentative: constant column sums") && !pL.get<bool>("tentative: calculate qr");
133
134 RCP<Matrix> EdgeMatrix = Get< RCP<Matrix> > (fineLevel, "A");
135 RCP<Matrix> D0 = Get< RCP<Matrix> > (fineLevel, "D0");
136 RCP<Matrix> NodeMatrix = Get< RCP<Matrix> > (fineLevel, "NodeAggMatrix");
137 RCP<Matrix> Pn = Get< RCP<Matrix> > (coarseLevel, "Pnodal");
138
139 const GO GO_INVALID = Teuchos::OrdinalTraits<GO>::invalid();
140 const LO LO_INVALID = Teuchos::OrdinalTraits<LO>::invalid();
141
142 // This needs to be an Operator because if NodeMatrix gets repartitioned away, we get an Operator on the level
143 RCP<Operator> CoarseNodeMatrix = Get< RCP<Operator> >(coarseLevel, "NodeAggMatrix");
144 int MyPID = EdgeMatrix.is_null()? -1 : EdgeMatrix->getRowMap()->getComm()->getRank();
145
146 // Matrix matrix params
147 RCP<ParameterList> mm_params = rcp(new ParameterList);;
148 if(pL.isSublist("matrixmatrix: kernel params"))
149 mm_params->sublist("matrixmatrix: kernel params") = pL.sublist("matrixmatrix: kernel params");
150
151
152 // Normalize P
153 if(!nodal_p_is_all_ones) {
154 // The parameters told us the nodal P isn't all ones, so we make a copy that is.
155 GetOStream(Runtime0) << "ReitzingerPFactory::BuildP(): Assuming Pn is not normalized" << std::endl;
156 RCP<Matrix> Pn_old = Pn;
157
158 Pn = Xpetra::MatrixFactory<SC,LO,GO,NO>::Build(Pn->getCrsGraph());
159 Pn->setAllToScalar(Teuchos::ScalarTraits<SC>::one());
160 Pn->fillComplete(Pn->getDomainMap(),Pn->getRangeMap());
161 }
162 else {
163 // The parameters claim P is all ones.
164 GetOStream(Runtime0) << "ReitzingerPFactory::BuildP(): Assuming Pn is normalized" << std::endl;
165 }
166
167 // TODO: We need to make sure Pn isn't normalized. Right now this has to be done explicitly by the user
168
169 // TODO: We need to look through and see which of these really need importers and which ones don't
170
171 /* Generate the D0 * Pn matrix and its transpose */
172 RCP<Matrix> D0_Pn, PnT_D0T, D0_Pn_nonghosted;
173 Teuchos::Array<int> D0_Pn_col_pids;
174 {
175 RCP<Matrix> dummy;
176 SubFactoryMonitor m2(*this, "Generate D0*Pn", coarseLevel);
177 D0_Pn = XMM::Multiply(*D0,false,*Pn,false,dummy,out0,true,true,"D0*Pn",mm_params);
178
179 // We don't want this guy getting accidently used later
180 if(!mm_params.is_null()) mm_params->remove("importer",false);
181
182 // Save this so we don't need to do the multiplication again later
183 D0_Pn_nonghosted = D0_Pn;
184
185 // Get owning PID information on columns for tie-breaking
186 if(!D0_Pn->getCrsGraph()->getImporter().is_null()) {
187 Xpetra::ImportUtils<LO,GO,NO> utils;
188 utils.getPids(*D0_Pn->getCrsGraph()->getImporter(),D0_Pn_col_pids,false);
189 }
190 else {
191 D0_Pn_col_pids.resize(D0_Pn->getCrsGraph()->getColMap()->getLocalNumElements(),MyPID);
192 }
193 }
194
195
196 {
197 // Get the transpose
198 SubFactoryMonitor m2(*this, "Transpose D0*Pn", coarseLevel);
199 PnT_D0T = Utilities::Transpose(*D0_Pn, true);
200 }
201
202
203 // We really need a ghosted version of D0_Pn here.
204 // The reason is that if there's only one fine edge between two coarse nodes, somebody is going
205 // to own the associated coarse edge. The sum/sign rule doesn't guarantee the fine owner is the coarse owner.
206 // So you can wind up with a situation that only guy who *can* register the coarse edge isn't the sum/sign
207 // owner. Adding more ghosting fixes that.
208 if(!PnT_D0T->getCrsGraph()->getImporter().is_null()) {
209 RCP<const Import> Importer = PnT_D0T->getCrsGraph()->getImporter();
210 RCP<const CrsMatrix> D0_Pn_crs = rcp_dynamic_cast<const CrsMatrixWrap>(D0_Pn)->getCrsMatrix();
211 RCP<Matrix> D0_Pn_new = rcp(new CrsMatrixWrap(CrsMatrixFactory::Build(D0_Pn_crs,*Importer,D0_Pn->getDomainMap(),Importer->getTargetMap())));
212 D0_Pn = D0_Pn_new;
213 // Get owning PID information on columns for tie-breaking
214 if(!D0_Pn->getCrsGraph()->getImporter().is_null()) {
215 Xpetra::ImportUtils<LO,GO,NO> utils;
216 utils.getPids(*D0_Pn->getCrsGraph()->getImporter(),D0_Pn_col_pids,false);
217 }
218 else {
219 D0_Pn_col_pids.resize(D0_Pn->getCrsGraph()->getColMap()->getLocalNumElements(),MyPID);
220 }
221 }
222
223
224 // FIXME: This is using deprecated interfaces
225 ArrayView<const LO> colind_E, colind_N;
226 ArrayView<const SC> values_E, values_N;
227
228 size_t Ne=EdgeMatrix->getLocalNumRows();
229 size_t Nn=NodeMatrix->getLocalNumRows();
230
231 // Upper bound on local number of coarse edges
232 size_t max_edges = (NodeMatrix->getLocalNumEntries() + Nn +1) / 2;
233 ArrayRCP<size_t> D0_rowptr(Ne+1);
234 ArrayRCP<LO> D0_colind(max_edges);
235 ArrayRCP<SC> D0_values(max_edges);
236 D0_rowptr[0] = 0;
237
238 LO current = 0;
239 LO Nnc = PnT_D0T->getRowMap()->getLocalNumElements();
240
241 for(LO i=0; i<(LO)Nnc; i++) {
242 // GO global_i = PnT_D0T->getRowMap()->getGlobalElement(i);
243
244 // FIXME: We don't really want an std::map here. This is just a first cut implementation
245 using value_type = bool;
246 std::map<LO, value_type> ce_map;
247
248 // FIXME: This is using deprecated interfaces
249 PnT_D0T->getLocalRowView(i,colind_E,values_E);
250
251 for(LO j=0; j<(LO)colind_E.size(); j++) {
252
253 // NOTE: Edges between procs will be via handled via the a version
254 // of ML's odd/even rule
255 // For this to function correctly, we make two assumptions:
256 // (a) The processor that owns a fine edge owns at least one of the attached nodes.
257 // (b) Aggregation is uncoupled.
258
259 // TODO: Add some debug code to check the assumptions
260
261 // Check to see if we own this edge and continue if we don't
262 GO edge_gid = PnT_D0T->getColMap()->getGlobalElement(colind_E[j]);
263 LO j_row = D0_Pn->getRowMap()->getLocalElement(edge_gid);
264 int pid0, pid1;
265 D0_Pn->getLocalRowView(j_row,colind_N,values_N);
266
267 // Skip incomplete rows
268 if(colind_N.size() != 2) continue;
269
270 pid0 = D0_Pn_col_pids[colind_N[0]];
271 pid1 = D0_Pn_col_pids[colind_N[1]];
272 // printf("[%d] Row %d considering edge (%d)%d -> (%d)%d\n",MyPID,global_i,colind_N[0],D0_Pn->getColMap()->getGlobalElement(colind_N[0]),colind_N[1],D0_Pn->getColMap()->getGlobalElement(colind_N[1]));
273
274 // Check to see who owns these nodes
275 // If the sum of owning procs is odd, the lower ranked proc gets it
276
277 bool zero_matches = pid0 == MyPID;
278 bool one_matches = pid1 == MyPID;
279 bool keep_shared_edge = false, own_both_nodes = false;
280 if(zero_matches && one_matches) {own_both_nodes=true;}
281 else {
282 int sum_is_odd = (pid0 + pid1) % 2;
283 int i_am_smaller = MyPID == std::min(pid0,pid1);
284 if(sum_is_odd && i_am_smaller) keep_shared_edge=true;
285 if(!sum_is_odd && !i_am_smaller) keep_shared_edge=true;
286 }
287 // printf("[%d] - matches %d/%d keep_shared = %d own_both = %d\n",MyPID,(int)zero_matches,(int)one_matches,(int)keep_shared_edge,(int)own_both_nodes);
288 if(!keep_shared_edge && !own_both_nodes) continue;
289
290
291 // We're doing this in GID space, but only because it allows us to explain
292 // the edge orientation as "always goes from lower GID to higher GID". This could
293 // be done entirely in local GIDs, but then the ordering is a little more confusing.
294 // This could be done in local indices later if we need the extra performance.
295 for(LO k=0; k<(LO)colind_N.size(); k++) {
296 LO my_colind = colind_N[k];
297 if(my_colind!=LO_INVALID && ((keep_shared_edge && my_colind != i) || (own_both_nodes && my_colind > i)) ) {
298 ce_map.emplace(std::make_pair(my_colind,true));
299 }
300 }//end for k < colind_N.size()
301 }// end for j < colind_E.size()
302
303
304 // std::map is sorted, so we'll just iterate through this
305 for(auto iter=ce_map.begin(); iter != ce_map.end(); iter++) {
306 LO col = iter->first;
307 // This shouldn't happen. But in case it did...
308 if(col == i) {
309 continue;
310 }
311
312 // ASSUMPTION: "i" is a valid local column id
313 D0_colind[current] = i;
314 D0_values[current] = -1;
315 current++;
316 D0_colind[current] = col;
317 D0_values[current] = 1;
318 current++;
319 D0_rowptr[current / 2] = current;
320 }
321
322 }// end for i < Nn
323
324 LO num_coarse_edges = current / 2;
325 D0_rowptr.resize(num_coarse_edges+1);
326 D0_colind.resize(current);
327 D0_values.resize(current);
328
329 // We're assuming that if the coarse NodeMatrix has no nodes on a rank, the coarse edge guy won't either.
330 // We check that here.
331 TEUCHOS_TEST_FOR_EXCEPTION( (num_coarse_edges > 0 && CoarseNodeMatrix.is_null()) ||
332 (num_coarse_edges == 0 && !CoarseNodeMatrix.is_null())
333 , Exceptions::RuntimeError, "MueLu::ReitzingerPFactory: Mismatched num_coarse_edges and NodeMatrix repartition.");
334
335
336 // Count the total number of edges
337 // NOTE: Since we solve the ownership issue above, this should do what we want
338 RCP<const Map> ownedCoarseEdgeMap = Xpetra::MapFactory<LO,GO,NO>::Build(EdgeMatrix->getRowMap()->lib(), GO_INVALID, num_coarse_edges,EdgeMatrix->getRowMap()->getIndexBase(),EdgeMatrix->getRowMap()->getComm());
339
340
341 // NOTE: This only works because of the assumptions above
342 RCP<const Map> ownedCoarseNodeMap = Pn->getDomainMap();
343 RCP<const Map> ownedPlusSharedCoarseNodeMap = D0_Pn->getCrsGraph()->getColMap();
344
345 // Create the coarse D0
346 RCP<CrsMatrix> D0_coarse;
347 {
348 SubFactoryMonitor m2(*this, "Build D0", coarseLevel);
349 // FIXME: We can be smarter with memory here
350 // TODO: Is there a smarter way to get this importer?
351 D0_coarse = CrsMatrixFactory::Build(ownedCoarseEdgeMap,ownedPlusSharedCoarseNodeMap,0);
352 TEUCHOS_TEST_FOR_EXCEPTION(D0_coarse.is_null(), Exceptions::RuntimeError, "MueLu::ReitzingerPFactory: CrsMatrixFatory failed.");
353
354 // FIXME: Deprecated code
355 ArrayRCP<size_t> ia;
356 ArrayRCP<LO> ja;
357 ArrayRCP<SC> val;
358 D0_coarse->allocateAllValues(current, ia, ja, val);
359 std::copy(D0_rowptr.begin(),D0_rowptr.end(),ia.begin());
360 std::copy(D0_colind.begin(),D0_colind.end(),ja.begin());
361 std::copy(D0_values.begin(),D0_values.end(),val.begin());
362 D0_coarse->setAllValues(ia, ja, val);
363
364#if 0
365 {
366 char fname[80];
367 printf("[%d] D0: ia.size() = %d ja.size() = %d\n",MyPID,(int)ia.size(),(int)ja.size());
368 printf("[%d] D0: ia :",MyPID);
369 for(int i=0; i<(int)ia.size(); i++)
370 printf("%d ",(int)ia[i]);
371 printf("\n[%d] D0: global ja :",MyPID);
372 for(int i=0; i<(int)ja.size(); i++)
373 printf("%d ",(int)ownedPlusSharedCoarseNodeMap->getGlobalElement(ja[i]));
374 printf("\n[%d] D0: local ja :",MyPID);
375 for(int i=0; i<(int)ja.size(); i++)
376 printf("%d ",(int)ja[i]);
377 printf("\n");
378
379 sprintf(fname,"D0_global_ja_%d_%d.dat",MyPID,fineLevel.GetLevelID());
380 FILE * f = fopen(fname,"w");
381 for(int i=0; i<(int)ja.size(); i++)
382 fprintf(f,"%d ",(int)ownedPlusSharedCoarseNodeMap->getGlobalElement(ja[i]));
383 fclose(f);
384
385 sprintf(fname,"D0_local_ja_%d_%d.dat",MyPID,fineLevel.GetLevelID());
386 f = fopen(fname,"w");
387 for(int i=0; i<(int)ja.size(); i++)
388 fprintf(f,"%d ",(int)ja[i]);
389 fclose(f);
390
391 }
392#endif
393 D0_coarse->expertStaticFillComplete(ownedCoarseNodeMap,ownedCoarseEdgeMap);
394 }
395 RCP<Matrix> D0_coarse_m = rcp(new CrsMatrixWrap(D0_coarse));
396 RCP<Teuchos::FancyOStream> fout = Teuchos::fancyOStream(Teuchos::rcpFromRef(std::cout));
397
398
399 // Create the Pe matrix, but with the extra entries. From ML's notes:
400 /* The general idea is that the matrix */
401 /* T_h P_n T_H^* */
402 /* is almost Pe. If we make sure that P_n contains 1's and -1's, the*/
403 /* matrix triple product will yield a matrix with +/- 1 and +/- 2's.*/
404 /* If we remove all the 1's and divide the 2's by 2. we arrive at Pe*/
405 RCP<Matrix> Pe;
406 {
407 SubFactoryMonitor m2(*this, "Generate Pe (pre-fix)", coarseLevel);
408
409 RCP<Matrix> dummy;
410 RCP<Matrix> Pn_D0cT = XMM::Multiply(*Pn,false,*D0_coarse_m,true,dummy,out0,true,true,"Pn*D0c'",mm_params);
411
412 // We don't want this guy getting accidently used later
413 if(!mm_params.is_null()) mm_params->remove("importer",false);
414
415 Pe = XMM::Multiply(*D0,false,*Pn_D0cT,false,dummy,out0,true,true,"D0*(Pn*D0c')",mm_params);
416
417 // TODO: Something like this *might* work. But this specifically, doesn't
418 // Pe = XMM::Multiply(*D0_Pn_nonghosted,false,*D0_coarse_m,true,dummy,out0,true,true,"(D0*Pn)*D0c'",mm_params);
419 }
420
421 /* Weed out the +/- entries, shrinking the matrix as we go */
422 {
423 SubFactoryMonitor m2(*this, "Generate Pe (post-fix)", coarseLevel);
424 Pe->resumeFill();
425 SC one = Teuchos::ScalarTraits<SC>::one();
426 MT two = 2*Teuchos::ScalarTraits<MT>::one();
427 SC zero = Teuchos::ScalarTraits<SC>::zero();
428 SC neg_one = - one;
429
430 RCP<const CrsMatrix> Pe_crs = rcp_dynamic_cast<const CrsMatrixWrap>(Pe)->getCrsMatrix();
431 TEUCHOS_TEST_FOR_EXCEPTION(Pe_crs.is_null(), Exceptions::RuntimeError, "MueLu::ReitzingerPFactory: Pe is not a crs matrix.");
432 ArrayRCP<const size_t > rowptr_const;
433 ArrayRCP<const LO> colind_const;
434 ArrayRCP<const SC> values_const;
435 Pe_crs->getAllValues(rowptr_const,colind_const,values_const);
436 ArrayRCP<size_t> rowptr = arcp_const_cast<size_t>(rowptr_const);
437 ArrayRCP<LO> colind = arcp_const_cast<LO>(colind_const);
438 ArrayRCP<SC> values = arcp_const_cast<SC>(values_const);
439 LO ct = 0;
440 LO lower = rowptr[0];
441 for(LO i=0; i<(LO)Ne; i++) {
442 for(size_t j=lower; j<rowptr[i+1]; j++) {
443 if (values[j] == one || values[j] == neg_one || values[j] == zero) {
444 // drop this guy
445 }
446 else {
447 colind[ct] = colind[j];
448 values[ct] = values[j] / two;
449 ct++;
450 }
451 }
452 lower = rowptr[i+1];
453 rowptr[i+1] = ct;
454 }
455 rowptr[Ne] = ct;
456 colind.resize(ct);
457 values.resize(ct);
458 rcp_const_cast<CrsMatrix>(Pe_crs)->setAllValues(rowptr,colind,values);
459
460 Pe->fillComplete(Pe->getDomainMap(),Pe->getRangeMap());
461 }
462
463 /* Check commuting property */
464 CheckCommutingProperty(*Pe,*D0_coarse_m,*D0,*Pn);
465
466 /* If we're repartitioning here, we need to cut down the communicators */
467 // NOTE: We need to do this *after* checking the commuting property, since
468 // that's going to need to fineLevel's communicators, not the repartitioned ones
469 if(update_communicators) {
470 //NOTE: We can only do D0 here. We have to do Ke_coarse=(Re Ke_fine Pe) in RebalanceAcFactory
471 RCP<const Teuchos::Comm<int> > newComm;
472 if(!CoarseNodeMatrix.is_null()) newComm = CoarseNodeMatrix->getDomainMap()->getComm();
473 RCP<const Map> newMap = Xpetra::MapFactory<LO,GO,NO>::copyMapWithNewComm(D0_coarse_m->getRowMap(),newComm);
474 D0_coarse_m->removeEmptyProcessesInPlace(newMap);
475
476 // The "in place" still leaves a dummy matrix here. That needs to go
477 if(newMap.is_null()) D0_coarse_m = Teuchos::null;
478
479 Set(coarseLevel,"InPlaceMap",newMap);
480 }
481
482 /* Set output on the level */
483 Set(coarseLevel,"P",Pe);
484 Set(coarseLevel,"Ptent",Pe);
485
486 Set(coarseLevel,"D0",D0_coarse_m);
487 coarseLevel.Set("D0",D0_coarse_m,NoFactory::get());
488 coarseLevel.AddKeepFlag("D0",NoFactory::get(), MueLu::Final);
489 coarseLevel.RemoveKeepFlag("D0",NoFactory::get(), MueLu::UserData);
490
491#if 0
492 {
493 int numProcs = Pe->getRowMap()->getComm()->getSize();
494 char fname[80];
495
496 sprintf(fname,"Pe_%d_%d.mat",numProcs,fineLevel.GetLevelID()); Xpetra::IO<SC,LO,GO,NO>::Write(fname,*Pe);
497 sprintf(fname,"Pn_%d_%d.mat",numProcs,fineLevel.GetLevelID()); Xpetra::IO<SC,LO,GO,NO>::Write(fname,*Pn);
498 sprintf(fname,"D0c_%d_%d.mat",numProcs,fineLevel.GetLevelID()); Xpetra::IO<SC,LO,GO,NO>::Write(fname,*D0_coarse_m);
499 sprintf(fname,"D0f_%d_%d.mat",numProcs,fineLevel.GetLevelID()); Xpetra::IO<SC,LO,GO,NO>::Write(fname,*D0);
500 }
501#endif
502
503 }// end Build
504
505 template <class Scalar,class LocalOrdinal, class GlobalOrdinal, class Node>
507 CheckCommutingProperty(const Matrix & Pe, const Matrix & D0_c, const Matrix& D0_f, const Matrix & Pn) const {
508 if(IsPrint(Statistics0)) {
509 using XMM = Xpetra::MatrixMatrix<SC,LO,GO,NO>;
510 using MT = typename Teuchos::ScalarTraits<SC>::magnitudeType;
511 SC one = Teuchos::ScalarTraits<SC>::one();
512 SC zero = Teuchos::ScalarTraits<SC>::zero();
513
514 RCP<Matrix> dummy;
515 Teuchos::FancyOStream &out0=GetBlackHole();
516 RCP<Matrix> left = XMM::Multiply(Pe,false,D0_c,false,dummy,out0);
517 RCP<Matrix> right = XMM::Multiply(D0_f,false,Pn,false,dummy,out0);
518
519 // We need a non-FC matrix for the add, sadly
520 RCP<CrsMatrix> sum_c = CrsMatrixFactory::Build(left->getRowMap(),left->getLocalMaxNumRowEntries()+right->getLocalMaxNumRowEntries());
521 RCP<Matrix> summation = rcp(new CrsMatrixWrap(sum_c));
522 XMM::TwoMatrixAdd(*left, false, one, *summation, zero);
523 XMM::TwoMatrixAdd(*right, false, -one, *summation, one);
524
525 MT norm = summation->getFrobeniusNorm();
526 GetOStream(Statistics0) << "CheckCommutingProperty: ||Pe D0_c - D0_f Pn || = "<<norm<<std::endl;
527
528 }
529
530 }//end CheckCommutingProperty
531
532
533
534} //namespace MueLu
535
536
537
538#define MUELU_REITZINGERPFACTORY_SHORT
539#endif // MUELU_REITZINGERPFACTORY_DEF_HPP
#define SET_VALID_ENTRY(name)
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 RemoveKeepFlag(const std::string &ename, const FactoryBase *factory, KeepType keep=MueLu::All)
int GetLevelID() const
Return level number.
Definition: MueLu_Level.cpp:76
void AddKeepFlag(const std::string &ename, const FactoryBase *factory=NoFactory::get(), KeepType keep=MueLu::Keep)
void Set(const std::string &ename, const T &entry, const FactoryBase *factory=NoFactory::get())
static const NoFactory * get()
void BuildP(Level &fineLevel, Level &coarseLevel) const
Abstract Build method.
void Build(Level &fineLevel, Level &coarseLevel) const
Build an object with this factory.
void CheckCommutingProperty(const Matrix &Pe, const Matrix &D0_c, const Matrix &D0_f, const Matrix &Pn) const
Utility method.
RCP< const ParameterList > GetValidParameterList() const
Return a const parameter list of valid parameters that setParameterList() will accept.
void DeclareInput(Level &fineLevel, Level &coarseLevel) const
Input.
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)
Namespace for MueLu classes and methods.
@ Final
Keep data only for this run. Used to keep data useful for Hierarchy::Iterate(). Data will be deleted ...
@ UserData
User data are always kept. This flag is set automatically when Level::Set("data", data) is used....
@ Runtime0
One-liner description of what is happening.
@ Statistics0
Print statistics that do not involve significant additional computation.