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