EpetraExt_BlockAdjacencyGraph.cpp

Go to the documentation of this file.
00001 //@HEADER
00002 // ***********************************************************************
00003 // 
00004 //     EpetraExt: Epetra Extended - Linear Algebra Services Package
00005 //                 Copyright (2001) Sandia Corporation
00006 // 
00007 // Under terms of Contract DE-AC04-94AL85000, there is a non-exclusive
00008 // license for use of this work by or on behalf of the U.S. Government.
00009 // 
00010 // This library is free software; you can redistribute it and/or modify
00011 // it under the terms of the GNU Lesser General Public License as
00012 // published by the Free Software Foundation; either version 2.1 of the
00013 // License, or (at your option) any later version.
00014 //  
00015 // This library is distributed in the hope that it will be useful, but
00016 // WITHOUT ANY WARRANTY; without even the implied warranty of
00017 // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
00018 // Lesser General Public License for more details.
00019 //  
00020 // You should have received a copy of the GNU Lesser General Public
00021 // License along with this library; if not, write to the Free Software
00022 // Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307
00023 // USA
00024 // Questions? Contact Michael A. Heroux (maherou@sandia.gov) 
00025 // 
00026 // ***********************************************************************
00027 //@HEADER
00028   
00029 #include <EpetraExt_BlockAdjacencyGraph.h>
00030 
00031 #include <Epetra_CrsMatrix.h>
00032 #include <Epetra_CrsGraph.h>
00033 #include <Epetra_Map.h>
00034 #include <Epetra_Comm.h>
00035 
00036 #include <math.h>
00037 #include <stdlib.h>
00038 
00039 namespace EpetraExt {
00040 
00041   int compare_ints(const void *a, const void *b)
00042   {
00043     return(*((int *) a) - *((int *) b));
00044   }
00045 
00046   int ceil31log2(int n)
00047   { // Given 1 <= n < 2^31, find l such that 2^(l-1) < n <= 2^(l)
00048    int l=0, m = 1;
00049    while( n > m & l < 31 )
00050    {
00051       m = 2*m;
00052       ++l;
00053    }
00054    return(l);
00055   }
00056 //  Purpose: Compute the block connectivity graph of a matrix.
00057 //  An nrr by nrr sparse matrix admits a (Dulmage-Mendelsohn)
00058 //  permutation to block triangular form with nbrr blocks.
00059 //  Abstractly, the block triangular form corresponds to a partition
00060 //  of the set of integers {0,...,n-1} into nbrr disjoint sets.
00061 //  The graph of the sparse matrix, with nrr vertices, may be compressed
00062 //  into the graph of the blocks, a graph with nbrr vertices, that is
00063 //  called here the block connectivity graph.
00064 //     The partition of the rows and columns of B is represented by
00065 //  r(0:nbrr),  0 = r(0) < r(1) < .. < r(nbrr) = nrr, 
00066 //  The graph (Mp,Mj) of the nbrr x nbrr matrix is represened by
00067 //  a sparse matrix in sparse coordinate format.
00068 //  Mp: row indices, dimension determined here (nzM).
00069 //  Mj: column indices, dimension determined here (nzM).
00070 //  The integer vector, weights, of block sizes  (dimension nbrr) is also
00071 //  computed here.  
00072 //     The case of nbrr proportional to nrr is critical.  One must
00073 //  efficiently look up the column indices of B in the partition.
00074 //  This is done here using a binary search tree, so that the
00075 //  look up cost is nzB*log2(nbrr).  
00076 
00077   Teuchos::RCP<Epetra_CrsGraph> BlockAdjacencyGraph::compute( Epetra_CrsGraph& B, int nbrr, std::vector<int>&r, std::vector<double>& weights, bool verbose)
00078   {
00079     // Check if the graph is on one processor.
00080     int myMatProc = -1, matProc = -1;
00081     int myPID = B.Comm().MyPID();
00082     for (int proc=0; proc<B.Comm().NumProc(); proc++)
00083       {
00084   if (B.NumGlobalEntries() == B.NumMyEntries())
00085     myMatProc = myPID;
00086       }
00087     B.Comm().MaxAll( &myMatProc, &matProc, 1 );
00088     
00089     if( matProc == -1)
00090       { cout << "FAIL for Global!  All CrsGraph entries must be on one processor!\n"; abort(); }
00091     
00092     int i= 0, j = 0, k, l = 0, p, pm, q = -1, ns;
00093     int tree_height;
00094     int error = -1;    /* error detected, possibly a problem with the input */
00095     int nrr;           /* number of rows in B */
00096     int nzM = 0;       /* number of edges in graph */
00097     int m = 0;         /* maximum number of nonzeros in any block row of B */
00098     int* colstack = 0; /* stack used to process each block row */
00099     int* bstree = 0;   /* binary search tree */
00100     std::vector<int> Mi, Mj, Mnum(nbrr+1,0);
00101     nrr = B.NumMyRows();
00102     if ( matProc == myPID && verbose )
00103       std::printf(" Matrix Size = %d      Number of Blocks = %d\n",nrr, nbrr);
00104     else
00105       nrr = -1;     /* Prevent processor from doing any computations */
00106     bstree = csr_bst(nbrr);  /* 0 : nbrr-1 */
00107     tree_height = ceil31log2(nbrr) + 1;
00108     error = -1;
00109 
00110     l = 0; j = 0; m = 0;
00111     for( i = 0; i < nrr; i++ ){
00112       if( i >= r[l+1] ){
00113   ++l;                 /* new block row */
00114   m = EPETRA_MAX(m,j) ;   /* nonzeros in block row */
00115   j = B.NumGlobalIndices(i);
00116       }else{
00117   j += B.NumGlobalIndices(i);
00118       }
00119     }
00120     /* one more time for the final block */
00121      m = EPETRA_MAX(m,j) ;   /* nonzeros in block row */
00122 
00123     colstack = (int*) malloc( EPETRA_MAX(m,1) * sizeof(int) );
00124     // The compressed graph is actually computed twice,
00125     // due to concerns about memory limitations.  First, 
00126     // without memory allocation, just nzM is computed.  
00127     // Next Mj is allocated. Then, the second time, the
00128     // arrays are actually populated.
00129     nzM = 0; q = -1; l = 0;
00130     int * indices;
00131     int numEntries;
00132     for( i = 0; i <= nrr; i++ ){
00133       if( i >= r[l+1] ){
00134   if( q > 0 ) qsort(colstack,q+1,sizeof(int),compare_ints); /* sort stack */
00135   if( q >= 0 ) ns = 1; /* l, colstack[0] M */
00136   for( j=1; j<=q ; j++ ){ /* delete copies */
00137     if( colstack[j] > colstack[j-1] ) ++ns;
00138   }
00139   nzM += ns; /*M->p[l+1] = M->p[l] + ns;*/
00140   ++l;
00141   q = -1;
00142       }
00143       if( i < nrr ){
00144   B.ExtractMyRowView( i, numEntries, indices );
00145   for( k = 0; k < numEntries; k++){
00146     j = indices[k];  ns = 0; p = 0;
00147     while( (r[bstree[p]] > j)  ||  (j >= r[bstree[p]+1])  ){
00148       if( r[bstree[p]] > j){
00149         p = 2*p+1;
00150       }else{
00151         if( r[bstree[p]+1] <= j) p = 2*p+2;
00152       }
00153       ++ns;
00154       if( p > nbrr || ns > tree_height ) {
00155         error = j;
00156         std::printf("error: p %d  nbrr %d  ns %d %d\n",p,nbrr,ns,j); break;
00157       }
00158     }
00159     colstack[++q] = bstree[p];
00160   }
00161   //if( error >-1 ){ std::printf("%d\n",error); break; }
00162         // p > nbrr is a fatal error that is ignored
00163       }
00164     }
00165     
00166     if ( matProc == myPID && verbose )
00167       std::printf("nzM =  %d \n", nzM );
00168     Mi.resize( nzM );
00169     Mj.resize( nzM );
00170     q = -1; l = 0; pm = -1;
00171     for( i = 0; i <= nrr; i++ ){
00172       if( i >= r[l+1] ){
00173   if( q > 0 ) qsort(colstack,q+1,sizeof(colstack[0]),compare_ints); /* sort stack */
00174   if( q >= 0 ){
00175     Mi[++pm] = l;
00176     Mj[pm] = colstack[0];
00177   }
00178   for( j=1; j<=q ; j++ ){ /* delete copies */
00179     if( colstack[j] > colstack[j-1] ){ /* l, colstack[j] */
00180       Mi[++pm] = l;
00181       Mj[pm] = colstack[j];
00182     }
00183   }
00184   ++l;
00185   Mnum[l] = pm + 1;
00186   
00187   /* sparse row format: M->p[l+1] = M->p[l] + ns; */
00188   q = -1;
00189       }
00190       if( i < nrr ){
00191   B.ExtractMyRowView( i, numEntries, indices );
00192   for( k = 0; k < numEntries; k++){
00193     j = indices[k]; ns = 0; p = 0;
00194     while( (r[bstree[p]] > j)  ||  (j >= r[bstree[p]+1])  ){
00195       if( r[bstree[p]] > j){
00196         p = 2*p+1;
00197       }else{
00198         if( r[bstree[p]+1] <= j) p = 2*p+2;
00199       }
00200       ++ns;
00201     }
00202     colstack[++q] = bstree[p];
00203   }
00204       }
00205     }
00206     if ( bstree ) free ( bstree );
00207     if ( colstack ) free( colstack );
00208     
00209     // Compute weights as number of rows in each block.
00210     weights.resize( nbrr );
00211     for( l=0; l<nbrr; l++) weights[l] = r[l+1] - r[l];
00212     
00213     // Compute Epetra_CrsGraph and return
00214     Teuchos::RCP<Epetra_Map> newMap;
00215     if ( matProc == myPID )
00216       newMap = Teuchos::rcp( new Epetra_Map(nbrr, nbrr, 0, B.Comm() ) );
00217     else
00218       newMap = Teuchos::rcp( new Epetra_Map( nbrr, 0, 0, B.Comm() ) );
00219     Teuchos::RCP<Epetra_CrsGraph> newGraph = Teuchos::rcp( new Epetra_CrsGraph( Copy, *newMap, 0 ) );
00220     for( l=0; l<newGraph->NumMyRows(); l++) {
00221       newGraph->InsertGlobalIndices( l, Mnum[l+1]-Mnum[l], &Mj[Mnum[l]] );
00222     }
00223     newGraph->FillComplete();
00224     
00225     return (newGraph);  
00226   }
00227   
00228   /*
00229    * bst(n) returns the complete binary tree, stored in an integer array of dimension n.
00230    * index i has children 2i+1 and 2i+2, root index 0.
00231    * A binary search of a sorted n-array: find array(k)
00232    * i=0; k = tree(i);
00233    * if < array( k ), i = 2i+1;
00234    * elif > array( k ), i = 2i+2;
00235    * otherwise exit
00236    */
00237   int* BlockAdjacencyGraph::csr_bst( int n )
00238   {
00239     int i, l=1, m, nstack = 0, nexp=0, os, *array, *stack;
00240     int max_nstack = 0;
00241     if( n == 0 ) return(NULL);
00242     while( l <= n ){
00243       l = 2*l;
00244       ++nexp;
00245     } /* n < l = 2^nexp */
00246     array = (int *) malloc( n * sizeof(int) );
00247     stack = (int *) malloc(3*nexp * sizeof(int) );
00248     stack[3*nstack] = 0; stack[3*nstack+1] = 0; stack[3*nstack+2] = n;
00249     ++nstack;
00250     /*if( debug ) std::printf("stack : %d %d %d\n", stack[0] , stack[1], stack[2] );*/
00251     while( nstack > 0 ){
00252       --nstack;
00253       i = stack[3*nstack]; os = stack[3*nstack+1]; m = stack[3*nstack+2];
00254       array[i] = csr_bstrootindex(m) + os; /* 5 */
00255       if( 2*i+2 < n){   /* right child                4     3      1                      3 - 5 - 1     */
00256   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 */
00257   ++nstack;
00258       }
00259       if( 2*i+1 < n){   /* left  child */
00260   stack[3*nstack] = 2*i+1; stack[3*nstack+1] = os; stack[3*nstack+2] = array[i] - os; /* 5 4 5 */
00261   ++nstack;
00262       }
00263       if( nstack > max_nstack ) max_nstack =  nstack;
00264     }
00265     free( stack );
00266     return(array);
00267   }
00268 
00269   /*
00270     Given a complete binary search tree with n nodes, return
00271     the index of the root node.  A nodeless tree, n=0, has no
00272     root node (-1).  Nodes are indexed from 0 to n-1.
00273   */
00274   int BlockAdjacencyGraph::csr_bstrootindex( int n )
00275   {
00276     int l = 1, nexp = 0, i;
00277     if( n == 0) return(-1);
00278     while( l <= n ){
00279       l = 2*l;
00280       ++nexp;
00281     } /* n < l = 2^nexp */
00282     i = l/2 - 1;
00283     if(n<4) return(i);
00284     if (n-i < l/4)
00285       return( n - l/4  );
00286     else
00287       return(i);  
00288   }
00289   
00290 } //namespace EpetraExt
00291 

Generated on Tue Jul 13 09:23:05 2010 for EpetraExt by  doxygen 1.4.7