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  } else {
106  // dof status is the dirichlet information of unsmooshed/unamalgamated A (fine level)
107  dofStatus = Teuchos::Array<char>(unamalgA->getRowMap()->getNodeNumElements() /*amalgP->getRowMap()->getNodeNumElements() * maxDofPerNode*/,'s');
108 
109  bool bHasZeroDiagonal = false;
110  Teuchos::ArrayRCP<const bool> dirOrNot = MueLu::Utilities<Scalar, LocalOrdinal, GlobalOrdinal, Node>::DetectDirichletRowsExt(*unamalgA,bHasZeroDiagonal,STS::magnitude(0.5));
111 
112  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());
113  for(decltype(dirOrNot.size()) i = 0; i < dirOrNot.size(); ++i) {
114  if(dirOrNot[i] == true) dofStatus[i] = 'p';
115  }
116  }
117 
118  // TODO: TAW the following check is invalid for SA-AMG based input prolongators
119  //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)");
120 
121  // extract CRS information from amalgamated prolongation operator
122  Teuchos::ArrayRCP<const size_t> amalgRowPtr(amalgP->getNodeNumRows());
123  Teuchos::ArrayRCP<const LocalOrdinal> amalgCols(amalgP->getNodeNumEntries());
124  Teuchos::ArrayRCP<const Scalar> amalgVals(amalgP->getNodeNumEntries());
125  Teuchos::RCP<CrsMatrixWrap> amalgPwrap = Teuchos::rcp_dynamic_cast<CrsMatrixWrap>(amalgP);
126  Teuchos::RCP<CrsMatrix> amalgPcrs = amalgPwrap->getCrsMatrix();
127  amalgPcrs->getAllValues(amalgRowPtr, amalgCols, amalgVals);
128 
129  // calculate number of dof rows for new prolongator
130  size_t paddedNrows = amalgP->getRowMap()->getNodeNumElements() * Teuchos::as<size_t>(maxDofPerNode);
131 
132  // reserve CSR arrays for new prolongation operator
133  Teuchos::ArrayRCP<size_t> newPRowPtr(paddedNrows+1);
134  Teuchos::ArrayRCP<LocalOrdinal> newPCols(amalgP->getNodeNumEntries() * maxDofPerNode);
135  Teuchos::ArrayRCP<Scalar> newPVals(amalgP->getNodeNumEntries() * maxDofPerNode);
136 
137  size_t rowCount = 0; // actual number of (local) in unamalgamated prolongator
138  if(fineIsPadded == true || fineLevel.GetLevelID() > 0) {
139 
140  // build prolongation operator for padded fine level matrices.
141  // Note: padded fine level dofs are transferred by injection.
142  // That is, these interpolation stencils do not take averages of
143  // coarse level variables. Further, fine level Dirichlet points
144  // also use injection.
145 
146  size_t cnt = 0; // local id counter
147  for (decltype(amalgRowPtr.size()) i = 0; i < amalgRowPtr.size() - 1; i++) {
148  // determine number of entries in amalgamated dof row i
149  size_t rowLength = amalgRowPtr[i+1] - amalgRowPtr[i];
150 
151  // loop over dofs per node (unamalgamation)
152  for(int j = 0; j < maxDofPerNode; j++) {
153  newPRowPtr[i*maxDofPerNode+j] = cnt;
154  if (dofStatus[i*maxDofPerNode+j] == 's') { // add only "standard" dofs to unamalgamated prolongator
155  // loop over column entries in amalgamated P
156  for (size_t k = 0; k < rowLength; k++) {
157  newPCols[cnt ] = amalgCols[k+amalgRowPtr[i]] * maxDofPerNode + j;
158  newPVals[cnt++] = amalgVals[k+amalgRowPtr[i]];
159  }
160 
161  }
162  }
163  }
164 
165  newPRowPtr[paddedNrows] = cnt; // close row CSR array
166  rowCount = paddedNrows;
167  } else {
168  // Build prolongation operator for non-padded fine level matrices.
169  // Need to map from non-padded dofs to padded dofs. For this, look
170  // at the status array and skip padded dofs.
171 
172  size_t cnt = 0; // local id counter
173 
174  for (decltype(amalgRowPtr.size()) i = 0; i < amalgRowPtr.size() - 1; i++) {
175  // determine number of entries in amalgamated dof row i
176  size_t rowLength = amalgRowPtr[i+1] - amalgRowPtr[i];
177 
178  // loop over dofs per node (unamalgamation)
179  for(int j = 0; j < maxDofPerNode; j++) {
180  // no interpolation for padded fine dofs as they do not exist
181 
182  if (dofStatus[i*maxDofPerNode+j] == 's') { // add only "standard" dofs to unamalgamated prolongator
183  newPRowPtr[rowCount++] = cnt;
184  // loop over column entries in amalgamated P
185  for (size_t k = 0; k < rowLength; k++) {
186  newPCols[cnt ] = amalgCols[k+amalgRowPtr[i]] * maxDofPerNode + j;
187  newPVals[cnt++] = amalgVals[k+amalgRowPtr[i]];
188  }
189 
190  }
191  if (dofStatus[i*maxDofPerNode+j] == 'd') { // Dirichlet handling
192  newPRowPtr[rowCount++] = cnt;
193  }
194  }
195  }
196  newPRowPtr[rowCount] = cnt; // close row CSR array
197  } // fineIsPadded == false
198 
199  // generate coarse domain map
200  // So far no support for gid offset or strided maps. This information
201  // could be gathered easily from the unamalgamated fine level operator A.
202  std::vector<size_t> stridingInfo(1, maxDofPerNode);
203 
204  GlobalOrdinal nCoarseDofs = amalgP->getDomainMap()->getNodeNumElements() * maxDofPerNode;
205  GlobalOrdinal indexBase = amalgP->getDomainMap()->getIndexBase();
206  RCP<const Map> coarseDomainMap = StridedMapFactory::Build(amalgP->getDomainMap()->lib(),
207  Teuchos::OrdinalTraits<Xpetra::global_size_t>::invalid(),
208  nCoarseDofs,
209  indexBase,
210  stridingInfo,
211  amalgP->getDomainMap()->getComm(),
212  -1 /* stridedBlockId */,
213  0 /*domainGidOffset */);
214 
215  size_t nColCoarseDofs = Teuchos::as<size_t>(amalgP->getColMap()->getNodeNumElements() * maxDofPerNode);
216  Teuchos::Array<GlobalOrdinal> unsmooshColMapGIDs(nColCoarseDofs);
217  for(size_t c = 0; c < amalgP->getColMap()->getNodeNumElements(); ++c) {
218  GlobalOrdinal gid = (amalgP->getColMap()->getGlobalElement(c)-indexBase) * maxDofPerNode + indexBase;
219 
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(),
225  Teuchos::OrdinalTraits<Xpetra::global_size_t>::invalid(),
226  unsmooshColMapGIDs(), //View,
227  indexBase,
228  amalgP->getDomainMap()->getComm());
229 
230  // Assemble unamalgamated P
231  Teuchos::RCP<CrsMatrix> unamalgPCrs = CrsMatrixFactory::Build(unamalgA->getRowMap(),
232  coarseColMap,
233  maxDofPerNode*amalgP->getNodeMaxNumRowEntries());
234  for (size_t i = 0; i < rowCount; i++) {
235  unamalgPCrs->insertLocalValues(i,
236  newPCols.view(newPRowPtr[i], newPRowPtr[i+1] - newPRowPtr[i]),
237  newPVals.view(newPRowPtr[i], newPRowPtr[i+1] - newPRowPtr[i]));
238  }
239  unamalgPCrs->fillComplete(coarseDomainMap, unamalgA->getRowMap());
240 
241  Teuchos::RCP<Matrix> unamalgP = Teuchos::rcp(new CrsMatrixWrap(unamalgPCrs));
242 
243  Set(coarseLevel,"P",unamalgP);
244  }
245 
246 
247 } /* MueLu */
248 
249 
250 #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.
Timer to be used in factories. Similar to Monitor but with additional timers.
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
MueLu::DefaultGlobalOrdinal GlobalOrdinal
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.