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

Generated on Thu Sep 18 12:31:43 2008 for EpetraExt by doxygen 1.3.9.1