AbstractLinAlgPack: C++ Interfaces For Vectors, Matrices And Related Linear Algebra Objects Version of the Day
AbstractLinAlgPack_TestMatrixSymSecant.cpp
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 // Redistribution and use in source and binary forms, with or without
00011 // modification, are permitted provided that the following conditions are
00012 // met:
00013 //
00014 // 1. Redistributions of source code must retain the above copyright
00015 // notice, this list of conditions and the following disclaimer.
00016 //
00017 // 2. Redistributions in binary form must reproduce the above copyright
00018 // notice, this list of conditions and the following disclaimer in the
00019 // documentation and/or other materials provided with the distribution.
00020 //
00021 // 3. Neither the name of the Corporation nor the names of the
00022 // contributors may be used to endorse or promote products derived from
00023 // this software without specific prior written permission.
00024 //
00025 // THIS SOFTWARE IS PROVIDED BY SANDIA CORPORATION "AS IS" AND ANY
00026 // EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
00027 // IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR
00028 // PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL SANDIA CORPORATION OR THE
00029 // CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL,
00030 // EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO,
00031 // PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR
00032 // PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF
00033 // LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING
00034 // NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
00035 // SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
00036 //
00037 // Questions? Contact Roscoe A. Bartlett (rabartl@sandia.gov) 
00038 // 
00039 // ***********************************************************************
00040 // @HEADER
00041 
00042 #include <ostream>
00043 #include <iomanip>
00044 
00045 #include <math.h>
00046 
00047 #include "AbstractLinAlgPack_TestMatrixSymSecant.hpp"
00048 #include "AbstractLinAlgPack_MatrixOp.hpp"
00049 #include "AbstractLinAlgPack_MatrixNonsing.hpp"
00050 #include "AbstractLinAlgPack_VectorSpace.hpp"
00051 #include "AbstractLinAlgPack_VectorMutable.hpp"
00052 #include "AbstractLinAlgPack_VectorOut.hpp"
00053 #include "AbstractLinAlgPack_VectorStdOps.hpp"
00054 #include "AbstractLinAlgPack_LinAlgOpPack.hpp"
00055 #include "AbstractLinAlgPack_assert_print_nan_inf.hpp"
00056 
00057 bool AbstractLinAlgPack::TestMatrixSymSecant(
00058   const MatrixOp        &B
00059   ,const Vector       &s
00060   ,const Vector       &y
00061   ,value_type               warning_tol
00062   ,value_type               error_tol
00063   ,bool                     print_all_warnings
00064   ,std::ostream             *out
00065   ,bool                     trase
00066   )
00067 {
00068   using std::setw;
00069   using std::endl;
00070   using AbstractLinAlgPack::sum;
00071   using AbstractLinAlgPack::V_InvMtV;
00072   using AbstractLinAlgPack::assert_print_nan_inf;
00073   using LinAlgOpPack::V_MtV;
00074 
00075   bool success = true;
00076   const char
00077     sep_line[] = "\n-----------------------------------------------------------------\n";
00078 
00079   // Check the secant property (B * s = y)
00080   {
00081     if( out && trase )
00082       *out
00083         << sep_line
00084         << "\n|(sum(B*s)-sum(y))/||y||inf| = ";
00085     VectorSpace::vec_mut_ptr_t
00086       Bs = y.space().create_member();
00087     V_MtV( Bs.get(), B, BLAS_Cpp::no_trans, s );
00088     const value_type
00089       sum_Bs = sum(*Bs),
00090       sum_y  = sum(y),
00091       nrm_y  = y.norm_inf(),
00092       err    = ::fabs( (sum_Bs - sum_y) / nrm_y );
00093     if( out && trase )
00094       *out
00095         <<"|("<<sum_Bs<<"-"<<sum_y<<")/"<<nrm_y<<"| = " << err << std::endl;
00096     if( err >= error_tol ) {
00097       if( out && trase )
00098         *out
00099           << "Error, above error = " << err << " >= error_tol = " << error_tol
00100           << "\nThe test has failed!\n";
00101       if(out && print_all_warnings) {
00102         *out
00103           << "\ns =\n" << s 
00104           << "\ny =\n" << y 
00105           << "\nBs =\n" << *Bs
00106           << endl;
00107       }
00108       success = false;
00109     }
00110     else if( err >= warning_tol ) {
00111       if( out && trase )
00112         *out
00113           << "Warning!, above error = " << err << " >= warning_tol = " << warning_tol << std::endl;
00114     }
00115   }
00116   // Check the secant property (s = inv(B)*y)
00117   const MatrixNonsing
00118     *B_nonsing = dynamic_cast<const MatrixNonsing*>(&B);
00119   if( B_nonsing ) {
00120     if( out && trase )
00121       *out
00122         << sep_line
00123         << "\n|(sum(inv(B)*y)-sum(s))/||s||inf| = ";
00124     VectorSpace::vec_mut_ptr_t
00125       InvBy = s.space().create_member();
00126     V_InvMtV( InvBy.get(), *B_nonsing, BLAS_Cpp::no_trans, y );
00127     const value_type
00128       sum_InvBy = sum(*InvBy),
00129       sum_s     = sum(s),
00130       nrm_s     = s.norm_inf(),
00131       err       = ::fabs( (sum_InvBy - sum_s) / nrm_s );
00132     if( out && trase )
00133       *out
00134         <<"|("<<sum_InvBy<<"-"<<sum_s<<")/"<<nrm_s<<"| = " << err << std::endl;
00135     if( err >= error_tol ) {
00136       if( out && trase )
00137         *out
00138           << "Error, above error = " << err << " >= error_tol = " << error_tol
00139           << "\nThe test has failed!\n";
00140       if(out && print_all_warnings) {
00141         *out
00142           << "\ns =\n" << s 
00143           << "\ny =\n" << y 
00144           << "\ninv(B)*y =\n" << *InvBy
00145           << endl;
00146       }
00147       success = false;
00148     }
00149     else if( err >= warning_tol ) {
00150       if( out && trase )
00151         *out
00152           << "Warning!, above error = " << err << " >= warning_tol = " << warning_tol << std::endl;
00153     }
00154   }
00155   if( out && trase )
00156     *out << sep_line;
00157   return success;
00158 }
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends