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 // This library is free software; you can redistribute it and/or modify
00013 // it under the terms of the GNU Lesser General Public License as
00014 // published by the Free Software Foundation; either version 2.1 of the
00015 // License, or (at your option) any later version.
00016 //  
00017 // This library is distributed in the hope that it will be useful, but
00018 // WITHOUT ANY WARRANTY; without even the implied warranty of
00019 // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
00020 // Lesser General Public License for more details.
00021 //  
00022 // You should have received a copy of the GNU Lesser General Public
00023 // License along with this library; if not, write to the Free Software
00024 // Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307
00025 // USA
00026 // Questions? Contact Roscoe A. Bartlett (rabartl@sandia.gov) 
00027 // 
00028 // ***********************************************************************
00029 // @HEADER
00030 
00031 #include <sstream>
00032 #include <typeinfo>
00033 #include <iomanip>
00034 
00035 #include "MoochoPack_ReducedHessianExactStd_Step.hpp"
00036 #include "MoochoPack_moocho_algo_conversion.hpp"
00037 #include "AbstractLinAlgPack_MatrixSymDenseInitialize.hpp"
00038 #include "IterationPack_print_algorithm_step.hpp"
00039 #include "ConstrainedOptPack/src/VectorWithNorms.h"
00040 #include "NLPInterfacePack_NLPSecondOrder.hpp"
00041 #include "AbstractLinAlgPack/src/AbstractLinAlgPack_MatrixSymOp.hpp"
00042 #include "DenseLinAlgPack_LinAlgOpPack.hpp"
00043 #include "DenseLinAlgPack_DMatrixAsTriSym.hpp"
00044 #include "DenseLinAlgPack_DMatrixOut.hpp"
00045 #include "DenseLinAlgPack_DVectorClass.hpp"
00046 #include "DenseLinAlgPack_DVectorOp.hpp"
00047 #include "DenseLinAlgPack_DVectorOut.hpp"
00048 #include "Midynamic_cast_verbose.h"
00049 
00050 namespace MoochoPack {
00051 
00052 bool ReducedHessianExactStd_Step::do_step(
00053     Algorithm& _algo, poss_type step_poss, IterationPack::EDoStepType type
00054   , poss_type assoc_step_poss)
00055 {
00056   using Teuchos::dyn_cast;
00057   using DenseLinAlgPack::nonconst_sym;
00058   using AbstractLinAlgPack::Mp_StMtMtM;
00059   typedef AbstractLinAlgPack::MatrixSymDenseInitialize  MatrixSymDenseInitialize;
00060   typedef AbstractLinAlgPack::MatrixSymOp     MatrixSymOp;
00061   using ConstrainedOptPack::NLPSecondOrder;
00062 
00063   NLPAlgo &algo = rsqp_algo(_algo);
00064   NLPAlgoState  &s    = algo.rsqp_state();
00065   NLPSecondOrder
00066 #ifdef _WINDOWS
00067         &nlp  = dynamic_cast<NLPSecondOrder&>(algo.nlp());
00068 #else
00069         &nlp  = dyn_cast<NLPSecondOrder>(algo.nlp());
00070 #endif
00071   MatrixSymOp
00072     *HL_sym_op = dynamic_cast<MatrixSymOp*>(&s.HL().get_k(0));
00073 
00074   EJournalOutputLevel olevel = algo.algo_cntr().journal_output_level();
00075   std::ostream& out = algo.track().journal_out();
00076 
00077   // print step header.
00078   if( static_cast<int>(olevel) >= static_cast<int>(PRINT_ALGORITHM_STEPS) ) {
00079     using IterationPack::print_algorithm_step;
00080     print_algorithm_step( algo, step_poss, type, assoc_step_poss, out );
00081   }
00082 
00083   // problem size
00084   size_type n   = nlp.n(),
00085         r   = nlp.r(),
00086         nind  = n - r;
00087 
00088   // Compute HL first (You may want to move this into its own step later?)
00089 
00090   if( !s.lambda().updated_k(-1) ) {
00091     if( static_cast<int>(olevel) >= static_cast<int>(PRINT_ALGORITHM_STEPS) ) {
00092       out << "Initializing lambda_km1 = nlp.get_lambda_init ... \n";
00093     }
00094     nlp.get_init_lagrange_mult( &s.lambda().set_k(-1).v(), NULL );
00095     if( static_cast<int>(olevel) >= static_cast<int>(PRINT_ALGORITHM_STEPS) ) {
00096       out << "||lambda_km1||inf = " << s.lambda().get_k(-1).norm_inf() << std::endl;
00097     }
00098     if( static_cast<int>(olevel) >= static_cast<int>(PRINT_VECTORS) ) {
00099       out << "lambda_km1 = \n" << s.lambda().get_k(-1)();
00100     }
00101   }
00102 
00103   nlp.set_HL( HL_sym_op );
00104   nlp.calc_HL( s.x().get_k(0)(), s.lambda().get_k(-1)(), false );
00105 
00106   if( static_cast<int>(olevel) >= static_cast<int>(PRINT_ITERATION_QUANTITIES) ) {
00107     s.HL().get_k(0).output( out << "\nHL_k = \n" );
00108   }
00109 
00110   // If rHL has already been updated for this iteration then just leave it.
00111   if( !s.rHL().updated_k(0) ) {
00112 
00113     if( !HL_sym_op ) {
00114       std::ostringstream omsg;
00115       omsg
00116         << "ReducedHessianExactStd_Step::do_step(...) : Error, "
00117         << "The matrix HL with the concrete type "
00118         << typeName(s.HL().get_k(0)) << " does not support the "
00119         << "MatrixSymOp iterface";
00120       throw std::logic_error( omsg.str() );
00121     }   
00122 
00123     MatrixSymDenseInitialize
00124       *rHL_sym_init = dynamic_cast<MatrixSymDenseInitialize*>(&s.rHL().set_k(0));
00125     if( !rHL_sym_init ) {
00126       std::ostringstream omsg;
00127       omsg
00128         << "ReducedHessianExactStd_Step::do_step(...) : Error, "
00129         << "The matrix rHL with the concrete type "
00130         << typeName(s.rHL().get_k(0)) << " does not support the "
00131         << "MatrixSymDenseInitialize iterface";
00132       throw std::logic_error( omsg.str() );
00133     }   
00134 
00135     // Compute the dense reduced Hessian
00136     DMatrix rHL_sym_store(nind,nind);
00137     DMatrixSliceSym rHL_sym(rHL_sym_store(),BLAS_Cpp::lower);
00138     Mp_StMtMtM( &rHL_sym, 1.0, MatrixSymOp::DUMMY_ARG, *HL_sym_op
00139           , s.Z().get_k(0), BLAS_Cpp::no_trans, 0.0 );
00140 
00141     if( (int)olevel >= (int)PRINT_ITERATION_QUANTITIES ) {
00142       out << "\nLower triangular partion of dense reduced Hessian (ignore nonzeros above diagonal):\nrHL_dense = \n" << rHL_sym_store(); 
00143     }
00144   
00145     // Set the reduced Hessain
00146     rHL_sym_init->initialize( rHL_sym );
00147 
00148     if( (int)olevel >= (int)PRINT_ITERATION_QUANTITIES ) {
00149       s.rHL().get_k(0).output( out << "\nrHL_k = \n" );
00150     }
00151   }
00152 
00153   return true;
00154 }
00155 
00156 void ReducedHessianExactStd_Step::print_step(
00157     const Algorithm& algo, poss_type step_poss, IterationPack::EDoStepType type
00158   , poss_type assoc_step_poss, std::ostream& out, const std::string& L ) const
00159 {
00160   out
00161     << L << "*** Calculate the exact reduced Hessian of the Lagrangian\n"
00162     << L << "if lambda_km1 is not updated then\n"
00163     << L << "    lambda_km1 = nlp.get_lambda_init\n"
00164     << L << "end\n"
00165     << L << "HL_k = HL(x_k,lambda_km1) <: R^(n+m) -> R^(n x n)\n"
00166     << L << "if rHL_k is not updated then\n"
00167     << L << "    rHL_dense = Z_k' * HL_k * Z_k  (MatrixSymOp interface for HL_k)\n"
00168     << L << "    rHL_k = rHL_dense (MatrixSymDenseInitialize interface for rHL_k)\n"
00169     << L << "end\n";
00170 }
00171 
00172 } // end namespace MoochoPack 
00173 
00174 #endif // 0

Generated on Wed May 12 21:52:31 2010 for MOOCHO (Single Doxygen Collection) by  doxygen 1.4.7