EpetraExt Development
EpetraExt_BTF_CrsGraph.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 #include <EpetraExt_BTF_CrsGraph.h>
00029 
00030 #include <Epetra_Import.h>
00031 #include <Epetra_CrsGraph.h>
00032 #include <Epetra_Map.h>
00033 
00034 #include <vector>
00035 
00036 using std::vector;
00037 
00038 #define MATTRANS_F77 F77_FUNC(mattrans,MATTRANS)
00039 #define GENBTF_F77   F77_FUNC(genbtf,GENBTF)
00040 extern "C" {
00041 extern void MATTRANS_F77( int*, int*, int*, int*, int*, int* );
00042 extern void GENBTF_F77( int*, int*, int*, int*, int*, int*, int*, int*, int*,
00043                      int*, int*, int*, int*, int*, int*, int*, int*, int*,
00044                      int*, int*, int* );
00045 }
00046 
00047 namespace EpetraExt {
00048 
00049 CrsGraph_BTF::
00050 ~CrsGraph_BTF()
00051 {
00052   if( NewRowMap_ ) delete NewRowMap_;
00053   if( NewDomainMap_ ) delete NewDomainMap_;
00054 
00055   if( NewGraph_ ) delete NewGraph_;
00056 }
00057 
00058 CrsGraph_BTF::NewTypeRef
00059 CrsGraph_BTF::
00060 operator()( OriginalTypeRef orig )
00061 {
00062   origObj_ = &orig;
00063 
00064   if( orig.RowMap().DistributedGlobal() )
00065   { cout << "FAIL for Global!\n"; abort(); }
00066   if( orig.IndicesAreGlobal() )
00067   { cout << "FAIL for Global Indices!\n"; abort(); }
00068 
00069   int n = orig.NumMyRows();
00070   int nnz = orig.NumMyNonzeros();
00071 
00072   //create std CRS format
00073   vector<int> ia(n+1,0);
00074   vector<int> ja(nnz);
00075   int cnt;
00076   for( int i = 0; i < n; ++i )
00077   {
00078     int * tmpP = &ja[ia[i]];
00079     orig.ExtractMyRowCopy( i, nnz-ia[i], cnt, tmpP );
00080     ia[i+1] = ia[i] + cnt;
00081   }
00082 
00083   //convert to Fortran indexing
00084   for( int i = 0; i < n+1; ++i ) ++ia[i];
00085   for( int i = 0; i < nnz; ++i ) ++ja[i];
00086 
00087 #ifdef BTF_VERBOSE
00088   cout << "-----------------------------------------\n";
00089   cout << "CRS Format Graph\n";
00090   cout << "-----------------------------------------\n";
00091   for( int i = 0; i < n; ++i )
00092   {
00093     cout << i << ": " << ia[i+1] << ": ";
00094     for( int j = ia[i]-1; j<ia[i+1]-1; ++j )
00095       cout << " " << ja[j];
00096     cout << endl;
00097   }
00098   cout << "-----------------------------------------\n";
00099 #endif
00100 
00101   vector<int> iat(n+1);
00102   vector<int> jat(nnz);
00103   int * jaf = &ja[0];
00104   int * iaf = &ia[0];
00105   int * jatf = &jat[0];
00106   int * iatf = &iat[0];
00107   MATTRANS_F77( &n, &n, jaf, iaf, jatf, iatf );
00108     
00109 #ifdef BTF_VERBOSE
00110   cout << "-----------------------------------------\n";
00111   cout << "CCS Format Graph\n";
00112   cout << "-----------------------------------------\n";
00113   for( int i = 0; i < n; ++i )
00114   {
00115     cout << i << ": " << iat[i+1] << ": ";
00116     for( int j = iat[i]-1; j<iat[i+1]-1; ++j )
00117       cout << " " << jat[j];
00118     cout << endl;
00119   }
00120   cout << "-----------------------------------------\n";
00121 #endif
00122 
00123   vector<int> w(10*n);
00124 
00125   vector<int> rowperm(n);
00126   vector<int> colperm(n);
00127 
00128   //horizontal block
00129   int nhrows, nhcols, hrzcmp;
00130   //square block
00131   int nsrows, sqcmpn;
00132   //vertial block
00133   int nvrows, nvcols, vrtcmp;
00134 
00135   vector<int> rcmstr(n+1);
00136   vector<int> ccmstr(n+1);
00137 
00138   int msglvl = 0;
00139   int output = 6;
00140 
00141   GENBTF_F77( &n, &n, &iat[0], &jat[0], &ia[0], &ja[0], &w[0],
00142           &rowperm[0], &colperm[0], &nhrows, &nhcols,
00143           &hrzcmp, &nsrows, &sqcmpn, &nvrows, &nvcols, &vrtcmp,
00144           &rcmstr[0], &ccmstr[0], &msglvl, &output );
00145 
00146   //convert back to C indexing
00147   for( int i = 0; i < n; ++i )
00148   {
00149     --rowperm[i];
00150     --colperm[i];
00151   }
00152   for( int i = 0; (i<n+1) && (rcmstr[i]!=n+1); ++i )
00153   {
00154     --rcmstr[i];
00155     --ccmstr[i];
00156   }
00157 
00158 #ifdef BTF_VERBOSE
00159   cout << "-----------------------------------------\n";
00160   cout << "BTF Output\n";
00161   cout << "-----------------------------------------\n";
00162   cout << "RowPerm and ColPerm\n";
00163   for( int i = 0; i<n; ++i )
00164     cout << rowperm[i] << "\t" << colperm[i] << endl;
00165   if( hrzcmp )
00166   {
00167     cout << "Num Horizontal: Rows, Cols, Comps\n";
00168     cout << nhrows << "\t" << nhcols << "\t" << hrzcmp << endl;
00169   }
00170   cout << "Num Square: Rows, Comps\n";
00171   cout << nsrows << "\t" << sqcmpn << endl;
00172   if( vrtcmp )
00173   {
00174     cout << "Num Vertical: Rows, Cols, Comps\n";
00175     cout << nvrows << "\t" << nvcols << "\t" << vrtcmp << endl;
00176   }
00177   cout << "Row, Col of upper left pt in blocks\n";
00178   for( int i = 0; (i<n+1) && (rcmstr[i]!=n+1); ++i )
00179     cout << i << " " << rcmstr[i] << " " << ccmstr[i] << endl;
00180   cout << "-----------------------------------------\n";
00181 #endif
00182 
00183   if( hrzcmp || vrtcmp )
00184   { cout << "FAILED! hrz cmp's:" << hrzcmp << " vrtcmp: " << vrtcmp << endl;
00185     exit(0); }
00186 
00187   //convert rowperm to OLD->NEW
00188   //reverse ordering of permutation to get upper triangular
00189   vector<int> rowperm_t( n );
00190   vector<int> colperm_t( n );
00191   for( int i = 0; i < n; ++i )
00192   {
00193 //    rowperm_t[ rowperm[i] ] = n-i;
00194 //    colperm[i] = n-colperm[i];
00195     rowperm_t[i] = rowperm[(n-1)-i];
00196     colperm_t[i] = colperm[(n-1)-i];
00197   }
00198 
00199   //Generate New Domain and Range Maps
00200   //for now, assume they start out as identical
00201   const Epetra_BlockMap & OldMap = orig.RowMap();
00202   vector<int> myElements( n );
00203   OldMap.MyGlobalElements( &myElements[0] );
00204 
00205   vector<int> newDomainElements( n );
00206   vector<int> newRangeElements( n );
00207   for( int i = 0; i < n; ++i )
00208   {
00209     newDomainElements[ i ] = myElements[ rowperm_t[i] ];
00210     newRangeElements[ i ] = myElements[ colperm_t[i] ];
00211 cout << i << "\t" << rowperm_t[i] << "\t" << colperm[i] << "\t" << myElements[i] << endl;
00212   }
00213 
00214   NewRowMap_ = new Epetra_Map( n, n, &newDomainElements[0], OldMap.IndexBase(), OldMap.Comm() );
00215   NewDomainMap_ = new Epetra_Map( n, n, &newRangeElements[0], OldMap.IndexBase(), OldMap.Comm() );
00216 
00217 #ifdef BTF_VERBOSE
00218   cout << "New Row Map\n";
00219   cout << *RowMap << endl;
00220   cout << "New Domain Map\n";
00221   cout << *DomainMap << endl;
00222 #endif
00223 
00224   //Generate New Graph
00225   NewGraph_ = new Epetra_CrsGraph( Copy, *NewRowMap_, 0 );
00226   Epetra_Import Importer( *NewRowMap_, OldMap );
00227   NewGraph_->Import( orig, Importer, Insert );
00228   NewGraph_->FillComplete( *NewDomainMap_, *NewRowMap_ );
00229 
00230 #ifdef BTF_VERBOSE
00231   cout << "New CrsGraph\n";
00232   cout << *NewGraph_ << endl;
00233 #endif
00234 
00235   newObj_ = NewGraph_;
00236 
00237   return *NewGraph_;
00238 }
00239 
00240 } //namespace EpetraExt
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Defines