00001
00002
00003
00004
00005
00006
00007
00008
00009
00010
00011
00012 #ifndef DIAGONAL_PRECONDITIONER_H
00013 #define DIAGONAL_PRECONDITIONER_H
00014
00015 #include "Epetra_Operator.h"
00016 #include "Epetra_Vector.h"
00017 #include "Epetra_CrsMatrix.h"
00018
00019 class Epetra_MultiVector;
00020 class Epetra_Map;
00021 class Epetra_Comm;
00022
00026 class DiagonalPreconditioner : public Epetra_Operator {
00027 public:
00028 DiagonalPreconditioner(const Epetra_CrsMatrix& mat)
00029 : inv_diag_(new Epetra_Vector(mat.RowMap(), false)),
00030 domain_map_(mat.OperatorDomainMap()),
00031 range_map_(mat.OperatorRangeMap()),
00032 comm_(mat.Comm()),
00033 use_transpose_(false)
00034 {
00035 mat.ExtractDiagonalCopy(*inv_diag_);
00036 for (int i = 0; i < inv_diag_->Map().NumMyElements(); i ++)
00037 (*inv_diag_)[i] = 1.0 / (*inv_diag_)[i];
00038 }
00039
00040 ~DiagonalPreconditioner() {
00041 delete inv_diag_;
00042 }
00043
00044 int SetUseTranspose(bool use_transpose) {
00045 use_transpose_ = use_transpose;
00046 return 0;
00047 }
00048
00049 int Apply(const Epetra_MultiVector& X, Epetra_MultiVector& Y) const {
00050 if (!X.Map().SameAs(OperatorDomainMap()))
00051 EPETRA_CHK_ERR(-1);
00052 if (!Y.Map().SameAs(OperatorRangeMap()))
00053 EPETRA_CHK_ERR(-2);
00054 if (Y.NumVectors() != X.NumVectors())
00055 EPETRA_CHK_ERR(-3);
00056 if (!domain_map_.SameAs(range_map_))
00057 EPETRA_CHK_ERR(-4);
00058
00059 return Y.Multiply(1.0, *inv_diag_, X, 0.0);
00060 }
00061
00062 int ApplyInverse(const Epetra_MultiVector& X, Epetra_MultiVector& Y) const {
00063 return Apply(X, Y);
00064 }
00065
00066 double NormInf() const {
00067 return 0.0;
00068 }
00069
00070 const char* Label() const {
00071 return "Identity operator";
00072 }
00073
00074 bool UseTranspose() const {
00075 return use_transpose_;
00076 }
00077
00078 bool HasNormInf() const {
00079 return false;
00080 }
00081
00082 const Epetra_Comm& Comm() const {
00083 return comm_;
00084 }
00085
00086 const Epetra_Map & OperatorDomainMap() const {
00087 return domain_map_;
00088 }
00089
00090 const Epetra_Map & OperatorRangeMap() const {
00091 return range_map_;
00092 }
00093
00094 private:
00095 Epetra_Vector* inv_diag_;
00096 const Epetra_Map& domain_map_;
00097 const Epetra_Map& range_map_;
00098 const Epetra_Comm& comm_;
00099 bool use_transpose_;
00100 };
00101
00102 #endif