src/operatortest.cpp

Go to the documentation of this file.
00001 //
00002 // C++ Implementation: operatortest
00003 //
00004 // Description:
00005 //
00006 //
00007 // Author: Roman Geus <roman.geus@psi.ch>, (C) 2004
00008 //
00009 // Copyright: See COPYING file that comes with this distribution
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         // s1 = x^T*OP*y
00089         double s1;
00090         EPETRA_CHK(apply(y, t));
00091         EPETRA_CHK(t.Dot(x, &s1));
00092 
00093         // s2 = y^T*OP*x
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     // z1 = OP * (x + y)
00118     t1 = x;
00119     EPETRA_CHK(t1.Update(1.0, y, 1.0));
00120     EPETRA_CHK(apply(t1, z1));
00121 
00122     // z2 = OP*x + OP*y
00123     EPETRA_CHK(apply(x, z2));
00124     EPETRA_CHK(apply(y, t2));
00125     EPETRA_CHK(z2.Update(1.0, t2, 1.0));
00126 
00127     // test z1 == z2
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     // test if x has not changed
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 }

Generated on Fri Oct 26 13:35:12 2007 for FEMAXX (Finite Element Maxwell Eigensolver) by  doxygen 1.4.7