00001
00002
00003
00004
00005
00006
00007
00008
00009
00010
00011
00012
00013 #include "Epetra_MultiVector.h"
00014 #include "Epetra_Vector.h"
00015 #include "Epetra_Operator.h"
00016 #include "Epetra_Map.h"
00017 #include "Epetra_Comm.h"
00018
00019 #include "operatortest.h"
00020
00021 #define EPETRA_CHK(EXPR) { \
00022 int ierr = EXPR; \
00023 if (op_->Comm().MyPID() == 0 && ierr != 0) { \
00024 std::cerr << "EPETRA error caught in " << __FILE__ << ", line "<< __LINE__ << " : " \
00025 << #EXPR << " == " << ierr << std::endl; \
00026 } \
00027 }
00028
00029 using std::endl;
00030
00031 const double OperatorTest::symmetric_tol_ = 1e-13;
00032 const double OperatorTest::linear_tol_ = 1e-11;
00033
00034 OperatorTest::OperatorTest(const Epetra_Operator* op)
00035 : op_(op),
00036 use_inverse_(false)
00037 {
00038 verbose_ = 1;
00039 nrows_ = op->OperatorRangeMap().NumGlobalElements();
00040 ncols_ = op->OperatorDomainMap().NumGlobalElements();
00041 }
00042
00043 void OperatorTest::set_operator(const Epetra_Operator* op, bool use_inverse) {
00044 op_ = op;
00045 use_inverse_ = use_inverse;
00046 nrows_ = op->OperatorRangeMap().NumGlobalElements();
00047 ncols_ = op->OperatorDomainMap().NumGlobalElements();
00048 }
00049
00050 void OperatorTest::dump_stats(std::ostream& ostr, std::string operator_name) {
00051 bool symmetric = is_symmetric();
00052 bool linear = is_linear();
00053 bool non_destructive = is_non_destructive();
00054
00055 if (op_->Comm().MyPID() == 0) {
00056 ostr << "Statistics for operator \"" << operator_name << "\"" << endl;
00057 if (symmetric) {
00058 ostr << " Appears to be symmetric" << endl;
00059 } else {
00060 ostr << " Is non-symmetric" << endl;
00061 }
00062 if (linear) {
00063 ostr << " Appears to be linear" << endl;
00064 } else {
00065 ostr << " Is non-linear" << endl;
00066 }
00067 if (non_destructive) {
00068 ostr << " Appears to be behaving nicely" << endl;
00069 } else {
00070 ostr << " Overwrites the input vector!!!" << endl;
00071 }
00072 ostr << endl;
00073 }
00074 }
00075
00076
00077 bool OperatorTest::is_symmetric() {
00078 if (ncols_ != nrows_) {
00079 return false;
00080 } else {
00081 Epetra_Vector x(op_->OperatorDomainMap());
00082 Epetra_Vector y(op_->OperatorDomainMap());
00083 Epetra_Vector t(op_->OperatorRangeMap());
00084
00085 x.Random();
00086 y.Random();
00087
00088
00089 double s1;
00090 EPETRA_CHK(apply(y, t));
00091 EPETRA_CHK(t.Dot(x, &s1));
00092
00093
00094 double s2;
00095 EPETRA_CHK(apply(x, t));
00096 EPETRA_CHK(t.Dot(y, &s2));
00097
00098 double err = fabs(s1 - s2)/fabs(s1);
00099
00100 if (verbose_ > 0 && op_->Comm().MyPID() == 0)
00101 std::cout << "OperatorTest::is_symmetric: error=" << err << endl;
00102 return fabs(s1 - s2)/s1 < symmetric_tol_;
00103 }
00104 }
00105
00106 bool OperatorTest::is_linear() {
00107 Epetra_Vector x(op_->OperatorDomainMap());
00108 Epetra_Vector y(op_->OperatorDomainMap());
00109 Epetra_Vector t1(op_->OperatorDomainMap());
00110 Epetra_Vector z1(op_->OperatorRangeMap());
00111 Epetra_Vector z2(op_->OperatorRangeMap());
00112 Epetra_Vector t2(op_->OperatorRangeMap());
00113
00114 x.Random();
00115 y.Random();
00116
00117
00118 t1 = x;
00119 EPETRA_CHK(t1.Update(1.0, y, 1.0));
00120 EPETRA_CHK(apply(t1, z1));
00121
00122
00123 EPETRA_CHK(apply(x, z2));
00124 EPETRA_CHK(apply(y, t2));
00125 EPETRA_CHK(z2.Update(1.0, t2, 1.0));
00126
00127
00128 EPETRA_CHK(z1.Update(-1.0, z2, 1.0));
00129 double z1_nrm2;
00130 EPETRA_CHK(z1.Norm2(&z1_nrm2));
00131
00132 if (verbose_ > 0 && op_->Comm().MyPID() == 0)
00133 std::cout << "OperatorTest::is_linear: error=" << z1_nrm2 << endl;
00134 return z1_nrm2 < linear_tol_;
00135 }
00136
00137 bool OperatorTest::is_non_destructive() {
00138 Epetra_Vector x(op_->OperatorDomainMap());
00139 Epetra_Vector y(op_->OperatorRangeMap());
00140 Epetra_Vector t(x);
00141
00142 EPETRA_CHK(apply(x, y));
00143
00144
00145 EPETRA_CHK(t.Update(-1.0, x, 1.0));
00146 double t_nrm2;
00147 EPETRA_CHK(t.Norm2(&t_nrm2));
00148 if (verbose_ > 0 && op_->Comm().MyPID() == 0)
00149 std::cout << "OperatorTest::is_non_destructive: error=" << t_nrm2 << endl;
00150 return t_nrm2 == 0.0;
00151 }
00152
00153
00154 int OperatorTest::apply(Epetra_MultiVector& X, Epetra_MultiVector& Y) {
00155 if (use_inverse_)
00156 return op_->ApplyInverse(X, Y);
00157 else
00158 return op_->Apply(X, Y);
00159 }