EpetraExt Package Browser (Single Doxygen Collection)  Development
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
EpetraExt_BlockAdjacencyGraph.cpp
Go to the documentation of this file.
1 //@HEADER
2 // ***********************************************************************
3 //
4 // EpetraExt: Epetra Extended - Linear Algebra Services Package
5 // Copyright (2011) Sandia Corporation
6 //
7 // Under the terms of Contract DE-AC04-94AL85000 with Sandia Corporation,
8 // the U.S. Government retains certain rights in this software.
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 
43 
44 #include <Epetra_CrsMatrix.h>
45 #include <Epetra_CrsGraph.h>
46 #include <Epetra_Map.h>
47 #include <Epetra_Comm.h>
48 
49 #include <math.h>
50 #include <cstdlib>
51 
52 namespace EpetraExt {
53 
54  int compare_ints(const void *a, const void *b)
55  {
56  return(*((int *) a) - *((int *) b));
57  }
58 
59  int ceil31log2(int n)
60  { // Given 1 <= n < 2^31, find l such that 2^(l-1) < n <= 2^(l)
61  int l=0, m = 1;
62  while ((n > m) && (l < 31)) {
63  m = 2*m;
64  ++l;
65  }
66  return(l);
67  }
68 // Purpose: Compute the block connectivity graph of a matrix.
69 // An nrr by nrr sparse matrix admits a (Dulmage-Mendelsohn)
70 // permutation to block triangular form with nbrr blocks.
71 // Abstractly, the block triangular form corresponds to a partition
72 // of the set of integers {0,...,n-1} into nbrr disjoint sets.
73 // The graph of the sparse matrix, with nrr vertices, may be compressed
74 // into the graph of the blocks, a graph with nbrr vertices, that is
75 // called here the block connectivity graph.
76 // The partition of the rows and columns of B is represented by
77 // r(0:nbrr), 0 = r(0) < r(1) < .. < r(nbrr) = nrr,
78 // The graph (Mp,Mj) of the nbrr x nbrr matrix is represened by
79 // a sparse matrix in sparse coordinate format.
80 // Mp: row indices, dimension determined here (nzM).
81 // Mj: column indices, dimension determined here (nzM).
82 // The integer vector, weights, of block sizes (dimension nbrr) is also
83 // computed here.
84 // The case of nbrr proportional to nrr is critical. One must
85 // efficiently look up the column indices of B in the partition.
86 // This is done here using a binary search tree, so that the
87 // look up cost is nzB*log2(nbrr).
88 
89  template<typename int_type>
90  Teuchos::RCP<Epetra_CrsGraph> BlockAdjacencyGraph::compute( Epetra_CrsGraph& B, int nbrr, std::vector<int>&r, std::vector<double>& weights, bool verbose)
91  {
92  // Check if the graph is on one processor.
93  int myMatProc = -1, matProc = -1;
94  int myPID = B.Comm().MyPID();
95  for (int proc=0; proc<B.Comm().NumProc(); proc++)
96  {
97  if (B.NumGlobalEntries64() == B.NumMyEntries())
98  myMatProc = myPID;
99  }
100  B.Comm().MaxAll( &myMatProc, &matProc, 1 );
101 
102  if( matProc == -1)
103  { std::cout << "FAIL for Global! All CrsGraph entries must be on one processor!\n"; abort(); }
104 
105  int i= 0, j = 0, k, l = 0, p, pm, q = -1, ns = 0;
106  int tree_height;
107  int error = -1; /* error detected, possibly a problem with the input */
108  (void) error; // silence "set but not used" warning
109  int nrr; /* number of rows in B */
110  int nzM = 0; /* number of edges in graph */
111  int m = 0; /* maximum number of nonzeros in any block row of B */
112  int* colstack = 0; /* stack used to process each block row */
113  int* bstree = 0; /* binary search tree */
114  std::vector<int_type> Mi, Mj, Mnum(nbrr+1,0);
115  nrr = B.NumMyRows();
116  if ( matProc == myPID )
117  {
118  if ( verbose ) { std::printf(" Matrix Size = %d Number of Blocks = %d\n",nrr, nbrr); }
119  }
120  else
121  {
122  nrr = -1; /* Prevent processor from doing any computations */
123  }
124  bstree = csr_bst(nbrr); /* 0 : nbrr-1 */
125  tree_height = ceil31log2(nbrr) + 1;
126  error = -1;
127 
128  l = 0; j = 0; m = 0;
129  for( i = 0; i < nrr; i++ ){
130  if( i >= r[l+1] ){
131  ++l; /* new block row */
132  m = EPETRA_MAX(m,j) ; /* nonzeros in block row */
133  j = B.NumGlobalIndices(i);
134  }else{
135  j += B.NumGlobalIndices(i);
136  }
137  }
138  /* one more time for the final block */
139  m = EPETRA_MAX(m,j) ; /* nonzeros in block row */
140 
141  colstack = (int*) malloc( EPETRA_MAX(m,1) * sizeof(int) );
142  // The compressed graph is actually computed twice,
143  // due to concerns about memory limitations. First,
144  // without memory allocation, just nzM is computed.
145  // Next Mj is allocated. Then, the second time, the
146  // arrays are actually populated.
147  nzM = 0; q = -1; l = 0;
148  int * indices;
149  int numEntries;
150  for( i = 0; i <= nrr; i++ ){
151  if( i >= r[l+1] ){
152  if( q > 0 ) std::qsort(colstack,q+1,sizeof(int),compare_ints); /* sort stack */
153  if( q >= 0 ) ns = 1; /* l, colstack[0] M */
154  for( j=1; j<=q ; j++ ){ /* delete copies */
155  if( colstack[j] > colstack[j-1] ) ++ns;
156  }
157  nzM += ns; /*M->p[l+1] = M->p[l] + ns;*/
158  ++l;
159  q = -1;
160  }
161  if( i < nrr ){
162  B.ExtractMyRowView( i, numEntries, indices );
163  for( k = 0; k < numEntries; k++){
164  j = indices[k]; ns = 0; p = 0;
165  while( (r[bstree[p]] > j) || (j >= r[bstree[p]+1]) ){
166  if( r[bstree[p]] > j){
167  p = 2*p+1;
168  }else{
169  if( r[bstree[p]+1] <= j) p = 2*p+2;
170  }
171  ++ns;
172  if( p > nbrr || ns > tree_height ) {
173  error = j;
174  std::printf("error: p %d nbrr %d ns %d %d\n",p,nbrr,ns,j); break;
175  }
176  }
177  colstack[++q] = bstree[p];
178  }
179  //if( error >-1 ){ std::printf("%d\n",error); break; }
180  // p > nbrr is a fatal error that is ignored
181  }
182  }
183 
184  if ( matProc == myPID && verbose )
185  std::printf("nzM = %d \n", nzM );
186  Mi.resize( nzM );
187  Mj.resize( nzM );
188  q = -1; l = 0; pm = -1;
189  for( i = 0; i <= nrr; i++ ){
190  if( i >= r[l+1] ){
191  if( q > 0 ) std::qsort(colstack,q+1,sizeof(colstack[0]),compare_ints); /* sort stack */
192  if( q >= 0 ){
193  Mi[++pm] = l;
194  Mj[pm] = colstack[0];
195  }
196  for( j=1; j<=q ; j++ ){ /* delete copies */
197  if( colstack[j] > colstack[j-1] ){ /* l, colstack[j] */
198  Mi[++pm] = l;
199  Mj[pm] = colstack[j];
200  }
201  }
202  ++l;
203  Mnum[l] = pm + 1;
204 
205  /* sparse row format: M->p[l+1] = M->p[l] + ns; */
206  q = -1;
207  }
208  if( i < nrr ){
209  B.ExtractMyRowView( i, numEntries, indices );
210  for( k = 0; k < numEntries; k++){
211  j = indices[k]; ns = 0; p = 0;
212  while( (r[bstree[p]] > j) || (j >= r[bstree[p]+1]) ){
213  if( r[bstree[p]] > j){
214  p = 2*p+1;
215  }else{
216  if( r[bstree[p]+1] <= j) p = 2*p+2;
217  }
218  ++ns;
219  }
220  colstack[++q] = bstree[p];
221  }
222  }
223  }
224  if ( bstree ) free ( bstree );
225  if ( colstack ) free( colstack );
226 
227  // Compute weights as number of rows in each block.
228  weights.resize( nbrr );
229  for( l=0; l<nbrr; l++) weights[l] = r[l+1] - r[l];
230 
231  // Compute Epetra_CrsGraph and return
233  if ( matProc == myPID )
234  newMap = Teuchos::rcp( new Epetra_Map(nbrr, nbrr, 0, B.Comm() ) );
235  else
236  newMap = Teuchos::rcp( new Epetra_Map( nbrr, 0, 0, B.Comm() ) );
237  Teuchos::RCP<Epetra_CrsGraph> newGraph = Teuchos::rcp( new Epetra_CrsGraph( Copy, *newMap, 0 ) );
238  for( l=0; l<newGraph->NumMyRows(); l++) {
239  newGraph->InsertGlobalIndices( l, Mnum[l+1]-Mnum[l], &Mj[Mnum[l]] );
240  }
241  newGraph->FillComplete();
242 
243  return (newGraph);
244  }
245 
246  Teuchos::RCP<Epetra_CrsGraph> BlockAdjacencyGraph::compute( Epetra_CrsGraph& B, int nbrr, std::vector<int>&r, std::vector<double>& weights, bool verbose) {
247 #ifndef EPETRA_NO_32BIT_GLOBAL_INDICES
248  if(B.RowMap().GlobalIndicesInt()) {
249  return compute<int>(B, nbrr, r, weights, verbose);
250  }
251  else
252 #endif
253 #ifndef EPETRA_NO_64BIT_GLOBAL_INDICES
254  if(B.RowMap().GlobalIndicesLongLong()) {
255  return compute<long long>(B, nbrr, r, weights, verbose);
256  }
257  else
258 #endif
259  throw "EpetraExt::BlockAdjacencyGraph::compute: GlobalIndices type unknown";
260  }
261 
262  /*
263  * bst(n) returns the complete binary tree, stored in an integer array of dimension n.
264  * index i has children 2i+1 and 2i+2, root index 0.
265  * A binary search of a sorted n-array: find array(k)
266  * i=0; k = tree(i);
267  * if < array( k ), i = 2i+1;
268  * elif > array( k ), i = 2i+2;
269  * otherwise exit
270  */
272  {
273  int i, l=1, m, nstack = 0, nexp=0, os, *array, *stack;
274  int max_nstack = 0;
275  if( n == 0 ) return(NULL);
276  while( l <= n ){
277  l = 2*l;
278  ++nexp;
279  } /* n < l = 2^nexp */
280  array = (int *) malloc( n * sizeof(int) );
281  stack = (int *) malloc(3*nexp * sizeof(int) );
282  stack[3*nstack] = 0; stack[3*nstack+1] = 0; stack[3*nstack+2] = n;
283  ++nstack;
284  /*if( debug ) std::printf("stack : %d %d %d\n", stack[0] , stack[1], stack[2] );*/
285  while( nstack > 0 ){
286  --nstack;
287  i = stack[3*nstack]; os = stack[3*nstack+1]; m = stack[3*nstack+2];
288  array[i] = csr_bstrootindex(m) + os; /* 5 */
289  if( 2*i+2 < n){ /* right child 4 3 1 3 - 5 - 1 */
290  stack[3*nstack] = 2*i+2; stack[3*nstack+1] = array[i] + 1 ; stack[3*nstack+2] = m-array[i]-1+os; /* 6 10 -3 */
291  ++nstack;
292  }
293  if( 2*i+1 < n){ /* left child */
294  stack[3*nstack] = 2*i+1; stack[3*nstack+1] = os; stack[3*nstack+2] = array[i] - os; /* 5 4 5 */
295  ++nstack;
296  }
297  if( nstack > max_nstack ) max_nstack = nstack;
298  }
299  free( stack );
300  return(array);
301  }
302 
303  /*
304  Given a complete binary search tree with n nodes, return
305  the index of the root node. A nodeless tree, n=0, has no
306  root node (-1). Nodes are indexed from 0 to n-1.
307  */
309  {
310  int l = 1, nexp = 0, i;
311  if( n == 0) return(-1);
312  while( l <= n ){
313  l = 2*l;
314  ++nexp;
315  } /* n < l = 2^nexp */
316  i = l/2 - 1;
317  if(n<4) return(i);
318  if (n-i < l/4)
319  return( n - l/4 );
320  else
321  return(i);
322  }
323 
324 } //namespace EpetraExt
325 
const Epetra_Comm & Comm() const
int NumGlobalIndices(long long Row) const
bool GlobalIndicesLongLong() const
int NumMyEntries() const
Teuchos::RCP< Epetra_CrsGraph > compute(Epetra_CrsGraph &B, int nbrr, std::vector< int > &r, std::vector< double > &weights, bool verbose=false)
Constructs an adjacency graph representing the block connectivity of the input graph, where nbrr is the number of block rows in B and r contains the row index where each block begins.
int NumMyRows() const
bool GlobalIndicesInt() const
int InsertGlobalIndices(int_type GlobalRow, int NumIndices, int_type *Indices)
virtual int MyPID() const =0
virtual int MaxAll(double *PartialMaxs, double *GlobalMaxs, int Count) const =0
TEUCHOS_DEPRECATED RCP< T > rcp(T *p, Dealloc_T dealloc, bool owns_mem)
int compare_ints(const void *a, const void *b)
Given an Epetra_CrsGraph that has block structure, an adjacency graph is constructed representing the...
const Epetra_BlockMap & RowMap() const
std::string error
int ExtractMyRowView(int LocalRow, int &NumIndices, int *&Indices) const
virtual int NumProc() const =0
long long NumGlobalEntries64() const
#define EPETRA_MAX(x, y)