IFPACK Development
Ifpack_SingletonFilter.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 "Ifpack_SingletonFilter.h"
00032 #include "Epetra_ConfigDefs.h"
00033 #include "Epetra_RowMatrix.h"
00034 #include "Epetra_Comm.h"
00035 #include "Epetra_Map.h"
00036 #include "Epetra_MultiVector.h"
00037 #include "Epetra_Vector.h"
00038 #include "Epetra_Import.h"
00039 
00040 //==============================================================================
00041 Ifpack_SingletonFilter::Ifpack_SingletonFilter(const Teuchos::RefCountPtr<Epetra_RowMatrix>& Matrix) :
00042   A_(Matrix),
00043   NumSingletons_(0),
00044   NumRows_(0),
00045   NumRowsA_(0),
00046   MaxNumEntries_(0),
00047   MaxNumEntriesA_(0),
00048   NumNonzeros_(0)
00049 {
00050   // use this filter only on serial matrices
00051   if (A_->Comm().NumProc() != 1) {
00052     cerr << "Ifpack_SingletonFilter can be used with Comm().NumProc() == 1" << endl;
00053     cerr << "only. This class is a tool for Ifpack_AdditiveSchwarz," << endl;
00054     cerr << "and it is not meant to be used otherwise." << endl;
00055     exit(EXIT_FAILURE);
00056   }
00057   
00058   if ((A_->NumMyRows() != A_->NumGlobalRows()) ||
00059      (A_->NumMyRows() != A_->NumMyCols()))
00060     IFPACK_CHK_ERRV(-1);
00061   
00062   NumRowsA_ = (A_->NumMyRows());
00063   MaxNumEntriesA_ = A_->MaxNumEntries();
00064 
00065   Indices_.resize(MaxNumEntriesA_);
00066   Values_.resize(MaxNumEntriesA_);
00067   Reorder_.resize(A_->NumMyRows());
00068 
00069   for (int i = 0 ; i < NumRowsA_ ; ++i)
00070     Reorder_[i] = -1;
00071 
00072   // first check how may singletons I do have
00073   for (int i = 0 ; i < NumRowsA_ ; ++i) {
00074     int Nnz;
00075     IFPACK_CHK_ERRV(A_->ExtractMyRowCopy(i,MaxNumEntriesA_,Nnz,
00076                     &Values_[0], &Indices_[0]));
00077     if (Nnz != 1) {
00078       Reorder_[i] = NumRows_++;
00079     }
00080     else {
00081       NumSingletons_++;
00082     }
00083   }
00084 
00085   InvReorder_.resize(NumRows_);
00086   for (int i = 0 ; i < NumRowsA_ ; ++i) {
00087     if (Reorder_[i] < 0)
00088       continue;
00089     InvReorder_[Reorder_[i]] = i;
00090   }
00091  
00092   NumEntries_.resize(NumRows_);
00093   SingletonIndex_.resize(NumSingletons_);
00094 
00095   // now compute the nonzeros per row
00096   int count = 0;
00097   for (int i = 0 ; i < A_->NumMyRows() ; ++i) {
00098 
00099     int Nnz;
00100     IFPACK_CHK_ERRV(A_->ExtractMyRowCopy(i,MaxNumEntriesA_,Nnz,
00101                       &Values_[0], &Indices_[0]));
00102     int ii = Reorder_[i];
00103     if (ii >= 0) {
00104       assert (Nnz != 1);
00105 
00106       NumEntries_[ii] = Nnz;
00107       NumNonzeros_ += Nnz;
00108       if (Nnz > MaxNumEntries_)
00109     MaxNumEntries_ = Nnz;
00110     }
00111     else {
00112       SingletonIndex_[count] = i;
00113       count++;
00114     }
00115   }
00116 
00117   Map_ = Teuchos::rcp( new Epetra_Map(NumRows_,0,Comm()) );
00118 
00119   // and finish up with the diagonal entry
00120   Diagonal_ = Teuchos::rcp( new Epetra_Vector(*Map_) );
00121 
00122   Epetra_Vector Diagonal(A_->Map());
00123   A_->ExtractDiagonalCopy(Diagonal);
00124   for (int i = 0 ; i < NumRows_ ; ++i) {
00125     int ii = InvReorder_[i];
00126     (*Diagonal_)[i] = Diagonal[ii];
00127   }
00128   
00129 }
00130 
00131 //==============================================================================
00132 int Ifpack_SingletonFilter::
00133 ExtractMyRowCopy(int MyRow, int Length, int & NumEntries, 
00134          double *Values, int * Indices) const
00135 {
00136   int Nnz;
00137 
00138   if (Length < NumEntries_[MyRow])
00139     IFPACK_CHK_ERR(-1);
00140 
00141   int Row = InvReorder_[MyRow];
00142   IFPACK_CHK_ERR(A_->ExtractMyRowCopy(Row,MaxNumEntriesA_,Nnz,
00143                      &Values_[0],&Indices_[0]));
00144   NumEntries = 0;
00145   for (int i = 0 ; i < Nnz ; ++i) {
00146     int ii = Reorder_[Indices_[i]];
00147     if ( ii >= 0) {
00148       Indices[NumEntries] = ii;
00149       Values[NumEntries] = Values_[i];
00150       NumEntries++;
00151     }
00152   }
00153   return(0);   
00154 }
00155 
00156 //==============================================================================
00157 int Ifpack_SingletonFilter::
00158 ExtractDiagonalCopy(Epetra_Vector & Diagonal) const
00159 {
00160   Diagonal = *Diagonal_;
00161   return(0);
00162 }
00163 
00164 //==============================================================================
00165 int Ifpack_SingletonFilter::
00166 Multiply(bool TransA, const Epetra_MultiVector& X, 
00167      Epetra_MultiVector& Y) const
00168 {
00169 
00170   int NumVectors = X.NumVectors();
00171   if (NumVectors != Y.NumVectors())
00172     IFPACK_CHK_ERR(-1);
00173 
00174   Y.PutScalar(0.0);
00175 
00176   vector<int> Indices(MaxNumEntries_);
00177   vector<double> Values(MaxNumEntries_);
00178 
00179   // cycle over all rows of the original matrix
00180   for (int i = 0 ; i < A_->NumMyRows() ; ++i) {
00181 
00182     if (Reorder_[i] < 0)
00183       continue; // skip singleton rows
00184     
00185     int Nnz;
00186     A_->ExtractMyRowCopy(i,MaxNumEntriesA_,Nnz,
00187              &Values[0], &Indices[0]);
00188     if (!TransA) {
00189       // no transpose first
00190       for (int j = 0 ; j < NumVectors ; ++j) {
00191     for (int k = 0 ; k < Nnz ; ++k) {
00192       if (Reorder_[i] >= 0)
00193         Y[j][i] += Values[k] * X[j][Reorder_[Indices[k]]];
00194     }
00195       }
00196     }
00197     else {
00198       // transpose here
00199       for (int j = 0 ; j < NumVectors ; ++j) {
00200     for (int k = 0 ; k < Nnz ; ++k) {
00201       if (Reorder_[i] >= 0)
00202         Y[j][Reorder_[Indices[k]]] += Values[k] * X[j][i];
00203     }
00204       }
00205     }
00206   }
00207 
00208   return(0);
00209 }
00210 
00211 //==============================================================================
00212 int Ifpack_SingletonFilter::
00213 Solve(bool Upper, bool Trans, bool UnitDiagonal, 
00214       const Epetra_MultiVector& X, Epetra_MultiVector& Y) const
00215 {
00216   return(-1);
00217 }
00218 
00219 //==============================================================================
00220 int Ifpack_SingletonFilter::
00221 Apply(const Epetra_MultiVector& X, Epetra_MultiVector& Y) const
00222 {
00223   IFPACK_CHK_ERR(Multiply(false,X,Y));
00224   return(0);
00225 }
00226 
00227 //==============================================================================
00228 int Ifpack_SingletonFilter::
00229 ApplyInverse(const Epetra_MultiVector& X, Epetra_MultiVector& Y) const
00230 {
00231   return(-1); // NOT IMPLEMENTED AT THIS STAGE
00232 }
00233 
00234 //==============================================================================
00235 int Ifpack_SingletonFilter::
00236 SolveSingletons(const Epetra_MultiVector& RHS, 
00237         Epetra_MultiVector& LHS)
00238 {
00239   for (int i = 0 ; i < NumSingletons_ ; ++i) {
00240     int ii = SingletonIndex_[i];
00241     // get the diagonal value for the singleton
00242     int Nnz;
00243     A_->ExtractMyRowCopy(ii,MaxNumEntriesA_,Nnz,
00244             &Values_[0], &Indices_[0]);
00245     for (int j = 0 ; j < Nnz ; ++j) {
00246       if (Indices_[j] == ii) {
00247     for (int k = 0 ; k < LHS.NumVectors() ; ++k)
00248       LHS[k][ii] = RHS[k][ii] / Values_[j];
00249       }
00250     }
00251   }
00252 
00253   return(0);
00254 }
00255 
00256 //==============================================================================
00257 int Ifpack_SingletonFilter::
00258 CreateReducedRHS(const Epetra_MultiVector& LHS,
00259          const Epetra_MultiVector& RHS, 
00260          Epetra_MultiVector& ReducedRHS)
00261 {
00262   int NumVectors = LHS.NumVectors();
00263 
00264   for (int i = 0 ; i < NumRows_ ; ++i)
00265     for (int k = 0 ; k < NumVectors ; ++k)
00266       ReducedRHS[k][i] = RHS[k][InvReorder_[i]];
00267 
00268   for (int i = 0 ; i < NumRows_ ; ++i) {
00269     int ii = InvReorder_[i];
00270     int Nnz;
00271     IFPACK_CHK_ERR(A_->ExtractMyRowCopy(ii,MaxNumEntriesA_,Nnz,
00272                     &Values_[0], &Indices_[0]));
00273 
00274     for (int j = 0 ; j < Nnz ; ++j) {
00275       if (Reorder_[Indices_[j]] == -1) {
00276     for (int k = 0 ; k < NumVectors ; ++k)
00277       ReducedRHS[k][i] -= Values_[j] * LHS[k][Indices_[j]];
00278       }
00279     }
00280   }
00281   return(0);
00282 }
00283 
00284 //==============================================================================
00285 int Ifpack_SingletonFilter::
00286 UpdateLHS(const Epetra_MultiVector& ReducedLHS,
00287       Epetra_MultiVector& LHS)
00288 {
00289   for (int i = 0 ; i < NumRows_ ; ++i)
00290     for (int k = 0 ; k < LHS.NumVectors() ; ++k)
00291       LHS[k][InvReorder_[i]] = ReducedLHS[k][i];
00292 
00293   return(0);
00294 }
 All Classes Files Functions Variables Enumerations Friends