EpetraExt Development
EpetraExt_BTF_CrsMatrix.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_CrsMatrix.h>
00042 
00043 #include <EpetraExt_Transpose_CrsGraph.h>
00044 
00045 #include <Epetra_Import.h>
00046 #include <Epetra_CrsMatrix.h>
00047 #include <Epetra_CrsGraph.h>
00048 #include <Epetra_Map.h>
00049 
00050 #include <vector>
00051 
00052 using std::vector;
00053 using std::cout;
00054 using std::endl;
00055 
00056 #define MATTRANS_F77 F77_FUNC(mattrans,MATTRANS)
00057 #define GENBTF_F77   F77_FUNC(genbtf,GENBTF)
00058                                                                                                   
00059 extern "C" {
00060 extern void MATTRANS_F77( int*, int*, int*, int*, int*, int* );
00061 extern void GENBTF_F77( int*, int*, int*, int*, int*, int*, int*, int*, int*,
00062                      int*, int*, int*, int*, int*, int*, int*, int*, int*,
00063                      int*, int*, int* );
00064 }
00065 
00066 namespace EpetraExt {
00067 
00068 CrsMatrix_BTF::
00069 ~CrsMatrix_BTF()
00070 {
00071   if( NewRowMap_ ) delete NewRowMap_;
00072   if( NewColMap_ ) delete NewColMap_;
00073 
00074   if( Importer_ ) delete Importer_;
00075 
00076   if( NewMatrix_ ) delete NewMatrix_;
00077   if( NewGraph_ ) delete NewGraph_;
00078 }
00079 
00080 CrsMatrix_BTF::NewTypeRef
00081 CrsMatrix_BTF::
00082 operator()( OriginalTypeRef orig )
00083 {
00084   origObj_ = &orig;
00085 
00086   if( orig.RowMap().DistributedGlobal() )
00087   { cout << "FAIL for Global!\n"; abort(); }
00088   if( orig.IndicesAreGlobal() )
00089   { cout << "FAIL for Global Indices!\n"; abort(); }
00090 
00091   int n = orig.NumMyRows();
00092   int nnz = orig.NumMyNonzeros();
00093 
00094   if( verbose_ )
00095   {
00096     cout << "Orig Matrix:\n";
00097     cout << orig << endl;
00098   }
00099 
00100   //create std CRS format
00101   //also create graph without zero elements
00102   vector<int> ia(n+1,0);
00103   int maxEntries = orig.MaxNumEntries();
00104   vector<int> ja_tmp(maxEntries);
00105   vector<double> jva_tmp(maxEntries);
00106   vector<int> ja(nnz);
00107   int cnt;
00108 
00109   const Epetra_BlockMap & OldRowMap = orig.RowMap();
00110   const Epetra_BlockMap & OldColMap = orig.ColMap();
00111   Epetra_CrsGraph strippedGraph( Copy, OldRowMap, OldColMap, 0 );
00112 
00113   for( int i = 0; i < n; ++i )
00114   {
00115     orig.ExtractMyRowCopy( i, maxEntries, cnt, &jva_tmp[0], &ja_tmp[0] );
00116     ia[i+1] = ia[i];
00117     for( int j = 0; j < cnt; ++j )
00118       if( fabs(jva_tmp[j]) > threshold_ )
00119         ja[ ia[i+1]++ ] = ja_tmp[j];
00120 
00121     int new_cnt = ia[i+1] - ia[i];
00122     strippedGraph.InsertMyIndices( i, new_cnt, &ja[ ia[i] ] );
00123   }
00124   nnz = ia[n];
00125   strippedGraph.FillComplete();
00126 
00127   if( verbose_ )
00128   {
00129     cout << "Stripped Graph\n";
00130     cout << strippedGraph;
00131   }
00132 
00133   vector<int> iat(n+1,0);
00134   vector<int> jat(nnz);
00135   for( int i = 0; i < n; ++i )
00136     for( int j = ia[i]; j < ia[i+1]; ++j )
00137       ++iat[ ja[j]+1 ];
00138   for( int i = 0; i < n; ++i )
00139     iat[i+1] += iat[i];
00140   for( int i = 0; i < n; ++i )
00141     for( int j = ia[i]; j < ia[i+1]; ++j )
00142       jat[ iat[ ja[j] ]++ ] = i;
00143   for( int i = 0; i < n; ++i )
00144     iat[n-i] = iat[n-i-1];
00145   iat[0] = 0;
00146 
00147   //convert to Fortran indexing
00148   for( int i = 0; i < n+1; ++i ) ++ia[i];
00149   for( int i = 0; i < nnz; ++i ) ++ja[i];
00150   for( int i = 0; i < n+1; ++i ) ++iat[i];
00151   for( int i = 0; i < nnz; ++i ) ++jat[i];
00152 
00153   if( verbose_ )
00154   {
00155     cout << "-----------------------------------------\n";
00156     cout << "CRS Format Graph\n";
00157     cout << "-----------------------------------------\n";
00158     for( int i = 0; i < n; ++i )
00159     {
00160       cout << i+1 << ": " << ia[i+1] << ": ";
00161       for( int j = ia[i]-1; j<ia[i+1]-1; ++j )
00162         cout << " " << ja[j];
00163       cout << endl;
00164     }
00165     cout << "-----------------------------------------\n";
00166   }
00167 
00168 /*
00169   vector<int> iat(n+1);
00170   vector<int> jat(nnz);
00171   int * jaf = &ja[0];
00172   int * iaf = &ia[0];
00173   int * jatf = &jat[0];
00174   int * iatf = &iat[0];
00175   MATTRANS_F77( &n, &n, &ja[0], &ia[0], &jat[0], &iat[0] );
00176 */
00177     
00178   if( verbose_ )
00179   {
00180     cout << "-----------------------------------------\n";
00181     cout << "CCS Format Graph\n";
00182     cout << "-----------------------------------------\n";
00183     for( int i = 0; i < n; ++i )
00184     {
00185       cout << i+1 << ": " << iat[i+1] << ": ";
00186       for( int j = iat[i]-1; j<iat[i+1]-1; ++j )
00187         cout << " " << jat[j];
00188       cout << endl;
00189     }
00190     cout << "-----------------------------------------\n";
00191   }
00192 
00193   vector<int> w(10*n);
00194 
00195   vector<int> rowperm(n);
00196   vector<int> colperm(n);
00197 
00198   //horizontal block
00199   int nhrows, nhcols, hrzcmp;
00200   //square block
00201   int nsrows, sqcmpn;
00202   //vertial block
00203   int nvrows, nvcols, vrtcmp;
00204 
00205   vector<int> rcmstr(n+1);
00206   vector<int> ccmstr(n+1);
00207 
00208   int msglvl = 0;
00209   int output = 6;
00210 
00211   GENBTF_F77( &n, &n, &iat[0], &jat[0], &ia[0], &ja[0], &w[0],
00212           &rowperm[0], &colperm[0], &nhrows, &nhcols,
00213           &hrzcmp, &nsrows, &sqcmpn, &nvrows, &nvcols, &vrtcmp,
00214           &rcmstr[0], &ccmstr[0], &msglvl, &output );
00215 
00216   //convert back to C indexing
00217   for( int i = 0; i < n; ++i )
00218   {
00219     --rowperm[i];
00220     --colperm[i];
00221   }
00222   for( int i = 0; (i<n+1) && (rcmstr[i]!=n+1); ++i )
00223   {
00224     --rcmstr[i];
00225     --ccmstr[i];
00226   }
00227 
00228   if( verbose_ )
00229   {
00230     cout << "-----------------------------------------\n";
00231     cout << "BTF Output\n";
00232     cout << "-----------------------------------------\n";
00233     cout << "RowPerm and ColPerm\n";
00234     for( int i = 0; i<n; ++i )
00235       cout << rowperm[i] << "\t" << colperm[i] << endl;
00236     if( hrzcmp )
00237     {
00238       cout << "Num Horizontal: Rows, Cols, Comps\n";
00239       cout << nhrows << "\t" << nhcols << "\t" << hrzcmp << endl;
00240     }
00241     cout << "Num Square: Rows, Comps\n";
00242     cout << nsrows << "\t" << sqcmpn << endl;
00243     if( vrtcmp )
00244     {
00245       cout << "Num Vertical: Rows, Cols, Comps\n";
00246       cout << nvrows << "\t" << nvcols << "\t" << vrtcmp << endl;
00247     }
00248     cout << "Row, Col of upper left pt in blocks\n";
00249     for( int i = 0; (i<n+1) && (rcmstr[i]!=n+1); ++i )
00250       cout << i << " " << rcmstr[i] << " " << ccmstr[i] << endl;
00251     cout << "-----------------------------------------\n";
00252   }
00253 
00254   if( hrzcmp || vrtcmp )
00255   {
00256     cout << "FAILED! hrz cmp's:" << hrzcmp << " vrtcmp: " << vrtcmp << endl;
00257     exit(0);
00258   }
00259 
00260   //convert rowperm to OLD->NEW
00261   //reverse ordering of permutation to get upper triangular
00262   vector<int> rowperm_t( n );
00263   vector<int> colperm_t( n );
00264   for( int i = 0; i < n; ++i )
00265   {
00266 //    rowperm_t[ rowperm[i] ] = i;
00267     rowperm_t[i] = rowperm[i];
00268     colperm_t[i] = colperm[i];
00269   }
00270 
00271   //Generate New Domain and Range Maps
00272   //for now, assume they start out as identical
00273   vector<int> myElements( n );
00274   OldRowMap.MyGlobalElements( &myElements[0] );
00275 
00276   vector<int> newDomainElements( n );
00277   vector<int> newRangeElements( n );
00278   for( int i = 0; i < n; ++i )
00279   {
00280     newDomainElements[ i ] = myElements[ rowperm_t[i] ];
00281     newRangeElements[ i ] = myElements[ colperm_t[i] ];
00282   }
00283 
00284   NewRowMap_ = new Epetra_Map( n, n, &newDomainElements[0], OldRowMap.IndexBase(), OldRowMap.Comm() );
00285   NewColMap_ = new Epetra_Map( n, n, &newRangeElements[0], OldColMap.IndexBase(), OldColMap.Comm() );
00286 
00287   if( verbose_ )
00288   {
00289     cout << "New Row Map\n";
00290     cout << *NewRowMap_ << endl;
00291     cout << "New ColMap\n";
00292     cout << *NewColMap_ << endl;
00293   }
00294 
00295   //Generate New Graph
00296   NewGraph_ = new Epetra_CrsGraph( Copy, *NewRowMap_, *NewColMap_, 0 );
00297   Importer_ = new Epetra_Import( *NewRowMap_, OldRowMap );
00298   NewGraph_->Import( strippedGraph, *Importer_, Insert );
00299   NewGraph_->FillComplete();
00300   if( verbose_ )
00301   {
00302     cout << "NewGraph\n";
00303     cout << *NewGraph_;
00304   }
00305 
00306   NewMatrix_ = new Epetra_CrsMatrix( Copy, *NewGraph_ );
00307   NewMatrix_->Import( orig, *Importer_, Insert );
00308   NewMatrix_->FillComplete();
00309 
00310   if( verbose_ )
00311   {
00312     cout << "New CrsMatrix\n";
00313     cout << *NewMatrix_ << endl;
00314   }
00315 
00316   newObj_ = NewMatrix_;
00317 
00318   return *NewMatrix_;
00319 }
00320 
00321 bool
00322 CrsMatrix_BTF::
00323 fwd()
00324 {
00325   int ret = NewMatrix_->Import( *origObj_, *Importer_, Insert );
00326   if (ret<0) return false;
00327   return true;
00328 }
00329 
00330 bool
00331 CrsMatrix_BTF::
00332 rvs()
00333 {
00334   int ret = origObj_->Export( *NewMatrix_, *Importer_, Insert );
00335   if (ret<0) return false;
00336   return true;
00337 }
00338 
00339 } //namespace EpetraExt
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Defines