MueLu Version of the Day
Loading...
Searching...
No Matches
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
54namespace 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()->getLocalNumElements() /*amalgP->getRowMap()->getLocalNumElements() * 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->getLocalNumRows());
123 Teuchos::ArrayRCP<const LocalOrdinal> amalgCols(amalgP->getLocalNumEntries());
124 Teuchos::ArrayRCP<const Scalar> amalgVals(amalgP->getLocalNumEntries());
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()->getLocalNumElements() * 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->getLocalNumEntries() * maxDofPerNode);
135 Teuchos::ArrayRCP<Scalar> newPVals(amalgP->getLocalNumEntries() * 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()->getLocalNumElements() * 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()->getLocalNumElements() * maxDofPerNode);
216 Teuchos::Array<GlobalOrdinal> unsmooshColMapGIDs(nColCoarseDofs);
217 for(size_t c = 0; c < amalgP->getColMap()->getLocalNumElements(); ++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->getLocalMaxNumRowEntries());
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_ */
MueLu::DefaultGlobalOrdinal GlobalOrdinal
Exception throws to report errors in the internal logical of the program.
Timer to be used in factories. Similar to Monitor but with additional timers.
Class that holds all level-specific information.
int GetLevelID() const
Return level number.
void Build(Level &fineLevel, Level &coarseLevel) const
Build an object with this factory.
void DeclareInput(Level &fineLevel, Level &coarseLevel) const
Input.
RCP< const ParameterList > GetValidParameterList() const
Return a const parameter list of valid parameters that setParameterList() will accept.
static Teuchos::ArrayRCP< const bool > DetectDirichletRowsExt(const Xpetra::Matrix< Scalar, LocalOrdinal, GlobalOrdinal, Node > &A, bool &bHasZeroDiagonal, const Magnitude &tol=Teuchos::ScalarTraits< Scalar >::zero())
Namespace for MueLu classes and methods.