Ifpack2 Templated Preconditioning Package  Version 1.0
Ifpack2_LinePartitioner_def.hpp
1 /*@HEADER
2 // ***********************************************************************
3 //
4 // Ifpack2: Tempated 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 
49 namespace Ifpack2 {
50 
51 template<class T>
52 inline 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
58 template<class GraphType,class Scalar>
60 LinePartitioner (const Teuchos::RCP<const row_graph_type>& graph) :
61  OverlappingPartitioner<GraphType> (graph), NumEqns_(1), threshold_(Teuchos::ScalarTraits<MT>::zero())
62 {}
63 
64 
65 template<class GraphType,class Scalar>
67 
68 
69 template<class GraphType,class Scalar>
70 void
72 setPartitionParameters(Teuchos::ParameterList& List) {
73  threshold_ = List.get("partitioner: line detection threshold",threshold_);
74  TEUCHOS_TEST_FOR_EXCEPTION(threshold_ < Teuchos::ScalarTraits<MT>::zero() || threshold_ > Teuchos::ScalarTraits<MT>::one(),
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 
85 template<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_->getNodeNumRows(),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_->getNodeNumRows(); 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 // ============================================================================
109 template<class GraphType,class Scalar>
110 int LinePartitioner<GraphType,Scalar>::Compute_Blocks_AutoLine(Teuchos::ArrayView<local_ordinal_type> blockIndices) const {
111  typedef local_ordinal_type LO;
112  const LO invalid = Teuchos::OrdinalTraits<LO>::invalid();
113  const Scalar zero = Teuchos::ScalarTraits<Scalar>::zero();
114  const MT mzero = Teuchos::ScalarTraits<MT>::zero();
115 
116  Teuchos::ArrayRCP<const Scalar> xvalsRCP, yvalsRCP, zvalsRCP;
117  Teuchos::ArrayView<const Scalar> 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  MT tol = threshold_;
123  size_t N = this->Graph_->getNodeNumRows();
124  size_t allocated_space = this->Graph_->getNodeMaxNumRowEntries();
125 
126  Teuchos::Array<LO> 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  Scalar x0 = (!xvals.is_null()) ? xvals[i/NumEqns_] : zero;
143  Scalar y0 = (!yvals.is_null()) ? yvals[i/NumEqns_] : zero;
144  Scalar 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 = mzero;
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<Scalar>(x0 - xvals[nn]);
152  if(!yvals.is_null()) mydist += square<Scalar>(y0 - yvals[nn]);
153  if(!zvals.is_null()) mydist += square<Scalar>(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 // ============================================================================
180 template<class GraphType,class Scalar>
181 void 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, MT tol, Teuchos::Array<local_ordinal_type> itemp, Teuchos::Array<MT> dtemp) const {
182  typedef local_ordinal_type LO;
183  const LO invalid = Teuchos::OrdinalTraits<LO>::invalid();
184  const Scalar zero = Teuchos::ScalarTraits<Scalar>::zero();
185  const MT mzero = Teuchos::ScalarTraits<MT>::zero();
186 
187  Teuchos::ArrayRCP<const Scalar> xvalsRCP, yvalsRCP, zvalsRCP;
188  Teuchos::ArrayView<const Scalar> 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_->getNodeNumRows();
194  size_t allocated_space = this->Graph_->getNodeMaxNumRowEntries();
195  Teuchos::ArrayView<LO> cols = itemp();
196  Teuchos::ArrayView<LO> indices = itemp.view(allocated_space,allocated_space);
197  Teuchos::ArrayView<MT> dist= dtemp();
198 
199  while (blockIndices[next] == invalid) {
200  // Get the next row
201  size_t nz=0;
202  LO neighbors_in_line=0;
203 
204  this->Graph_->getLocalRowCopy(next,cols(),nz);
205  Scalar x0 = (!xvals.is_null()) ? xvals[next/NumEqns_] : zero;
206  Scalar y0 = (!yvals.is_null()) ? yvals[next/NumEqns_] : zero;
207  Scalar z0 = (!zvals.is_null()) ? zvals[next/NumEqns_] : zero;
208 
209  // Calculate neighbor distances & sort
210  LO neighbor_len=0;
211  for(size_t i=0; i<nz; i+=NumEqns) {
212  MT mydist = mzero;
213  if(cols[i] >=(LO)N) continue; // Check for off-proc entries
214  LO nn = cols[i] / NumEqns;
215  if(blockIndices[nn]==LineID) neighbors_in_line++;
216  if(!xvals.is_null()) mydist += square<Scalar>(x0 - xvals[nn]);
217  if(!yvals.is_null()) mydist += square<Scalar>(y0 - yvals[nn]);
218  if(!zvals.is_null()) mydist += square<Scalar>(z0 - zvals[nn]);
219  dist[neighbor_len] = Teuchos::ScalarTraits<MT>::squareroot(mydist);
220  indices[neighbor_len]=cols[i];
221  neighbor_len++;
222  }
223  // If more than one of my neighbors is already in this line. I
224  // can't be because I'd create a cycle
225  if(neighbors_in_line > 1) break;
226 
227  // Otherwise add me to the line
228  for(LO k=0; k<NumEqns; k++)
229  blockIndices[next + k] = LineID;
230 
231  // Try to find the next guy in the line (only check the closest two that aren't element 0 (diagonal))
232  Teuchos::ArrayView<MT> dist_view = dist(0,neighbor_len);
233  Tpetra::sort2(dist_view.begin(),dist_view.end(),indices.begin());
234 
235  if(neighbor_len > 2 && indices[1] != last && blockIndices[indices[1]] == -1 && dist[1]/dist[neighbor_len-1] < tol) {
236  last=next;
237  next=indices[1];
238  }
239  else if(neighbor_len > 3 && indices[2] != last && blockIndices[indices[2]] == -1 && dist[2]/dist[neighbor_len-1] < tol) {
240  last=next;
241  next=indices[2];
242  }
243  else {
244  // I have no further neighbors in this line
245  break;
246  }
247  }
248 }
249 
250 
251 
252 }// namespace Ifpack2
253 
254 #define IFPACK2_LINEPARTITIONER_INSTANT(S,LO,GO,N) \
255  template class Ifpack2::LinePartitioner<Tpetra::CrsGraph< LO, GO, N >,S >; \
256  template class Ifpack2::LinePartitioner<Tpetra::RowGraph< LO, GO, N >,S >;
257 
258 #endif // IFPACK2_LINEPARTITIONER_DEF_HPP
LinePartitioner(const Teuchos::RCP< const row_graph_type > &graph)
Constructor.
Definition: Ifpack2_LinePartitioner_def.hpp:60
void computePartitions()
Compute the partitions.
Definition: Ifpack2_LinePartitioner_def.hpp:86
virtual ~LinePartitioner()
Destructor.
Definition: Ifpack2_LinePartitioner_def.hpp:66
Definition: Ifpack2_Container.hpp:761
Create overlapping partitions of a local graph.
Definition: Ifpack2_OverlappingPartitioner_decl.hpp:78
void setPartitionParameters(Teuchos::ParameterList &List)
Set the partitioner&#39;s parameters (none for linear partitioning).
Definition: Ifpack2_LinePartitioner_def.hpp:72
Preconditioners and smoothers for Tpetra sparse matrices.
Definition: Ifpack2_AdditiveSchwarz_decl.hpp:72
Ifpack2::LinePartitioner: A class to define partitions into a set of lines.
Definition: Ifpack2_LinePartitioner_decl.hpp:77