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