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 /*
00047  * Compute the block connectivity graph of an Epetra_CrsGraph.
00048  * Given a sparse graph B, and a partition of the rows and
00049  * columns of B,  r(0:nbrr),
00050  * 0 = r(0) < r(1) < .. < r(nbrr) = B->n,
00051  * return the graph of the corresponding nbrr x nbrr matrix.
00052  * Weights has dimension nbrr and is determined here too.
00053  *
00054  */
00055 
00056   Teuchos::RCP<Epetra_CrsGraph> BlockAdjacencyGraph::compute( Epetra_CrsGraph& B, int nbrr, std::vector<int>&r, std::vector<double>& weights)
00057   {
00058     // Check if the graph is on one processor.
00059     int myMatProc = -1, matProc = -1;
00060     int myPID = B.Comm().MyPID();
00061     for (int proc=0; proc<B.Comm().NumProc(); proc++)
00062       {
00063   if (B.NumGlobalEntries() == B.NumMyEntries())
00064     myMatProc = myPID;
00065       }
00066     B.Comm().MaxAll( &myMatProc, &matProc, 1 );
00067     
00068     if( matProc == -1)
00069       { cout << "FAIL for Global!  All CrsGraph entries must be on one processor!\n"; abort(); }
00070     
00071     int i= 0, j = 0, k, l = 0, p, pm, q = -1, ns, log2nbrr = 20;
00072     int error = -1;    /* error detected, possibly a problem with the input */
00073     int nrr;           /* number of rows in B */
00074     int nzM = 0;       /* number of edges in graph */
00075     int m = 0;         /* maximum number of nonzeros in any block row of B */
00076     int* colstack = 0; /* stack used to process each block row */
00077     int* bstree = 0;   /* binary search tree */
00078     std::vector<int> Mi, Mj, Mnum(nbrr+1,0);
00079     nrr = B.NumMyRows();
00080     if ( matProc == myPID )
00081       std::printf(" nrr = %d      nbrr = %d\n",nrr, nbrr);
00082     else
00083       nrr = -1;     /* Prevent processor from doing any computations */
00084     bstree = csr_bst(nbrr);  /* 0 : nbrr-1 */
00085     error = -1;
00086     l = 0; j = 0; m = 0;
00087     for( i = 0; i < nrr; i++ ){
00088       if( i >= r[l+1] ){
00089   ++l;                 /* new block row */
00090   m = EPETRA_MAX(m,j) ;   /* nonzeros in block row */
00091   j = B.NumGlobalIndices(i);
00092       }else{
00093   j += B.NumGlobalIndices(i);
00094       }
00095     }
00096     /* one more time for the final block */
00097      m = EPETRA_MAX(m,j) ;   /* nonzeros in block row */
00098 
00099     colstack = (int*) malloc( EPETRA_MAX(m,1) * sizeof(int) );
00100     nzM = 0; q = -1; l = 0;
00101     int * indices;
00102     int numEntries;
00103     for( i = 0; i <= nrr; i++ ){
00104       if( i >= r[l+1] ){
00105   if( q > 0 ) qsort(colstack,q+1,sizeof(int),compare_ints); /* sort stack */
00106   if( q >= 0 ) ns = 1; /* l, colstack[0] M */
00107   for( j=1; j<=q ; j++ ){ /* delete copies */
00108     if( colstack[j] > colstack[j-1] ) ++ns;
00109   }
00110   nzM += ns; /*M->p[l+1] = M->p[l] + ns;*/
00111   ++l;
00112   q = -1;
00113       }
00114       if( i < nrr ){
00115   B.ExtractMyRowView( i, numEntries, indices );
00116   for( k = 0; k < numEntries; k++){
00117     j = indices[k];  ns = 0; p = 0;
00118     while( (r[bstree[p]] > j)  ||  (j >= r[bstree[p]+1])  ){
00119       if( r[bstree[p]] > j){
00120         p = 2*p+1;
00121       }else{
00122         if( r[bstree[p]+1] <= j) p = 2*p+2;
00123       }
00124       ++ns;
00125       if( p > nbrr || ns > 20 ) {
00126         error = j;
00127         std::printf("error: p %d  nbrr %d  ns %d %d\n",p,nbrr,ns,j); break;
00128       }
00129     }
00130     colstack[++q] = bstree[p];
00131   }
00132   if( error >-1 ){ std::printf("%d\n",error); break; }
00133       }
00134     }
00135     
00136     if ( matProc == myPID )
00137       std::printf("nzM =  %d \n", nzM );
00138     Mi.resize( nzM );
00139     Mj.resize( nzM );
00140     q = -1; l = 0; pm = -1;
00141     for( i = 0; i <= nrr; i++ ){
00142       if( i >= r[l+1] ){
00143   if( q > 0 ) qsort(colstack,q+1,sizeof(colstack[0]),compare_ints); /* sort stack */
00144   if( q >= 0 ){
00145     Mi[++pm] = l;
00146     Mj[pm] = colstack[0];
00147   }
00148   for( j=1; j<=q ; j++ ){ /* delete copies */
00149     if( colstack[j] > colstack[j-1] ){ /* l, colstack[j] */
00150       Mi[++pm] = l;
00151       Mj[pm] = colstack[j];
00152     }
00153   }
00154   ++l;
00155   Mnum[l] = pm + 1;
00156   
00157   /* sparse row format: M->p[l+1] = M->p[l] + ns; */
00158   q = -1;
00159       }
00160       if( i < nrr ){
00161   B.ExtractMyRowView( i, numEntries, indices );
00162   for( k = 0; k < numEntries; k++){
00163     j = indices[k]; ns = 0; p = 0;
00164     while( (r[bstree[p]] > j)  ||  (j >= r[bstree[p]+1])  ){
00165       if( r[bstree[p]] > j){
00166         p = 2*p+1;
00167       }else{
00168         if( r[bstree[p]+1] <= j) p = 2*p+2;
00169       }
00170       ++ns;
00171       if( p > nbrr || ns > 14 ) {
00172         error = j;
00173         std::printf("error: p %d  nbrr %d  ns %d %d\n",p,nbrr,ns,j); break;
00174       }
00175     }
00176     colstack[++q] = bstree[p];
00177   }
00178   if( error >-1 ){ std::printf("%d\n",error); break; }
00179       }
00180     }
00181     if ( bstree ) free ( bstree );
00182     if ( colstack ) free( colstack );
00183     
00184     // Compute weights as number of rows in each block.
00185     weights.resize( nbrr );
00186     for( l=0; l<nbrr; l++) weights[l] = r[l+1] - r[l];
00187     
00188     // Compute Epetra_CrsGraph and return
00189     Teuchos::RCP<Epetra_Map> newMap;
00190     if ( matProc == myPID )
00191       newMap = Teuchos::rcp( new Epetra_Map(nbrr, nbrr, 0, B.Comm() ) );
00192     else
00193       newMap = Teuchos::rcp( new Epetra_Map( nbrr, 0, 0, B.Comm() ) );
00194     Teuchos::RCP<Epetra_CrsGraph> newGraph = Teuchos::rcp( new Epetra_CrsGraph( Copy, *newMap, 0 ) );
00195     for( l=0; l<newGraph->NumMyRows(); l++) {
00196       newGraph->InsertGlobalIndices( l, Mnum[l+1]-Mnum[l], &Mj[Mnum[l]] );
00197     }
00198     newGraph->FillComplete();
00199     
00200     return (newGraph);  
00201   }
00202   
00203   /*
00204    * bst(n) returns the complete binary tree, stored in an integer array of dimension n.
00205    * index i has children 2i+1 and 2i+2, root index 0.
00206    * A binary search of a sorted n-array: find array(k)
00207    * i=0; k = tree(i);
00208    * if < array( k ), i = 2i+1;
00209    * elif > array( k ), i = 2i+2;
00210    * otherwise exit
00211    */
00212   int* BlockAdjacencyGraph::csr_bst( int n )
00213   {
00214     int i, l=1, m, nstack = 0, nexp=0, os, *array, *stack;
00215     int max_nstack = 0;
00216     if( n == 0 ) return(NULL);
00217     while( l <= n ){
00218       l = 2*l;
00219       ++nexp;
00220     } /* n < l = 2^nexp */
00221     array = (int *) malloc( n * sizeof(int) );
00222     stack = (int *) malloc(3*nexp * sizeof(int) );
00223     stack[3*nstack] = 0; stack[3*nstack+1] = 0; stack[3*nstack+2] = n;
00224     ++nstack;
00225     /*if( debug ) std::printf("stack : %d %d %d\n", stack[0] , stack[1], stack[2] );*/
00226     while( nstack > 0 ){
00227       --nstack;
00228       i = stack[3*nstack]; os = stack[3*nstack+1]; m = stack[3*nstack+2];
00229       array[i] = csr_bstrootindex(m) + os; /* 5 */
00230       if( 2*i+2 < n){   /* right child                4     3      1                      3 - 5 - 1     */
00231   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 */
00232   ++nstack;
00233       }
00234       if( 2*i+1 < n){   /* left  child */
00235   stack[3*nstack] = 2*i+1; stack[3*nstack+1] = os; stack[3*nstack+2] = array[i] - os; /* 5 4 5 */
00236   ++nstack;
00237       }
00238       if( nstack > max_nstack ) max_nstack =  nstack;
00239     }
00240     free( stack );
00241     return(array);
00242   }
00243 
00244   /*
00245     Given a complete binary search tree with n nodes, return
00246     the index of the root node.  A nodeless tree, n=0, has no
00247     root node (-1).  Nodes are indexed from 0 to n-1.
00248   */
00249   int BlockAdjacencyGraph::csr_bstrootindex( int n )
00250   {
00251     int l = 1, nexp = 0, i;
00252     if( n == 0) return(-1);
00253     while( l <= n ){
00254       l = 2*l;
00255       ++nexp;
00256     } /* n < l = 2^nexp */
00257     i = l/2 - 1;
00258     if(n<4) return(i);
00259     if (n-i < l/4)
00260       return( n - l/4  );
00261     else
00262       return(i);  
00263   }
00264   
00265 } //namespace EpetraExt
00266 

Generated on Wed May 12 21:40:37 2010 for EpetraExt by  doxygen 1.4.7