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