MueLu Version of the Day
Loading...
Searching...
No Matches
MueLu_ReplicatePFactory_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_REPLICATEPFACTORY_DEF_HPP
47#define MUELU_REPLICATEPFACTORY_DEF_HPP
48
49#include <stdlib.h>
50#include <iomanip>
51
52
53// #include <Teuchos_LAPACK.hpp>
54#include <Teuchos_SerialDenseMatrix.hpp>
55#include <Teuchos_SerialDenseVector.hpp>
56#include <Teuchos_SerialDenseSolver.hpp>
57
58#include <Xpetra_CrsMatrixWrap.hpp>
59#include <Xpetra_ImportFactory.hpp>
60#include <Xpetra_Matrix.hpp>
61#include <Xpetra_MapFactory.hpp>
62#include <Xpetra_MultiVectorFactory.hpp>
63#include <Xpetra_VectorFactory.hpp>
64
65#include <Xpetra_IO.hpp>
66
69
70#include "MueLu_MasterList.hpp"
71#include "MueLu_Monitor.hpp"
72
73namespace MueLu {
74
75 template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
77 RCP<ParameterList> validParamList = rcp(new ParameterList());
78 validParamList->setEntry("replicate: npdes",ParameterEntry(1));
79
80 return validParamList;
81 }
82
83 template <class Scalar,class LocalOrdinal, class GlobalOrdinal, class Node>
85 DeclareInput(Level& fineLevel, Level& /* coarseLevel */) const {
86// Input(fineLevel, "Psubblock");
87 }
88
89 template <class Scalar,class LocalOrdinal, class GlobalOrdinal, class Node>
91 Level& coarseLevel) const {
92 return BuildP(fineLevel, coarseLevel);
93 }
94
95 template <class Scalar,class LocalOrdinal, class GlobalOrdinal, class Node>
97 Level& coarseLevel) const {
98 FactoryMonitor m(*this, "Build", coarseLevel);
99
100 RCP<Matrix> Psubblock = coarseLevel.Get< RCP<Matrix> >("Psubblock", NoFactory::get());
101 const ParameterList& pL = GetParameterList();
102 const LO dofPerNode = as<LO>(pL.get<int> ("replicate: npdes"));
103
104 Teuchos::ArrayRCP<const size_t> amalgRowPtr(Psubblock->getLocalNumRows());
105 Teuchos::ArrayRCP<const LocalOrdinal> amalgCols(Psubblock->getLocalNumEntries());
106 Teuchos::ArrayRCP<const Scalar> amalgVals(Psubblock->getLocalNumEntries());
107 Teuchos::RCP<CrsMatrixWrap> Psubblockwrap = Teuchos::rcp_dynamic_cast<CrsMatrixWrap>(Psubblock);
108 Teuchos::RCP<CrsMatrix> Psubblockcrs = Psubblockwrap->getCrsMatrix();
109 Psubblockcrs->getAllValues(amalgRowPtr, amalgCols, amalgVals);
110
111 size_t paddedNrows = Psubblock->getRowMap()->getLocalNumElements() * Teuchos::as<size_t>(dofPerNode);
112 Teuchos::ArrayRCP<size_t> newPRowPtr(paddedNrows+1);
113 Teuchos::ArrayRCP<LocalOrdinal> newPCols(Psubblock->getLocalNumEntries() * dofPerNode);
114 Teuchos::ArrayRCP<Scalar> newPVals(Psubblock->getLocalNumEntries() * dofPerNode);
115 size_t cnt = 0; // local id counter
116 for (decltype(amalgRowPtr.size()) i = 0; i < amalgRowPtr.size() - 1; i++) {
117 size_t rowLength = amalgRowPtr[i+1] - amalgRowPtr[i];
118 for(int j = 0; j < dofPerNode; j++) {
119 newPRowPtr[i*dofPerNode+j] = cnt;
120 for (size_t k = 0; k < rowLength; k++) {
121 newPCols[cnt ] = amalgCols[k+amalgRowPtr[i]] * dofPerNode + j;
122 newPVals[cnt++] = amalgVals[k+amalgRowPtr[i]];
123 }
124 }
125 }
126
127 newPRowPtr[paddedNrows] = cnt; // close row CSR array
128 std::vector<size_t> stridingInfo(1, dofPerNode);
129
130 GlobalOrdinal nCoarseDofs = Psubblock->getDomainMap()->getLocalNumElements() * dofPerNode;
131 GlobalOrdinal indexBase = Psubblock->getDomainMap()->getIndexBase();
132 RCP<const Map> coarseDomainMap = StridedMapFactory::Build(Psubblock->getDomainMap()->lib(),
133 Teuchos::OrdinalTraits<Xpetra::global_size_t>::invalid(),
134 nCoarseDofs,
135 indexBase,
136 stridingInfo,
137 Psubblock->getDomainMap()->getComm(),
138 -1 /* stridedBlockId */,
139 0 /*domainGidOffset */);
140
141 size_t nColCoarseDofs = Teuchos::as<size_t>(Psubblock->getColMap()->getLocalNumElements() * dofPerNode);
142 Teuchos::Array<GlobalOrdinal> unsmooshColMapGIDs(nColCoarseDofs);
143 for(size_t c = 0; c < Psubblock->getColMap()->getLocalNumElements(); ++c) {
144 GlobalOrdinal gid = (Psubblock->getColMap()->getGlobalElement(c)-indexBase) * dofPerNode + indexBase;
145
146 for(int i = 0; i < dofPerNode; ++i) {
147 unsmooshColMapGIDs[c * dofPerNode + i] = gid + i;
148 }
149 }
150 Teuchos::RCP<Map> coarseColMap = MapFactory::Build(Psubblock->getDomainMap()->lib(),
151 Teuchos::OrdinalTraits<Xpetra::global_size_t>::invalid(),
152 unsmooshColMapGIDs(), //View,
153 indexBase,
154 Psubblock->getDomainMap()->getComm());
155
156 Teuchos::Array<GlobalOrdinal> unsmooshRowMapGIDs(paddedNrows);
157 for(size_t c = 0; c < Psubblock->getRowMap()->getLocalNumElements(); ++c) {
158 GlobalOrdinal gid = (Psubblock->getRowMap()->getGlobalElement(c)-indexBase) * dofPerNode + indexBase;
159
160 for(int i = 0; i < dofPerNode; ++i) {
161 unsmooshRowMapGIDs[c * dofPerNode + i] = gid + i;
162 }
163 }
164 Teuchos::RCP<Map> fineRowMap = MapFactory::Build(Psubblock->getDomainMap()->lib(),
165 Teuchos::OrdinalTraits<Xpetra::global_size_t>::invalid(),
166 unsmooshRowMapGIDs(), //View,
167 indexBase,
168 Psubblock->getDomainMap()->getComm());
169
170 Teuchos::RCP<CrsMatrix> bigPCrs = CrsMatrixFactory::Build(fineRowMap, coarseColMap,
171 dofPerNode*Psubblock->getLocalMaxNumRowEntries());
172 for (size_t i = 0; i < paddedNrows; i++) {
173 bigPCrs->insertLocalValues(i,
174 newPCols.view(newPRowPtr[i], newPRowPtr[i+1] - newPRowPtr[i]),
175 newPVals.view(newPRowPtr[i], newPRowPtr[i+1] - newPRowPtr[i]));
176 }
177 bigPCrs->fillComplete(coarseDomainMap, fineRowMap);
178
179 Teuchos::RCP<Matrix> bigP = Teuchos::rcp(new CrsMatrixWrap(bigPCrs));
180
181 Set(coarseLevel, "P", bigP);
182 }
183
184
185} //namespace MueLu
186
187#define MUELU_REPLICATEPFACTORY_SHORT
188#endif // MUELU_REPLICATEPFACTORY_DEF_HPP
MueLu::DefaultGlobalOrdinal GlobalOrdinal
Timer to be used in factories. Similar to Monitor but with additional timers.
Class that holds all level-specific information.
T & Get(const std::string &ename, const FactoryBase *factory=NoFactory::get())
Get data without decrementing associated storage counter (i.e., read-only access)....
static const NoFactory * get()
void DeclareInput(Level &fineLevel, Level &coarseLevel) const
Input.
void BuildP(Level &fineLevel, Level &coarseLevel) const
Abstract Build method.
void Build(Level &fineLevel, Level &coarseLevel) const
Build an object with this factory.
RCP< const ParameterList > GetValidParameterList() const
Return a const parameter list of valid parameters that setParameterList() will accept.
Namespace for MueLu classes and methods.