EpetraExt Development
EpetraExt_BlockAdjacencyGraph.cpp
Go to the documentation of this file.
00001 //@HEADER
00002 // ***********************************************************************
00003 //
00004 //     EpetraExt: Epetra Extended - Linear Algebra Services Package
00005 //                 Copyright (2011) Sandia Corporation
00006 //
00007 // Under the terms of Contract DE-AC04-94AL85000 with Sandia Corporation,
00008 // the U.S. Government retains certain rights in this software.
00009 //
00010 // Redistribution and use in source and binary forms, with or without
00011 // modification, are permitted provided that the following conditions are
00012 // met:
00013 //
00014 // 1. Redistributions of source code must retain the above copyright
00015 // notice, this list of conditions and the following disclaimer.
00016 //
00017 // 2. Redistributions in binary form must reproduce the above copyright
00018 // notice, this list of conditions and the following disclaimer in the
00019 // documentation and/or other materials provided with the distribution.
00020 //
00021 // 3. Neither the name of the Corporation nor the names of the
00022 // contributors may be used to endorse or promote products derived from
00023 // this software without specific prior written permission.
00024 //
00025 // THIS SOFTWARE IS PROVIDED BY SANDIA CORPORATION "AS IS" AND ANY
00026 // EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
00027 // IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR
00028 // PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL SANDIA CORPORATION OR THE
00029 // CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL,
00030 // EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO,
00031 // PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR
00032 // PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF
00033 // LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING
00034 // NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
00035 // SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
00036 //
00037 // Questions? Contact Michael A. Heroux (maherou@sandia.gov)
00038 //
00039 // ***********************************************************************
00040 //@HEADER
00041 
00042 #include <EpetraExt_BlockAdjacencyGraph.h>
00043 
00044 #include <Epetra_CrsMatrix.h>
00045 #include <Epetra_CrsGraph.h>
00046 #include <Epetra_Map.h>
00047 #include <Epetra_Comm.h>
00048 
00049 #include <math.h>
00050 #include <cstdlib>
00051 
00052 namespace EpetraExt {
00053 
00054   int compare_ints(const void *a, const void *b)
00055   {
00056     return(*((int *) a) - *((int *) b));
00057   }
00058 
00059   int ceil31log2(int n)
00060   { // Given 1 <= n < 2^31, find l such that 2^(l-1) < n <= 2^(l)
00061     int l=0, m = 1;
00062     while ((n > m) && (l < 31)) {
00063       m = 2*m;
00064       ++l;
00065     }
00066     return(l);
00067   }
00068 //  Purpose: Compute the block connectivity graph of a matrix.
00069 //  An nrr by nrr sparse matrix admits a (Dulmage-Mendelsohn)
00070 //  permutation to block triangular form with nbrr blocks.
00071 //  Abstractly, the block triangular form corresponds to a partition
00072 //  of the set of integers {0,...,n-1} into nbrr disjoint sets.
00073 //  The graph of the sparse matrix, with nrr vertices, may be compressed
00074 //  into the graph of the blocks, a graph with nbrr vertices, that is
00075 //  called here the block connectivity graph.
00076 //     The partition of the rows and columns of B is represented by
00077 //  r(0:nbrr),  0 = r(0) < r(1) < .. < r(nbrr) = nrr,
00078 //  The graph (Mp,Mj) of the nbrr x nbrr matrix is represened by
00079 //  a sparse matrix in sparse coordinate format.
00080 //  Mp: row indices, dimension determined here (nzM).
00081 //  Mj: column indices, dimension determined here (nzM).
00082 //  The integer vector, weights, of block sizes  (dimension nbrr) is also
00083 //  computed here.
00084 //     The case of nbrr proportional to nrr is critical.  One must
00085 //  efficiently look up the column indices of B in the partition.
00086 //  This is done here using a binary search tree, so that the
00087 //  look up cost is nzB*log2(nbrr).
00088 
00089   template<typename int_type>
00090   Teuchos::RCP<Epetra_CrsGraph> BlockAdjacencyGraph::compute( Epetra_CrsGraph& B, int nbrr, std::vector<int>&r, std::vector<double>& weights, bool verbose)
00091   {
00092     // Check if the graph is on one processor.
00093     int myMatProc = -1, matProc = -1;
00094     int myPID = B.Comm().MyPID();
00095     for (int proc=0; proc<B.Comm().NumProc(); proc++)
00096       {
00097         if (B.NumGlobalEntries64() == B.NumMyEntries())
00098           myMatProc = myPID;
00099       }
00100     B.Comm().MaxAll( &myMatProc, &matProc, 1 );
00101 
00102     if( matProc == -1)
00103       { std::cout << "FAIL for Global!  All CrsGraph entries must be on one processor!\n"; abort(); }
00104 
00105     int i= 0, j = 0, k, l = 0, p, pm, q = -1, ns;
00106     int tree_height;
00107     int error = -1;    /* error detected, possibly a problem with the input */
00108     (void) error; // silence "set but not used" warning
00109     int nrr;           /* number of rows in B */
00110     int nzM = 0;       /* number of edges in graph */
00111     int m = 0;         /* maximum number of nonzeros in any block row of B */
00112     int* colstack = 0; /* stack used to process each block row */
00113     int* bstree = 0;   /* binary search tree */
00114     std::vector<int_type> Mi, Mj, Mnum(nbrr+1,0);
00115     nrr = B.NumMyRows();
00116     if ( matProc == myPID && verbose )
00117       std::printf(" Matrix Size = %d      Number of Blocks = %d\n",nrr, nbrr);
00118     else
00119       nrr = -1;     /* Prevent processor from doing any computations */
00120     bstree = csr_bst(nbrr);  /* 0 : nbrr-1 */
00121     tree_height = ceil31log2(nbrr) + 1;
00122     error = -1;
00123 
00124     l = 0; j = 0; m = 0;
00125     for( i = 0; i < nrr; i++ ){
00126       if( i >= r[l+1] ){
00127         ++l;                 /* new block row */
00128         m = EPETRA_MAX(m,j) ;   /* nonzeros in block row */
00129         j = B.NumGlobalIndices(i);
00130       }else{
00131         j += B.NumGlobalIndices(i);
00132       }
00133     }
00134     /* one more time for the final block */
00135      m = EPETRA_MAX(m,j) ;   /* nonzeros in block row */
00136 
00137     colstack = (int*) malloc( EPETRA_MAX(m,1) * sizeof(int) );
00138     // The compressed graph is actually computed twice,
00139     // due to concerns about memory limitations.  First,
00140     // without memory allocation, just nzM is computed.
00141     // Next Mj is allocated. Then, the second time, the
00142     // arrays are actually populated.
00143     nzM = 0; q = -1; l = 0;
00144     int * indices;
00145     int numEntries;
00146     for( i = 0; i <= nrr; i++ ){
00147       if( i >= r[l+1] ){
00148         if( q > 0 ) std::qsort(colstack,q+1,sizeof(int),compare_ints); /* sort stack */
00149         if( q >= 0 ) ns = 1; /* l, colstack[0] M */
00150         for( j=1; j<=q ; j++ ){ /* delete copies */
00151           if( colstack[j] > colstack[j-1] ) ++ns;
00152         }
00153         nzM += ns; /*M->p[l+1] = M->p[l] + ns;*/
00154         ++l;
00155         q = -1;
00156       }
00157       if( i < nrr ){
00158         B.ExtractMyRowView( i, numEntries, indices );
00159         for( k = 0; k < numEntries; k++){
00160           j = indices[k];  ns = 0; p = 0;
00161           while( (r[bstree[p]] > j)  ||  (j >= r[bstree[p]+1])  ){
00162             if( r[bstree[p]] > j){
00163               p = 2*p+1;
00164             }else{
00165               if( r[bstree[p]+1] <= j) p = 2*p+2;
00166             }
00167             ++ns;
00168             if( p > nbrr || ns > tree_height ) {
00169               error = j;
00170               std::printf("error: p %d  nbrr %d  ns %d %d\n",p,nbrr,ns,j); break;
00171             }
00172           }
00173           colstack[++q] = bstree[p];
00174         }
00175         //if( error >-1 ){ std::printf("%d\n",error); break; }
00176         // p > nbrr is a fatal error that is ignored
00177       }
00178     }
00179 
00180     if ( matProc == myPID && verbose )
00181       std::printf("nzM =  %d \n", nzM );
00182     Mi.resize( nzM );
00183     Mj.resize( nzM );
00184     q = -1; l = 0; pm = -1;
00185     for( i = 0; i <= nrr; i++ ){
00186       if( i >= r[l+1] ){
00187         if( q > 0 ) std::qsort(colstack,q+1,sizeof(colstack[0]),compare_ints); /* sort stack */
00188         if( q >= 0 ){
00189           Mi[++pm] = l;
00190           Mj[pm] = colstack[0];
00191         }
00192         for( j=1; j<=q ; j++ ){ /* delete copies */
00193           if( colstack[j] > colstack[j-1] ){ /* l, colstack[j] */
00194             Mi[++pm] = l;
00195             Mj[pm] = colstack[j];
00196           }
00197         }
00198         ++l;
00199         Mnum[l] = pm + 1;
00200 
00201         /* sparse row format: M->p[l+1] = M->p[l] + ns; */
00202         q = -1;
00203       }
00204       if( i < nrr ){
00205         B.ExtractMyRowView( i, numEntries, indices );
00206         for( k = 0; k < numEntries; k++){
00207           j = indices[k]; ns = 0; p = 0;
00208           while( (r[bstree[p]] > j)  ||  (j >= r[bstree[p]+1])  ){
00209             if( r[bstree[p]] > j){
00210               p = 2*p+1;
00211             }else{
00212               if( r[bstree[p]+1] <= j) p = 2*p+2;
00213             }
00214             ++ns;
00215           }
00216           colstack[++q] = bstree[p];
00217         }
00218       }
00219     }
00220     if ( bstree ) free ( bstree );
00221     if ( colstack ) free( colstack );
00222 
00223     // Compute weights as number of rows in each block.
00224     weights.resize( nbrr );
00225     for( l=0; l<nbrr; l++) weights[l] = r[l+1] - r[l];
00226 
00227     // Compute Epetra_CrsGraph and return
00228     Teuchos::RCP<Epetra_Map> newMap;
00229     if ( matProc == myPID )
00230       newMap = Teuchos::rcp( new Epetra_Map(nbrr, nbrr, 0, B.Comm() ) );
00231     else
00232       newMap = Teuchos::rcp( new Epetra_Map( nbrr, 0, 0, B.Comm() ) );
00233     Teuchos::RCP<Epetra_CrsGraph> newGraph = Teuchos::rcp( new Epetra_CrsGraph( Copy, *newMap, 0 ) );
00234     for( l=0; l<newGraph->NumMyRows(); l++) {
00235       newGraph->InsertGlobalIndices( l, Mnum[l+1]-Mnum[l], &Mj[Mnum[l]] );
00236     }
00237     newGraph->FillComplete();
00238 
00239     return (newGraph);
00240   }
00241 
00242   Teuchos::RCP<Epetra_CrsGraph> BlockAdjacencyGraph::compute( Epetra_CrsGraph& B, int nbrr, std::vector<int>&r, std::vector<double>& weights, bool verbose) {
00243 #ifndef EPETRA_NO_32BIT_GLOBAL_INDICES
00244   if(B.RowMap().GlobalIndicesInt()) {
00245     return compute<int>(B, nbrr, r, weights, verbose);
00246   }
00247   else
00248 #endif
00249 #ifndef EPETRA_NO_64BIT_GLOBAL_INDICES
00250   if(B.RowMap().GlobalIndicesLongLong()) {
00251     return compute<long long>(B, nbrr, r, weights, verbose);
00252   }
00253   else
00254 #endif
00255     throw "EpetraExt::BlockAdjacencyGraph::compute: GlobalIndices type unknown";
00256   }
00257 
00258   /*
00259    * bst(n) returns the complete binary tree, stored in an integer array of dimension n.
00260    * index i has children 2i+1 and 2i+2, root index 0.
00261    * A binary search of a sorted n-array: find array(k)
00262    * i=0; k = tree(i);
00263    * if < array( k ), i = 2i+1;
00264    * elif > array( k ), i = 2i+2;
00265    * otherwise exit
00266    */
00267   int* BlockAdjacencyGraph::csr_bst( int n )
00268   {
00269     int i, l=1, m, nstack = 0, nexp=0, os, *array, *stack;
00270     int max_nstack = 0;
00271     if( n == 0 ) return(NULL);
00272     while( l <= n ){
00273       l = 2*l;
00274       ++nexp;
00275     } /* n < l = 2^nexp */
00276     array = (int *) malloc( n * sizeof(int) );
00277     stack = (int *) malloc(3*nexp * sizeof(int) );
00278     stack[3*nstack] = 0; stack[3*nstack+1] = 0; stack[3*nstack+2] = n;
00279     ++nstack;
00280     /*if( debug ) std::printf("stack : %d %d %d\n", stack[0] , stack[1], stack[2] );*/
00281     while( nstack > 0 ){
00282       --nstack;
00283       i = stack[3*nstack]; os = stack[3*nstack+1]; m = stack[3*nstack+2];
00284       array[i] = csr_bstrootindex(m) + os; /* 5 */
00285       if( 2*i+2 < n){   /* right child                4     3      1                      3 - 5 - 1     */
00286         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 */
00287         ++nstack;
00288       }
00289       if( 2*i+1 < n){   /* left  child */
00290         stack[3*nstack] = 2*i+1; stack[3*nstack+1] = os; stack[3*nstack+2] = array[i] - os; /* 5 4 5 */
00291         ++nstack;
00292       }
00293       if( nstack > max_nstack ) max_nstack =  nstack;
00294     }
00295     free( stack );
00296     return(array);
00297   }
00298 
00299   /*
00300     Given a complete binary search tree with n nodes, return
00301     the index of the root node.  A nodeless tree, n=0, has no
00302     root node (-1).  Nodes are indexed from 0 to n-1.
00303   */
00304   int BlockAdjacencyGraph::csr_bstrootindex( int n )
00305   {
00306     int l = 1, nexp = 0, i;
00307     if( n == 0) return(-1);
00308     while( l <= n ){
00309       l = 2*l;
00310       ++nexp;
00311     } /* n < l = 2^nexp */
00312     i = l/2 - 1;
00313     if(n<4) return(i);
00314     if (n-i < l/4)
00315       return( n - l/4  );
00316     else
00317       return(i);
00318   }
00319 
00320 } //namespace EpetraExt
00321 
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Defines