IFPACK  Development
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 
50 // ============================================================================
51 inline void Ifpack_LinePartitioner::local_automatic_line_search(int NumEqns, int * blockIndices, int last, int next, int LineID, double tol, int *itemp, double * dtemp) const {
52  double *xvals=xcoord_, *yvals=ycoord_, *zvals=zcoord_;
53 
54  int N = NumMyRows();
55 
56  int allocated_space = MaxNumEntries();
57  int * cols = itemp;
58  int * indices = &itemp[allocated_space];
59  double * dist = dtemp;
60 
61 
62  while (blockIndices[next] == -1) {
63  // Get the next row
64  int n=0;
65  int neighbors_in_line=0;
66 
67  Graph_->ExtractMyRowCopy(next,allocated_space,n,cols);
68  double x0 = (xvals) ? xvals[next/NumEqns] : 0.0;
69  double y0 = (yvals) ? yvals[next/NumEqns] : 0.0;
70  double z0 = (zvals) ? zvals[next/NumEqns] : 0.0;
71 
72  // Calculate neighbor distances & sort
73  int neighbor_len=0;
74  for(int i=0; i<n; i+=NumEqns) {
75  double mydist = 0.0;
76  if(cols[i] >=N) continue; // Check for off-proc entries
77  int nn = cols[i] / NumEqns;
78  if(blockIndices[nn]==LineID) neighbors_in_line++;
79  if(xvals) mydist += (x0 - xvals[nn]) * (x0 - xvals[nn]);
80  if(yvals) mydist += (y0 - yvals[nn]) * (y0 - yvals[nn]);
81  if(zvals) mydist += (z0 - zvals[nn]) * (z0 - zvals[nn]);
82  dist[neighbor_len] = sqrt(mydist);
83  indices[neighbor_len]=cols[i];
84  neighbor_len++;
85  }
86  // If more than one of my neighbors is already in this line. I
87  // can't be because I'd create a cycle
88  if(neighbors_in_line > 1) break;
89 
90  // Otherwise add me to the line
91  for(int k=0; k<NumEqns; k++)
92  blockIndices[next + k] = LineID;
93 
94  // Try to find the next guy in the line (only check the closest two that aren't element 0 (diagonal))
95  Epetra_Util::Sort(true,neighbor_len,dist,0,0,1,&indices,0,0);
96 
97  if(neighbor_len > 2 && indices[1] != last && blockIndices[indices[1]] == -1 && dist[1]/dist[neighbor_len-1] < tol) {
98  last=next;
99  next=indices[1];
100  }
101  else if(neighbor_len > 3 && indices[2] != last && blockIndices[indices[2]] == -1 && dist[2]/dist[neighbor_len-1] < tol) {
102  last=next;
103  next=indices[2];
104  }
105  else {
106  // I have no further neighbors in this line
107  break;
108  }
109  }
110 }
111 
112 // ============================================================================
113 int Ifpack_LinePartitioner::Compute_Blocks_AutoLine(int * blockIndices) const {
114  double *xvals=xcoord_, *yvals=ycoord_, *zvals=zcoord_;
115  int NumEqns = NumEqns_;
116  double tol = threshold_;
117  int N = NumMyRows();
118  int allocated_space = MaxNumEntries();
119 
120  int * cols = new int[2*allocated_space];
121  int * indices = &cols[allocated_space];
122  double * dist = new double[allocated_space];
123 
124  int * itemp = new int[2*allocated_space];
125  double *dtemp = new double[allocated_space];
126 
127  int num_lines = 0;
128 
129  for(int i=0; i<N; i+=NumEqns) {
130  int nz=0;
131  // Short circuit if I've already been blocked
132  if(blockIndices[i] !=-1) continue;
133 
134  // Get neighbors and sort by distance
135  Graph_->ExtractMyRowCopy(i,allocated_space,nz,cols);
136  double x0 = (xvals) ? xvals[i/NumEqns] : 0.0;
137  double y0 = (yvals) ? yvals[i/NumEqns] : 0.0;
138  double z0 = (zvals) ? zvals[i/NumEqns] : 0.0;
139 
140  int neighbor_len=0;
141  for(int j=0; j<nz; j+=NumEqns) {
142  double mydist = 0.0;
143  int nn = cols[j] / NumEqns;
144  if(cols[j] >=N) continue; // Check for off-proc entries
145  if(xvals) mydist += (x0 - xvals[nn]) * (x0 - xvals[nn]);
146  if(yvals) mydist += (y0 - yvals[nn]) * (y0 - yvals[nn]);
147  if(zvals) mydist += (z0 - zvals[nn]) * (z0 - zvals[nn]);
148  dist[neighbor_len] = sqrt(mydist);
149  indices[neighbor_len]=cols[j];
150  neighbor_len++;
151  }
152 
153  Epetra_Util::Sort(true,neighbor_len,dist,0,0,1,&indices,0,0);
154 
155  // Number myself
156  for(int k=0; k<NumEqns; k++)
157  blockIndices[i + k] = num_lines;
158 
159  // Fire off a neighbor line search (nearest neighbor)
160  if(neighbor_len > 2 && dist[1]/dist[neighbor_len-1] < tol) {
161  local_automatic_line_search(NumEqns,blockIndices,i,indices[1],num_lines,tol,itemp,dtemp);
162  }
163  // Fire off a neighbor line search (second nearest neighbor)
164  if(neighbor_len > 3 && dist[2]/dist[neighbor_len-1] < tol) {
165  local_automatic_line_search(NumEqns,blockIndices,i,indices[2],num_lines,tol,itemp,dtemp);
166  }
167 
168  num_lines++;
169  }
170 
171  // Cleanup
172  delete [] cols;
173  delete [] dist;
174  delete [] itemp;
175  delete [] dtemp;
176 
177  return num_lines;
178 }
179 //==============================================================================
181 {
182  // Sanity Checks
183  if(!xcoord_ && !ycoord_ && !zcoord_) IFPACK_CHK_ERR(-1);
184  if((int)Partition_.size() != NumMyRows()) IFPACK_CHK_ERR(-2);
185 
186  // Short circuit
187  if(Partition_.size() == 0) {NumLocalParts_ = 0; return 0;}
188 
189  // Set partitions to -1 to initialize algorithm
190  for(int i=0; i<NumMyRows(); i++)
191  Partition_[i] = -1;
192 
193  // Use the auto partitioner
194  NumLocalParts_ = Compute_Blocks_AutoLine(&Partition_[0]);
195 
196  // Resize Parts_
197  Parts_.resize(NumLocalParts_);
198  return(0);
199 }
virtual int ExtractMyRowCopy(int MyRow, int LenOfIndices, int &NumIndices, int *Indices) const =0
Extracts a copy of input local row.
int MaxNumEntries() const
Returns the max number of local entries in a row.
int NumMyRows() const
Returns the number of local rows.
std::vector< int > Partition_
Partition_[i] contains the ID of non-overlapping part it belongs to.
int ComputePartitions()
Computes the partitions. Returns 0 if successful.
const Ifpack_Graph * Graph_
Reference to the graph to be partitioned.
int NumLocalParts_
Number of local subgraphs.
std::vector< std::vector< int > > Parts_
Parts_[i][j] is the ID of the j-th row contained in the (overlapping)