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 namespace LinAlgOpPack {
00048   using AbstractLinAlgPack::Vp_StMtV;
00049 }
00050 
00051 namespace MoochoPack {
00052 
00053 bool TangentialStepWithoutBounds_Step::do_step(
00054   Algorithm& _algo, poss_type step_poss, IterationPack::EDoStepType type
00055   ,poss_type assoc_step_poss
00056   )
00057 {
00058   using BLAS_Cpp::no_trans;
00059   using Teuchos::dyn_cast;
00060   using AbstractLinAlgPack::assert_print_nan_inf;
00061   using AbstractLinAlgPack::Vt_S;
00062   using AbstractLinAlgPack::Vp_StV;
00063   using AbstractLinAlgPack::V_InvMtV;
00064   using LinAlgOpPack::V_StV;
00065   using LinAlgOpPack::V_MtV;
00066 
00067   NLPAlgo &algo = rsqp_algo(_algo);
00068   NLPAlgoState  &s    = algo.rsqp_state();
00069 
00070   EJournalOutputLevel olevel = algo.algo_cntr().journal_output_level();
00071   EJournalOutputLevel ns_olevel = algo.algo_cntr().null_space_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( algo.nlp().num_bounded_x() )
00081     TEST_FOR_EXCEPTION(
00082       true, std::logic_error
00083       ,"TangentialStepWithoutBounds_Step::do_step(...): Error, "
00084       "can't solve for pz for NLP with undecomposed constraints or "
00085       "has bounds on the variables");
00086 
00087   // Comupte qp_grad which is an approximation to rGf + Z' * HL * Y * py
00088 
00089   // qp_grad_k = rGf_k
00090   VectorMutable &qp_grad_k = s.qp_grad().set_k(0) = s.rGf().get_k(0);
00091 
00092   IterQuantityAccess<value_type>               &zeta_iq       = s.zeta();
00093   IterQuantityAccess<VectorMutable>      &w_iq          = s.w();
00094   if( w_iq.updated_k(0) && zeta_iq.updated_k(0) ) {
00095     // qp_grad += zeta * w
00096     Vp_StV( &qp_grad_k, zeta_iq.get_k(0), w_iq.get_k(0) );
00097   }
00098 
00099   // Solve the system pz = - inv(rHL) * qp_grad
00100   VectorMutable               &pz_k  = s.pz().set_k(0);
00101   const MatrixSymOpNonsing  &rHL_k = dyn_cast<MatrixSymOpNonsing>(s.rHL().get_k(0));
00102   V_InvMtV( &pz_k, rHL_k, no_trans, qp_grad_k );
00103   Vt_S( &pz_k, -1.0 );
00104 
00105   // nu = 0.0
00106   s.nu().set_k(0) = 0.0;
00107 
00108   // Zpz = Z * pz
00109   V_MtV( &s.Zpz().set_k(0), s.Z().get_k(0), no_trans, pz_k );
00110 
00111   if( static_cast<int>(olevel) >= static_cast<int>(PRINT_ALGORITHM_STEPS) ) {
00112     out << "\n||pz_k||inf   = " << s.pz().get_k(0).norm_inf()
00113       << "\n||Zpz_k||2    = " << s.Zpz().get_k(0).norm_2()  << std::endl;
00114   }
00115 
00116   if( static_cast<int>(ns_olevel) >= static_cast<int>(PRINT_VECTORS) ) {
00117     out << "\npz_k = \n" << s.pz().get_k(0);
00118     out << std::endl;
00119   }
00120 
00121   if( static_cast<int>(olevel) >= static_cast<int>(PRINT_VECTORS) ) {
00122     out << "\nnu_k = \n" << s.nu().get_k(0);
00123     out << "\nZpz_k = \n" << s.Zpz().get_k(0);
00124     out << std::endl;
00125   }
00126 
00127   if(algo.algo_cntr().check_results()) {
00128     assert_print_nan_inf(s.pz().get_k(0),  "pz_k",true,&out);
00129     assert_print_nan_inf(s.Zpz().get_k(0), "Zpz_k",true,&out);
00130   }
00131 
00132   return true;
00133 }
00134 
00135 void TangentialStepWithoutBounds_Step::print_step( const Algorithm& algo
00136   , poss_type step_poss, IterationPack::EDoStepType type, poss_type assoc_step_poss
00137   , std::ostream& out, const std::string& L ) const
00138 {
00139   out
00140     << L << "*** Calculate the null space step by solving an unconstrainted QP\n"
00141     << L << "qp_grad_k = rGf_k + zeta_k * w_k\n"
00142     << L << "solve:\n"
00143     << L << "    min     qp_grad_k' * pz_k + 1/2 * pz_k' * rHL_k * pz_k\n"
00144     << L << "    pz_k <: R^(n-r)\n"
00145     << L << "Zpz_k = Z_k * pz_k\n"
00146     << L << "nu_k = 0\n";
00147 }
00148 
00149 } // end namespace MoochoPack

Generated on Thu Sep 18 12:35:18 2008 for MOOCHO (Single Doxygen Collection) by doxygen 1.3.9.1