MoochoPack : Framework for Large-Scale Optimization Algorithms Version of the Day
MoochoPack_DampenCrossTermStd_Step.cpp
00001 #if 0
00002 
00003 // @HEADER
00004 // ***********************************************************************
00005 // 
00006 // Moocho: Multi-functional Object-Oriented arCHitecture for Optimization
00007 //                  Copyright (2003) Sandia Corporation
00008 // 
00009 // Under terms of Contract DE-AC04-94AL85000, there is a non-exclusive
00010 // license for use of this work by or on behalf of the U.S. Government.
00011 // 
00012 // Redistribution and use in source and binary forms, with or without
00013 // modification, are permitted provided that the following conditions are
00014 // met:
00015 //
00016 // 1. Redistributions of source code must retain the above copyright
00017 // notice, this list of conditions and the following disclaimer.
00018 //
00019 // 2. Redistributions in binary form must reproduce the above copyright
00020 // notice, this list of conditions and the following disclaimer in the
00021 // documentation and/or other materials provided with the distribution.
00022 //
00023 // 3. Neither the name of the Corporation nor the names of the
00024 // contributors may be used to endorse or promote products derived from
00025 // this software without specific prior written permission.
00026 //
00027 // THIS SOFTWARE IS PROVIDED BY SANDIA CORPORATION "AS IS" AND ANY
00028 // EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
00029 // IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR
00030 // PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL SANDIA CORPORATION OR THE
00031 // CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL,
00032 // EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO,
00033 // PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR
00034 // PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF
00035 // LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING
00036 // NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
00037 // SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
00038 //
00039 // Questions? Contact Roscoe A. Bartlett (rabartl@sandia.gov) 
00040 // 
00041 // ***********************************************************************
00042 // @HEADER
00043 
00044 #include <ostream>
00045 #include <sstream>
00046 #include <limits>
00047 
00048 #include "MoochoPack_DampenCrossTermStd_Step.hpp"
00049 #include "MoochoPack_moocho_algo_conversion.hpp"
00050 #include "IterationPack_print_algorithm_step.hpp"
00051 #include "ConstrainedOptPack/src/VectorWithNorms.h"
00052 #include "AbstractLinAlgPack/src/MatrixWithOpFactorized.hpp"
00053 #include "DenseLinAlgPack_LinAlgOpPack.hpp"
00054 #include "DenseLinAlgPack_DVectorClass.hpp"
00055 #include "DenseLinAlgPack_DVectorOut.hpp"
00056 
00057 MoochoPack::DampenCrossTermStd_Step::DampenCrossTermStd_Step(const value_type& frac_descent)
00058   : frac_descent_(frac_descent)
00059 {}
00060 
00061 bool MoochoPack::DampenCrossTermStd_Step::do_step(Algorithm& _algo
00062   , poss_type step_poss, IterationPack::EDoStepType type, poss_type assoc_step_poss)
00063 {
00064   using AbstractLinAlgPack::V_InvMtV;
00065   using DenseLinAlgPack::norm_inf;
00066   using DenseLinAlgPack::dot;
00067 
00068   NLPAlgo &algo = rsqp_algo(_algo);
00069   NLPAlgoState  &s    = algo.rsqp_state();
00070 
00071   EJournalOutputLevel olevel = algo.algo_cntr().journal_output_level();
00072   std::ostream& out = algo.track().journal_out();
00073 
00074   // print step header.
00075   if( static_cast<int>(olevel) >= static_cast<int>(PRINT_ALGORITHM_STEPS) ) {
00076     using IterationPack::print_algorithm_step;
00077     print_algorithm_step( _algo, step_poss, type, assoc_step_poss, out );
00078   }
00079 
00080   if( s.w().updated_k(0) ) {
00081 
00082     // inv(rHL_k) * rGf_k
00083     const DVectorSlice rGf_k = s.rGf().get_k(0)();
00084     DVector Inv_rHL_rGf;
00085     V_InvMtV( &Inv_rHL_rGf, dynamic_cast<MatrixWithOpFactorized&>(s.rHL().get_k(0))
00086       , BLAS_Cpp::no_trans, rGf_k );
00087     
00088     const value_type
00089       small_num     = 1e-20,
00090       rGfT_Inv_rHL_rGf  = dot( Inv_rHL_rGf(), rGf_k ),          // rGf_k'*inv(rHL_k)*rGf_k
00091       rGfT_Inv_rHL_w    = dot( Inv_rHL_rGf(), s.w().get_k(0)() ),   // rGf_k'*inv(rHL_k)*w_k
00092       term        = -(1.0-frac_descent()) * (rGfT_Inv_rHL_rGf + 2*small_num)
00093                   / (rGfT_Inv_rHL_w + small_num);
00094 
00095     if( rGfT_Inv_rHL_w >= 0.0 ) {
00096       // We know that the descent property will be satisfied for all zeta_k > 0
00097       // so set zeta_k = 1
00098       s.zeta().set_k(0) = 1.0;
00099     }
00100     else {
00101       // For some zeta_k > 0 the descent property will be violated so we may have to
00102       // cut zeta_k back from 1.
00103       s.zeta().set_k(0) = std::_MIN( term, 1.0 );
00104     }
00105 
00106     if( static_cast<int>(olevel) >= static_cast<int>(PRINT_ALGORITHM_STEPS) ) {
00107       out << "\nterm1 = rGf_k'*inv(rHL_k)*rGf_k                = "  <<  rGfT_Inv_rHL_rGf;
00108       out << "\nterm2 = rGf_k'*inv(rHL_k)*w_k                  = "  <<  rGfT_Inv_rHL_w;
00109       out << "\n(1-frac_descent)*(term1+2*small)/(term2+small) = "  <<  term;
00110       out << "\nzeta_k                                         = "  <<  s.zeta().get_k(0)
00111         << std::endl;
00112     }
00113 
00114     if( static_cast<int>(olevel) >= static_cast<int>(PRINT_VECTORS) ) {
00115       out << "\ninv(rHL_k)*rGf_k = "  <<  Inv_rHL_rGf();
00116     }
00117 
00118     if( rGfT_Inv_rHL_rGf < 0.0 ) {
00119       std::ostringstream omsg;
00120       omsg
00121         << "Error, rGf_k'*inv(rHL_k)*rGf_k = " << rGfT_Inv_rHL_rGf << " < 0.0 and therefore "
00122         << "the reduced Hessian rHL_k can not be positive definite";
00123       if( (int)(olevel) >= (int)(PRINT_ALGORITHM_STEPS) ) {
00124         out << omsg.str();
00125       }
00126       throw std::runtime_error( std::string("DampenCrossTermStd_Step::do_step(...) : ")
00127                     + omsg.str() );
00128     }
00129   }
00130 
00131   return true;
00132 }
00133 
00134 void MoochoPack::DampenCrossTermStd_Step::print_step( const Algorithm& algo
00135   , poss_type step_poss, IterationPack::EDoStepType type, poss_type assoc_step_poss
00136   , std::ostream& out, const std::string& L ) const
00137 {
00138   out
00139     << L << "*** Compute the dampening parameter for the reduced QP cross term w_k\n"
00140     << L << "default: frac_descent = " << frac_descent() << std::endl
00141     << L << "if w_k is update then\n"
00142     << L << "  find zeta_k s.t.\n"
00143     << L << "    Gf_k'*Z_k*pz_k approx\n"
00144     << L << "       - zeta_k * rGf_k'*inv(rHL_k)*w_k - rGf_k'*inv(rHL_k)*rGf_k\n"
00145     << L << "       <= - frac_descent * rGf_k'*inv(rHL_k)*rGf_k\n"
00146     << L << "end\n";
00147 }
00148 
00149 #endif // 0
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends