IFPACK  Development
 All Classes Namespaces Files Functions Variables Enumerations Friends Pages
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 
52 using Teuchos::RefCountPtr;
53 using Teuchos::rcp;
54 
55 
56 #ifdef HAVE_IFPACK_EPETRAEXT
57 #include "EpetraExt_MatrixMatrix.h"
58 
59 
60 Ifpack_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 
77 void Ifpack_IHSS::Destroy(){
78 }
79 
80 
81 
82 int 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 
93 int Ifpack_IHSS::SetParameters(Teuchos::ParameterList& parameterlist){
94  List_=parameterlist;
95  return 0;
96 }
97 
98 
99 int 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 
185 int 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 
237 std::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 
248 double 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 
256 int 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 GlobalIndicesLongLong() const
double ElapsedTime(void) const
virtual const Epetra_Map & OperatorDomainMap() const =0
bool GlobalIndicesInt() const
int FillComplete(bool OptimizeDataStorage=true)
int Initialize()
Initialize L and U with values from user matrix A.
Definition: Ifpack_ICT.cpp:138
const Epetra_Map & RowMap() const
virtual int Apply(const Epetra_MultiVector &X, Epetra_MultiVector &Y) const =0
virtual const Epetra_Map & OperatorRangeMap() const =0
virtual int SetParameters(Teuchos::ParameterList &List)=0
Sets all parameters for the preconditioner.
Ifpack_Preconditioner: basic class for preconditioning in Ifpack.
int NumMyRows() const
virtual int InsertGlobalValues(int GlobalRow, int NumEntries, const double *Values, const int *Indices)
Ifpack: a function class to define Ifpack preconditioners.
Definition: Ifpack.h:138
int GRID(int LRID_in) const
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
virtual int Compute()=0
Computes all it is necessary to apply the preconditioner.
void ResetStartTime(void)