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 (NumMyRows_ == 0)
00127     IFPACK_CHK_ERR(-1); // strange graph this one
00128   
00129   if ((RootNode_ < 0) || (RootNode_ >= NumMyRows_))
00130     RootNode_ = 0;
00131     
00132   Reorder_.resize(NumMyRows_);
00133 
00134   for (int i = 0 ; i < NumMyRows_ ; ++i)
00135     Reorder_[i] = -1;
00136 
00137   vector<int> tmp;
00138   tmp.push_back(RootNode_);
00139 
00140   int count = NumMyRows_ - 1;
00141   int Length = Graph.MaxMyNumEntries();
00142   vector<int> Indices(Length);
00143   
00144   Reorder_[RootNode_] = count;
00145   count--;
00146 
00147   // stop when no nodes have been added in the previous level
00148 
00149   while (tmp.size()) {
00150 
00151     vector<int> tmp2;
00152 
00153     // for each node in the previous level, look for non-marked
00154     // neighbors. 
00155     for (int i = 0 ; i < (int)tmp.size() ; ++i) {
00156       int NumEntries;
00157       IFPACK_CHK_ERR(Graph.ExtractMyRowCopy(tmp[i], Length,
00158                          NumEntries, &Indices[0]));
00159 
00160       if (Length > 1)
00161     sort(Indices.begin(), Indices.begin() + Length);
00162 
00163       for (int j = 0 ; j < NumEntries ; ++j) {
00164     int col = Indices[j];
00165     if (col >= NumMyRows_) 
00166       continue;
00167 
00168     if (Reorder_[col] == -1) {
00169       Reorder_[col] = count;
00170       count--;
00171       if (col != tmp[i]) {
00172         tmp2.push_back(col);
00173       }
00174     }
00175       }
00176     }
00177 
00178     // if no nodes have been found but we still have
00179     // rows to walk through, to localize the next -1 
00180     // and restart.
00181     // FIXME: I can replace with STL
00182     if ((tmp2.size() == 0) && (count != -1)) {
00183       for (int i = 0 ; i < NumMyRows_ ; ++i)
00184     if (Reorder_[i] == -1) {
00185       tmp2.push_back(i);
00186       Reorder_[i] = count--;
00187       break;
00188     }
00189     }
00190 
00191     // prepare for the next level
00192     tmp = tmp2;
00193   }
00194 
00195   // check nothing went wrong
00196   for (int i = 0 ; i < NumMyRows_ ; ++i) {
00197     if (Reorder_[i] == -1)
00198       IFPACK_CHK_ERR(-1);
00199   }
00200   
00201   // build inverse reorder (will be used by ExtractMyRowCopy() 
00202   InvReorder_.resize(NumMyRows_);
00203 
00204   for (int i = 0 ; i < NumMyRows_ ; ++i)
00205     InvReorder_[i] = -1;
00206 
00207   for (int i = 0 ; i < NumMyRows_ ; ++i)
00208     InvReorder_[Reorder_[i]] = i;
00209 
00210   for (int i = 0 ; i < NumMyRows_ ; ++i) {
00211     if (InvReorder_[i] == -1)
00212       IFPACK_CHK_ERR(-1);
00213   }
00214 
00215   IsComputed_ = true;
00216   return(0);
00217 }
00218 
00219 //==============================================================================
00220 int Ifpack_RCMReordering::Reorder(const int i) const
00221 {
00222 #ifdef IFPACK_ABC
00223   if (!IsComputed())
00224     IFPACK_CHK_ERR(-1);
00225   if ((i < 0) || (i >= NumMyRows_))
00226     IFPACK_CHK_ERR(-1);
00227 #endif
00228 
00229   return(Reorder_[i]);
00230 }
00231 
00232 //==============================================================================
00233 int Ifpack_RCMReordering::InvReorder(const int i) const
00234 {
00235 #ifdef IFPACK_ABC
00236   if (!IsComputed())
00237     IFPACK_CHK_ERR(-1);
00238   if ((i < 0) || (i >= NumMyRows_))
00239     IFPACK_CHK_ERR(-1);
00240 #endif
00241 
00242   return(InvReorder_[i]);
00243 }
00244 //==============================================================================
00245 int Ifpack_RCMReordering::P(const Epetra_MultiVector& Xorig,
00246                 Epetra_MultiVector& X) const
00247 {  
00248   int NumVectors = X.NumVectors();
00249 
00250   for (int j = 0 ; j < NumVectors ; ++j) {
00251     for (int i = 0 ; i < NumMyRows_ ; ++i) {
00252       int np = Reorder_[i];
00253       X[j][np] = Xorig[j][i];
00254     }
00255   }
00256 
00257   return(0);
00258 }
00259 
00260 //==============================================================================
00261 int Ifpack_RCMReordering::Pinv(const Epetra_MultiVector& Xorig,
00262                    Epetra_MultiVector& X) const
00263 {
00264   int NumVectors = X.NumVectors();
00265 
00266   for (int j = 0 ; j < NumVectors ; ++j) {
00267     for (int i = 0 ; i < NumMyRows_ ; ++i) {
00268       int np = Reorder_[i];
00269       X[j][i] = Xorig[j][np];
00270     }
00271   }
00272 
00273   return(0);
00274 }
00275 
00276 //==============================================================================
00277 ostream& Ifpack_RCMReordering::Print(std::ostream& os) const
00278 {
00279   os << "*** Ifpack_RCMReordering" << endl << endl;
00280   if (!IsComputed())
00281     os << "*** Reordering not yet computed." << endl;
00282   
00283   os << "*** Number of local rows = " << NumMyRows_ << endl;
00284   os << "*** Root node = " << RootNode_ << endl;
00285   os << endl;
00286   os << "Local Row\tReorder[i]\tInvReorder[i]" << endl;
00287   for (int i = 0 ; i < NumMyRows_ ; ++i) {
00288     os << '\t' << i << "\t\t" << Reorder_[i] << "\t\t" << InvReorder_[i] << endl;
00289   }
00290    
00291   return(os);
00292 }

Generated on Tue Oct 20 12:48:54 2009 for IFPACK by doxygen 1.4.7