EpetraExt Package Browser (Single Doxygen Collection) Development
EpetraExt_AMD_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 
00042 #include "EpetraExt_ConfigDefs.h"
00043 #ifdef HAVE_EXPERIMENTAL
00044 #ifdef HAVE_GRAPH_REORDERINGS
00045 
00046 #include <EpetraExt_AMD_CrsGraph.h>
00047 
00048 #include <Epetra_Import.h>
00049 #include <Epetra_CrsGraph.h>
00050 #include <Epetra_Map.h>
00051 
00052 #include <vector>
00053 
00054 extern "C" {
00055 #include <amd.h>
00056 }
00057 
00058 namespace EpetraExt {
00059 
00060 CrsGraph_AMD::
00061 ~CrsGraph_AMD()
00062 {
00063   if( NewMap_ ) delete NewMap_;
00064 
00065   if( NewGraph_ ) delete NewGraph_;
00066 }
00067 
00068 CrsGraph_AMD::NewTypeRef
00069 CrsGraph_AMD::
00070 operator()( OriginalTypeRef orig )
00071 {
00072   origObj_ = &orig;
00073 
00074   int n = orig.NumMyRows();
00075   int nnz = orig.NumMyNonzeros();
00076 
00077   //create std CRS format
00078   std::vector<int> ia(n+1,0);
00079   std::vector<int> ja(nnz);
00080   int cnt;
00081   for( int i = 0; i < n; ++i )
00082   {
00083     int * tmpP = &ja[ia[i]];
00084     orig.ExtractMyRowCopy( i, nnz-ia[i], cnt, tmpP );
00085     ia[i+1] = ia[i] + cnt;
00086   }
00087 
00088   //trim down to local only
00089   std::vector<int> iat(n+1);
00090   std::vector<int> jat(nnz);
00091   int loc = 0;
00092   for( int i = 0; i < n; ++i )
00093   {
00094     iat[i] = loc;
00095     for( int j = ia[i]; j < ia[i+1]; ++j )
00096     {
00097       if( ja[j] < n )
00098         jat[loc++] = ja[j];
00099       else
00100   break;
00101     }
00102   }
00103   iat[n] = loc;
00104 
00105 
00106   if( verbose_ )
00107   {
00108     std::cout << "Orig Graph\n";
00109     std::cout << orig << std::endl;
00110     std::cout << "-----------------------------------------\n";
00111     std::cout << "CRS Format Graph\n";
00112     std::cout << "-----------------------------------------\n";
00113     for( int i = 0; i < n; ++i )
00114     {
00115       std::cout << i << ": " << iat[i+1] << ": ";
00116       for( int j = iat[i]; j<iat[i+1]; ++j )
00117         std::cout << " " << jat[j];
00118       std::cout << std::endl;
00119     }
00120     std::cout << "-----------------------------------------\n";
00121   }
00122 
00123   std::vector<int> perm(n);
00124   std::vector<double> info(AMD_INFO);
00125 
00126   amd_order( n, &iat[0], &jat[0], &perm[0], NULL, &info[0] ); 
00127 
00128   if( info[AMD_STATUS] == AMD_INVALID )
00129     std::cout << "AMD ORDERING: Invalid!!!!\n";
00130 
00131   if( verbose_ )
00132   {
00133     std::cout << "-----------------------------------------\n";
00134     std::cout << "AMD Output\n";
00135     std::cout << "-----------------------------------------\n";
00136     std::cout << "STATUS: " << info[AMD_STATUS] << std::endl;
00137     std::cout << "SYMM: " << info[AMD_SYMMETRY] << std::endl;
00138     std::cout << "N: " << info[AMD_N] << std::endl;
00139     std::cout << "NZ: " << info[AMD_NZ] << std::endl;
00140     std::cout << "SYMM: " << info[AMD_SYMMETRY] << std::endl;
00141     std::cout << "NZDIAG: " << info[AMD_NZDIAG] << std::endl;
00142     std::cout << "NZ A+At: " << info[AMD_NZ_A_PLUS_AT] << std::endl;
00143     std::cout << "NDENSE: " << info[AMD_SYMMETRY] << std::endl;
00144     std::cout << "Perm\n";
00145     for( int i = 0; i<n; ++i )
00146       std::cout << perm[i] << std::endl;
00147     std::cout << "-----------------------------------------\n";
00148   }
00149 
00150   //Generate New Domain and Range Maps
00151   //for now, assume they start out as identical
00152   const Epetra_BlockMap & OldMap = orig.RowMap();
00153   int nG = orig.NumGlobalRows();
00154 
00155   std::vector<int> newElements( n );
00156   for( int i = 0; i < n; ++i )
00157     newElements[i] = OldMap.GID( perm[i] );
00158 
00159   NewMap_ = new Epetra_Map( nG, n, &newElements[0], OldMap.IndexBase(), OldMap.Comm() );
00160 
00161   if( verbose_ )
00162   {
00163     std::cout << "Old Map\n";
00164     std::cout << OldMap << std::endl;
00165     std::cout << "New Map\n";
00166     std::cout << *NewMap_ << std::endl;
00167   }
00168 
00169   //Generate New Graph
00170   NewGraph_ = new Epetra_CrsGraph( Copy, *NewMap_, 0 );
00171   Epetra_Import Importer( *NewMap_, OldMap );
00172   NewGraph_->Import( orig, Importer, Insert );
00173   NewGraph_->FillComplete();
00174 
00175   if( verbose_ )
00176   {
00177     std::cout << "New CrsGraph\n";
00178     std::cout << *NewGraph_ << std::endl;
00179   }
00180 
00181   newObj_ = NewGraph_;
00182 
00183   return *NewGraph_;
00184 }
00185 
00186 } //namespace EpetraExt
00187 
00188 #endif //HAVE_GRAPH_REORDERINGS
00189 #endif //HAVE_EXPERIMENTAL
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Defines