MueLu Version of the Day
Loading...
Searching...
No Matches
MueLu_SegregatedAFactory_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// Tobias Wiesner (tawiesn@sandia.gov)
43//
44// ***********************************************************************
45//
46// @HEADER
47#ifndef MUELU_SEGREGATEDAFACTORY_DEF_HPP
48#define MUELU_SEGREGATEDAFACTORY_DEF_HPP
49
50#include <Xpetra_Matrix.hpp>
51#include <Xpetra_MatrixFactory.hpp>
52
54
55#include "MueLu_FactoryManager.hpp"
56#include "MueLu_Level.hpp"
57//#include "MueLu_MasterList.hpp"
58#include "MueLu_Monitor.hpp"
59
60namespace MueLu {
61
62 template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
64 RCP<ParameterList> validParamList = rcp(new ParameterList());
65
66#define SET_VALID_ENTRY(name) validParamList->setEntry(name, MasterList::getEntry(name))
67#undef SET_VALID_ENTRY
68
69 validParamList->set< RCP<const FactoryBase> >("A", Teuchos::null, "Generating factory of the matrix A used for filtering");
70
71 validParamList->set< std::string > ("map: name", "", "Name of map (Xpetra::Map) provided by user containing the special DOFs.");
72 validParamList->set< std::string > ("map: factory", "null", "Name of generating factory for 'map: name'");
73
74 return validParamList;
75 }
76
77 template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
79 Input(currentLevel, "A");
80
81 const ParameterList& pL = GetParameterList();
82 std::string mapName = pL.get<std::string> ("map: name");
83 std::string mapFact = pL.get<std::string> ("map: factory");
84 if(mapFact == "NoFactory") {
85 currentLevel.DeclareInput(mapName,MueLu::NoFactory::get(),this);
86 } else if (mapFact != "null") {
87 Teuchos::RCP<const FactoryBase> fact = currentLevel.GetFactoryManager()->GetFactory(mapFact);
88 currentLevel.DeclareInput(mapName,fact.get(),this);
89 } else {
90 currentLevel.DeclareInput(mapName,NULL,this);
91 }
92 }
93
94 template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
96 FactoryMonitor m(*this, "Matrix filtering (segregation)", currentLevel);
97
98
99 RCP<Matrix> Ain = Get< RCP<Matrix> >(currentLevel, "A");
100
101 const ParameterList& pL = GetParameterList();
102 std::string mapName = pL.get<std::string> ("map: name");
103 std::string mapFact = pL.get<std::string> ("map: factory");
104
105 Teuchos::RCP<const FactoryBase> mapFactory = Teuchos::null;
106 if(mapFact == "NoFactory") {
107 mapFactory = MueLu::NoFactory::getRCP();
108 } else if (mapFact != "null") {
109 mapFactory = currentLevel.GetFactoryManager()->GetFactory(mapFact);
110 }
111
112 Teuchos::RCP<Map> map = currentLevel.Get< Teuchos::RCP<Map> >(mapName,mapFactory.get());
113
114 // create new empty Operator
115 Teuchos::RCP<Matrix> Aout = MatrixFactory::Build(Ain->getRowMap(), Ain->getGlobalMaxNumRowEntries());
116
117 size_t numLocalRows = Ain->getLocalNumRows();
118 for(size_t row=0; row<numLocalRows; row++) { // how can i replace this by a parallel for?
119 GlobalOrdinal grid = Ain->getRowMap()->getGlobalElement(row); // global row id
120 bool isInMap = map->isNodeGlobalElement(grid);
121
122 // extract row information from input matrix
123 Teuchos::ArrayView<const LocalOrdinal> indices;
124 Teuchos::ArrayView<const Scalar> vals;
125 Ain->getLocalRowView(row, indices, vals);
126
127 // just copy all values in output
128 Teuchos::ArrayRCP<GlobalOrdinal> indout(indices.size(),Teuchos::ScalarTraits<GlobalOrdinal>::zero());
129 Teuchos::ArrayRCP<Scalar> valout(indices.size(),Teuchos::ScalarTraits<Scalar>::zero());
130
131 size_t nNonzeros = 0;
132 for(size_t i=0; i<(size_t)indices.size(); i++) { // or can be parallelize this loop?
133 GlobalOrdinal gcid = Ain->getColMap()->getGlobalElement(indices[i]); // global column id
134 bool isInMap2 = map->isNodeGlobalElement(gcid);
135
136 if (isInMap == isInMap2) {
137 indout [nNonzeros] = gcid;
138 valout [nNonzeros] = vals[i];
139 nNonzeros++;
140 }
141 }
142 indout.resize(nNonzeros);
143 valout.resize(nNonzeros);
144
145 Aout->insertGlobalValues(Ain->getRowMap()->getGlobalElement(row), indout.view(0,indout.size()), valout.view(0,valout.size()));
146 }
147
148 Aout->fillComplete(Ain->getDomainMap(), Ain->getRangeMap());
149
150 // copy block size information
151 Aout->SetFixedBlockSize(Ain->GetFixedBlockSize());
152
153 GetOStream(Statistics0, 0) << "Nonzeros in A (input): " << Ain->getGlobalNumEntries() << ", Nonzeros after filtering A: " << Aout->getGlobalNumEntries() << std::endl;
154
155 Set(currentLevel, "A", Aout);
156 }
157
158} //namespace MueLu
159
160#endif // MUELU_SEGREGATEDAFACTORY_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.
void DeclareInput(const std::string &ename, const FactoryBase *factory, const FactoryBase *requestedBy=NoFactory::get())
Callback from FactoryBase::CallDeclareInput() and FactoryBase::DeclareInput()
const RCP< const FactoryManagerBase > GetFactoryManager()
returns the current factory manager
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 RCP< const NoFactory > getRCP()
Static Get() functions.
static const NoFactory * get()
void DeclareInput(Level &currentLevel) const
Input.
RCP< const ParameterList > GetValidParameterList() const
Return a const parameter list of valid parameters that setParameterList() will accept.
void Build(Level &currentLevel) const
Build method.
Namespace for MueLu classes and methods.
@ Statistics0
Print statistics that do not involve significant additional computation.