IFPACK Development
Loading...
Searching...
No Matches
Ifpack_LinePartitioner.cpp
1/*@HEADER
2// ***********************************************************************
3//
4// Ifpack: Object-Oriented Algebraic Preconditioner Package
5// Copyright (2002) 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#include "Ifpack_ConfigDefs.h"
44#include "Ifpack_Partitioner.h"
45#include "Ifpack_OverlappingPartitioner.h"
46#include "Ifpack_LinePartitioner.h"
47#include "Ifpack_Graph.h"
48#include "Epetra_Util.h"
49#include <cfloat>
50
51// ============================================================================
52
53inline double compute_distance_coordinates(double x0, double y0, double z0, int nn, const double*xvals,const double*yvals,const double*zvals) {
54 double mydist=0.0;
55 if(xvals) mydist += (x0 - xvals[nn]) * (x0 - xvals[nn]);
56 if(yvals) mydist += (y0 - yvals[nn]) * (y0 - yvals[nn]);
57 if(zvals) mydist += (z0 - zvals[nn]) * (z0 - zvals[nn]);
58 return sqrt(mydist);
59}
60
61inline double compute_distance_matrix_entries(const double *vals, int id, int NumEqns) {
62 double sum=0.0;
63 for(int i=0; i<NumEqns; i++)
64 sum+=std::abs(vals[id+i]);
65 return sum;
66}
67
68
69// ============================================================================
70inline void Ifpack_LinePartitioner::local_automatic_line_search(int NumEqns, int * blockIndices, int last, int next, int LineID, double tol, int *itemp, double * dtemp) const {
71 double *xvals=xcoord_, *yvals=ycoord_, *zvals=zcoord_;
72
73 int N = NumMyRows();
74
75 int allocated_space = MaxNumEntries();
76 int * cols = itemp;
77 int * indices = &itemp[allocated_space];
78 double * merit= dtemp;
79 double * vals = &dtemp[allocated_space];
80
81 while (blockIndices[next] == -1) {
82 // Get the next row
83 int n=0;
84 int neighbors_in_line=0;
85
86 if(mode_==MATRIX_ENTRIES) Matrix_->ExtractMyRowCopy(next,allocated_space,n,vals,cols);
87 else Graph_->ExtractMyRowCopy(next,allocated_space,n,cols);
88
89 // Coordinate distance info
90 double x0 = (xvals) ? xvals[next/NumEqns] : 0.0;
91 double y0 = (yvals) ? yvals[next/NumEqns] : 0.0;
92 double z0 = (zvals) ? zvals[next/NumEqns] : 0.0;
93
94 // Calculate neighbor distances & sort
95 int neighbor_len=0;
96 for(int i=0; i<n; i+=NumEqns) {
97 if(cols[i] >=N) continue; // Check for off-proc entries
98 int nn = cols[i] / NumEqns;
99 if(blockIndices[nn]==LineID) neighbors_in_line++;
100 if(mode_==COORDINATES) merit[neighbor_len] = compute_distance_coordinates(x0,y0,z0,nn,xvals,yvals,zvals);
101 else {
102 merit[neighbor_len] = - compute_distance_matrix_entries(vals,i,NumEqns); // Make this negative here, so we can use the same if tests at coordinates
103 // Boost the diagonal here to ensure it goes first
104 if(cols[i]==next) merit[neighbor_len] = -DBL_MAX;
105 }
106
107 indices[neighbor_len]=cols[i];
108 neighbor_len++;
109 }
110 // If more than one of my neighbors is already in this line. I
111 // can't be because I'd create a cycle
112 if(neighbors_in_line > 1) break;
113
114 // Otherwise add me to the line
115 for(int k=0; k<NumEqns; k++)
116 blockIndices[next + k] = LineID;
117
118 // Try to find the next guy in the line (only check the closest two that aren't element 0 (diagonal))
119 Epetra_Util::Sort(true,neighbor_len,merit,0,0,1,&indices,0,0);
120
121 if(neighbor_len > 2 && indices[1] != last && blockIndices[indices[1]] == -1 && merit[1] < tol*merit[neighbor_len-1]) {
122 last=next;
123 next=indices[1];
124 }
125 else if(neighbor_len > 3 && indices[2] != last && blockIndices[indices[2]] == -1 && merit[2] < tol*merit[neighbor_len-1]) {
126 last=next;
127 next=indices[2];
128 }
129 else {
130 // I have no further neighbors in this line
131 break;
132 }
133 }
134}
135
136// ============================================================================
137int Ifpack_LinePartitioner::Compute_Blocks_AutoLine(int * blockIndices) const {
138 double *xvals=xcoord_, *yvals=ycoord_, *zvals=zcoord_;
139 int NumEqns = NumEqns_;
140 double tol = threshold_;
141 int N = NumMyRows();
142 int allocated_space = MaxNumEntries();
143
144 int * cols = new int[2*allocated_space];
145 int * indices = &cols[allocated_space];
146 double * merit = new double[2*allocated_space];
147 double * vals = &merit[allocated_space];
148
149 int * itemp = new int[2*allocated_space];
150 double *dtemp = new double[2*allocated_space];
151
152
153 int num_lines = 0;
154
155 for(int i=0; i<N; i+=NumEqns) {
156 int nz=0;
157 // Short circuit if I've already been blocked
158 if(blockIndices[i] !=-1) continue;
159
160 // Get neighbors and sort by distance
161 if(mode_==MATRIX_ENTRIES) Matrix_->ExtractMyRowCopy(i,allocated_space,nz,vals,cols);
162 else Graph_->ExtractMyRowCopy(i,allocated_space,nz,cols);
163
164 double x0 = (xvals) ? xvals[i/NumEqns] : 0.0;
165 double y0 = (yvals) ? yvals[i/NumEqns] : 0.0;
166 double z0 = (zvals) ? zvals[i/NumEqns] : 0.0;
167
168 int neighbor_len=0;
169 for(int j=0; j<nz; j+=NumEqns) {
170 int nn = cols[j] / NumEqns;
171 if(cols[j] >=N) continue; // Check for off-proc entries
172 if(mode_==COORDINATES) merit[neighbor_len] = compute_distance_coordinates(x0,y0,z0,nn,xvals,yvals,zvals);
173 else {
174 merit[neighbor_len] = - compute_distance_matrix_entries(vals,j,NumEqns); // Make this negative here, so we can use the same if tests at coordinates
175 // Boost the diagonal here to ensure it goes first
176 if(cols[j]==i) merit[neighbor_len] = -DBL_MAX;
177 }
178 indices[neighbor_len] = cols[j];
179 neighbor_len++;
180 }
181
182 Epetra_Util::Sort(true,neighbor_len,merit,0,0,1,&indices,0,0);
183
184 // Number myself
185 for(int k=0; k<NumEqns; k++)
186 blockIndices[i + k] = num_lines;
187
188 // Fire off a neighbor line search (nearest neighbor)
189 if(neighbor_len > 2 && merit[1] < tol*merit[neighbor_len-1]) {
190 local_automatic_line_search(NumEqns,blockIndices,i,indices[1],num_lines,tol,itemp,dtemp);
191 }
192 // Fire off a neighbor line search (second nearest neighbor)
193 if(neighbor_len > 3 && merit[2] < tol*merit[neighbor_len-1]) {
194 local_automatic_line_search(NumEqns,blockIndices,i,indices[2],num_lines,tol,itemp,dtemp);
195 }
196
197 num_lines++;
198 }
199
200 // Cleanup
201 delete [] cols;
202 delete [] merit;
203 delete [] itemp;
204 delete [] dtemp;
205
206 return num_lines;
207}
208//==============================================================================
210{
211 // Sanity Checks
212 if(mode_==COORDINATES && !xcoord_ && !ycoord_ && !zcoord_) IFPACK_CHK_ERR(-1);
213 if((int)Partition_.size() != NumMyRows()) IFPACK_CHK_ERR(-2);
214
215 // Short circuit
216 if(Partition_.size() == 0) {NumLocalParts_ = 0; return 0;}
217
218 // Set partitions to -1 to initialize algorithm
219 for(int i=0; i<NumMyRows(); i++)
220 Partition_[i] = -1;
221
222 // Use the auto partitioner
223 NumLocalParts_ = Compute_Blocks_AutoLine(&Partition_[0]);
224
225 // Resize Parts_
226 Parts_.resize(NumLocalParts_);
227 return(0);
228}
virtual int ExtractMyRowCopy(int MyRow, int Length, int &NumEntries, double *Values, int *Indices) const=0
static void EPETRA_LIB_DLL_EXPORT Sort(bool SortAscending, int NumKeys, T *Keys, int NumDoubleCompanions, double **DoubleCompanions, int NumIntCompanions, int **IntCompanions, int NumLongLongCompanions, long long **LongLongCompanions)
virtual int ExtractMyRowCopy(int MyRow, int LenOfIndices, int &NumIndices, int *Indices) const =0
Extracts a copy of input local row.
int ComputePartitions()
Computes the partitions. Returns 0 if successful.
std::vector< std::vector< int > > Parts_
Parts_[i][j] is the ID of the j-th row contained in the (overlapping)
int NumLocalParts_
Number of local subgraphs.
int NumMyRows() const
Returns the number of local rows.
int MaxNumEntries() const
Returns the max number of local entries in a row.
std::vector< int > Partition_
Partition_[i] contains the ID of non-overlapping part it belongs to.
const Ifpack_Graph * Graph_
Reference to the graph to be partitioned.