Ifpack2 Templated Preconditioning Package Version 1.0
Loading...
Searching...
No Matches
Ifpack2_LinePartitioner_def.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_LINE_PARTITIONER_DEF_HPP
44#define IFPACK2_LINE_PARTITIONER_DEF_HPP
45
46#include "Tpetra_CrsGraph.hpp"
47#include "Tpetra_Util.hpp"
48
49namespace Ifpack2 {
50
51template<class T>
52inline typename Teuchos::ScalarTraits<T>::magnitudeType square(T x) {
53 return Teuchos::ScalarTraits<T>::magnitude(x) * Teuchos::ScalarTraits<T>::magnitude(x);
54}
55
56//==============================================================================
57// Constructor
58template<class GraphType,class Scalar>
60LinePartitioner (const Teuchos::RCP<const row_graph_type>& graph) :
61 OverlappingPartitioner<GraphType> (graph), NumEqns_(1), threshold_(0.0)
62{}
63
64
65template<class GraphType,class Scalar>
67
68
69template<class GraphType,class Scalar>
70void
72setPartitionParameters(Teuchos::ParameterList& List) {
73 threshold_ = List.get("partitioner: line detection threshold",threshold_);
74 TEUCHOS_TEST_FOR_EXCEPTION(threshold_ < 0.0 || threshold_ > 1.0,
75 std::runtime_error,"Ifpack2::LinePartitioner: threshold not valid");
76
77 NumEqns_ = List.get("partitioner: PDE equations",NumEqns_);
78 TEUCHOS_TEST_FOR_EXCEPTION(NumEqns_<1,std::runtime_error,"Ifpack2::LinePartitioner: NumEqns not valid");
79
80 coord_ = List.get("partitioner: coordinates",coord_);
81 TEUCHOS_TEST_FOR_EXCEPTION(coord_.is_null(),std::runtime_error,"Ifpack2::LinePartitioner: coordinates not defined");
82}
83
84
85template<class GraphType,class Scalar>
87 const local_ordinal_type invalid = Teuchos::OrdinalTraits<local_ordinal_type>::invalid();
88
89 // Sanity Checks
90 TEUCHOS_TEST_FOR_EXCEPTION(coord_.is_null(),std::runtime_error,"Ifpack2::LinePartitioner: coordinates not defined");
91 TEUCHOS_TEST_FOR_EXCEPTION((size_t)this->Partition_.size() != this->Graph_->getLocalNumRows(),std::runtime_error,"Ifpack2::LinePartitioner: partition size error");
92
93 // Short circuit
94 if(this->Partition_.size() == 0) {this->NumLocalParts_ = 0; return;}
95
96 // Set partitions to invalid to initialize algorithm
97 for(size_t i=0; i<this->Graph_->getLocalNumRows(); i++)
98 this->Partition_[i] = invalid;
99
100 // Use the auto partitioner
101 this->NumLocalParts_ = this->Compute_Blocks_AutoLine(this->Partition_());
102
103 // Resize Parts_
104 this->Parts_.resize(this->NumLocalParts_);
105}
106
107
108// ============================================================================
109template<class GraphType,class Scalar>
110int LinePartitioner<GraphType,Scalar>::Compute_Blocks_AutoLine(Teuchos::ArrayView<local_ordinal_type> blockIndices) const {
111 typedef local_ordinal_type LO;
112 typedef magnitude_type MT;
113 const LO invalid = Teuchos::OrdinalTraits<LO>::invalid();
114 const MT zero = Teuchos::ScalarTraits<MT>::zero();
115
116 Teuchos::ArrayRCP<const MT> xvalsRCP, yvalsRCP, zvalsRCP;
117 Teuchos::ArrayView<const MT> xvals, yvals, zvals;
118 xvalsRCP = coord_->getData(0); xvals = xvalsRCP();
119 if(coord_->getNumVectors() > 1) { yvalsRCP = coord_->getData(1); yvals = yvalsRCP(); }
120 if(coord_->getNumVectors() > 2) { zvalsRCP = coord_->getData(2); zvals = zvalsRCP(); }
121
122 double tol = threshold_;
123 size_t N = this->Graph_->getLocalNumRows();
124 size_t allocated_space = this->Graph_->getLocalMaxNumRowEntries();
125
126 nonconst_local_inds_host_view_type cols("cols",allocated_space);
127 Teuchos::Array<LO> indices(allocated_space);
128 Teuchos::Array<MT> dist(allocated_space);
129
130 Teuchos::Array<LO> itemp(2*allocated_space);
131 Teuchos::Array<MT> dtemp(allocated_space);
132
133 LO num_lines = 0;
134
135 for(LO i=0; i<(LO)N; i+=NumEqns_) {
136 size_t nz=0;
137 // Short circuit if I've already been blocked
138 if(blockIndices[i] != invalid) continue;
139
140 // Get neighbors and sort by distance
141 this->Graph_->getLocalRowCopy(i,cols,nz);
142 MT x0 = (!xvals.is_null()) ? xvals[i/NumEqns_] : zero;
143 MT y0 = (!yvals.is_null()) ? yvals[i/NumEqns_] : zero;
144 MT z0 = (!zvals.is_null()) ? zvals[i/NumEqns_] : zero;
145
146 LO neighbor_len=0;
147 for(size_t j=0; j<nz; j+=NumEqns_) {
148 MT mydist = zero;
149 LO nn = cols[j] / NumEqns_;
150 if(cols[j] >=(LO)N) continue; // Check for off-proc entries
151 if(!xvals.is_null()) mydist += square<MT>(x0 - xvals[nn]);
152 if(!yvals.is_null()) mydist += square<MT>(y0 - yvals[nn]);
153 if(!zvals.is_null()) mydist += square<MT>(z0 - zvals[nn]);
154 dist[neighbor_len] = Teuchos::ScalarTraits<MT>::squareroot(mydist);
155 indices[neighbor_len]=cols[j];
156 neighbor_len++;
157 }
158
159 Teuchos::ArrayView<MT> dist_view = dist(0,neighbor_len);
160 Tpetra::sort2(dist_view.begin(),dist_view.end(),indices.begin());
161
162 // Number myself
163 for(LO k=0; k<NumEqns_; k++)
164 blockIndices[i + k] = num_lines;
165
166 // Fire off a neighbor line search (nearest neighbor)
167 if(neighbor_len > 2 && dist[1]/dist[neighbor_len-1] < tol) {
168 local_automatic_line_search(NumEqns_,blockIndices,i,indices[1],num_lines,tol,itemp,dtemp);
169 }
170 // Fire off a neighbor line search (second nearest neighbor)
171 if(neighbor_len > 3 && dist[2]/dist[neighbor_len-1] < tol) {
172 local_automatic_line_search(NumEqns_,blockIndices,i,indices[2],num_lines,tol,itemp,dtemp);
173 }
174
175 num_lines++;
176 }
177 return num_lines;
178}
179// ============================================================================
180template<class GraphType,class Scalar>
181void LinePartitioner<GraphType,Scalar>::local_automatic_line_search(int NumEqns, Teuchos::ArrayView <local_ordinal_type> blockIndices, local_ordinal_type last, local_ordinal_type next, local_ordinal_type LineID, double tol, Teuchos::Array<local_ordinal_type> itemp, Teuchos::Array<magnitude_type> dtemp) const {
182 typedef local_ordinal_type LO;
183 typedef magnitude_type MT;
184 const LO invalid = Teuchos::OrdinalTraits<LO>::invalid();
185 const MT zero = Teuchos::ScalarTraits<MT>::zero();
186
187 Teuchos::ArrayRCP<const MT> xvalsRCP, yvalsRCP, zvalsRCP;
188 Teuchos::ArrayView<const MT> xvals, yvals, zvals;
189 xvalsRCP = coord_->getData(0); xvals = xvalsRCP();
190 if(coord_->getNumVectors() > 1) { yvalsRCP = coord_->getData(1); yvals = yvalsRCP(); }
191 if(coord_->getNumVectors() > 2) { zvalsRCP = coord_->getData(2); zvals = zvalsRCP(); }
192
193 size_t N = this->Graph_->getLocalNumRows();
194 size_t allocated_space = this->Graph_->getLocalMaxNumRowEntries();
195
196 nonconst_local_inds_host_view_type cols(itemp.data(),allocated_space);
197 Teuchos::ArrayView<LO> indices = itemp.view(allocated_space,allocated_space);
198 Teuchos::ArrayView<MT> dist= dtemp();
199
200 while (blockIndices[next] == invalid) {
201 // Get the next row
202 size_t nz=0;
203 LO neighbors_in_line=0;
204
205 this->Graph_->getLocalRowCopy(next,cols,nz);
206 MT x0 = (!xvals.is_null()) ? xvals[next/NumEqns_] : zero;
207 MT y0 = (!yvals.is_null()) ? yvals[next/NumEqns_] : zero;
208 MT z0 = (!zvals.is_null()) ? zvals[next/NumEqns_] : zero;
209
210 // Calculate neighbor distances & sort
211 LO neighbor_len=0;
212 for(size_t i=0; i<nz; i+=NumEqns) {
213 MT mydist = zero;
214 if(cols[i] >=(LO)N) continue; // Check for off-proc entries
215 LO nn = cols[i] / NumEqns;
216 if(blockIndices[nn]==LineID) neighbors_in_line++;
217 if(!xvals.is_null()) mydist += square<MT>(x0 - xvals[nn]);
218 if(!yvals.is_null()) mydist += square<MT>(y0 - yvals[nn]);
219 if(!zvals.is_null()) mydist += square<MT>(z0 - zvals[nn]);
220 dist[neighbor_len] = Teuchos::ScalarTraits<MT>::squareroot(mydist);
221 indices[neighbor_len]=cols[i];
222 neighbor_len++;
223 }
224 // If more than one of my neighbors is already in this line. I
225 // can't be because I'd create a cycle
226 if(neighbors_in_line > 1) break;
227
228 // Otherwise add me to the line
229 for(LO k=0; k<NumEqns; k++)
230 blockIndices[next + k] = LineID;
231
232 // Try to find the next guy in the line (only check the closest two that aren't element 0 (diagonal))
233 Teuchos::ArrayView<MT> dist_view = dist(0,neighbor_len);
234 Tpetra::sort2(dist_view.begin(),dist_view.end(),indices.begin());
235
236 if(neighbor_len > 2 && indices[1] != last && blockIndices[indices[1]] == -1 && dist[1]/dist[neighbor_len-1] < tol) {
237 last=next;
238 next=indices[1];
239 }
240 else if(neighbor_len > 3 && indices[2] != last && blockIndices[indices[2]] == -1 && dist[2]/dist[neighbor_len-1] < tol) {
241 last=next;
242 next=indices[2];
243 }
244 else {
245 // I have no further neighbors in this line
246 break;
247 }
248 }
249}
250
251
252
253}// namespace Ifpack2
254
255#define IFPACK2_LINEPARTITIONER_INSTANT(S,LO,GO,N) \
256 template class Ifpack2::LinePartitioner<Tpetra::CrsGraph< LO, GO, N >,S >; \
257 template class Ifpack2::LinePartitioner<Tpetra::RowGraph< LO, GO, N >,S >;
258
259#endif // IFPACK2_LINEPARTITIONER_DEF_HPP
Ifpack2::LinePartitioner: A class to define partitions into a set of lines.
Definition: Ifpack2_LinePartitioner_decl.hpp:78
LinePartitioner(const Teuchos::RCP< const row_graph_type > &graph)
Constructor.
Definition: Ifpack2_LinePartitioner_def.hpp:60
void setPartitionParameters(Teuchos::ParameterList &List)
Set the partitioner's parameters (none for linear partitioning).
Definition: Ifpack2_LinePartitioner_def.hpp:72
virtual ~LinePartitioner()
Destructor.
Definition: Ifpack2_LinePartitioner_def.hpp:66
void computePartitions()
Compute the partitions.
Definition: Ifpack2_LinePartitioner_def.hpp:86
Create overlapping partitions of a local graph.
Definition: Ifpack2_OverlappingPartitioner_decl.hpp:78
Preconditioners and smoothers for Tpetra sparse matrices.
Definition: Ifpack2_AdditiveSchwarz_decl.hpp:74