44 #include "Ifpack_IHSS.h" 47 #include "Teuchos_ParameterList.hpp" 48 #include "Teuchos_RefCountPtr.hpp" 52 using Teuchos::RefCountPtr;
56 #ifdef HAVE_IFPACK_EPETRAEXT 57 #include "EpetraExt_MatrixMatrix.h" 60 Ifpack_IHSS::Ifpack_IHSS(Epetra_RowMatrix* A):
61 IsInitialized_(false),
72 Epetra_CrsMatrix *Acrs=
dynamic_cast<Epetra_CrsMatrix*
>(A);
73 TEUCHOS_TEST_FOR_EXCEPT(!Acrs)
77 void Ifpack_IHSS::Destroy(){
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_);
93 int Ifpack_IHSS::SetParameters(Teuchos::ParameterList& parameterlist){
99 int Ifpack_IHSS::Compute(){
104 Epetra_CrsMatrix *Askew=0,*Aherm=0;
106 Time_.ResetStartTime();
109 rv=EpetraExt::MatrixMatrix::Add(*A_,
false,.5,*A_,
true,.5,Aherm);
110 Aherm->FillComplete();
111 if(rv) IFPACK_CHK_ERR(-1);
114 Epetra_Vector avec(Aherm->RowMap());
115 IFPACK_CHK_ERR(Aherm->ExtractDiagonalCopy(avec));
123 avec.MaxValue(&Alpha_);
126 for(
int i=0;i<Aherm->NumMyRows();i++) avec[i]+=Alpha_;
127 IFPACK_CHK_ERR(Aherm->ReplaceDiagonalValues(avec));
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);
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);
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);
153 throw "Ifpack_IHSS::Compute: Unable to deduce GlobalIndices type";
155 Askew->FillComplete();
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);
163 IFPACK_CHK_ERR(Pherm->Compute());
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);
171 IFPACK_CHK_ERR(Pskew->
Compute());
175 sprintf(Label_,
"IFPACK IHSS (H,S)=(%s/%s)",htype.c_str(),stype.c_str());
180 ComputeTime_ += Time_.ElapsedTime();
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;
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) );
198 initial_guess_is_zero=
true;
201 Xcopy = Teuchos::rcp( &X,
false );
203 Epetra_MultiVector T1(Y),T2(*Xcopy);
215 for(
int i=0;i<NumSweeps_;i++){
217 if(!initial_guess_is_zero || i >0 ){
219 T2.Update(2*Alpha_,Y,-1,T1,1);
221 Pherm_->ApplyInverse(T2,Y);
225 T2.Scale(1.0,*Xcopy);
226 T2.Update(2*Alpha_,Y,-1,T1,1.0);
227 Pskew_->ApplyInverse(T2,Y);
232 ApplyInverseTime_ += Time_.ElapsedTime();
237 std::ostream& Ifpack_IHSS::Print(std::ostream& os)
const{
240 os<<
"Ifpack_IHSS"<<endl;
241 os<<
"-Hermetian preconditioner"<<endl;
242 os<<
"-Skew Hermetian preconditioner"<<endl;
248 double Ifpack_IHSS::Condest(
const Ifpack_CondestType CT,
251 Epetra_RowMatrix* Matrix_in){
256 int Ifpack_IHSS::PowerMethod(Epetra_Operator * Op,
const int MaximumIterations,
double& lambda_max)
261 double RQ_top, RQ_bottom, norm;
262 Epetra_Vector x(Op->OperatorDomainMap());
263 Epetra_Vector y(Op->OperatorRangeMap());
266 if (norm == 0.0) IFPACK_CHK_ERR(-1);
270 for (
int iter = 0; iter < MaximumIterations; ++iter)
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));
int Initialize()
Initialize L and U with values from user matrix A.
virtual int SetParameters(Teuchos::ParameterList &List)=0
Sets all parameters for the preconditioner.
Ifpack_Preconditioner: basic class for preconditioning in Ifpack.
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...
virtual int Compute()=0
Computes all it is necessary to apply the preconditioner.