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

Generated on Tue Jul 13 09:30:53 2010 for MOOCHO (Single Doxygen Collection) by  doxygen 1.4.7