AbstractLinAlgPack_TestMatrixSymSecant.cpp

Go to the documentation of this file.
00001 // @HEADER
00002 // ***********************************************************************
00003 // 
00004 // Moocho: Multi-functional Object-Oriented arCHitecture for Optimization
00005 //                  Copyright (2003) Sandia Corporation
00006 // 
00007 // Under terms of Contract DE-AC04-94AL85000, there is a non-exclusive
00008 // license for use of this work by or on behalf of the U.S. Government.
00009 // 
00010 // This library is free software; you can redistribute it and/or modify
00011 // it under the terms of the GNU Lesser General Public License as
00012 // published by the Free Software Foundation; either version 2.1 of the
00013 // License, or (at your option) any later version.
00014 //  
00015 // This library is distributed in the hope that it will be useful, but
00016 // WITHOUT ANY WARRANTY; without even the implied warranty of
00017 // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
00018 // Lesser General Public License for more details.
00019 //  
00020 // You should have received a copy of the GNU Lesser General Public
00021 // License along with this library; if not, write to the Free Software
00022 // Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307
00023 // USA
00024 // Questions? Contact Roscoe A. Bartlett (rabartl@sandia.gov) 
00025 // 
00026 // ***********************************************************************
00027 // @HEADER
00028 
00029 #include <ostream>
00030 #include <iomanip>
00031 
00032 #include <math.h>
00033 
00034 #include "AbstractLinAlgPack_TestMatrixSymSecant.hpp"
00035 #include "AbstractLinAlgPack_MatrixOp.hpp"
00036 #include "AbstractLinAlgPack_MatrixNonsing.hpp"
00037 #include "AbstractLinAlgPack_VectorSpace.hpp"
00038 #include "AbstractLinAlgPack_VectorMutable.hpp"
00039 #include "AbstractLinAlgPack_VectorOut.hpp"
00040 #include "AbstractLinAlgPack_VectorStdOps.hpp"
00041 #include "AbstractLinAlgPack_LinAlgOpPack.hpp"
00042 #include "AbstractLinAlgPack_assert_print_nan_inf.hpp"
00043 
00044 bool AbstractLinAlgPack::TestMatrixSymSecant(
00045   const MatrixOp        &B
00046   ,const Vector       &s
00047   ,const Vector       &y
00048   ,value_type               warning_tol
00049   ,value_type               error_tol
00050   ,bool                     print_all_warnings
00051   ,std::ostream             *out
00052   ,bool                     trase
00053   )
00054 {
00055   using std::setw;
00056   using std::endl;
00057   using AbstractLinAlgPack::sum;
00058   using AbstractLinAlgPack::V_InvMtV;
00059   using AbstractLinAlgPack::assert_print_nan_inf;
00060   using LinAlgOpPack::V_MtV;
00061 
00062   bool success = true;
00063   const char
00064     sep_line[] = "\n-----------------------------------------------------------------\n";
00065 
00066   // Check the secant property (B * s = y)
00067   {
00068     if( out && trase )
00069       *out
00070         << sep_line
00071         << "\n|(sum(B*s)-sum(y))/||y||inf| = ";
00072     VectorSpace::vec_mut_ptr_t
00073       Bs = y.space().create_member();
00074     V_MtV( Bs.get(), B, BLAS_Cpp::no_trans, s );
00075     const value_type
00076       sum_Bs = sum(*Bs),
00077       sum_y  = sum(y),
00078       nrm_y  = y.norm_inf(),
00079       err    = ::fabs( (sum_Bs - sum_y) / nrm_y );
00080     if( out && trase )
00081       *out
00082         <<"|("<<sum_Bs<<"-"<<sum_y<<")/"<<nrm_y<<"| = " << err << std::endl;
00083     if( err >= error_tol ) {
00084       if( out && trase )
00085         *out
00086           << "Error, above error = " << err << " >= error_tol = " << error_tol
00087           << "\nThe test has failed!\n";
00088       if(out && print_all_warnings) {
00089         *out
00090           << "\ns =\n" << s 
00091           << "\ny =\n" << y 
00092           << "\nBs =\n" << *Bs
00093           << endl;
00094       }
00095       success = false;
00096     }
00097     else if( err >= warning_tol ) {
00098       if( out && trase )
00099         *out
00100           << "Warning!, above error = " << err << " >= warning_tol = " << warning_tol << std::endl;
00101     }
00102   }
00103   // Check the secant property (s = inv(B)*y)
00104   const MatrixNonsing
00105     *B_nonsing = dynamic_cast<const MatrixNonsing*>(&B);
00106   if( B_nonsing ) {
00107     if( out && trase )
00108       *out
00109         << sep_line
00110         << "\n|(sum(inv(B)*y)-sum(s))/||s||inf| = ";
00111     VectorSpace::vec_mut_ptr_t
00112       InvBy = s.space().create_member();
00113     V_InvMtV( InvBy.get(), *B_nonsing, BLAS_Cpp::no_trans, y );
00114     const value_type
00115       sum_InvBy = sum(*InvBy),
00116       sum_s     = sum(s),
00117       nrm_s     = s.norm_inf(),
00118       err       = ::fabs( (sum_InvBy - sum_s) / nrm_s );
00119     if( out && trase )
00120       *out
00121         <<"|("<<sum_InvBy<<"-"<<sum_s<<")/"<<nrm_s<<"| = " << err << std::endl;
00122     if( err >= error_tol ) {
00123       if( out && trase )
00124         *out
00125           << "Error, above error = " << err << " >= error_tol = " << error_tol
00126           << "\nThe test has failed!\n";
00127       if(out && print_all_warnings) {
00128         *out
00129           << "\ns =\n" << s 
00130           << "\ny =\n" << y 
00131           << "\ninv(B)*y =\n" << *InvBy
00132           << endl;
00133       }
00134       success = false;
00135     }
00136     else if( err >= warning_tol ) {
00137       if( out && trase )
00138         *out
00139           << "Warning!, above error = " << err << " >= warning_tol = " << warning_tol << std::endl;
00140     }
00141   }
00142   if( out && trase )
00143     *out << sep_line;
00144   return success;
00145 }

Generated on Tue Oct 20 12:51:44 2009 for MOOCHO (Single Doxygen Collection) by doxygen 1.4.7