43 #ifndef IFPACK2_LINE_PARTITIONER_DEF_HPP 44 #define IFPACK2_LINE_PARTITIONER_DEF_HPP 46 #include "Tpetra_CrsGraph.hpp" 47 #include "Tpetra_Util.hpp" 52 inline typename Teuchos::ScalarTraits<T>::magnitudeType square(T x) {
53 return Teuchos::ScalarTraits<T>::magnitude(x) * Teuchos::ScalarTraits<T>::magnitude(x);
58 template<
class GraphType,
class Scalar>
65 template<
class GraphType,
class Scalar>
69 template<
class GraphType,
class Scalar>
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");
77 NumEqns_ = List.get(
"partitioner: PDE equations",NumEqns_);
78 TEUCHOS_TEST_FOR_EXCEPTION(NumEqns_<1,std::runtime_error,
"Ifpack2::LinePartitioner: NumEqns not valid");
80 coord_ = List.get(
"partitioner: coordinates",coord_);
81 TEUCHOS_TEST_FOR_EXCEPTION(coord_.is_null(),std::runtime_error,
"Ifpack2::LinePartitioner: coordinates not defined");
85 template<
class GraphType,
class Scalar>
87 const local_ordinal_type invalid = Teuchos::OrdinalTraits<local_ordinal_type>::invalid();
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");
94 if(this->Partition_.size() == 0) {this->NumLocalParts_ = 0;
return;}
97 for(
size_t i=0; i<this->Graph_->getNodeNumRows(); i++)
98 this->Partition_[i] = invalid;
101 this->NumLocalParts_ = this->Compute_Blocks_AutoLine(this->Partition_());
104 this->Parts_.resize(this->NumLocalParts_);
109 template<
class GraphType,
class Scalar>
111 typedef local_ordinal_type LO;
112 const LO invalid = Teuchos::OrdinalTraits<LO>::invalid();
113 const double zero = Teuchos::ScalarTraits<double>::zero();
115 Teuchos::ArrayRCP<const double> xvalsRCP, yvalsRCP, zvalsRCP;
116 Teuchos::ArrayView<const double> xvals, yvals, zvals;
117 xvalsRCP = coord_->getData(0); xvals = xvalsRCP();
118 if(coord_->getNumVectors() > 1) { yvalsRCP = coord_->getData(1); yvals = yvalsRCP(); }
119 if(coord_->getNumVectors() > 2) { zvalsRCP = coord_->getData(2); zvals = zvalsRCP(); }
121 double tol = threshold_;
122 size_t N = this->Graph_->getNodeNumRows();
123 size_t allocated_space = this->Graph_->getNodeMaxNumRowEntries();
125 Teuchos::Array<LO> cols(allocated_space);
126 Teuchos::Array<LO> indices(allocated_space);
127 Teuchos::Array<double> dist(allocated_space);
129 Teuchos::Array<LO> itemp(2*allocated_space);
130 Teuchos::Array<double> dtemp(allocated_space);
134 for(LO i=0; i<(LO)N; i+=NumEqns_) {
137 if(blockIndices[i] != invalid)
continue;
140 this->Graph_->getLocalRowCopy(i,cols(),nz);
141 double x0 = (!xvals.is_null()) ? xvals[i/NumEqns_] : zero;
142 double y0 = (!yvals.is_null()) ? yvals[i/NumEqns_] : zero;
143 double z0 = (!zvals.is_null()) ? zvals[i/NumEqns_] : zero;
146 for(
size_t j=0; j<nz; j+=NumEqns_) {
147 double mydist = zero;
148 LO nn = cols[j] / NumEqns_;
149 if(cols[j] >=(LO)N)
continue;
150 if(!xvals.is_null()) mydist += square<double>(x0 - xvals[nn]);
151 if(!yvals.is_null()) mydist += square<double>(y0 - yvals[nn]);
152 if(!zvals.is_null()) mydist += square<double>(z0 - zvals[nn]);
153 dist[neighbor_len] = Teuchos::ScalarTraits<double>::squareroot(mydist);
154 indices[neighbor_len]=cols[j];
158 Teuchos::ArrayView<double> dist_view = dist(0,neighbor_len);
159 Tpetra::sort2(dist_view.begin(),dist_view.end(),indices.begin());
162 for(LO k=0; k<NumEqns_; k++)
163 blockIndices[i + k] = num_lines;
166 if(neighbor_len > 2 && dist[1]/dist[neighbor_len-1] < tol) {
167 local_automatic_line_search(NumEqns_,blockIndices,i,indices[1],num_lines,tol,itemp,dtemp);
170 if(neighbor_len > 3 && dist[2]/dist[neighbor_len-1] < tol) {
171 local_automatic_line_search(NumEqns_,blockIndices,i,indices[2],num_lines,tol,itemp,dtemp);
179 template<
class GraphType,
class Scalar>
180 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,
double tol, Teuchos::Array<local_ordinal_type> itemp, Teuchos::Array<double> dtemp)
const {
181 typedef local_ordinal_type LO;
182 const LO invalid = Teuchos::OrdinalTraits<LO>::invalid();
183 const double zero = Teuchos::ScalarTraits<double>::zero();
185 Teuchos::ArrayRCP<const double> xvalsRCP, yvalsRCP, zvalsRCP;
186 Teuchos::ArrayView<const double> xvals, yvals, zvals;
187 xvalsRCP = coord_->getData(0); xvals = xvalsRCP();
188 if(coord_->getNumVectors() > 1) { yvalsRCP = coord_->getData(1); yvals = yvalsRCP(); }
189 if(coord_->getNumVectors() > 2) { zvalsRCP = coord_->getData(2); zvals = zvalsRCP(); }
191 size_t N = this->Graph_->getNodeNumRows();
192 size_t allocated_space = this->Graph_->getNodeMaxNumRowEntries();
193 Teuchos::ArrayView<LO> cols = itemp();
194 Teuchos::ArrayView<LO> indices = itemp.view(allocated_space,allocated_space);
195 Teuchos::ArrayView<double> dist= dtemp();
197 while (blockIndices[next] == invalid) {
200 LO neighbors_in_line=0;
202 this->Graph_->getLocalRowCopy(next,cols(),nz);
203 double x0 = (!xvals.is_null()) ? xvals[next/NumEqns_] : zero;
204 double y0 = (!yvals.is_null()) ? yvals[next/NumEqns_] : zero;
205 double z0 = (!zvals.is_null()) ? zvals[next/NumEqns_] : zero;
209 for(
size_t i=0; i<nz; i+=NumEqns) {
210 double mydist = zero;
211 if(cols[i] >=(LO)N)
continue;
212 LO nn = cols[i] / NumEqns;
213 if(blockIndices[nn]==LineID) neighbors_in_line++;
214 if(!xvals.is_null()) mydist += square<double>(x0 - xvals[nn]);
215 if(!yvals.is_null()) mydist += square<double>(y0 - yvals[nn]);
216 if(!zvals.is_null()) mydist += square<double>(z0 - zvals[nn]);
217 dist[neighbor_len] = Teuchos::ScalarTraits<double>::squareroot(mydist);
218 indices[neighbor_len]=cols[i];
223 if(neighbors_in_line > 1)
break;
226 for(LO k=0; k<NumEqns; k++)
227 blockIndices[next + k] = LineID;
230 Teuchos::ArrayView<double> dist_view = dist(0,neighbor_len);
231 Tpetra::sort2(dist_view.begin(),dist_view.end(),indices.begin());
233 if(neighbor_len > 2 && indices[1] != last && blockIndices[indices[1]] == -1 && dist[1]/dist[neighbor_len-1] < tol) {
237 else if(neighbor_len > 3 && indices[2] != last && blockIndices[indices[2]] == -1 && dist[2]/dist[neighbor_len-1] < tol) {
252 #define IFPACK2_LINEPARTITIONER_INSTANT(S,LO,GO,N) \ 253 template class Ifpack2::LinePartitioner<Tpetra::CrsGraph< LO, GO, N >,S >; \ 254 template class Ifpack2::LinePartitioner<Tpetra::RowGraph< LO, GO, N >,S >; 256 #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
Create overlapping partitions of a local graph.
Definition: Ifpack2_OverlappingPartitioner_decl.hpp:78
void setPartitionParameters(Teuchos::ParameterList &List)
Set the partitioner'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