Ifpack2 Templated Preconditioning Package Version 1.0
Loading...
Searching...
No Matches
Ifpack2_CreateOverlapGraph.hpp
1/*@HEADER
2// ***********************************************************************
3//
4// Ifpack2: Templated Object-Oriented Algebraic Preconditioner Package
5// Copyright (2009) Sandia Corporation
6//
7// Under terms of Contract DE-AC04-94AL85000, there is a non-exclusive
8// license for use of this work by or on behalf of the U.S. Government.
9//
10// Redistribution and use in source and binary forms, with or without
11// modification, are permitted provided that the following conditions are
12// met:
13//
14// 1. Redistributions of source code must retain the above copyright
15// notice, this list of conditions and the following disclaimer.
16//
17// 2. Redistributions in binary form must reproduce the above copyright
18// notice, this list of conditions and the following disclaimer in the
19// documentation and/or other materials provided with the distribution.
20//
21// 3. Neither the name of the Corporation nor the names of the
22// contributors may be used to endorse or promote products derived from
23// this software without specific prior written permission.
24//
25// THIS SOFTWARE IS PROVIDED BY SANDIA CORPORATION "AS IS" AND ANY
26// EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
27// IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR
28// PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL SANDIA CORPORATION OR THE
29// CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL,
30// EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO,
31// PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR
32// PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF
33// LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING
34// NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
35// SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
36//
37// Questions? Contact Michael A. Heroux (maherou@sandia.gov)
38//
39// ***********************************************************************
40//@HEADER
41*/
42
43#ifndef IFPACK2_CREATEOVERLAPGRAPH_HPP
44#define IFPACK2_CREATEOVERLAPGRAPH_HPP
45
46#include "Ifpack2_ConfigDefs.hpp"
47#include "Tpetra_CrsGraph.hpp"
48#include "Tpetra_CrsMatrix.hpp"
49#include "Tpetra_Import.hpp"
50#include "Teuchos_RefCountPtr.hpp"
51
52
53namespace Ifpack2 {
54
71template<class GraphType>
72Teuchos::RCP<const GraphType>
73createOverlapGraph (const Teuchos::RCP<const GraphType>& inputGraph,
74 const int overlapLevel)
75{
76 using Teuchos::RCP;
77 using Teuchos::rcp;
78 typedef Tpetra::Map<typename GraphType::local_ordinal_type,
79 typename GraphType::global_ordinal_type,
80 typename GraphType::node_type> map_type;
81 typedef Tpetra::Import<typename GraphType::local_ordinal_type,
82 typename GraphType::global_ordinal_type,
83 typename GraphType::node_type> import_type;
84 TEUCHOS_TEST_FOR_EXCEPTION(
85 overlapLevel < 0, std::invalid_argument,
86 "Ifpack2::createOverlapGraph: overlapLevel must be >= 0, "
87 "but you specified overlapLevel = " << overlapLevel << ".");
88
89 const int numProcs = inputGraph->getMap ()->getComm ()->getSize ();
90 if (overlapLevel == 0 || numProcs < 2) {
91 return inputGraph;
92 }
93
94 RCP<const map_type> overlapRowMap = inputGraph->getRowMap ();
95 RCP<const map_type> domainMap = inputGraph->getDomainMap ();
96 RCP<const map_type> rangeMap = inputGraph->getRangeMap ();
97
98 RCP<GraphType> overlapGraph;
99 RCP<const GraphType> oldGraph;
100 RCP<const map_type> oldRowMap;
101 for (int level = 0; level < overlapLevel; ++level) {
102 oldGraph = overlapGraph;
103 oldRowMap = overlapRowMap;
104
105 RCP<const import_type> overlapImporter;
106 if (level == 0) {
107 overlapImporter = inputGraph->getImporter ();
108 } else {
109 overlapImporter = oldGraph->getImporter ();
110 }
111
112 overlapRowMap = overlapImporter->getTargetMap ();
113 if (level < overlapLevel - 1) {
114 overlapGraph = rcp (new GraphType (overlapRowMap, 0));
115 }
116 else {
117 // On last iteration, we want to filter out all columns except those that
118 // correspond to rows in the graph. This ensures that our graph is square
119 overlapGraph = rcp (new GraphType (overlapRowMap, overlapRowMap, 0));
120 }
121
122 overlapGraph->doImport (*inputGraph, *overlapImporter, Tpetra::INSERT);
123 overlapGraph->fillComplete (domainMap, rangeMap);
124 }
125
126 return overlapGraph;
127}
128
138template<class MatrixType>
139Teuchos::RCP<const MatrixType>
140createOverlapMatrix (const Teuchos::RCP<const MatrixType>& inputMatrix,
141 const int overlapLevel)
142{
143 using Teuchos::RCP;
144 using Teuchos::rcp;
145 typedef typename MatrixType::map_type map_type;
146 typedef Tpetra::Import<typename MatrixType::local_ordinal_type,
147 typename MatrixType::global_ordinal_type,
148 typename MatrixType::node_type> import_type;
149
150 TEUCHOS_TEST_FOR_EXCEPTION(
151 overlapLevel < 0, std::invalid_argument,
152 "Ifpack2::createOverlapMatrix: overlapLevel must be >= 0, "
153 "but you specified overlapLevel = " << overlapLevel << ".");
154
155 const int numProcs = inputMatrix->getMap ()->getComm ()->getSize ();
156 if (overlapLevel == 0 || numProcs < 2) {
157 return inputMatrix;
158 }
159
160 RCP<const map_type> overlapRowMap = inputMatrix->getRowMap ();
161 RCP<const map_type> domainMap = inputMatrix->getDomainMap ();
162 RCP<const map_type> rangeMap = inputMatrix->getRangeMap ();
163
164 RCP<MatrixType> overlapMatrix;
165 RCP<const MatrixType> oldMatrix;
166 RCP<const map_type> oldRowMap;
167 for (int level = 0; level < overlapLevel; ++level) {
168 oldMatrix = overlapMatrix;
169 oldRowMap = overlapRowMap;
170
171 RCP<const import_type> overlapImporter;
172 if (level == 0) {
173 overlapImporter = inputMatrix->getGraph ()->getImporter ();
174 } else {
175 overlapImporter = oldMatrix->getGraph ()->getImporter ();
176 }
177
178 overlapRowMap = overlapImporter->getTargetMap ();
179 if (level < overlapLevel - 1) {
180 overlapMatrix = rcp (new MatrixType (overlapRowMap, 0));
181 }
182 else {
183 // On last iteration, we want to filter out all columns except those that
184 // correspond to rows in the matrix. This assures that our matrix is square
185 overlapMatrix = rcp (new MatrixType (overlapRowMap, overlapRowMap, 0));
186 }
187
188 overlapMatrix->doImport (*inputMatrix, *overlapImporter, Tpetra::INSERT);
189 overlapMatrix->fillComplete (domainMap, rangeMap);
190 }
191
192 return overlapMatrix;
193}
194
195} // namespace Ifpack2
196
197#endif // IFPACK2_CREATEOVERLAPGRAPH_HPP
Preconditioners and smoothers for Tpetra sparse matrices.
Definition: Ifpack2_AdditiveSchwarz_decl.hpp:74
Teuchos::RCP< const MatrixType > createOverlapMatrix(const Teuchos::RCP< const MatrixType > &inputMatrix, const int overlapLevel)
Construct an overlapped matrix for use with Ifpack2 preconditioners.
Definition: Ifpack2_CreateOverlapGraph.hpp:140
Teuchos::RCP< const GraphType > createOverlapGraph(const Teuchos::RCP< const GraphType > &inputGraph, const int overlapLevel)
Construct an overlapped graph for use with Ifpack2 preconditioners.
Definition: Ifpack2_CreateOverlapGraph.hpp:73