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

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