Ifpack_RCMReordering.cpp

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_RCMReordering.h"
00038 
00039 //==============================================================================
00040 Ifpack_RCMReordering::
00041 Ifpack_RCMReordering() :
00042   RootNode_(0),
00043   NumMyRows_(0),
00044   IsComputed_(false)
00045 {
00046 }
00047 
00048 //==============================================================================
00049 Ifpack_RCMReordering::
00050 Ifpack_RCMReordering(const Ifpack_RCMReordering& RHS) :
00051   RootNode_(RHS.RootNode()),
00052   NumMyRows_(RHS.NumMyRows()),
00053   IsComputed_(RHS.IsComputed())
00054 {
00055   Reorder_.resize(NumMyRows());
00056   InvReorder_.resize(NumMyRows());
00057   for (int i = 0 ; i < NumMyRows() ; ++i) {
00058     Reorder_[i] = RHS.Reorder(i);
00059     InvReorder_[i] = RHS.InvReorder(i);
00060   }
00061 }
00062 
00063 //==============================================================================
00064 Ifpack_RCMReordering& Ifpack_RCMReordering::
00065 operator=(const Ifpack_RCMReordering& RHS)
00066 {
00067   if (this == &RHS) {
00068     return (*this);
00069   }
00070 
00071   NumMyRows_ = RHS.NumMyRows(); // set number of local rows
00072   RootNode_ = RHS.RootNode(); // set root node
00073   IsComputed_ = RHS.IsComputed();
00074   // resize vectors, and copy values from RHS
00075   Reorder_.resize(NumMyRows()); 
00076   InvReorder_.resize(NumMyRows());
00077   if (IsComputed()) {
00078     for (int i = 0 ; i < NumMyRows_ ; ++i) {
00079       Reorder_[i] = RHS.Reorder(i);
00080       InvReorder_[i] = RHS.InvReorder(i);
00081     }
00082   }
00083   return (*this);
00084 }
00085 
00086 //==============================================================================
00087 int Ifpack_RCMReordering::
00088 SetParameter(const string Name, const int Value)
00089 {
00090   if (Name == "reorder: root node")
00091     RootNode_ = Value;
00092   return(0);
00093 }
00094 
00095 //==============================================================================
00096 int Ifpack_RCMReordering::
00097 SetParameter(const string Name, const double Value)
00098 {
00099   return(0);
00100 }
00101 
00102 //==============================================================================
00103 int Ifpack_RCMReordering::
00104 SetParameters(Teuchos::ParameterList& List)
00105 {
00106   RootNode_ = List.get("reorder: root node", RootNode_);
00107   return(0);
00108 }
00109 
00110 //==============================================================================
00111 int Ifpack_RCMReordering::Compute(const Epetra_RowMatrix& Matrix)
00112 {
00113   Ifpack_Graph_Epetra_RowMatrix Graph(Teuchos::rcp(&Matrix,false));
00114 
00115   IFPACK_CHK_ERR(Compute(Graph));
00116 
00117   return(0);
00118 }
00119 
00120 //==============================================================================
00121 int Ifpack_RCMReordering::Compute(const Ifpack_Graph& Graph)
00122 {
00123   IsComputed_ = false;
00124   NumMyRows_ = Graph.NumMyRows();
00125   
00126   if ((RootNode_ < 0) || (RootNode_ >= NumMyRows_))
00127     RootNode_ = 0;
00128     
00129   Reorder_.resize(NumMyRows_);
00130 
00131   // the case where one processor holds no chunk of the graph happens...
00132   if (!NumMyRows_) {
00133     InvReorder_.resize(NumMyRows_);
00134     IsComputed_ = true;
00135     return(0);
00136   }
00137 
00138   for (int i = 0 ; i < NumMyRows_ ; ++i)
00139     Reorder_[i] = -1;
00140 
00141   vector<int> tmp;
00142   tmp.push_back(RootNode_);
00143 
00144   int count = NumMyRows_ - 1;
00145   int Length = Graph.MaxMyNumEntries();
00146   vector<int> Indices(Length);
00147   
00148   Reorder_[RootNode_] = count;
00149   count--;
00150 
00151   // stop when no nodes have been added in the previous level
00152 
00153   while (tmp.size()) {
00154 
00155     vector<int> tmp2;
00156 
00157     // for each node in the previous level, look for non-marked
00158     // neighbors. 
00159     for (int i = 0 ; i < (int)tmp.size() ; ++i) {
00160       int NumEntries;
00161       IFPACK_CHK_ERR(Graph.ExtractMyRowCopy(tmp[i], Length,
00162                          NumEntries, &Indices[0]));
00163 
00164       if (Length > 1)
00165     sort(Indices.begin(), Indices.begin() + Length);
00166 
00167       for (int j = 0 ; j < NumEntries ; ++j) {
00168     int col = Indices[j];
00169     if (col >= NumMyRows_) 
00170       continue;
00171 
00172     if (Reorder_[col] == -1) {
00173       Reorder_[col] = count;
00174       count--;
00175       if (col != tmp[i]) {
00176         tmp2.push_back(col);
00177       }
00178     }
00179       }
00180     }
00181 
00182     // if no nodes have been found but we still have
00183     // rows to walk through, to localize the next -1 
00184     // and restart.
00185     // FIXME: I can replace with STL
00186     if ((tmp2.size() == 0) && (count != -1)) {
00187       for (int i = 0 ; i < NumMyRows_ ; ++i)
00188     if (Reorder_[i] == -1) {
00189       tmp2.push_back(i);
00190       Reorder_[i] = count--;
00191       break;
00192     }
00193     }
00194 
00195     // prepare for the next level
00196     tmp = tmp2;
00197   }
00198 
00199   // check nothing went wrong
00200   for (int i = 0 ; i < NumMyRows_ ; ++i) {
00201     if (Reorder_[i] == -1)
00202       IFPACK_CHK_ERR(-1);
00203   }
00204   
00205   // build inverse reorder (will be used by ExtractMyRowCopy() 
00206   InvReorder_.resize(NumMyRows_);
00207 
00208   for (int i = 0 ; i < NumMyRows_ ; ++i)
00209     InvReorder_[i] = -1;
00210 
00211   for (int i = 0 ; i < NumMyRows_ ; ++i)
00212     InvReorder_[Reorder_[i]] = i;
00213 
00214   for (int i = 0 ; i < NumMyRows_ ; ++i) {
00215     if (InvReorder_[i] == -1)
00216       IFPACK_CHK_ERR(-1);
00217   }
00218 
00219   IsComputed_ = true;
00220   return(0);
00221 }
00222 
00223 //==============================================================================
00224 int Ifpack_RCMReordering::Reorder(const int i) const
00225 {
00226 #ifdef IFPACK_ABC
00227   if (!IsComputed())
00228     IFPACK_CHK_ERR(-1);
00229   if ((i < 0) || (i >= NumMyRows_))
00230     IFPACK_CHK_ERR(-1);
00231 #endif
00232 
00233   return(Reorder_[i]);
00234 }
00235 
00236 //==============================================================================
00237 int Ifpack_RCMReordering::InvReorder(const int i) const
00238 {
00239 #ifdef IFPACK_ABC
00240   if (!IsComputed())
00241     IFPACK_CHK_ERR(-1);
00242   if ((i < 0) || (i >= NumMyRows_))
00243     IFPACK_CHK_ERR(-1);
00244 #endif
00245 
00246   return(InvReorder_[i]);
00247 }
00248 //==============================================================================
00249 int Ifpack_RCMReordering::P(const Epetra_MultiVector& Xorig,
00250                 Epetra_MultiVector& X) const
00251 {  
00252   int NumVectors = X.NumVectors();
00253 
00254   for (int j = 0 ; j < NumVectors ; ++j) {
00255     for (int i = 0 ; i < NumMyRows_ ; ++i) {
00256       int np = Reorder_[i];
00257       X[j][np] = Xorig[j][i];
00258     }
00259   }
00260 
00261   return(0);
00262 }
00263 
00264 //==============================================================================
00265 int Ifpack_RCMReordering::Pinv(const Epetra_MultiVector& Xorig,
00266                    Epetra_MultiVector& X) const
00267 {
00268   int NumVectors = X.NumVectors();
00269 
00270   for (int j = 0 ; j < NumVectors ; ++j) {
00271     for (int i = 0 ; i < NumMyRows_ ; ++i) {
00272       int np = Reorder_[i];
00273       X[j][i] = Xorig[j][np];
00274     }
00275   }
00276 
00277   return(0);
00278 }
00279 
00280 //==============================================================================
00281 ostream& Ifpack_RCMReordering::Print(std::ostream& os) const
00282 {
00283   os << "*** Ifpack_RCMReordering" << endl << endl;
00284   if (!IsComputed())
00285     os << "*** Reordering not yet computed." << endl;
00286   
00287   os << "*** Number of local rows = " << NumMyRows_ << endl;
00288   os << "*** Root node = " << RootNode_ << endl;
00289   os << endl;
00290   os << "Local Row\tReorder[i]\tInvReorder[i]" << endl;
00291   for (int i = 0 ; i < NumMyRows_ ; ++i) {
00292     os << '\t' << i << "\t\t" << Reorder_[i] << "\t\t" << InvReorder_[i] << endl;
00293   }
00294    
00295   return(os);
00296 }
 All Classes Files Functions Variables Enumerations Friends
Generated on Wed Apr 13 10:05:23 2011 for IFPACK by  doxygen 1.6.3