MoochoPack_TangentialStepWithoutBounds_Step.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 
00031 #include "MoochoPack_TangentialStepWithoutBounds_Step.hpp"
00032 #include "MoochoPack_EvalNewPointTailoredApproach_Step.hpp"
00033 #include "MoochoPack_Exceptions.hpp"
00034 #include "MoochoPack_moocho_algo_conversion.hpp"
00035 #include "IterationPack_print_algorithm_step.hpp"
00036 #include "NLPInterfacePack_NLPDirect.hpp"
00037 #include "AbstractLinAlgPack_MatrixSymOpNonsing.hpp"
00038 #include "AbstractLinAlgPack_MatrixOpOut.hpp"
00039 #include "AbstractLinAlgPack_VectorMutable.hpp"
00040 #include "AbstractLinAlgPack_VectorStdOps.hpp"
00041 #include "AbstractLinAlgPack_VectorOut.hpp"
00042 #include "AbstractLinAlgPack_assert_print_nan_inf.hpp"
00043 #include "AbstractLinAlgPack_LinAlgOpPack.hpp"
00044 #include "Teuchos_dyn_cast.hpp"
00045 #include "Teuchos_TestForException.hpp"
00046 
00047 
00048 namespace LinAlgOpPack {
00049   using AbstractLinAlgPack::Vp_StMtV;
00050 }
00051 
00052 
00053 namespace MoochoPack {
00054 
00055 
00056 TangentialStepWithoutBounds_Step::TangentialStepWithoutBounds_Step()
00057   :max_pz_norm_(-1.0),
00058    num_pz_damp_iters_(0)
00059 {}
00060 
00061 
00062 bool TangentialStepWithoutBounds_Step::do_step(
00063   Algorithm& _algo, poss_type step_poss, IterationPack::EDoStepType type
00064   ,poss_type assoc_step_poss
00065   )
00066 {
00067 
00068   using std::endl;
00069   using BLAS_Cpp::no_trans;
00070   using Teuchos::dyn_cast;
00071   using AbstractLinAlgPack::assert_print_nan_inf;
00072   using AbstractLinAlgPack::Vt_S;
00073   using AbstractLinAlgPack::Vp_StV;
00074   using AbstractLinAlgPack::V_InvMtV;
00075   using LinAlgOpPack::V_StV;
00076   using LinAlgOpPack::V_MtV;
00077 
00078   NLPAlgo &algo = rsqp_algo(_algo);
00079   NLPAlgoState &s = algo.rsqp_state();
00080   
00081   EJournalOutputLevel olevel = algo.algo_cntr().journal_output_level();
00082   EJournalOutputLevel ns_olevel = algo.algo_cntr().null_space_journal_output_level();
00083   std::ostream& out = algo.track().journal_out();
00084 
00085   // print step header.
00086   if( static_cast<int>(olevel) >= static_cast<int>(PRINT_ALGORITHM_STEPS) ) {
00087     using IterationPack::print_algorithm_step;
00088     print_algorithm_step( algo, step_poss, type, assoc_step_poss, out );
00089   }
00090 
00091   if( algo.nlp().num_bounded_x() )
00092     TEST_FOR_EXCEPTION(
00093       true, std::logic_error
00094       ,"TangentialStepWithoutBounds_Step::do_step(...): Error, "
00095       "can't solve for pz for NLP with undecomposed constraints or "
00096       "has bounds on the variables");
00097 
00098   // Comupte qp_grad which is an approximation to rGf + Z' * HL * Y * py
00099 
00100   // qp_grad_k = rGf_k
00101   VectorMutable &qp_grad_k = s.qp_grad().set_k(0) = s.rGf().get_k(0);
00102 
00103   IterQuantityAccess<value_type> &zeta_iq = s.zeta();
00104   IterQuantityAccess<VectorMutable> &w_iq = s.w();
00105   if( w_iq.updated_k(0) && zeta_iq.updated_k(0) ) {
00106     // qp_grad += zeta * w
00107     Vp_StV( &qp_grad_k, zeta_iq.get_k(0), w_iq.get_k(0) );
00108   }
00109 
00110   // Solve the system pz = - inv(rHL) * qp_grad
00111   VectorMutable &pz_k = s.pz().set_k(0);
00112   const MatrixSymOpNonsing &rHL_k = dyn_cast<MatrixSymOpNonsing>(s.rHL().get_k(0));
00113   V_InvMtV( &pz_k, rHL_k, no_trans, qp_grad_k );
00114   Vt_S( &pz_k, -1.0 );
00115 
00116   // nu = 0.0
00117   s.nu().set_k(0) = 0.0;
00118 
00119   value_type pz_norm_inf = -1.0;
00120   if( static_cast<int>(olevel) >= static_cast<int>(PRINT_ALGORITHM_STEPS) ) {
00121     pz_norm_inf = pz_k.norm_inf();
00122     out << "\n||pz_k||inf   = " << pz_norm_inf << endl;
00123   }
00124   if( static_cast<int>(ns_olevel) >= static_cast<int>(PRINT_VECTORS) )
00125     out << "\npz_k = \n" << pz_k << std::endl;
00126 
00127   // Check to see if we need to dampen pz
00128   const bool dampen_pz = max_pz_norm() >= 0.0 && s.k() <= num_pz_damp_iters();
00129   if( static_cast<int>(olevel) >= static_cast<int>(PRINT_ALGORITHM_STEPS) ) {
00130     out
00131       << "\nChecking if we need to dampen pz:"
00132       << "  ( (max_pz_norm="<<max_pz_norm()<<") >= 0.0 )"
00133       << " && ( (k="<<s.k()<<") <= (num_pz_damp_iters="<<num_pz_damp_iters()<<") ) : "
00134       << ( dampen_pz ? "true, dampen pz ..." : "false, no dampen pz!" )
00135       << endl;
00136   }
00137 
00138   if ( dampen_pz ) {
00139     // pz_new = ( max_pz_norm / pz_norm_inf ) * pz
00140     if (pz_norm_inf < 0.0 )
00141       pz_norm_inf = pz_k.norm_inf();
00142     const value_type pz_damp_factor = ( max_pz_norm() / pz_norm_inf );
00143     if( static_cast<int>(olevel) >= static_cast<int>(PRINT_ALGORITHM_STEPS) ) {
00144       out
00145         << "\npz_damp_factor = max_pz_norm / ||pz||inf = "
00146         << max_pz_norm() << " / " << pz_norm_inf << " = "
00147         << pz_damp_factor;
00148     }
00149     Vt_S( &pz_k, pz_damp_factor );
00150   }
00151 
00152   // Zpz = Z * pz
00153   V_MtV( &s.Zpz().set_k(0), s.Z().get_k(0), no_trans, pz_k );
00154 
00155   if( static_cast<int>(olevel) >= static_cast<int>(PRINT_ALGORITHM_STEPS) ) {
00156     out << "\n||pz_k||inf   = " << s.pz().get_k(0).norm_inf()
00157         << "\n||Zpz_k||2    = " << s.Zpz().get_k(0).norm_2()  << std::endl;
00158   }
00159 
00160   if( static_cast<int>(ns_olevel) >= static_cast<int>(PRINT_VECTORS) ) {
00161     out << "\npz_k = \n" << s.pz().get_k(0);
00162     out << std::endl;
00163   }
00164 
00165   if( static_cast<int>(olevel) >= static_cast<int>(PRINT_VECTORS) ) {
00166     out << "\nnu_k = \n" << s.nu().get_k(0);
00167     out << "\nZpz_k = \n" << s.Zpz().get_k(0);
00168     out << std::endl;
00169   }
00170 
00171   if(algo.algo_cntr().check_results()) {
00172     assert_print_nan_inf(s.pz().get_k(0),  "pz_k",true,&out);
00173     assert_print_nan_inf(s.Zpz().get_k(0), "Zpz_k",true,&out);
00174   }
00175 
00176   return true;
00177 }
00178 
00179 
00180 void TangentialStepWithoutBounds_Step::print_step( const Algorithm& algo
00181   , poss_type step_poss, IterationPack::EDoStepType type, poss_type assoc_step_poss
00182   , std::ostream& out, const std::string& L ) const
00183 {
00184   out
00185     << L << "*** Calculate the null space step by solving an unconstrainted QP\n"
00186     << L << "qp_grad_k = rGf_k + zeta_k * w_k\n"
00187     << L << "solve:\n"
00188     << L << "    min     qp_grad_k' * pz_k + 1/2 * pz_k' * rHL_k * pz_k\n"
00189     << L << "    pz_k <: R^(n-r)\n"
00190     << L << "Zpz_k = Z_k * pz_k\n"
00191     << L << "nu_k = 0\n";
00192 }
00193 
00194 
00195 } // end namespace MoochoPack

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