Ifpack Package Browser (Single Doxygen Collection)  Development
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
Ifpack_LinePartitioner.cpp
Go to the documentation of this file.
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"
46 #include "Ifpack_LinePartitioner.h"
47 #include "Ifpack_Graph.h"
48 #include "Epetra_Util.h"
49 #include <cfloat>
50 
51 // ============================================================================
52 
53 inline 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 
61 inline 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 // ============================================================================
70 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 {
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 // ============================================================================
137 int 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
224 
225  // Resize Parts_
226  Parts_.resize(NumLocalParts_);
227  return(0);
228 }
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.
int Compute_Blocks_AutoLine(int *blockIndices) const
double compute_distance_coordinates(double x0, double y0, double z0, int nn, const double *xvals, const double *yvals, const double *zvals)
const double tol
const int N
const Ifpack_Graph * Graph_
Reference to the graph to be partitioned.
void local_automatic_line_search(int NumEqns, int *blockIndices, int last, int next, int LineID, double tol, int *itemp, double *dtemp) const
const Epetra_RowMatrix * Matrix_
int NumLocalParts_
Number of local subgraphs.
double compute_distance_matrix_entries(const double *vals, int id, int NumEqns)
virtual int ExtractMyRowCopy(int MyRow, int Length, int &NumEntries, double *Values, int *Indices) const =0
static void Sort(bool SortAscending, int NumKeys, T *Keys, int NumDoubleCompanions, double **DoubleCompanions, int NumIntCompanions, int **IntCompanions, int NumLongLongCompanions, long long **LongLongCompanions)
std::vector< std::vector< int > > Parts_
Parts_[i][j] is the ID of the j-th row contained in the (overlapping)
#define IFPACK_CHK_ERR(ifpack_err)
int n