IFPACK Development
Ifpack_DiagonalFilter.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 // Redistribution and use in source and binary forms, with or without
00011 // modification, are permitted provided that the following conditions are
00012 // met:
00013 //
00014 // 1. Redistributions of source code must retain the above copyright
00015 // notice, this list of conditions and the following disclaimer.
00016 //
00017 // 2. Redistributions in binary form must reproduce the above copyright
00018 // notice, this list of conditions and the following disclaimer in the
00019 // documentation and/or other materials provided with the distribution.
00020 //
00021 // 3. Neither the name of the Corporation nor the names of the
00022 // contributors may be used to endorse or promote products derived from
00023 // this software without specific prior written permission.
00024 //
00025 // THIS SOFTWARE IS PROVIDED BY SANDIA CORPORATION "AS IS" AND ANY
00026 // EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
00027 // IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR
00028 // PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL SANDIA CORPORATION OR THE
00029 // CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL,
00030 // EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO,
00031 // PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR
00032 // PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF
00033 // LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING
00034 // NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
00035 // SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
00036 //
00037 // Questions? Contact Michael A. Heroux (maherou@sandia.gov)
00038 //
00039 // ***********************************************************************
00041 */
00042
00043 #include "Ifpack_ConfigDefs.h"
00044 #include "Ifpack_DiagonalFilter.h"
00045 #include "Epetra_ConfigDefs.h"
00046 #include "Epetra_RowMatrix.h"
00047 #include "Epetra_Comm.h"
00048 #include "Epetra_Map.h"
00049 #include "Epetra_MultiVector.h"
00050 #include "Epetra_Vector.h"
00051
00052 //==============================================================================
00053 Ifpack_DiagonalFilter::Ifpack_DiagonalFilter(const Teuchos::RefCountPtr<Epetra_RowMatrix>& Matrix,
00054                          double AbsoluteThreshold,
00055                          double RelativeThreshold) :
00056   A_(Matrix),
00057   AbsoluteThreshold_(AbsoluteThreshold),
00058   RelativeThreshold_(RelativeThreshold)
00059 {
00060   Epetra_Time Time(Comm());
00061
00062   pos_.resize(NumMyRows());
00063   val_.resize(NumMyRows());
00064
00065   std::vector<int> Indices(MaxNumEntries());
00066   std::vector<double> Values(MaxNumEntries());
00067   int NumEntries;
00068
00069   for (int MyRow = 0 ; MyRow < NumMyRows() ; ++MyRow) {
00070
00071     pos_[MyRow] = -1;
00072     val_[MyRow] = 0.0;
00073     int ierr = A_->ExtractMyRowCopy(MyRow, MaxNumEntries(), NumEntries,
00074                     &Values[0], &Indices[0]);
00075     assert (ierr == 0);
00076
00077     for (int i = 0 ; i < NumEntries ; ++i) {
00078       if (Indices[i] == MyRow) {
00079     pos_[MyRow] = i;
00080     val_[MyRow] = Values[i] * (RelativeThreshold_ - 1) +
00081       AbsoluteThreshold_ * EPETRA_SGN(Values[i]);
00082       }
00083       break;
00084     }
00085   }
00086   cout << "TIME = " << Time.ElapsedTime() << endl;
00087 }
00088
00089 //==============================================================================
00090 int Ifpack_DiagonalFilter::
00091 ExtractMyRowCopy(int MyRow, int Length, int& NumEntries,
00092          double* Values, int* Indices) const
00093 {
00094
00095   IFPACK_CHK_ERR(A_->ExtractMyRowCopy(MyRow, Length, NumEntries,
00096                      Values,Indices));
00097
00098   if (pos_[MyRow] != -1)
00099     Values[pos_[MyRow]] += val_[MyRow];
00100
00101   return(0);
00102 }
00103
00104 //==============================================================================
00105 int Ifpack_DiagonalFilter::
00106 Multiply(bool TransA, const Epetra_MultiVector& X,
00107      Epetra_MultiVector& Y) const
00108 {
00109
00110   if (X.NumVectors() != Y.NumVectors())
00111     IFPACK_CHK_ERR(-2);
00112
00113   IFPACK_CHK_ERR(A_->Multiply(TransA, X, Y));
00114
00115   for (int v = 0 ; v < X.NumVectors() ; ++v)
00116     for (int i = 0 ; i < NumMyRows() ; ++i)
00117       Y[v][i] += val_[i] * X[v][i];
00118
00119
00120   return(0);
00121 }