Ifpack_DropFilter.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_DropFilter.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 
00039 //==============================================================================
00040 Ifpack_DropFilter::Ifpack_DropFilter(const Teuchos::RefCountPtr<Epetra_RowMatrix>& Matrix,
00041                      double DropTol) :
00042   A_(Matrix),
00043   DropTol_(DropTol),
00044   MaxNumEntries_(0),
00045   MaxNumEntriesA_(0),
00046   NumNonzeros_(0)
00047 {
00048   // use this filter only on serial matrices
00049   if (A_->Comm().NumProc() != 1) {
00050     cerr << "Ifpack_DropFilter can be used with Comm().NumProc() == 1" << endl;
00051     cerr << "only. This class is a tool for Ifpack_AdditiveSchwarz," << endl;
00052     cerr << "and it is not meant to be used otherwise." << endl;
00053     exit(EXIT_FAILURE);
00054   }
00055   
00056   if ((A_->NumMyRows() != A_->NumGlobalRows()) ||
00057       (A_->NumMyRows() != A_->NumMyCols()))
00058     IFPACK_CHK_ERRV(-2);
00059 
00060   NumRows_ = A_->NumMyRows();
00061   MaxNumEntriesA_ = A_->MaxNumEntries();
00062 
00063   NumEntries_.resize(NumRows_);
00064   Indices_.resize(MaxNumEntriesA_);
00065   Values_.resize(MaxNumEntriesA_);
00066 
00067   vector<int>    Ind(MaxNumEntriesA_);
00068   vector<double> Val(MaxNumEntriesA_);
00069 
00070   for (int i = 0 ; i < NumRows_ ; ++i) {
00071     NumEntries_[i] = MaxNumEntriesA_;
00072     int Nnz;
00073     IFPACK_CHK_ERRV(ExtractMyRowCopy(i,MaxNumEntriesA_,Nnz,
00074                      &Val[0], &Ind[0]));
00075  
00076     NumEntries_[i] = Nnz;
00077     NumNonzeros_ += Nnz;
00078     if (Nnz > MaxNumEntries_)
00079       MaxNumEntries_ = Nnz;
00080   }
00081 
00082 }
00083 
00084 //==============================================================================
00085 int Ifpack_DropFilter::
00086 ExtractMyRowCopy(int MyRow, int Length, int & NumEntries, 
00087          double *Values, int * Indices) const
00088 {
00089   if (Length < NumEntries_[MyRow])
00090     IFPACK_CHK_ERR(-1);
00091 
00092   int Nnz;
00093 
00094   IFPACK_CHK_ERR(A_->ExtractMyRowCopy(MyRow,MaxNumEntriesA_,Nnz,
00095                      &Values_[0],&Indices_[0]));
00096 
00097   // loop over all nonzero elements of row MyRow,
00098   // and drop elements below specified threshold.
00099   // Also, add value to zero diagonal
00100   int count = 0;
00101   for (int i = 0 ; i < Nnz ; ++i) {
00102 
00103     // if element is above specified tol, add to the
00104     // user's defined arrays. Check that we are not
00105     // exceeding allocated space. Do not drop any diagonal entry.
00106     if ((Indices_[i] == MyRow) || (IFPACK_ABS(Values_[i]) >= DropTol_)) {
00107       if (count == Length)
00108     IFPACK_CHK_ERR(-1);
00109       Values[count] = Values_[i];
00110       Indices[count] = Indices_[i];
00111       count++;
00112     }
00113   }
00114 
00115   NumEntries = count;
00116   return(0);
00117 }
00118 
00119 //==============================================================================
00120 int Ifpack_DropFilter::
00121 ExtractDiagonalCopy(Epetra_Vector & Diagonal) const
00122 {
00123   IFPACK_CHK_ERR(A_->ExtractDiagonalCopy(Diagonal));
00124   return(0);
00125 }
00126 
00127 //==============================================================================
00128 int Ifpack_DropFilter::
00129 Multiply(bool TransA, const Epetra_MultiVector& X, 
00130      Epetra_MultiVector& Y) const
00131 {
00132   // NOTE: I suppose that the matrix has been localized,
00133   // hence all maps are trivial.
00134   int NumVectors = X.NumVectors();
00135   if (NumVectors != Y.NumVectors())
00136     IFPACK_CHK_ERR(-1);
00137 
00138   Y.PutScalar(0.0);
00139 
00140   vector<int> Indices(MaxNumEntries_);
00141   vector<double> Values(MaxNumEntries_);
00142 
00143   for (int i = 0 ; i < NumRows_ ; ++i) {
00144 
00145     int Nnz;
00146     ExtractMyRowCopy(i,MaxNumEntries_,Nnz,
00147              &Values[0], &Indices[0]);
00148     if (!TransA) {
00149       // no transpose first
00150       for (int j = 0 ; j < NumVectors ; ++j) {
00151     for (int k = 0 ; k < Nnz ; ++k) {
00152       Y[j][i] += Values[k] * X[j][Indices[k]];
00153     }
00154       }
00155     }
00156     else {
00157       // transpose here
00158       for (int j = 0 ; j < NumVectors ; ++j) {
00159     for (int k = 0 ; k < Nnz ; ++k) {
00160       Y[j][Indices[k]] += Values[k] * X[j][i];
00161     }
00162       }
00163     }
00164   }
00165 
00166   return(0);
00167 }
00168 
00169 //==============================================================================
00170 int Ifpack_DropFilter::
00171 Solve(bool Upper, bool Trans, bool UnitDiagonal, 
00172       const Epetra_MultiVector& X, Epetra_MultiVector& Y) const
00173 {
00174   IFPACK_CHK_ERR(-99);
00175 }
00176 
00177 //==============================================================================
00178 int Ifpack_DropFilter::
00179 Apply(const Epetra_MultiVector& X, Epetra_MultiVector& Y) const
00180 {
00181   IFPACK_RETURN(Multiply(UseTranspose(),X,Y));
00182 }
00183 
00184 //==============================================================================
00185 int Ifpack_DropFilter::
00186 ApplyInverse(const Epetra_MultiVector& X, Epetra_MultiVector& Y) const
00187 {
00188   IFPACK_CHK_ERR(-99);
00189 }
00190 
00191 //==============================================================================
00192 int Ifpack_DropFilter::InvRowSums(Epetra_Vector& x) const
00193 {
00194   IFPACK_CHK_ERR(-1);
00195 }
00196 
00197 //==============================================================================
00198 int Ifpack_DropFilter::InvColSums(Epetra_Vector& x) const
00199 {
00200   IFPACK_CHK_ERR(-1);
00201 }
 All Classes Files Functions Variables Enumerations Friends
Generated on Wed Apr 13 10:05:23 2011 for IFPACK by  doxygen 1.6.3