Ifpack Package Browser (Single Doxygen Collection) Development
Loading...
Searching...
No Matches
Ifpack_DiagonalFilter.cpp
Go to the documentation of this file.
1/*@HEADER
2// ***********************************************************************
3//
4// Ifpack: Object-Oriented Algebraic Preconditioner Package
5// Copyright (2002) Sandia Corporation
6//
7// Under terms of Contract DE-AC04-94AL85000, there is a non-exclusive
8// license for use of this work by or on behalf of the U.S. Government.
9//
10// Redistribution and use in source and binary forms, with or without
11// modification, are permitted provided that the following conditions are
12// met:
13//
14// 1. Redistributions of source code must retain the above copyright
15// notice, this list of conditions and the following disclaimer.
16//
17// 2. Redistributions in binary form must reproduce the above copyright
18// notice, this list of conditions and the following disclaimer in the
19// documentation and/or other materials provided with the distribution.
20//
21// 3. Neither the name of the Corporation nor the names of the
22// contributors may be used to endorse or promote products derived from
23// this software without specific prior written permission.
24//
25// THIS SOFTWARE IS PROVIDED BY SANDIA CORPORATION "AS IS" AND ANY
26// EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
27// IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR
28// PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL SANDIA CORPORATION OR THE
29// CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL,
30// EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO,
31// PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR
32// PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF
33// LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING
34// NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
35// SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
36//
37// Questions? Contact Michael A. Heroux (maherou@sandia.gov)
38//
39// ***********************************************************************
40//@HEADER
41*/
42
43#include "Ifpack_ConfigDefs.h"
45#include "Epetra_ConfigDefs.h"
46#include "Epetra_RowMatrix.h"
47#include "Epetra_Comm.h"
48#include "Epetra_Map.h"
49#include "Epetra_MultiVector.h"
50#include "Epetra_Vector.h"
51
52//==============================================================================
53Ifpack_DiagonalFilter::Ifpack_DiagonalFilter(const Teuchos::RefCountPtr<Epetra_RowMatrix>& Matrix,
54 double AbsoluteThreshold,
55 double RelativeThreshold) :
56 A_(Matrix),
57 AbsoluteThreshold_(AbsoluteThreshold),
58 RelativeThreshold_(RelativeThreshold)
59{
60 using std::cout;
61 using std::endl;
62
63 Epetra_Time Time(Comm());
64
65 pos_.resize(NumMyRows());
66 val_.resize(NumMyRows());
67
68 std::vector<int> Indices(MaxNumEntries());
69 std::vector<double> Values(MaxNumEntries());
70 int NumEntries;
71
72 for (int MyRow = 0 ; MyRow < NumMyRows() ; ++MyRow) {
73
74 pos_[MyRow] = -1;
75 val_[MyRow] = 0.0;
76 int ierr = A_->ExtractMyRowCopy(MyRow, MaxNumEntries(), NumEntries,
77 &Values[0], &Indices[0]);
78 assert (ierr == 0);
79 (void) ierr; // avoid build warning when assert() is defined out
80
81 for (int i = 0 ; i < NumEntries ; ++i) {
82 if (Indices[i] == MyRow) {
83 pos_[MyRow] = i;
84 val_[MyRow] = Values[i] * (RelativeThreshold_ - 1) +
85 AbsoluteThreshold_ * EPETRA_SGN(Values[i]);
86 }
87 break;
88 }
89 }
90 cout << "TIME = " << Time.ElapsedTime() << endl;
91}
92
93//==============================================================================
95ExtractMyRowCopy(int MyRow, int Length, int& NumEntries,
96 double* Values, int* Indices) const
97{
98
99 IFPACK_CHK_ERR(A_->ExtractMyRowCopy(MyRow, Length, NumEntries,
100 Values,Indices));
101
102 if (pos_[MyRow] != -1)
103 Values[pos_[MyRow]] += val_[MyRow];
104
105 return(0);
106}
107
108//==============================================================================
110Multiply(bool TransA, const Epetra_MultiVector& X,
111 Epetra_MultiVector& Y) const
112{
113
114 if (X.NumVectors() != Y.NumVectors())
115 IFPACK_CHK_ERR(-2);
116
117 IFPACK_CHK_ERR(A_->Multiply(TransA, X, Y));
118
119 for (int v = 0 ; v < X.NumVectors() ; ++v)
120 for (int i = 0 ; i < NumMyRows() ; ++i)
121 Y[v][i] += val_[i] * X[v][i];
122
123
124 return(0);
125}
#define EPETRA_SGN(x)
#define IFPACK_CHK_ERR(ifpack_err)
int NumVectors() const
double ElapsedTime(void) const
std::vector< double > val_
Stores as additional diagonal contribution due to the filter.
double RelativeThreshold_
Multiplies A(i,i) by this value.
const Epetra_Comm & Comm() const
std::vector< int > pos_
Stores the position of the diagonal element, or -1 if not present.
Ifpack_DiagonalFilter(const Teuchos::RefCountPtr< Epetra_RowMatrix > &Matrix, double AbsoluteThreshold, double RelativeThreshold)
Constructor.
virtual int Multiply(bool TransA, const Epetra_MultiVector &X, Epetra_MultiVector &Y) const
Teuchos::RefCountPtr< Epetra_RowMatrix > A_
Pointer to the matrix to be filtered.
virtual int MaxNumEntries() const
Returns the maximum number of entries.
double AbsoluteThreshold_
This value (times the sgn(A(i,i)) is added to the diagonal elements.
virtual int ExtractMyRowCopy(int MyRow, int Length, int &NumEntries, double *Values, int *Indices) const
virtual int NumMyRows() const