MoochoPack : Framework for Large-Scale Optimization Algorithms Version of the Day
MoochoPack_TangentialStepIP_Step.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 <iostream>
00044 
00045 #include "MoochoPack_TangentialStepIP_Step.hpp"
00046 #include "MoochoPack_EvalNewPointTailoredApproach_Step.hpp"
00047 #include "MoochoPack_Exceptions.hpp"
00048 #include "MoochoPack_IpState.hpp"
00049 #include "MoochoPack_moocho_algo_conversion.hpp"
00050 #include "IterationPack_print_algorithm_step.hpp"
00051 #include "NLPInterfacePack_NLPDirect.hpp"
00052 #include "AbstractLinAlgPack_MatrixSymDiagStd.hpp"
00053 #include "AbstractLinAlgPack_MatrixSymOpNonsing.hpp"
00054 #include "AbstractLinAlgPack_MatrixOpOut.hpp"
00055 #include "AbstractLinAlgPack_VectorMutable.hpp"
00056 #include "AbstractLinAlgPack_VectorStdOps.hpp"
00057 #include "AbstractLinAlgPack_VectorOut.hpp"
00058 #include "AbstractLinAlgPack_assert_print_nan_inf.hpp"
00059 #include "AbstractLinAlgPack_LinAlgOpPack.hpp"
00060 #include "Teuchos_dyn_cast.hpp"
00061 #include "Teuchos_Assert.hpp"
00062 
00063 namespace MoochoPack {
00064 
00065 bool TangentialStepIP_Step::do_step(
00066   Algorithm& _algo, poss_type step_poss, IterationPack::EDoStepType type
00067   ,poss_type assoc_step_poss
00068   )
00069   {
00070   using BLAS_Cpp::no_trans;
00071   using Teuchos::dyn_cast;
00072   using AbstractLinAlgPack::assert_print_nan_inf;
00073   using LinAlgOpPack::Vt_S;
00074   using LinAlgOpPack::Vp_StV;
00075   using LinAlgOpPack::V_StV;
00076   using LinAlgOpPack::V_MtV;
00077   using LinAlgOpPack::V_InvMtV;
00078    using LinAlgOpPack::M_StM;
00079   using LinAlgOpPack::Mp_StM;
00080   using LinAlgOpPack::assign;
00081 
00082   NLPAlgo &algo = rsqp_algo(_algo);
00083   IpState     &s      = dyn_cast<IpState>(_algo.state());
00084 
00085   EJournalOutputLevel olevel = algo.algo_cntr().journal_output_level();
00086   std::ostream& out = algo.track().journal_out();
00087 
00088   // print step header.
00089   if( static_cast<int>(olevel) >= static_cast<int>(PRINT_ALGORITHM_STEPS) ) {
00090     using IterationPack::print_algorithm_step;
00091     print_algorithm_step( algo, step_poss, type, assoc_step_poss, out );
00092   }
00093 
00094   // Compute qp_grad which is an approximation to rGf + Z'*(mu*(invXu*e-invXl*e) + no_cross_term
00095   // minimize round off error by calc'ing Z'*(Gf + mu*(invXu*e-invXl*e))
00096 
00097   // qp_grad_k = Z'*(Gf + mu*(invXu*e-invXl*e))
00098   const MatrixSymDiagStd  &invXu = s.invXu().get_k(0);
00099   const MatrixSymDiagStd  &invXl = s.invXl().get_k(0);
00100   const value_type            &mu    = s.barrier_parameter().get_k(0);
00101   const MatrixOp          &Z_k   = s.Z().get_k(0);
00102 
00103   Teuchos::RCP<VectorMutable> rhs = s.Gf().get_k(0).clone();
00104   Vp_StV( rhs.get(), mu,      invXu.diag() );
00105   Vp_StV( rhs.get(), -1.0*mu, invXl.diag() );
00106   
00107   if( (int)olevel >= (int)PRINT_ALGORITHM_STEPS ) 
00108     {
00109     out << "\n||Gf_k + mu_k*(invXu_k-invXl_k)||inf = " << rhs->norm_inf() << std::endl;
00110     }
00111   if( (int)olevel >= (int)PRINT_VECTORS)
00112     {
00113     out << "\nGf_k + mu_k*(invXu_k-invXl_k) =\n" << *rhs;
00114     }
00115 
00116   VectorMutable &qp_grad_k = s.qp_grad().set_k(0);
00117   V_MtV(&qp_grad_k, Z_k, BLAS_Cpp::trans, *rhs);
00118   
00119   if( (int)olevel >= (int)PRINT_ALGORITHM_STEPS ) 
00120     {
00121     out << "\n||Z_k'*(Gf_k + mu_k*(invXu_k-invXl_k))||inf = " << qp_grad_k.norm_inf() << std::endl;
00122     }
00123   if( (int)olevel >= (int)PRINT_VECTORS )
00124     {
00125     out << "\nZ_k'*(Gf_k + mu_k*(invXu_k-invXl_k)) =\n" << qp_grad_k;
00126     }
00127 
00128   // error check for cross term
00129   value_type         &zeta    = s.zeta().set_k(0);
00130   const Vector &w_sigma = s.w_sigma().get_k(0);
00131   
00132   // need code to calculate damping parameter
00133   zeta = 1.0;
00134 
00135   Vp_StV(&qp_grad_k, zeta, w_sigma);
00136 
00137   if( (int)olevel >= (int)PRINT_ALGORITHM_STEPS ) 
00138     {
00139     out << "\n||qp_grad_k||inf = " << qp_grad_k.norm_inf() << std::endl;
00140     }
00141   if( (int)olevel >= (int)PRINT_VECTORS ) 
00142     {
00143     out << "\nqp_grad_k =\n" << qp_grad_k;
00144     }
00145 
00146   // build the "Hessian" term B = rHL + rHB
00147   // should this be MatrixSymOpNonsing
00148   const MatrixSymOp      &rHL_k = s.rHL().get_k(0);
00149   const MatrixSymOp      &rHB_k = s.rHB().get_k(0);
00150   MatrixSymOpNonsing &B_k   = dyn_cast<MatrixSymOpNonsing>(s.B().set_k(0));
00151   if (B_k.cols() != Z_k.cols())
00152     {
00153     // Initialize space in rHB
00154     dyn_cast<MatrixSymInitDiag>(B_k).init_identity(Z_k.space_rows(), 0.0);
00155     }
00156 
00157   //  M_StM(&B_k, 1.0, rHL_k, no_trans);
00158   assign(&B_k, rHL_k, BLAS_Cpp::no_trans);
00159   if( (int)olevel >= (int)PRINT_VECTORS ) 
00160     {
00161     out << "\nB_k = rHL_k =\n" << B_k;
00162     }
00163   Mp_StM(&B_k, 1.0, rHB_k, BLAS_Cpp::no_trans);
00164   if( (int)olevel >= (int)PRINT_VECTORS ) 
00165     {
00166     out << "\nB_k = rHL_k + rHB_k =\n" << B_k;
00167     }
00168 
00169   // Solve the system pz = - inv(rHL) * qp_grad
00170   VectorMutable   &pz_k  = s.pz().set_k(0);
00171   V_InvMtV( &pz_k, B_k, no_trans, qp_grad_k );
00172   Vt_S( &pz_k, -1.0 );
00173 
00174   // Zpz = Z * pz
00175   V_MtV( &s.Zpz().set_k(0), s.Z().get_k(0), no_trans, pz_k );
00176 
00177   if( (int)olevel >= (int)PRINT_ALGORITHM_STEPS )
00178     {
00179     out << "\n||pz||inf   = " << s.pz().get_k(0).norm_inf()
00180       << "\nsum(Zpz)    = " << AbstractLinAlgPack::sum(s.Zpz().get_k(0))  << std::endl;
00181     }
00182 
00183   if( (int)olevel >= (int)PRINT_VECTORS )
00184     {
00185     out << "\npz_k = \n" << s.pz().get_k(0);
00186     out << "\nnu_k = \n" << s.nu().get_k(0);
00187     out << "\nZpz_k = \n" << s.Zpz().get_k(0);
00188     out << std::endl;
00189     }
00190 
00191   if(algo.algo_cntr().check_results())
00192     {
00193     assert_print_nan_inf(s.pz().get_k(0),  "pz_k",true,&out);
00194     assert_print_nan_inf(s.Zpz().get_k(0), "Zpz_k",true,&out);
00195     }
00196 
00197   return true;
00198   }
00199 
00200 void TangentialStepIP_Step::print_step( const Algorithm& algo
00201   , poss_type step_poss, IterationPack::EDoStepType type, poss_type assoc_step_poss
00202   , std::ostream& out, const std::string& L ) const
00203   {
00204   out
00205     << L << "*** Calculate the null space step by solving an unconstrainted QP\n"
00206     << L << "zeta_k = 1.0\n"
00207     << L << "qp_grad_k = Z_k'*(Gf_k + mu_k*(invXu_k-invXl_k)) + zeta_k*w_sigma_k\n"
00208     << L << "B_k = rHL_k + rHB_k\n"
00209     << L << "pz_k = -inv(B_k)*qp_grad_k\n"
00210     << L << "Zpz_k = Z_k*pz_k\n"
00211     ;
00212   }
00213 
00214 } // end namespace MoochoPack
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends