MueLu  Version of the Day
MueLu_UnsmooshFactory_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 // Tobias Wiesner (tawiesn@sandia.gov)
42 // Ray Tuminaro (rstumin@sandia.gov)
43 //
44 // ***********************************************************************
45 //
46 // @HEADER
47 #ifndef PACKAGES_MUELU_SRC_GRAPH_MUELU_UNSMOOSHFACTORY_DEF_HPP_
48 #define PACKAGES_MUELU_SRC_GRAPH_MUELU_UNSMOOSHFACTORY_DEF_HPP_
49 
50 #include "MueLu_Monitor.hpp"
51 
53 
54 namespace MueLu {
55 
56  template <class Scalar,class LocalOrdinal, class GlobalOrdinal, class Node>
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 for unamalgamated matrix. Row map of (unamalgamted) output prolongation operator should match row map of this A.");
63  validParamList->set< RCP<const FactoryBase> >("P", Teuchos::null, "Generating factory of the (amalgamated) prolongator P");
64  validParamList->set< RCP<const FactoryBase> >("DofStatus", Teuchos::null, "Generating factory for dofStatus array (usually the VariableDofLaplacdianFactory)");
65 
66  validParamList->set< int > ("maxDofPerNode", 1, "Maximum number of DOFs per node");
67  validParamList->set< bool > ("fineIsPadded" , false, "true if finest level input matrix is padded");
68 
69  return validParamList;
70  }
71 
72  template <class Scalar,class LocalOrdinal, class GlobalOrdinal, class Node>
74  //const ParameterList& pL = GetParameterList();
75  Input(fineLevel, "A");
76  Input(coarseLevel, "P");
77 
78  // DofStatus only provided on the finest level (by user)
79  // On the coarser levels it is auto-generated using the DBC information from the unamalgamated matrix A
80  if(fineLevel.GetLevelID() == 0)
81  Input(fineLevel, "DofStatus");
82  }
83 
84  template <class Scalar,class LocalOrdinal, class GlobalOrdinal, class Node>
86  FactoryMonitor m(*this, "Build", coarseLevel);
87  typedef Teuchos::ScalarTraits<SC> STS;
88 
89  const ParameterList & pL = GetParameterList();
90 
91  // extract matrices (unamalgamated A and amalgamated P)
92  RCP<Matrix> unamalgA = Get< RCP<Matrix> >(fineLevel, "A");
93  RCP<Matrix> amalgP = Get< RCP<Matrix> >(coarseLevel, "P");
94 
95  // extract user parameters
96  int maxDofPerNode = pL.get<int> ("maxDofPerNode");
97  bool fineIsPadded = pL.get<bool>("fineIsPadded");
98 
99  // get dofStatus information
100  // On the finest level it is provided by the user. On the coarser levels it is constructed
101  // using the DBC information of the matrix A
102  Teuchos::Array<char> dofStatus;
103  if(fineLevel.GetLevelID() == 0) {
104  dofStatus = Get<Teuchos::Array<char> >(fineLevel, "DofStatus");
105  TEUCHOS_TEST_FOR_EXCEPTION(Teuchos::as<size_t>(dofStatus.size()) == Teuchos::as<size_t>(unamalgA->getRowMap()->getNodeNumElements()), MueLu::Exceptions::RuntimeError,"MueLu::UnsmooshFactory::Build: User provided dofStatus on level 0 does not fit to size of unamalgamted A");
106  } else {
107  // dof status is the dirichlet information of unsmooshed/unamalgamated A (fine level)
108  dofStatus = Teuchos::Array<char>(unamalgA->getRowMap()->getNodeNumElements() /*amalgP->getRowMap()->getNodeNumElements() * maxDofPerNode*/,'s');
109 
110  bool bHasZeroDiagonal = false;
112 
113  TEUCHOS_TEST_FOR_EXCEPTION(dirOrNot.size() != dofStatus.size(), MueLu::Exceptions::RuntimeError,"MueLu::UnsmooshFactory::Build: inconsistent number of coarse DBC array and dofStatus array. dirOrNot.size() = " << dirOrNot.size() << " dofStatus.size() = " << dofStatus.size());
114  for(decltype(dirOrNot.size()) i = 0; i < dirOrNot.size(); ++i) {
115  if(dirOrNot[i] == true) dofStatus[i] = 'p';
116  }
117  }
118 
119  // TODO: TAW the following check is invalid for SA-AMG based input prolongators
120  //TEUCHOS_TEST_FOR_EXCEPTION(amalgP->getDomainMap()->isSameAs(*amalgP->getColMap()) == false, MueLu::Exceptions::RuntimeError,"MueLu::UnsmooshFactory::Build: only support for non-overlapping aggregates. (column map of Ptent must be the same as domain map of Ptent)");
121 
122  // extract CRS information from amalgamated prolongation operator
123  Teuchos::ArrayRCP<const size_t> amalgRowPtr(amalgP->getNodeNumRows());
124  Teuchos::ArrayRCP<const LocalOrdinal> amalgCols(amalgP->getNodeNumEntries());
125  Teuchos::ArrayRCP<const Scalar> amalgVals(amalgP->getNodeNumEntries());
126  Teuchos::RCP<CrsMatrixWrap> amalgPwrap = Teuchos::rcp_dynamic_cast<CrsMatrixWrap>(amalgP);
127  Teuchos::RCP<CrsMatrix> amalgPcrs = amalgPwrap->getCrsMatrix();
128  amalgPcrs->getAllValues(amalgRowPtr, amalgCols, amalgVals);
129 
130  // calculate number of dof rows for new prolongator
131  size_t paddedNrows = amalgP->getRowMap()->getNodeNumElements() * Teuchos::as<size_t>(maxDofPerNode);
132 
133  // reserve CSR arrays for new prolongation operator
134  Teuchos::ArrayRCP<size_t> newPRowPtr(paddedNrows+1);
135  Teuchos::ArrayRCP<LocalOrdinal> newPCols(amalgP->getNodeNumEntries() * maxDofPerNode);
136  Teuchos::ArrayRCP<Scalar> newPVals(amalgP->getNodeNumEntries() * maxDofPerNode);
137 
138  size_t rowCount = 0; // actual number of (local) in unamalgamated prolongator
139  if(fineIsPadded == true || fineLevel.GetLevelID() > 0) {
140 
141  // build prolongation operator for padded fine level matrices.
142  // Note: padded fine level dofs are transferred by injection.
143  // That is, these interpolation stencils do not take averages of
144  // coarse level variables. Further, fine level Dirichlet points
145  // also use injection.
146 
147  size_t cnt = 0; // local id counter
148  for (decltype(amalgRowPtr.size()) i = 0; i < amalgRowPtr.size() - 1; i++) {
149  // determine number of entries in amalgamated dof row i
150  size_t rowLength = amalgRowPtr[i+1] - amalgRowPtr[i];
151 
152  // loop over dofs per node (unamalgamation)
153  for(int j = 0; j < maxDofPerNode; j++) {
154  newPRowPtr[i*maxDofPerNode+j] = cnt;
155  if (dofStatus[i*maxDofPerNode+j] == 's') { // add only "standard" dofs to unamalgamated prolongator
156  // loop over column entries in amalgamated P
157  for (size_t k = 0; k < rowLength; k++) {
158  newPCols[cnt ] = amalgCols[k+amalgRowPtr[i]] * maxDofPerNode + j;
159  newPVals[cnt++] = amalgVals[k+amalgRowPtr[i]];
160  }
161 
162  }
163  }
164  }
165 
166  newPRowPtr[paddedNrows] = cnt; // close row CSR array
167  rowCount = paddedNrows;
168  } else {
169  // Build prolongation operator for non-padded fine level matrices.
170  // Need to map from non-padded dofs to padded dofs. For this, look
171  // at the status array and skip padded dofs.
172 
173  size_t cnt = 0; // local id counter
174 
175  for (decltype(amalgRowPtr.size()) i = 0; i < amalgRowPtr.size() - 1; i++) {
176  // determine number of entries in amalgamated dof row i
177  size_t rowLength = amalgRowPtr[i+1] - amalgRowPtr[i];
178 
179  // loop over dofs per node (unamalgamation)
180  for(int j = 0; j < maxDofPerNode; j++) {
181  // no interpolation for padded fine dofs as they do not exist
182 
183  if (dofStatus[i*maxDofPerNode+j] == 's') { // add only "standard" dofs to unamalgamated prolongator
184  newPRowPtr[rowCount++] = cnt;
185  // loop over column entries in amalgamated P
186  for (size_t k = 0; k < rowLength; k++) {
187  newPCols[cnt ] = amalgCols[k+amalgRowPtr[i]] * maxDofPerNode + j;
188  newPVals[cnt++] = amalgVals[k+amalgRowPtr[i]];
189  }
190 
191  }
192  if (dofStatus[i*maxDofPerNode+j] == 'd') { // Dirichlet handling
193  newPRowPtr[rowCount++] = cnt;
194  }
195  }
196  }
197  newPRowPtr[rowCount] = cnt; // close row CSR array
198  } // fineIsPadded == false
199 
200  // generate coarse domain map
201  // So far no support for gid offset or strided maps. This information
202  // could be gathered easily from the unamalgamated fine level operator A.
203  std::vector<size_t> stridingInfo(1,maxDofPerNode);
204 
205  GlobalOrdinal nCoarseDofs = amalgP->getDomainMap()->getNodeNumElements() * maxDofPerNode;
206  GlobalOrdinal indexBase = amalgP->getDomainMap()->getIndexBase();
207  RCP<const Map> coarseDomainMap = StridedMapFactory::Build(amalgP->getDomainMap()->lib(),
209  nCoarseDofs,
210  indexBase,
211  stridingInfo,
212  amalgP->getDomainMap()->getComm(),
213  -1 /* stridedBlockId */,
214  0 /*domainGidOffset */);
215 
216  size_t nColCoarseDofs = Teuchos::as<size_t>(amalgP->getColMap()->getNodeNumElements() * maxDofPerNode);
217  Teuchos::Array<GlobalOrdinal> unsmooshColMapGIDs(nColCoarseDofs);
218  for(size_t c = 0; c < amalgP->getColMap()->getNodeNumElements(); ++c) {
219  GlobalOrdinal gid = amalgP->getColMap()->getGlobalElement(c) * maxDofPerNode;
220  for(int i = 0; i < maxDofPerNode; ++i) {
221  unsmooshColMapGIDs[c * maxDofPerNode + i] = gid + i;
222  }
223  }
224  Teuchos::RCP<Map> coarseColMap = MapFactory::Build(amalgP->getDomainMap()->lib(),
226  unsmooshColMapGIDs(), //View,
227  indexBase,
228  amalgP->getDomainMap()->getComm());
229 
230  // Assemble unamalgamated P
231  Teuchos::RCP<CrsMatrix> unamalgPCrs = CrsMatrixFactory::Build(unamalgA->getRowMap(),coarseColMap, 3);
232  for (decltype(rowCount) i = 0; i < rowCount; i++) {
233  unamalgPCrs->insertLocalValues(i, newPCols.view(newPRowPtr[i],newPRowPtr[i+1]-newPRowPtr[i]),
234  newPVals.view(newPRowPtr[i],newPRowPtr[i+1]-newPRowPtr[i]));
235  }
236  unamalgPCrs->fillComplete(coarseDomainMap,unamalgA->getRowMap());
237 
238  Teuchos::RCP<Matrix> unamalgP = Teuchos::rcp(new CrsMatrixWrap(unamalgPCrs));
239 
240  Set(coarseLevel,"P",unamalgP);
241  }
242 
243 
244 } /* MueLu */
245 
246 
247 #endif /* PACKAGES_MUELU_SRC_GRAPH_MUELU_UNSMOOSHFACTORY_DEF_HPP_ */
RCP< const ParameterList > GetValidParameterList() const
Return a const parameter list of valid parameters that setParameterList() will accept.
size_type size() const
T & get(const std::string &name, T def_value)
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)
Namespace for MueLu classes and methods.
void Build(Level &fineLevel, Level &coarseLevel) const
Build an object with this factory.
int GetLevelID() const
Return level number.
Definition: MueLu_Level.cpp:76
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 Teuchos::ArrayRCP< const bool > DetectDirichletRowsExt(const Xpetra::Matrix< Scalar, LocalOrdinal, GlobalOrdinal, Node > &A, bool &bHasZeroDiagonal, const Magnitude &tol=Teuchos::ScalarTraits< Scalar >::zero())
void DeclareInput(Level &fineLevel, Level &coarseLevel) const
Input.
Exception throws to report errors in the internal logical of the program.