MueLu  Version of the Day
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
MueLu_AmalgamationFactory_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_AMALGAMATIONFACTORY_DEF_HPP
47 #define MUELU_AMALGAMATIONFACTORY_DEF_HPP
48 
49 #include <Xpetra_Matrix.hpp>
50 
51 #include "MueLu_AmalgamationFactory.hpp"
52 
53 #include "MueLu_Level.hpp"
54 #include "MueLu_AmalgamationInfo.hpp"
55 #include "MueLu_Monitor.hpp"
56 
57 namespace MueLu {
58 
59  template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
61  RCP<ParameterList> validParamList = rcp(new ParameterList());
62  validParamList->set< RCP<const FactoryBase> >("A", Teuchos::null, "Generating factory of the matrix A");
63  return validParamList;
64  }
65 
66  template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
68  Input(currentLevel, "A"); // sub-block from blocked A
69  }
70 
71  template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
73  {
74  FactoryMonitor m(*this, "Build", currentLevel);
75 
76  RCP<Matrix> A = Get< RCP<Matrix> >(currentLevel, "A");
77 
78  LO fullblocksize = 1; // block dim for fixed size blocks
79  GO offset = 0; // global offset of dof gids
80  LO blockid = -1; // block id in strided map
81  LO nStridedOffset = 0; // DOF offset for strided block id "blockid" (default = 0)
82  LO stridedblocksize = fullblocksize; // size of strided block id "blockid" (default = fullblocksize, only if blockid!=-1 stridedblocksize <= fullblocksize)
83  // GO indexBase = A->getRowMap()->getIndexBase(); // index base for maps (unused)
84 
85  // 1) check for blocking/striding information
86 
87  if (A->IsView("stridedMaps") && Teuchos::rcp_dynamic_cast<const StridedMap>(A->getRowMap("stridedMaps")) != Teuchos::null) {
88  Xpetra::viewLabel_t oldView = A->SwitchToView("stridedMaps"); // NOTE: "stridedMaps are always non-overlapping (correspond to range and domain maps!)
89  RCP<const StridedMap> stridedRowMap = Teuchos::rcp_dynamic_cast<const StridedMap>(A->getRowMap());
90  TEUCHOS_TEST_FOR_EXCEPTION(stridedRowMap == Teuchos::null,Exceptions::BadCast,"MueLu::CoalesceFactory::Build: cast to strided row map failed.");
91  fullblocksize = stridedRowMap->getFixedBlockSize();
92  offset = stridedRowMap->getOffset();
93  blockid = stridedRowMap->getStridedBlockId();
94 
95  if (blockid > -1) {
96  std::vector<size_t> stridingInfo = stridedRowMap->getStridingData();
97  for (size_t j = 0; j < Teuchos::as<size_t>(blockid); j++)
98  nStridedOffset += stridingInfo[j];
99  stridedblocksize = Teuchos::as<LocalOrdinal>(stridingInfo[blockid]);
100 
101  } else {
102  stridedblocksize = fullblocksize;
103  }
104  oldView = A->SwitchToView(oldView);
105  GetOStream(Runtime1) << "AmalagamationFactory::Build():" << " found fullblocksize=" << fullblocksize << " and stridedblocksize=" << stridedblocksize << " from strided maps. offset=" << offset << std::endl;
106 
107  } else {
108  GetOStream(Warnings0) << "AmalagamationFactory::Build(): no striding information available. Use blockdim=1 with offset=0" << std::endl;
109  }
110 
111  // build node row map (uniqueMap) and node column map (nonUniqueMap)
112  // the arrays rowTranslation and colTranslation contain the local node id
113  // given a local dof id. They are only necessary for the CoalesceDropFactory if
114  // fullblocksize > 1
115  RCP<const Map> uniqueMap, nonUniqueMap;
116  RCP<AmalgamationInfo> amalgamationData;
117  RCP<Array<LO> > rowTranslation = Teuchos::null;
118  RCP<Array<LO> > colTranslation = Teuchos::null;
119 
120  if (fullblocksize > 1) {
121  // mfh 14 Apr 2015: These need to have different names than
122  // rowTranslation and colTranslation, in order to avoid
123  // shadowing warnings (-Wshadow with GCC). Alternately, it
124  // looks like you could just assign to the existing variables in
125  // this scope, rather than creating new ones.
126  RCP<Array<LO> > theRowTranslation = rcp(new Array<LO>);
127  RCP<Array<LO> > theColTranslation = rcp(new Array<LO>);
128  AmalgamateMap(*(A->getRowMap()), *A, uniqueMap, *theRowTranslation);
129  AmalgamateMap(*(A->getColMap()), *A, nonUniqueMap, *theColTranslation);
130 
131  amalgamationData = rcp(new AmalgamationInfo(theRowTranslation,
132  theColTranslation,
133  uniqueMap,
134  nonUniqueMap,
135  A->getColMap(),
136  fullblocksize,
137  offset,
138  blockid,
139  nStridedOffset,
140  stridedblocksize) );
141  } else {
142  amalgamationData = rcp(new AmalgamationInfo(rowTranslation, // Teuchos::null
143  colTranslation, // Teuchos::null
144  A->getRowMap(), // unique map of graph
145  A->getColMap(), // non-unique map of graph
146  A->getColMap(), // column map of A
147  fullblocksize,
148  offset,
149  blockid,
150  nStridedOffset,
151  stridedblocksize) );
152  }
153 
154  // store (un)amalgamation information on current level
155  Set(currentLevel, "UnAmalgamationInfo", amalgamationData);
156  }
157 
158  template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
159  void AmalgamationFactory<Scalar, LocalOrdinal, GlobalOrdinal, Node>::AmalgamateMap(const Map& sourceMap, const Matrix& A, RCP<const Map>& amalgamatedMap, Array<LO>& translation) {
160  typedef typename ArrayView<const GO>::size_type size_type;
161  typedef std::map<GO,size_type> container;
162 
163  GO indexBase = sourceMap.getIndexBase();
164  ArrayView<const GO> elementAList = sourceMap.getNodeElementList();
165  size_type numElements = elementAList.size();
166  container filter; // TODO: replace std::set with an object having faster lookup/insert, hashtable for instance
167 
168  GO offset = 0;
169  LO blkSize = A.GetFixedBlockSize();
170  if (A.IsView("stridedMaps") == true) {
171  Teuchos::RCP<const Map> myMap = A.getRowMap("stridedMaps");
172  Teuchos::RCP<const StridedMap> strMap = Teuchos::rcp_dynamic_cast<const StridedMap>(myMap);
173  TEUCHOS_TEST_FOR_EXCEPTION(strMap == null, Exceptions::RuntimeError, "Map is not of type StridedMap");
174  offset = strMap->getOffset();
175  blkSize = Teuchos::as<const LO>(strMap->getFixedBlockSize());
176  }
177 
178  Array<GO> elementList(numElements);
179  translation.resize(numElements);
180 
181  size_type numRows = 0;
182  for (size_type id = 0; id < numElements; id++) {
183  GO dofID = elementAList[id];
184  GO nodeID = AmalgamationFactory::DOFGid2NodeId(dofID, blkSize, offset, indexBase);
185 
186  typename container::iterator it = filter.find(nodeID);
187  if (it == filter.end()) {
188  filter[nodeID] = numRows;
189 
190  translation[id] = numRows;
191  elementList[numRows] = nodeID;
192 
193  numRows++;
194 
195  } else {
196  translation[id] = it->second;
197  }
198  }
199  elementList.resize(numRows);
200 
201  amalgamatedMap = MapFactory::Build(sourceMap.lib(), Teuchos::OrdinalTraits<Xpetra::global_size_t>::invalid(), elementList, indexBase, sourceMap.getComm());
202 
203  }
204 
205  template <class Scalar,class LocalOrdinal, class GlobalOrdinal, class Node>
206  const GlobalOrdinal AmalgamationFactory<Scalar, LocalOrdinal, GlobalOrdinal, Node>::DOFGid2NodeId(GlobalOrdinal gid, LocalOrdinal blockSize, const GlobalOrdinal offset, const GlobalOrdinal indexBase) {
207  // here, the assumption is, that the node map has the same indexBase as the dof map
208  GlobalOrdinal globalblockid = ((GlobalOrdinal) gid - offset - indexBase) / blockSize + indexBase;
209  return globalblockid;
210  }
211 
212 } //namespace MueLu
213 
214 #endif /* MUELU_SUBBLOCKUNAMALGAMATIONFACTORY_DEF_HPP */
215 
Important warning messages (one line)
Exception indicating invalid cast attempted.
GlobalOrdinal GO
ParameterList & set(std::string const &name, T const &value, std::string const &docString="", RCP< const ParameterEntryValidator > const &validator=null)
Timer to be used in factories. Similar to Monitor but with additional timers.
#define TEUCHOS_TEST_FOR_EXCEPTION(throw_exception_test, Exception, msg)
size_type size() const
LocalOrdinal LO
static const GlobalOrdinal DOFGid2NodeId(GlobalOrdinal gid, LocalOrdinal blockSize, const GlobalOrdinal offset, const GlobalOrdinal indexBase)
translate global (row/column) id to global amalgamation block id
void Build(Level &currentLevel) const
Build an object with this factory.
TEUCHOS_DEPRECATED RCP< T > rcp(T *p, Dealloc_T dealloc, bool owns_mem)
Class that holds all level-specific information.
Definition: MueLu_Level.hpp:99
static void AmalgamateMap(const Map &sourceMap, const Matrix &A, RCP< const Map > &amalgamatedMap, Array< LO > &translation)
Method to create merged map for systems of PDEs.
void resize(size_type new_size, const value_type &x=value_type())
RCP< const ParameterList > GetValidParameterList() const
Return a const parameter list of valid parameters that setParameterList() will accept.
std::string viewLabel_t
Exception throws to report errors in the internal logical of the program.
Description of what is happening (more verbose)
void DeclareInput(Level &currentLevel) const
Input.
minimal container class for storing amalgamation information