IFPACK Development
Loading...
Searching...
No Matches
Ifpack_IHSS.cpp
1/*
2//@HEADER
3// ***********************************************************************
4//
5// Ifpack: Object-Oriented Algebraic Preconditioner Package
6// Copyright (2002) Sandia Corporation
7//
8// Under terms of Contract DE-AC04-94AL85000, there is a non-exclusive
9// license for use of this work by or on behalf of the U.S. Government.
10//
11// Redistribution and use in source and binary forms, with or without
12// modification, are permitted provided that the following conditions are
13// met:
14//
15// 1. Redistributions of source code must retain the above copyright
16// notice, this list of conditions and the following disclaimer.
17//
18// 2. Redistributions in binary form must reproduce the above copyright
19// notice, this list of conditions and the following disclaimer in the
20// documentation and/or other materials provided with the distribution.
21//
22// 3. Neither the name of the Corporation nor the names of the
23// contributors may be used to endorse or promote products derived from
24// this software without specific prior written permission.
25//
26// THIS SOFTWARE IS PROVIDED BY SANDIA CORPORATION "AS IS" AND ANY
27// EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
28// IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR
29// PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL SANDIA CORPORATION OR THE
30// CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL,
31// EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO,
32// PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR
33// PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF
34// LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING
35// NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
36// SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
37//
38// Questions? Contact Michael A. Heroux (maherou@sandia.gov)
39//
40// ***********************************************************************
41//@HEADER
42*/
43
44#include "Ifpack_IHSS.h"
45#include "Ifpack.h"
46#include "Ifpack_Utils.h"
47#include "Teuchos_ParameterList.hpp"
48#include "Teuchos_RefCountPtr.hpp"
49
50
51
52using Teuchos::RefCountPtr;
53using Teuchos::rcp;
54
55
56#ifdef HAVE_IFPACK_EPETRAEXT
57#include "EpetraExt_MatrixMatrix.h"
58
59
60Ifpack_IHSS::Ifpack_IHSS(Epetra_RowMatrix* A):
61 IsInitialized_(false),
62 IsComputed_(false),
63 Label_(),
64 EigMaxIters_(10),
65 EigRatio_(30.0),
66 LambdaMax_(-1.0),
67 Alpha_(-1.0),
68 NumSweeps_(1),
69
70 Time_(A->Comm())
71{
72 Epetra_CrsMatrix *Acrs=dynamic_cast<Epetra_CrsMatrix*>(A);
73 TEUCHOS_TEST_FOR_EXCEPT(!Acrs)
74 A_=rcp(Acrs,false);
75}
76
77void Ifpack_IHSS::Destroy(){
78}
79
80
81
82int Ifpack_IHSS::Initialize(){
83 EigMaxIters_ = List_.get("ihss: eigenvalue max iterations",EigMaxIters_);
84 EigRatio_ = List_.get("ihss: ratio eigenvalue", EigRatio_);
85 NumSweeps_ = List_.get("ihss: sweeps",NumSweeps_);
86
87 // Counters
88 IsInitialized_=true;
89 NumInitialize_++;
90 return 0;
91}
92
93int Ifpack_IHSS::SetParameters(Teuchos::ParameterList& parameterlist){
94 List_=parameterlist;
95 return 0;
96}
97
98
99int Ifpack_IHSS::Compute(){
100 if(!IsInitialized_) Initialize();
101
102 int rv;
103 Ifpack Factory;
104 Epetra_CrsMatrix *Askew=0,*Aherm=0;
105 Ifpack_Preconditioner *Pskew=0, *Pherm=0;
106 Time_.ResetStartTime();
107
108 // Create Aherm (w/o diagonal)
109 rv=EpetraExt::MatrixMatrix::Add(*A_,false,.5,*A_,true,.5,Aherm);
110 Aherm->FillComplete();
111 if(rv) IFPACK_CHK_ERR(-1);
112
113 // Grab Aherm's diagonal
114 Epetra_Vector avec(Aherm->RowMap());
115 IFPACK_CHK_ERR(Aherm->ExtractDiagonalCopy(avec));
116
117
118 // Compute alpha using the Bai, Golub & Ng 2003 formula, not the more multigrid-appropriate Hamilton, Benzi and Haber 2007.
119 // PowerMethod(Aherm, EigMaxIters_,LambdaMax_);
120 // Alpha_=LambdaMax_ / sqrt(EigRatio_);
121
122 // Try something more Hamilton inspired, using the maximum diagonal value of Aherm.
123 avec.MaxValue(&Alpha_);
124
125 // Add alpha to the diagonal of Aherm
126 for(int i=0;i<Aherm->NumMyRows();i++) avec[i]+=Alpha_;
127 IFPACK_CHK_ERR(Aherm->ReplaceDiagonalValues(avec));
128 Aherm_=rcp(Aherm);
129
130 // Compute Askew (and add diagonal)
131 Askew=new Epetra_CrsMatrix(Copy,A_->RowMap(),0);
132 rv=EpetraExt::MatrixMatrix::Add(*A_,false,.5,*A_,true,-.5,Askew);
133 if(rv) IFPACK_CHK_ERR(-2);
134
135#ifndef EPETRA_NO_32BIT_GLOBAL_INDICES
136 if(Askew->RowMap().GlobalIndicesInt()) {
137 for(int i=0;i<Askew->NumMyRows();i++) {
138 int gid=Askew->GRID(i);
139 Askew->InsertGlobalValues(gid,1,&Alpha_,&gid);
140 }
141 }
142 else
143#endif
144#ifndef EPETRA_NO_64BIT_GLOBAL_INDICES
145 if(Askew->RowMap().GlobalIndicesLongLong()) {
146 for(int i=0;i<Askew->NumMyRows();i++) {
147 long long gid=Askew->GRID64(i);
148 Askew->InsertGlobalValues(gid,1,&Alpha_,&gid);
149 }
150 }
151 else
152#endif
153 throw "Ifpack_IHSS::Compute: Unable to deduce GlobalIndices type";
154
155 Askew->FillComplete();
156 Askew_=rcp(Askew);
157
158 // Compute preconditioner for Aherm
159 Teuchos::ParameterList PLh=List_.sublist("ihss: hermetian list");
160 std::string htype=List_.get("ihss: hermetian type","ILU");
161 Pherm= Factory.Create(htype, Aherm);
162 Pherm->SetParameters(PLh);
163 IFPACK_CHK_ERR(Pherm->Compute());
164 Pherm_=rcp(Pherm);
165
166 // Compute preconditoner for Askew
167 Teuchos::ParameterList PLs=List_.sublist("ihss: skew hermetian list");
168 std::string stype=List_.get("ihss: skew hermetian type","ILU");
169 Pskew= Factory.Create(stype, Askew);
170 Pskew->SetParameters(PLs);
171 IFPACK_CHK_ERR(Pskew->Compute());
172 Pskew_=rcp(Pskew);
173
174 // Label
175 sprintf(Label_, "IFPACK IHSS (H,S)=(%s/%s)",htype.c_str(),stype.c_str());
176
177 // Counters
178 IsComputed_=true;
179 NumCompute_++;
180 ComputeTime_ += Time_.ElapsedTime();
181 return 0;
182}
183
184
185int Ifpack_IHSS::ApplyInverse(const Epetra_MultiVector& X, Epetra_MultiVector& Y) const{
186 if(!IsComputed_) return -1;
187 Time_.ResetStartTime();
188 bool initial_guess_is_zero=false;
189
190 // AztecOO gives X and Y pointing to the same memory location,
191 // need to create an auxiliary vector, Xcopy
192 Teuchos::RefCountPtr<const Epetra_MultiVector> Xcopy;
193 Epetra_MultiVector Temp(X);
194 if (X.Pointers()[0] == Y.Pointers()[0]){
195 Xcopy = Teuchos::rcp( new Epetra_MultiVector(X) );
196 // Since the user didn't give us anything better, our initial guess is zero.
197 Y.Scale(0.0);
198 initial_guess_is_zero=true;
199 }
200 else
201 Xcopy = Teuchos::rcp( &X, false );
202
203 Epetra_MultiVector T1(Y),T2(*Xcopy);
204
205 // Note: Since Aherm and Askew are actually (aI+H) and (aI+S) accordingly (to save memory),
206 // the application thereof needs to be a little different.
207 // The published algorithm is:
208 // temp = (aI+H)^{-1} [ (aI-S) y + x ]
209 // y = (aI+S)^{-1} [ (aI-H) temp + x ]
210 //
211 // But we're doing:
212 // temp = (aI+H)^{-1} [ 2a y - Shat y + x ]
213 // y = (aI+S)^{-1} [ 2 a temp - Hhat temp + x ]
214
215 for(int i=0;i<NumSweeps_;i++){
216 // temp = (aI+H)^{-1} [ 2a y - Shat y + x ]
217 if(!initial_guess_is_zero || i >0 ){
218 Askew_->Apply(Y,T1);
219 T2.Update(2*Alpha_,Y,-1,T1,1);
220 }
221 Pherm_->ApplyInverse(T2,Y);
222
223 // y = (aI+S)^{-1} [ 2 a temp - Hhat temp + x ]
224 Aherm_->Apply(Y,T1);
225 T2.Scale(1.0,*Xcopy);
226 T2.Update(2*Alpha_,Y,-1,T1,1.0);
227 Pskew_->ApplyInverse(T2,Y);
228 }
229
230 // Counter update
231 NumApplyInverse_++;
232 ApplyInverseTime_ += Time_.ElapsedTime();
233 return 0;
234}
235
236
237std::ostream& Ifpack_IHSS::Print(std::ostream& os) const{
238 using std::endl;
239
240 os<<"Ifpack_IHSS"<<endl;
241 os<<"-Hermetian preconditioner"<<endl;
242 os<<"-Skew Hermetian preconditioner"<<endl;
243 os<<endl;
244 return os;
245}
246
247
248double Ifpack_IHSS::Condest(const Ifpack_CondestType /* CT */,
249 const int /* MaxIters */,
250 const double /* Tol */,
251 Epetra_RowMatrix* /* Matrix_in */){
252 return -1.0;
253}
254
255
256int Ifpack_IHSS::PowerMethod(Epetra_Operator * Op,const int MaximumIterations, double& lambda_max)
257{
258
259 // this is a simple power method
260 lambda_max = 0.0;
261 double RQ_top, RQ_bottom, norm;
264 x.Random();
265 x.Norm2(&norm);
266 if (norm == 0.0) IFPACK_CHK_ERR(-1);
267
268 x.Scale(1.0 / norm);
269
270 for (int iter = 0; iter < MaximumIterations; ++iter)
271 {
272 Op->Apply(x, y);
273 IFPACK_CHK_ERR(y.Dot(x, &RQ_top));
274 IFPACK_CHK_ERR(x.Dot(x, &RQ_bottom));
275 lambda_max = RQ_top / RQ_bottom;
276 IFPACK_CHK_ERR(y.Norm2(&norm));
277 if (norm == 0.0) IFPACK_CHK_ERR(-1);
278 IFPACK_CHK_ERR(x.Update(1.0 / norm, y, 0.0));
279 }
280
281 return(0);
282}
283
284#endif
bool GlobalIndicesInt() const
bool GlobalIndicesLongLong() const
const Epetra_Map & RowMap() const
virtual int InsertGlobalValues(int GlobalRow, int NumEntries, const double *Values, const int *Indices)
int FillComplete(bool OptimizeDataStorage=true)
int GRID(int LRID_in) const
int NumMyRows() const
int Scale(double ScalarValue)
double ** Pointers() const
virtual int Apply(const Epetra_MultiVector &X, Epetra_MultiVector &Y) const=0
virtual const Epetra_Map & OperatorDomainMap() const=0
virtual const Epetra_Map & OperatorRangeMap() const=0
Ifpack_Preconditioner: basic class for preconditioning in Ifpack.
virtual int Compute()=0
Computes all it is necessary to apply the preconditioner.
virtual int SetParameters(Teuchos::ParameterList &List)=0
Sets all parameters for the preconditioner.
Ifpack: a function class to define Ifpack preconditioners.
Definition Ifpack.h:138
static Ifpack_Preconditioner * Create(EPrecType PrecType, Epetra_RowMatrix *Matrix, const int overlap=0, bool overrideSerialDefault=false)
Creates an instance of Ifpack_Preconditioner given the enum value of the preconditioner type (can not...
Definition Ifpack.cpp:253