Ifpack Package Browser (Single Doxygen Collection) Development
Ifpack_AMDReordering.cpp
Go to the documentation of this file.
00001 /*@HEADER
00002 // ***********************************************************************
00003 //
00004 //       Ifpack: Object-Oriented Algebraic Preconditioner Package
00005 //                 Copyright (2002) 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 
00030 #include "Ifpack_ConfigDefs.h"
00031 #include "Teuchos_ParameterList.hpp"
00032 #include "Teuchos_RefCountPtr.hpp"
00033 #include "Epetra_MultiVector.h"
00034 #include "Ifpack_Graph.h"
00035 #include "Epetra_RowMatrix.h"
00036 #include "Ifpack_Graph_Epetra_RowMatrix.h"
00037 #include "Ifpack_AMDReordering.h"
00038 
00039 extern "C" {
00040 #include <amesos_amd.h>
00041 }
00042 
00043 //==============================================================================
00044 Ifpack_AMDReordering::
00045 Ifpack_AMDReordering() :
00046   NumMyRows_(0),
00047   IsComputed_(false)
00048 {
00049 }
00050 
00051 //==============================================================================
00052 Ifpack_AMDReordering::
00053 Ifpack_AMDReordering(const Ifpack_AMDReordering& RHS) :
00054   NumMyRows_(RHS.NumMyRows()),
00055   IsComputed_(RHS.IsComputed())
00056 {
00057   Reorder_.resize(NumMyRows());
00058   InvReorder_.resize(NumMyRows());
00059   for (int i = 0 ; i < NumMyRows() ; ++i) {
00060     Reorder_[i] = RHS.Reorder(i);
00061     InvReorder_[i] = RHS.InvReorder(i);
00062   }
00063 }
00064 
00065 //==============================================================================
00066 Ifpack_AMDReordering& Ifpack_AMDReordering::
00067 operator=(const Ifpack_AMDReordering& RHS)
00068 {
00069   if (this == &RHS) {
00070     return (*this);
00071   }
00072 
00073   NumMyRows_ = RHS.NumMyRows(); // set number of local rows
00074   IsComputed_ = RHS.IsComputed();
00075   // resize vectors, and copy values from RHS
00076   Reorder_.resize(NumMyRows()); 
00077   InvReorder_.resize(NumMyRows());
00078   if (IsComputed()) {
00079     for (int i = 0 ; i < NumMyRows_ ; ++i) {
00080       Reorder_[i] = RHS.Reorder(i);
00081       InvReorder_[i] = RHS.InvReorder(i);
00082     }
00083   }
00084   return (*this);
00085 }
00086 
00087 //==============================================================================
00088 int Ifpack_AMDReordering::
00089 SetParameter(const string Name, const int Value)
00090 {
00091   return(0);
00092 }
00093 
00094 //==============================================================================
00095 int Ifpack_AMDReordering::
00096 SetParameter(const string Name, const double Value)
00097 {
00098   return(0);
00099 }
00100 
00101 //==============================================================================
00102 int Ifpack_AMDReordering::
00103 SetParameters(Teuchos::ParameterList& List)
00104 {
00105   return(0);
00106 }
00107 
00108 //==============================================================================
00109 int Ifpack_AMDReordering::Compute(const Epetra_RowMatrix& Matrix)
00110 {
00111   Ifpack_Graph_Epetra_RowMatrix Graph(Teuchos::rcp(&Matrix,false));
00112 
00113   IFPACK_CHK_ERR(Compute(Graph));
00114 
00115   return(0);
00116 }
00117 
00118 //==============================================================================
00119 int Ifpack_AMDReordering::Compute(const Ifpack_Graph& Graph)
00120 {
00121   IsComputed_ = false;
00122   NumMyRows_ = Graph.NumMyRows();
00123   int NumNz = Graph.NumMyNonzeros();
00124   
00125   if (NumMyRows_ == 0)
00126     IFPACK_CHK_ERR(-1); // strange graph this one
00127   
00128   // Extract CRS format
00129   vector<int> ia(NumMyRows_+1,0);
00130   vector<int> ja(NumNz);
00131   int cnt;
00132   for( int i = 0; i < NumMyRows_; ++i )
00133   {
00134     int * tmpP = &ja[ia[i]];
00135     Graph.ExtractMyRowCopy( i, NumNz-ia[i], cnt, tmpP );
00136     ia[i+1] = ia[i] + cnt;
00137   }
00138 
00139   // Trim down to local only
00140   vector<int> iat(NumMyRows_+1);
00141   vector<int> jat(NumNz);
00142   int loc = 0;
00143   for( int i = 0; i < NumMyRows_; ++i )
00144   {
00145     iat[i] = loc;
00146     for( int j = ia[i]; j < ia[i+1]; ++j )
00147     {
00148       if( ja[j] < NumMyRows_ )
00149         jat[loc++] = ja[j];
00150       else
00151   break;
00152     }
00153   }
00154   iat[NumMyRows_] = loc;
00155 
00156   // Compute AMD permutation
00157   Reorder_.resize(NumMyRows_);
00158   vector<double> info(AMD_INFO);
00159 
00160   amesos_amd_order( NumMyRows_, &iat[0], &jat[0], &Reorder_[0], NULL, &info[0] );
00161 
00162   if( info[AMD_STATUS] == AMD_INVALID )
00163     cout << "AMD ORDERING: Invalid!!!!\n";
00164 
00165   // Build inverse reorder (will be used by ExtractMyRowCopy() 
00166   InvReorder_.resize(NumMyRows_);
00167 
00168   for (int i = 0 ; i < NumMyRows_ ; ++i)
00169     InvReorder_[i] = -1;
00170 
00171   for (int i = 0 ; i < NumMyRows_ ; ++i)
00172     InvReorder_[Reorder_[i]] = i;
00173 
00174   for (int i = 0 ; i < NumMyRows_ ; ++i) {
00175     if (InvReorder_[i] == -1)
00176       IFPACK_CHK_ERR(-1);
00177   }
00178 
00179   IsComputed_ = true;
00180   return(0);
00181 }
00182 
00183 //==============================================================================
00184 int Ifpack_AMDReordering::Reorder(const int i) const
00185 {
00186 #ifdef IFPACK_ABC
00187   if (!IsComputed())
00188     IFPACK_CHK_ERR(-1);
00189   if ((i < 0) || (i >= NumMyRows_))
00190     IFPACK_CHK_ERR(-1);
00191 #endif
00192 
00193   return(Reorder_[i]);
00194 }
00195 
00196 //==============================================================================
00197 int Ifpack_AMDReordering::InvReorder(const int i) const
00198 {
00199 #ifdef IFPACK_ABC
00200   if (!IsComputed())
00201     IFPACK_CHK_ERR(-1);
00202   if ((i < 0) || (i >= NumMyRows_))
00203     IFPACK_CHK_ERR(-1);
00204 #endif
00205 
00206   return(InvReorder_[i]);
00207 }
00208 //==============================================================================
00209 int Ifpack_AMDReordering::P(const Epetra_MultiVector& Xorig,
00210           Epetra_MultiVector& X) const
00211 {  
00212   int NumVectors = X.NumVectors();
00213 
00214   for (int j = 0 ; j < NumVectors ; ++j) {
00215     for (int i = 0 ; i < NumMyRows_ ; ++i) {
00216       int np = Reorder_[i];
00217       X[j][np] = Xorig[j][i];
00218     }
00219   }
00220 
00221   return(0);
00222 }
00223 
00224 //==============================================================================
00225 int Ifpack_AMDReordering::Pinv(const Epetra_MultiVector& Xorig,
00226              Epetra_MultiVector& X) const
00227 {
00228   int NumVectors = X.NumVectors();
00229 
00230   for (int j = 0 ; j < NumVectors ; ++j) {
00231     for (int i = 0 ; i < NumMyRows_ ; ++i) {
00232       int np = Reorder_[i];
00233       X[j][i] = Xorig[j][np];
00234     }
00235   }
00236 
00237   return(0);
00238 }
00239 
00240 //==============================================================================
00241 ostream& Ifpack_AMDReordering::Print(std::ostream& os) const
00242 {
00243   os << "*** Ifpack_AMDReordering" << endl << endl;
00244   if (!IsComputed())
00245     os << "*** Reordering not yet computed." << endl;
00246   
00247   os << "*** Number of local rows = " << NumMyRows_ << endl;
00248   os << endl;
00249   os << "Local Row\tReorder[i]\tInvReorder[i]" << endl;
00250   for (int i = 0 ; i < NumMyRows_ ; ++i) {
00251     os << '\t' << i << "\t\t" << Reorder_[i] << "\t\t" << InvReorder_[i] << endl;
00252   }
00253    
00254   return(os);
00255 }
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Defines