EpetraExt Package Browser (Single Doxygen Collection) Development
Loading...
Searching...
No Matches
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
52namespace 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
232 Teuchos::RCP<Epetra_Map> newMap;
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
#define EPETRA_MAX(x, y)
Copy
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,...
bool GlobalIndicesInt() const
bool GlobalIndicesLongLong() const
virtual int MaxAll(double *PartialMaxs, double *GlobalMaxs, int Count) const=0
virtual int MyPID() const=0
int ExtractMyRowView(int LocalRow, int &NumIndices, int *&Indices) const
const Epetra_BlockMap & RowMap() const
long long NumGlobalEntries64() const
const Epetra_Comm & Comm() const
int NumMyEntries() const
int NumMyRows() const
int NumGlobalIndices(long long Row) const
EpetraExt::BlockCrsMatrix: A class for constructing a distributed block matrix.
int compare_ints(const void *a, const void *b)
Given an Epetra_CrsGraph that has block structure, an adjacency graph is constructed representing the...