ShyLU  Version of the Day
shylu_util.cpp
Go to the documentation of this file.
1 
2 //@HEADER
3 // ************************************************************************
4 //
5 // ShyLU: Hybrid preconditioner package
6 // Copyright 2012 Sandia Corporation
7 //
8 // Under the terms of Contract DE-AC04-94AL85000 with Sandia Corporation,
9 // the U.S. Government retains certain rights in this software.
10 //
11 // Redistribution and use in source and binary forms, with or without
12 // modification, are permitted provided that the following conditions are
13 // met:
14 //
15 // 1. Redistributions of source code must retain the above copyright
16 // notice, this list of conditions and the following disclaimer.
17 //
18 // 2. Redistributions in binary form must reproduce the above copyright
19 // notice, this list of conditions and the following disclaimer in the
20 // documentation and/or other materials provided with the distribution.
21 //
22 // 3. Neither the name of the Corporation nor the names of the
23 // contributors may be used to endorse or promote products derived from
24 // this software without specific prior written permission.
25 //
26 // THIS SOFTWARE IS PROVIDED BY SANDIA CORPORATION "AS IS" AND ANY
27 // EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
28 // IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR
29 // PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL SANDIA CORPORATION OR THE
30 // CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL,
31 // EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO,
32 // PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR
33 // PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF
34 // LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING
35 // NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
36 // SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
37 //
38 // Questions? Contact Michael A. Heroux (maherou@sandia.gov)
39 //
40 // ************************************************************************
41 //@HEADER
42 
51 #include <assert.h>
52 #include <fstream>
53 
54 #include "shylu_util.h"
55 
56 #ifdef HAVE_SHYLUCORE_MPI
57 #include "Epetra_MpiComm.h"
58 #else
59 #include "Epetra_SerialComm.h"
60 #endif
61 #include "Epetra_CrsMatrix.h"
62 
63 // EpetraExt includes
64 #include "EpetraExt_RowMatrixOut.h"
65 
66 //Teuchos includes
68 
69 //Isorropia includes
70 #include "Isorropia_Epetra.hpp"
71 #include "Isorropia_EpetraRedistributor.hpp"
72 #include "Isorropia_EpetraPartitioner.hpp"
73 
74 
75 
76 using namespace std;
77 
78 
79 // Currently takes onle MpiComm
80 Epetra_CrsMatrix *balanceAndRedistribute(Epetra_CrsMatrix *A,
81  Teuchos::ParameterList isoList)
82 {
83  // int myPID = A->Comm().MyPID(); // unused
84 
85  // Debug [
86  Epetra_Map ARowMap = A->RowMap();
87  // int nrows = ARowMap.NumMyElements(); // unused
88  // int *rows = ARowMap.MyGlobalElements(); // unused
89  // ]
90 
91  // ==================== Symbolic factorization =========================
92  // 1. Partition and redistribute [
93  Isorropia::Epetra::Partitioner *partitioner = new
94  Isorropia::Epetra::Partitioner(A, isoList, false);
95  partitioner->partition();
96 
97  Isorropia::Epetra::Redistributor rd(partitioner);
98  Epetra_CrsMatrix *newA;
99  rd.redistribute(*A, newA);
100  // ]
101  EpetraExt::RowMatrixToMatlabFile("A.mat", *newA);
102 
103  delete partitioner;
104  return newA;
105 }
106 
107 /* TODO : Do this only for Debug ? */
108 void checkMaps(Epetra_CrsMatrix *A)
109 {
110  // Get column map
111  Epetra_Map AColMap = A->ColMap();
112  // int ncols = AColMap.NumMyElements(); // unused
113  // int *cols = AColMap.MyGlobalElements(); // unused
114 
115  // Get domain map
116  Epetra_Map ADomainMap = A->DomainMap();
117 
118  int nelems = ADomainMap.NumMyElements();
119 #ifndef NDEBUG
120  // mfh 25 May 2015: Only used in an assert() below.
121  // assert() is defined to nothing in a release build.
122  int *dom_cols = ADomainMap.MyGlobalElements();
123 #endif // NDEBUG
124 
125  // Get range map
126  Epetra_Map ARangeMap = A->RangeMap();
127  // int npts = ARangeMap.NumMyElements(); // unused
128 #ifndef NDEBUG
129  // mfh 25 May 2015: Only used in an assert() below.
130  // assert() is defined to nothing in a release build.
131  int *ran_cols = ARangeMap.MyGlobalElements();
132 #endif // NDEBUG
133 
134  // Get row map
135  Epetra_Map ARowMap = A->RowMap();
136  // int nrows = ARowMap.NumMyElements(); // unused
137 #ifndef NDEBUG
138  // mfh 25 May 2015: Only used in an assert() below.
139  // assert() is defined to nothing in a release build.
140  int *rows = ARowMap.MyGlobalElements();
141 #endif // NDEBUG
142 
143  //cout <<"In PID ="<< A->Comm().MyPID() <<" #cols="<< ncols << " #rows="<<
144  //nrows <<" #domain elems="<< nelems <<" #range elems="<< npts << endl;
145  // See if domain map == range map == row map
146  for (int i = 0; i < nelems ; i++)
147  {
148  // Will this always be the case ? We will find out if assertion fails !
149  assert(dom_cols[i] == ran_cols[i]);
150  assert(rows[i] == ran_cols[i]);
151  }
152 }
153 
154 // TODO: SNumGlobalCols never used
155 void findLocalColumns(Epetra_CrsMatrix *A, int *gvals, int &SNumGlobalCols)
156 {
157 
158  int n = A->NumGlobalRows();
159  // Get column map
160  Epetra_Map AColMap = A->ColMap();
161  int ncols = AColMap.NumMyElements();
162  int *cols = AColMap.MyGlobalElements(); // TODO : Indexing using GID !
163 
164  // 2. Find column permutation [
165  // Find all columns in this proc
166  int *vals = new int[n]; // vector of size n, not ncols
167  for (int i = 0; i < n ; i++)
168  {
169  vals[i] = 0;
170  gvals[i] = 0;
171  }
172 
173  // Null columns in A are not part of any proc
174  for (int i = 0; i < ncols ; i++)
175  {
176  vals[cols[i]] = 1; // Set to 1 for locally owned columns
177  }
178 
179  // Bottleneck?: Compute the column permutation
180  A->Comm().SumAll(vals, gvals, n);
181 
182  SNumGlobalCols = 0;
183  for (int i = 0; i < n ; i++)
184  {
185  //cout << gvals[i] ;
186  if (gvals[i] > 1)
187  SNumGlobalCols++;
188  }
189  //cout << endl;
190  //cout << "Snum Global cols=" << SNumGlobalCols << endl;
191 
192  delete[] vals;
193  return;
194 }
195 
196 // This function uses a very simple tie-breaking heuristic to find a
197 // "narrow" separator from a wide separator. The vertices in the proc with
198 // smaller procID will become part of the separator
199 // This is not a true narrow separator, which needs a vertex cover algorithm.
200 // This is like a medium separator !
201 // TODO : This assumes symmetry I guess, Check
202 void findNarrowSeparator(Epetra_CrsMatrix *A, int *gvals)
203 {
204  int nentries;
205  double *values;
206  int *indices;
207  int n = A->NumGlobalRows();
208 
209  int myPID = A->Comm().MyPID();
210 
211  // Get row map
212  Epetra_Map rMap = A->RowMap();
213  Epetra_Map cMap = A->ColMap();
214  int *rows = rMap.MyGlobalElements();
215  int relems = rMap.NumMyElements();
216 
217  int *vals = new int[n]; // vector of size n, not ncols
218  int *allGIDs = new int[n]; // vector of size n, not ncols
219  for (int i = 0; i < n ; i++) // initialize to zero
220  {
221  vals[i] = 0;
222  }
223 
224  // Rows are uniquely owned, so this will work
225  for (int i = 0; i < relems ; i++)
226  {
227  vals[rows[i]] = myPID; // I own relems[i]
228  }
229 
230 
231  // **************** Collective communication **************
232  // This will not scale well for very large number of nodes
233  // But on the node this should be fine
234  A->Comm().SumAll(vals, allGIDs, n);
235 
236  // At this point all procs know who owns what rows
237  for (int i = 0; i < n ; i++) // initialize to zero
238  vals[i] = 0;
239 
240  int gid, cgid;
241  for (int i = 0; i < relems; i++)
242  {
243  gid = rows[i];
244  //cout << "PID=" << myPID << " " << "rowid=" << gid ;
245  if (gvals[gid] != 1)
246  {
247  //cout << " in the sep ";
248  bool movetoBlockDiagonal = false;
249  // mfh 25 May 2015: This call used to assign its (int)
250  // return value to 'err'. I got rid of this, because
251  // 'err' was unused. This resulted in a "set but unused
252  // variable" warning.
253  (void) A->ExtractMyRowView(i, nentries, values, indices);
254  //cout << " with nentries= "<< nentries;
255 
256  assert(nentries != 0);
257  for (int j = 0; j < nentries; j++)
258  {
259  cgid = cMap.GID(indices[j]);
260  assert(cgid != -1);
261  if (gvals[cgid] == 1 || allGIDs[cgid] == myPID)
262  continue; // simplify the rest
263 
264  /*if (numProcs == 2)
265  {*/
266  if (allGIDs[cgid] < myPID)
267  {
268  // row cgid is owned by a proc with smaller PID
269  movetoBlockDiagonal = true;
270  //cout << "\t mving to diag because of column" << cgid;
271  }
272  else
273  {
274  // There is at least one edge from this vertex to a
275  // vertex in a proc with PID > myPID, cannot move
276  // to diagonal. This is too restrictive, but
277  // important for correctness, until we can use a
278  // vertex cover algorithm.
279  movetoBlockDiagonal = false;
280  break;
281  //cout << "\tNo problem with cgid=" << cgid << "in sep";
282  }
283  /*}
284  else
285  {
286  if (myPID == 0 && allGIDs[cgid] == numProcs-1)
287  {
288  // row cgid is owned by a proc with smaller PID
289  movetoBlockDiagonal = true;
290  cout << "\t I AM HERE mving to diag because of column" << cgid;
291  }
292  else if (myPID == numProcs-1 && allGIDs[cgid] == 0)
293  {
294  cout << "\t I AM HERE to continue " << cgid;
295  continue;
296  }
297  else if (allGIDs[cgid] < myPID)
298  {
299  // row cgid is owned by a proc with smaller PID
300  movetoBlockDiagonal = true;
301  cout << "\t mving to diag because of column" << cgid;
302  }
303  else
304  {
305  //cout << "\tNo problem with cgid=" << cgid << "in sep";
306  }
307  }*/
308  }
309  if (movetoBlockDiagonal)
310  {
311  //cout << "Moving to Diagonal";
312  vals[gid] = 1;
313  gvals[gid] = 1; // The smaller PIDs have to know about this
314  // change. Send the change using gvals.
315  }
316  }
317  else
318  {
319  // do nothing, in the diagonal block already
320  //cout << "In the diagonal block";
321  }
322  //cout << endl;
323  }
324 
325  // Reuse allGIDs to propagate the result of moving to diagonal
326  for (int i = 0; i < n ; i++) // initialize to zero
327  allGIDs[i] = 0;
328 
329  A->Comm().SumAll(vals, allGIDs, n);
330  for (int i = 0; i < n ; i++)
331  {
332  if (allGIDs[i] == 1)
333  {
334  // Some interface columns will have gvals[1] after this
335  // as the separator is narrow now.
336  gvals[i] = 1; // GIDs as indices assumption
337  }
338  }
339 
340  delete[] vals;
341  delete[] allGIDs;
342 
343 }
344 
345 void findBlockElems(Epetra_CrsMatrix *A, int nrows, int *rows, int *gvals,
346  int Lnr, int *LeftElems,
347  int Rnr, int *RightElems, string s1, string s2, bool cols)
348 {
349 
350  int gid;
351  int rcnt = 0; int lcnt = 0;
352  // Assemble ids in two arrays
353  ostringstream ssmsg1;
354  ostringstream ssmsg2;
355 
356 #ifdef DUMP_MATRICES
357  ostringstream fnamestr;
358  fnamestr << s1 << ".mat";
359  string fname = fnamestr.str();
360  ofstream outD(fname.c_str());
361 
362  ostringstream fnamestrR;
363  fnamestrR << s2 << ".mat";
364  string fnameR = fnamestrR.str();
365  ofstream outR(fnameR.c_str());
366 #endif
367 
368  ssmsg1 << s1;
369  ssmsg2 << s2;
370  for (int i = 0; i < nrows; i++)
371  {
372  gid = rows[i];
373  assert (gvals[gid] >= 1);
374  // If the row is local & row/column is not shared then row/column
375  // belongs to D (this is not true for R, which can have more columns
376  // than D)
377  //if (A->LRID(gid) != -1 && gvals[gid] == 1)
378  if (gvals[gid] == 1)
379  {
380  if (cols && A->LRID(gid) == -1) continue;
381  assert(lcnt < Lnr);
382  LeftElems[lcnt++] = gid;
383  ssmsg1 << gid << " ";
384 #ifdef DUMP_MATRICES
385  outD << gid << endl;
386 #endif
387  }
388  else
389  {
390  assert(rcnt < Rnr);
391  RightElems[rcnt++] = gid;
392  ssmsg2 << gid << " ";
393 #ifdef DUMP_MATRICES
394  outR << gid << endl;
395 #endif
396  }
397  }
398 
399 #ifdef DUMP_MATRICES
400  outD.close();
401  outR.close();
402 #endif
403 
404 #ifdef DEBUG
405  cout << ssmsg1.str() << endl;
406  cout << ssmsg2.str() << endl;
407 #endif
408  ssmsg1.clear(); ssmsg1.str("");
409  ssmsg2.clear(); ssmsg2.str("");
410 
411  assert(lcnt == Lnr);
412  assert(rcnt == Rnr);
413 }
Utilities for ShyLU.