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    {
00064       m = 2*m;
00065       ++l;
00066    }
00067    return(l);
00068   }
00069 //  Purpose: Compute the block connectivity graph of a matrix.
00070 //  An nrr by nrr sparse matrix admits a (Dulmage-Mendelsohn)
00071 //  permutation to block triangular form with nbrr blocks.
00072 //  Abstractly, the block triangular form corresponds to a partition
00073 //  of the set of integers {0,...,n-1} into nbrr disjoint sets.
00074 //  The graph of the sparse matrix, with nrr vertices, may be compressed
00075 //  into the graph of the blocks, a graph with nbrr vertices, that is
00076 //  called here the block connectivity graph.
00077 //     The partition of the rows and columns of B is represented by
00078 //  r(0:nbrr),  0 = r(0) < r(1) < .. < r(nbrr) = nrr, 
00079 //  The graph (Mp,Mj) of the nbrr x nbrr matrix is represened by
00080 //  a sparse matrix in sparse coordinate format.
00081 //  Mp: row indices, dimension determined here (nzM).
00082 //  Mj: column indices, dimension determined here (nzM).
00083 //  The integer vector, weights, of block sizes  (dimension nbrr) is also
00084 //  computed here.  
00085 //     The case of nbrr proportional to nrr is critical.  One must
00086 //  efficiently look up the column indices of B in the partition.
00087 //  This is done here using a binary search tree, so that the
00088 //  look up cost is nzB*log2(nbrr).  
00089 
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.NumGlobalEntries() == B.NumMyEntries())
00098     myMatProc = myPID;
00099       }
00100     B.Comm().MaxAll( &myMatProc, &matProc, 1 );
00101     
00102     if( matProc == -1)
00103       { 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     int nrr;           /* number of rows in B */
00109     int nzM = 0;       /* number of edges in graph */
00110     int m = 0;         /* maximum number of nonzeros in any block row of B */
00111     int* colstack = 0; /* stack used to process each block row */
00112     int* bstree = 0;   /* binary search tree */
00113     std::vector<int> Mi, Mj, Mnum(nbrr+1,0);
00114     nrr = B.NumMyRows();
00115     if ( matProc == myPID && verbose )
00116       std::printf(" Matrix Size = %d      Number of Blocks = %d\n",nrr, nbrr);
00117     else
00118       nrr = -1;     /* Prevent processor from doing any computations */
00119     bstree = csr_bst(nbrr);  /* 0 : nbrr-1 */
00120     tree_height = ceil31log2(nbrr) + 1;
00121     error = -1;
00122 
00123     l = 0; j = 0; m = 0;
00124     for( i = 0; i < nrr; i++ ){
00125       if( i >= r[l+1] ){
00126   ++l;                 /* new block row */
00127   m = EPETRA_MAX(m,j) ;   /* nonzeros in block row */
00128   j = B.NumGlobalIndices(i);
00129       }else{
00130   j += B.NumGlobalIndices(i);
00131       }
00132     }
00133     /* one more time for the final block */
00134      m = EPETRA_MAX(m,j) ;   /* nonzeros in block row */
00135 
00136     colstack = (int*) malloc( EPETRA_MAX(m,1) * sizeof(int) );
00137     // The compressed graph is actually computed twice,
00138     // due to concerns about memory limitations.  First, 
00139     // without memory allocation, just nzM is computed.  
00140     // Next Mj is allocated. Then, the second time, the
00141     // arrays are actually populated.
00142     nzM = 0; q = -1; l = 0;
00143     int * indices;
00144     int numEntries;
00145     for( i = 0; i <= nrr; i++ ){
00146       if( i >= r[l+1] ){
00147   if( q > 0 ) std::qsort(colstack,q+1,sizeof(int),compare_ints); /* sort stack */
00148   if( q >= 0 ) ns = 1; /* l, colstack[0] M */
00149   for( j=1; j<=q ; j++ ){ /* delete copies */
00150     if( colstack[j] > colstack[j-1] ) ++ns;
00151   }
00152   nzM += ns; /*M->p[l+1] = M->p[l] + ns;*/
00153   ++l;
00154   q = -1;
00155       }
00156       if( i < nrr ){
00157   B.ExtractMyRowView( i, numEntries, indices );
00158   for( k = 0; k < numEntries; k++){
00159     j = indices[k];  ns = 0; p = 0;
00160     while( (r[bstree[p]] > j)  ||  (j >= r[bstree[p]+1])  ){
00161       if( r[bstree[p]] > j){
00162         p = 2*p+1;
00163       }else{
00164         if( r[bstree[p]+1] <= j) p = 2*p+2;
00165       }
00166       ++ns;
00167       if( p > nbrr || ns > tree_height ) {
00168         error = j;
00169         std::printf("error: p %d  nbrr %d  ns %d %d\n",p,nbrr,ns,j); break;
00170       }
00171     }
00172     colstack[++q] = bstree[p];
00173   }
00174   //if( error >-1 ){ std::printf("%d\n",error); break; }
00175         // p > nbrr is a fatal error that is ignored
00176       }
00177     }
00178     
00179     if ( matProc == myPID && verbose )
00180       std::printf("nzM =  %d \n", nzM );
00181     Mi.resize( nzM );
00182     Mj.resize( nzM );
00183     q = -1; l = 0; pm = -1;
00184     for( i = 0; i <= nrr; i++ ){
00185       if( i >= r[l+1] ){
00186   if( q > 0 ) std::qsort(colstack,q+1,sizeof(colstack[0]),compare_ints); /* sort stack */
00187   if( q >= 0 ){
00188     Mi[++pm] = l;
00189     Mj[pm] = colstack[0];
00190   }
00191   for( j=1; j<=q ; j++ ){ /* delete copies */
00192     if( colstack[j] > colstack[j-1] ){ /* l, colstack[j] */
00193       Mi[++pm] = l;
00194       Mj[pm] = colstack[j];
00195     }
00196   }
00197   ++l;
00198   Mnum[l] = pm + 1;
00199   
00200   /* sparse row format: M->p[l+1] = M->p[l] + ns; */
00201   q = -1;
00202       }
00203       if( i < nrr ){
00204   B.ExtractMyRowView( i, numEntries, indices );
00205   for( k = 0; k < numEntries; k++){
00206     j = indices[k]; ns = 0; p = 0;
00207     while( (r[bstree[p]] > j)  ||  (j >= r[bstree[p]+1])  ){
00208       if( r[bstree[p]] > j){
00209         p = 2*p+1;
00210       }else{
00211         if( r[bstree[p]+1] <= j) p = 2*p+2;
00212       }
00213       ++ns;
00214     }
00215     colstack[++q] = bstree[p];
00216   }
00217       }
00218     }
00219     if ( bstree ) free ( bstree );
00220     if ( colstack ) free( colstack );
00221     
00222     // Compute weights as number of rows in each block.
00223     weights.resize( nbrr );
00224     for( l=0; l<nbrr; l++) weights[l] = r[l+1] - r[l];
00225     
00226     // Compute Epetra_CrsGraph and return
00227     Teuchos::RCP<Epetra_Map> newMap;
00228     if ( matProc == myPID )
00229       newMap = Teuchos::rcp( new Epetra_Map(nbrr, nbrr, 0, B.Comm() ) );
00230     else
00231       newMap = Teuchos::rcp( new Epetra_Map( nbrr, 0, 0, B.Comm() ) );
00232     Teuchos::RCP<Epetra_CrsGraph> newGraph = Teuchos::rcp( new Epetra_CrsGraph( Copy, *newMap, 0 ) );
00233     for( l=0; l<newGraph->NumMyRows(); l++) {
00234       newGraph->InsertGlobalIndices( l, Mnum[l+1]-Mnum[l], &Mj[Mnum[l]] );
00235     }
00236     newGraph->FillComplete();
00237     
00238     return (newGraph);  
00239   }
00240   
00241   /*
00242    * bst(n) returns the complete binary tree, stored in an integer array of dimension n.
00243    * index i has children 2i+1 and 2i+2, root index 0.
00244    * A binary search of a sorted n-array: find array(k)
00245    * i=0; k = tree(i);
00246    * if < array( k ), i = 2i+1;
00247    * elif > array( k ), i = 2i+2;
00248    * otherwise exit
00249    */
00250   int* BlockAdjacencyGraph::csr_bst( int n )
00251   {
00252     int i, l=1, m, nstack = 0, nexp=0, os, *array, *stack;
00253     int max_nstack = 0;
00254     if( n == 0 ) return(NULL);
00255     while( l <= n ){
00256       l = 2*l;
00257       ++nexp;
00258     } /* n < l = 2^nexp */
00259     array = (int *) malloc( n * sizeof(int) );
00260     stack = (int *) malloc(3*nexp * sizeof(int) );
00261     stack[3*nstack] = 0; stack[3*nstack+1] = 0; stack[3*nstack+2] = n;
00262     ++nstack;
00263     /*if( debug ) std::printf("stack : %d %d %d\n", stack[0] , stack[1], stack[2] );*/
00264     while( nstack > 0 ){
00265       --nstack;
00266       i = stack[3*nstack]; os = stack[3*nstack+1]; m = stack[3*nstack+2];
00267       array[i] = csr_bstrootindex(m) + os; /* 5 */
00268       if( 2*i+2 < n){   /* right child                4     3      1                      3 - 5 - 1     */
00269   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 */
00270   ++nstack;
00271       }
00272       if( 2*i+1 < n){   /* left  child */
00273   stack[3*nstack] = 2*i+1; stack[3*nstack+1] = os; stack[3*nstack+2] = array[i] - os; /* 5 4 5 */
00274   ++nstack;
00275       }
00276       if( nstack > max_nstack ) max_nstack =  nstack;
00277     }
00278     free( stack );
00279     return(array);
00280   }
00281 
00282   /*
00283     Given a complete binary search tree with n nodes, return
00284     the index of the root node.  A nodeless tree, n=0, has no
00285     root node (-1).  Nodes are indexed from 0 to n-1.
00286   */
00287   int BlockAdjacencyGraph::csr_bstrootindex( int n )
00288   {
00289     int l = 1, nexp = 0, i;
00290     if( n == 0) return(-1);
00291     while( l <= n ){
00292       l = 2*l;
00293       ++nexp;
00294     } /* n < l = 2^nexp */
00295     i = l/2 - 1;
00296     if(n<4) return(i);
00297     if (n-i < l/4)
00298       return( n - l/4  );
00299     else
00300       return(i);  
00301   }
00302   
00303 } //namespace EpetraExt
00304 
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Defines