MOOCHO (Single Doxygen Collection) Version of the Day
MoochoPack_ReducedHessianExactStd_Step.cpp
Go to the documentation of this file.
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 <sstream>
00045 #include <typeinfo>
00046 #include <iomanip>
00047 
00048 #include "MoochoPack_ReducedHessianExactStd_Step.hpp"
00049 #include "MoochoPack_moocho_algo_conversion.hpp"
00050 #include "AbstractLinAlgPack_MatrixSymDenseInitialize.hpp"
00051 #include "IterationPack_print_algorithm_step.hpp"
00052 #include "ConstrainedOptPack/src/VectorWithNorms.h"
00053 #include "NLPInterfacePack_NLPSecondOrder.hpp"
00054 #include "AbstractLinAlgPack/src/AbstractLinAlgPack_MatrixSymOp.hpp"
00055 #include "DenseLinAlgPack_LinAlgOpPack.hpp"
00056 #include "DenseLinAlgPack_DMatrixAsTriSym.hpp"
00057 #include "DenseLinAlgPack_DMatrixOut.hpp"
00058 #include "DenseLinAlgPack_DVectorClass.hpp"
00059 #include "DenseLinAlgPack_DVectorOp.hpp"
00060 #include "DenseLinAlgPack_DVectorOut.hpp"
00061 #include "Midynamic_cast_verbose.h"
00062 
00063 namespace MoochoPack {
00064 
00065 bool ReducedHessianExactStd_Step::do_step(
00066     Algorithm& _algo, poss_type step_poss, IterationPack::EDoStepType type
00067   , poss_type assoc_step_poss)
00068 {
00069   using Teuchos::dyn_cast;
00070   using DenseLinAlgPack::nonconst_sym;
00071   using AbstractLinAlgPack::Mp_StMtMtM;
00072   typedef AbstractLinAlgPack::MatrixSymDenseInitialize  MatrixSymDenseInitialize;
00073   typedef AbstractLinAlgPack::MatrixSymOp     MatrixSymOp;
00074   using ConstrainedOptPack::NLPSecondOrder;
00075 
00076   NLPAlgo &algo = rsqp_algo(_algo);
00077   NLPAlgoState  &s    = algo.rsqp_state();
00078   NLPSecondOrder
00079 #ifdef _WINDOWS
00080         &nlp  = dynamic_cast<NLPSecondOrder&>(algo.nlp());
00081 #else
00082         &nlp  = dyn_cast<NLPSecondOrder>(algo.nlp());
00083 #endif
00084   MatrixSymOp
00085     *HL_sym_op = dynamic_cast<MatrixSymOp*>(&s.HL().get_k(0));
00086 
00087   EJournalOutputLevel olevel = algo.algo_cntr().journal_output_level();
00088   std::ostream& out = algo.track().journal_out();
00089 
00090   // print step header.
00091   if( static_cast<int>(olevel) >= static_cast<int>(PRINT_ALGORITHM_STEPS) ) {
00092     using IterationPack::print_algorithm_step;
00093     print_algorithm_step( algo, step_poss, type, assoc_step_poss, out );
00094   }
00095 
00096   // problem size
00097   size_type n   = nlp.n(),
00098         r   = nlp.r(),
00099         nind  = n - r;
00100 
00101   // Compute HL first (You may want to move this into its own step later?)
00102 
00103   if( !s.lambda().updated_k(-1) ) {
00104     if( static_cast<int>(olevel) >= static_cast<int>(PRINT_ALGORITHM_STEPS) ) {
00105       out << "Initializing lambda_km1 = nlp.get_lambda_init ... \n";
00106     }
00107     nlp.get_init_lagrange_mult( &s.lambda().set_k(-1).v(), NULL );
00108     if( static_cast<int>(olevel) >= static_cast<int>(PRINT_ALGORITHM_STEPS) ) {
00109       out << "||lambda_km1||inf = " << s.lambda().get_k(-1).norm_inf() << std::endl;
00110     }
00111     if( static_cast<int>(olevel) >= static_cast<int>(PRINT_VECTORS) ) {
00112       out << "lambda_km1 = \n" << s.lambda().get_k(-1)();
00113     }
00114   }
00115 
00116   nlp.set_HL( HL_sym_op );
00117   nlp.calc_HL( s.x().get_k(0)(), s.lambda().get_k(-1)(), false );
00118 
00119   if( static_cast<int>(olevel) >= static_cast<int>(PRINT_ITERATION_QUANTITIES) ) {
00120     s.HL().get_k(0).output( out << "\nHL_k = \n" );
00121   }
00122 
00123   // If rHL has already been updated for this iteration then just leave it.
00124   if( !s.rHL().updated_k(0) ) {
00125 
00126     if( !HL_sym_op ) {
00127       std::ostringstream omsg;
00128       omsg
00129         << "ReducedHessianExactStd_Step::do_step(...) : Error, "
00130         << "The matrix HL with the concrete type "
00131         << typeName(s.HL().get_k(0)) << " does not support the "
00132         << "MatrixSymOp iterface";
00133       throw std::logic_error( omsg.str() );
00134     }   
00135 
00136     MatrixSymDenseInitialize
00137       *rHL_sym_init = dynamic_cast<MatrixSymDenseInitialize*>(&s.rHL().set_k(0));
00138     if( !rHL_sym_init ) {
00139       std::ostringstream omsg;
00140       omsg
00141         << "ReducedHessianExactStd_Step::do_step(...) : Error, "
00142         << "The matrix rHL with the concrete type "
00143         << typeName(s.rHL().get_k(0)) << " does not support the "
00144         << "MatrixSymDenseInitialize iterface";
00145       throw std::logic_error( omsg.str() );
00146     }   
00147 
00148     // Compute the dense reduced Hessian
00149     DMatrix rHL_sym_store(nind,nind);
00150     DMatrixSliceSym rHL_sym(rHL_sym_store(),BLAS_Cpp::lower);
00151     Mp_StMtMtM( &rHL_sym, 1.0, MatrixSymOp::DUMMY_ARG, *HL_sym_op
00152           , s.Z().get_k(0), BLAS_Cpp::no_trans, 0.0 );
00153 
00154     if( (int)olevel >= (int)PRINT_ITERATION_QUANTITIES ) {
00155       out << "\nLower triangular partion of dense reduced Hessian (ignore nonzeros above diagonal):\nrHL_dense = \n" << rHL_sym_store(); 
00156     }
00157   
00158     // Set the reduced Hessain
00159     rHL_sym_init->initialize( rHL_sym );
00160 
00161     if( (int)olevel >= (int)PRINT_ITERATION_QUANTITIES ) {
00162       s.rHL().get_k(0).output( out << "\nrHL_k = \n" );
00163     }
00164   }
00165 
00166   return true;
00167 }
00168 
00169 void ReducedHessianExactStd_Step::print_step(
00170     const Algorithm& algo, poss_type step_poss, IterationPack::EDoStepType type
00171   , poss_type assoc_step_poss, std::ostream& out, const std::string& L ) const
00172 {
00173   out
00174     << L << "*** Calculate the exact reduced Hessian of the Lagrangian\n"
00175     << L << "if lambda_km1 is not updated then\n"
00176     << L << "    lambda_km1 = nlp.get_lambda_init\n"
00177     << L << "end\n"
00178     << L << "HL_k = HL(x_k,lambda_km1) <: R^(n+m) -> R^(n x n)\n"
00179     << L << "if rHL_k is not updated then\n"
00180     << L << "    rHL_dense = Z_k' * HL_k * Z_k  (MatrixSymOp interface for HL_k)\n"
00181     << L << "    rHL_k = rHL_dense (MatrixSymDenseInitialize interface for rHL_k)\n"
00182     << L << "end\n";
00183 }
00184 
00185 } // end namespace MoochoPack 
00186 
00187 #endif // 0
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Defines