44 #include <Epetra_CrsMatrix.h> 45 #include <Epetra_CrsGraph.h> 46 #include <Epetra_Map.h> 47 #include <Epetra_Comm.h> 56 return(*((
int *) a) - *((
int *) b));
62 while ((n > m) && (l < 31)) {
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)
93 int myMatProc = -1, matProc = -1;
94 int myPID = B.Comm().MyPID();
95 for (
int proc=0; proc<B.Comm().NumProc(); proc++)
97 if (B.NumGlobalEntries64() == B.NumMyEntries())
100 B.Comm().MaxAll( &myMatProc, &matProc, 1 );
103 { std::cout <<
"FAIL for Global! All CrsGraph entries must be on one processor!\n"; abort(); }
105 int i= 0, j = 0, k, l = 0, p, pm, q = -1, ns = 0;
114 std::vector<int_type> Mi, Mj, Mnum(nbrr+1,0);
116 if ( matProc == myPID )
118 if ( verbose ) { std::printf(
" Matrix Size = %d Number of Blocks = %d\n",nrr, nbrr); }
124 bstree = csr_bst(nbrr);
129 for( i = 0; i < nrr; i++ ){
132 m = EPETRA_MAX(m,j) ;
133 j = B.NumGlobalIndices(i);
135 j += B.NumGlobalIndices(i);
139 m = EPETRA_MAX(m,j) ;
141 colstack = (
int*) malloc( EPETRA_MAX(m,1) *
sizeof(int) );
147 nzM = 0; q = -1; l = 0;
150 for( i = 0; i <= nrr; i++ ){
152 if( q > 0 ) std::qsort(colstack,q+1,
sizeof(
int),
compare_ints);
154 for( j=1; j<=q ; j++ ){
155 if( colstack[j] > colstack[j-1] ) ++ns;
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){
169 if( r[bstree[p]+1] <= j) p = 2*p+2;
172 if( p > nbrr || ns > tree_height ) {
174 std::printf(
"error: p %d nbrr %d ns %d %d\n",p,nbrr,ns,j);
break;
177 colstack[++q] = bstree[p];
184 if ( matProc == myPID && verbose )
185 std::printf(
"nzM = %d \n", nzM );
188 q = -1; l = 0; pm = -1;
189 for( i = 0; i <= nrr; i++ ){
191 if( q > 0 ) std::qsort(colstack,q+1,
sizeof(colstack[0]),
compare_ints);
194 Mj[pm] = colstack[0];
196 for( j=1; j<=q ; j++ ){
197 if( colstack[j] > colstack[j-1] ){
199 Mj[pm] = colstack[j];
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){
216 if( r[bstree[p]+1] <= j) p = 2*p+2;
220 colstack[++q] = bstree[p];
224 if ( bstree ) free ( bstree );
225 if ( colstack ) free( colstack );
228 weights.resize( nbrr );
229 for( l=0; l<nbrr; l++) weights[l] = r[l+1] - r[l];
232 Teuchos::RCP<Epetra_Map> newMap;
233 if ( matProc == myPID )
234 newMap = Teuchos::rcp(
new Epetra_Map(nbrr, nbrr, 0, B.Comm() ) );
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]] );
241 newGraph->FillComplete();
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);
253 #ifndef EPETRA_NO_64BIT_GLOBAL_INDICES 254 if(B.RowMap().GlobalIndicesLongLong()) {
255 return compute<long long>(B, nbrr, r, weights, verbose);
259 throw "EpetraExt::BlockAdjacencyGraph::compute: GlobalIndices type unknown";
271 int* BlockAdjacencyGraph::csr_bst(
int n )
273 int i, l=1, m, nstack = 0, nexp=0, os, *array, *stack;
275 if( n == 0 )
return(NULL);
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;
287 i = stack[3*nstack]; os = stack[3*nstack+1]; m = stack[3*nstack+2];
288 array[i] = csr_bstrootindex(m) + os;
290 stack[3*nstack] = 2*i+2; stack[3*nstack+1] = array[i] + 1 ; stack[3*nstack+2] = m-array[i]-1+os;
294 stack[3*nstack] = 2*i+1; stack[3*nstack+1] = os; stack[3*nstack+2] = array[i] - os;
297 if( nstack > max_nstack ) max_nstack = nstack;
308 int BlockAdjacencyGraph::csr_bstrootindex(
int n )
310 int l = 1, nexp = 0, i;
311 if( n == 0)
return(-1);
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.
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...