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;
111  Teuchos::ArrayRCP<const bool> dirOrNot = MueLu::Utilities<Scalar, LocalOrdinal, GlobalOrdinal, Node>::DetectDirichletRowsExt(*unamalgA,bHasZeroDiagonal,STS::magnitude(0.5));
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(),
208  Teuchos::OrdinalTraits<Xpetra::global_size_t>::invalid(),
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(),
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(),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.
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
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.