MueLu  Version of the Day
MueLu_BlockedCoordinatesTransferFactory_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_BLOCKEDCOORDINATESTRANSFER_FACTORY_DEF_HPP
47 #define MUELU_BLOCKEDCOORDINATESTRANSFER_FACTORY_DEF_HPP
48 
49 #include "Xpetra_ImportFactory.hpp"
50 #include "Xpetra_MultiVectorFactory.hpp"
51 #include "Xpetra_MapFactory.hpp"
52 #include "Xpetra_IO.hpp"
53 
54 #include "MueLu_CoarseMapFactory.hpp"
55 #include "MueLu_Aggregates.hpp"
57 
58 #include "MueLu_Level.hpp"
59 #include "MueLu_Monitor.hpp"
60 
61 namespace MueLu {
62 
63  template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
65  RCP<ParameterList> validParamList = rcp(new ParameterList());
66 
67  validParamList->set<RCP<const FactoryBase> >("Coordinates", Teuchos::null, "Factory for coordinates generation");
68  validParamList->set<RCP<const FactoryBase> >("CoarseMap", Teuchos::null, "Generating factory of the coarse map");
69  return validParamList;
70  }
71 
72  template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
74  Input(coarseLevel, "CoarseMap");
75 
76  // Make sure the Level knows I need these sub-Factories
77  const size_t numSubFactories = NumFactories();
78  for(size_t i=0; i<numSubFactories; i++) {
79  const RCP<const FactoryBase>& myFactory = subFactories_[i];
80  coarseLevel.DeclareInput("Coordinates", myFactory.getRawPtr(), this);
81  }
82 
83  // call DeclareInput of all user-given transfer factories
84  for (std::vector<RCP<const FactoryBase> >::const_iterator it = subFactories_.begin(); it != subFactories_.end(); ++it)
85  (*it)->CallDeclareInput(coarseLevel);
86  }
87 
88  template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
90  FactoryMonitor m(*this, "Build", coarseLevel);
91 
92  typedef Xpetra::MultiVector<typename Teuchos::ScalarTraits<Scalar>::coordinateType,LO,GO,NO> dMV;
93  typedef Xpetra::BlockedMultiVector<typename Teuchos::ScalarTraits<Scalar>::coordinateType,LO,GO,NO> dBV;
94 
95  GetOStream(Runtime0) << "Transferring (blocked) coordinates" << std::endl;
96 
97  const size_t numSubFactories = NumFactories();
98  std::vector<RCP<const Map> > subBlockMaps(numSubFactories);
99  std::vector<RCP<dMV> > subBlockCoords(numSubFactories);
100 
101  if (coarseLevel.IsAvailable("Coordinates", this)) {
102  GetOStream(Runtime0) << "Reusing coordinates" << std::endl;
103  return;
104  }
105 
106  // Get components
107  for(size_t i=0; i<numSubFactories; i++) {
108  GetOStream(Runtime1) << "Generating Coordinates for block " << i <<"/"<<numSubFactories <<std::endl;
109  const RCP<const FactoryBase>& myFactory = subFactories_[i];
110  myFactory->CallBuild(coarseLevel);
111  subBlockCoords[i] = coarseLevel.Get<RCP<dMV> >("Coordinates", myFactory.get());
112  subBlockMaps[i] = subBlockCoords[i]->getMap();
113  }
114 
115  // Blocked Map
116  RCP<const BlockedMap> coarseCoordMapBlocked;
117 
118  {
119  // coarseMap is being used to set up the domain map of tentative P, and therefore, the row map of Ac
120  // Therefore, if we amalgamate coarseMap, logical nodes in the coordinates vector would correspond to
121  // logical blocks in the matrix
122  RCP<const BlockedMap> coarseMap = Get< RCP<const BlockedMap> >(coarseLevel, "CoarseMap");
123  bool thyraMode = coarseMap->getThyraMode();
124 
125  ArrayView<const GO> elementAList = coarseMap->getFullMap()->getNodeElementList();
126 
127  LO blkSize = 1;
128  if (rcp_dynamic_cast<const StridedMap>(coarseMap->getMap(0, thyraMode)) != Teuchos::null)
129  blkSize = rcp_dynamic_cast<const StridedMap>(coarseMap->getMap(0, thyraMode))->getFixedBlockSize();
130 
131  for(size_t i=1; i<numSubFactories; i++) {
132  LO otherBlkSize = 1;
133  if (rcp_dynamic_cast<const StridedMap>(coarseMap->getMap(i, thyraMode)) != Teuchos::null)
134  otherBlkSize = rcp_dynamic_cast<const StridedMap>(coarseMap->getMap(i, thyraMode))->getFixedBlockSize();
135  TEUCHOS_TEST_FOR_EXCEPTION(otherBlkSize != blkSize, Exceptions::RuntimeError, "BlockedCoordinatesTransferFactory: Subblocks have different Block sizes. This is not yet supported.");
136  }
137 
138  GO indexBase = coarseMap->getFullMap()->getIndexBase();
139  size_t numElements = elementAList.size() / blkSize;
140  Array<GO> elementList(numElements);
141 
142  // Amalgamate the map
143  for (LO i = 0; i < Teuchos::as<LO>(numElements); i++)
144  elementList[i] = (elementAList[i*blkSize]-indexBase)/blkSize + indexBase;
145 
146  RCP<const Map> coarseCoordMap = MapFactory::Build(coarseMap->getFullMap()->lib(),
147  Teuchos::OrdinalTraits<Xpetra::global_size_t>::invalid(), elementList, indexBase, coarseMap->getFullMap()->getComm());
148 
149  coarseCoordMapBlocked = rcp(new BlockedMap(coarseCoordMap, subBlockMaps, thyraMode));
150  }
151 
152  // Build blocked coordinates vector
153  RCP<dBV> bcoarseCoords = rcp(new dBV(coarseCoordMapBlocked,subBlockCoords));
154 
155  // Turn the blocked coordinates vector into an unblocked one
156  RCP<dMV> coarseCoords = bcoarseCoords->Merge();
157  Set<RCP<dMV> >(coarseLevel, "Coordinates", coarseCoords);
158  }
159 
160  template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
162  subFactories_.push_back(factory);
163  }
164 
165 
166 
167 } // namespace MueLu
168 
169 #endif // MUELU_BLOCKEDCOORDINATESTRANSFER_FACTORY_DEF_HPP
T & Get(const std::string &ename, const FactoryBase *factory=NoFactory::get())
Get data without decrementing associated storage counter (i.e., read-only access). Usage: Level->Get< RCP<Matrix> >("A", factory) if factory == NULL => use default factory.
Timer to be used in factories. Similar to Monitor but with additional timers.
One-liner description of what is happening.
void DeclareInput(Level &finelevel, Level &coarseLevel) const
Specifies the data that this class needs, and the factories that generate that data.
Namespace for MueLu classes and methods.
void Build(Level &fineLevel, Level &coarseLevel) const
Build an object with this factory.
Class that holds all level-specific information.
Definition: MueLu_Level.hpp:99
void AddFactory(const RCP< const FactoryBase > &factory)
Add (sub) coords factory in the end of list of factories in BlockedCoordinatesTransferFactory.
bool IsAvailable(const std::string &ename, const FactoryBase *factory=NoFactory::get()) const
Test whether a need&#39;s value has been saved.
Exception throws to report errors in the internal logical of the program.
RCP< const ParameterList > GetValidParameterList() const
Return a const parameter list of valid parameters that setParameterList() will accept.
Description of what is happening (more verbose)
void DeclareInput(const std::string &ename, const FactoryBase *factory, const FactoryBase *requestedBy=NoFactory::get())
Callback from FactoryBase::CallDeclareInput() and FactoryBase::DeclareInput()