MoochoPack_ReducedHessianSerialization_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 // 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 <fstream>
00030 
00031 #include "MoochoPack_ReducedHessianSerialization_Step.hpp"
00032 #include "MoochoPack_moocho_algo_conversion.hpp"
00033 #include "MoochoPack_Exceptions.hpp"
00034 #include "IterationPack_print_algorithm_step.hpp"
00035 #include "AbstractLinAlgPack_VectorSpace.hpp"
00036 #include "AbstractLinAlgPack_VectorStdOps.hpp"
00037 #include "AbstractLinAlgPack_VectorMutable.hpp"
00038 #include "AbstractLinAlgPack_VectorOut.hpp"
00039 #include "AbstractLinAlgPack_MatrixSymOpNonsing.hpp"
00040 #include "AbstractLinAlgPack_MatrixSymInitDiag.hpp"
00041 #include "AbstractLinAlgPack_MatrixOpOut.hpp"
00042 #include "AbstractLinAlgPack_LinAlgOpPack.hpp"
00043 #include "SerializationPack_Serializable.hpp"
00044 #include "Teuchos_dyn_cast.hpp"
00045 
00046 namespace MoochoPack {
00047 
00048 ReducedHessianSerialization_Step::ReducedHessianSerialization_Step(
00049   const std::string    &reduced_hessian_input_file_name
00050   ,const std::string   &reduced_hessian_output_file_name
00051   )
00052   :reduced_hessian_input_file_name_(reduced_hessian_input_file_name)
00053   ,reduced_hessian_output_file_name_(reduced_hessian_output_file_name)
00054 {}
00055 
00056 // Overridden from AlgorithmStep
00057 
00058 bool ReducedHessianSerialization_Step::do_step(
00059   Algorithm& _algo, poss_type step_poss, IterationPack::EDoStepType type
00060   ,poss_type assoc_step_poss
00061   )
00062 {
00063   using Teuchos::dyn_cast;
00064   using SerializationPack::Serializable;
00065 
00066   NLPAlgo       &algo = rsqp_algo(_algo);
00067   NLPAlgoState  &s    = algo.rsqp_state();
00068   
00069   EJournalOutputLevel olevel = algo.algo_cntr().journal_output_level();
00070   EJournalOutputLevel ns_olevel = algo.algo_cntr().null_space_journal_output_level();
00071   std::ostream& out = algo.track().journal_out();
00072 
00073   // print step header.
00074   if( static_cast<int>(olevel) >= static_cast<int>(PRINT_ALGORITHM_STEPS) ) {
00075     using IterationPack::print_algorithm_step;
00076     print_algorithm_step( _algo, step_poss, type, assoc_step_poss, out );
00077   }
00078 
00079   IterQuantityAccess<MatrixSymOp>  &rHL_iq = s.rHL();
00080 
00081   if( !rHL_iq.updated_k(0) && reduced_hessian_input_file_name().length() ) {
00082     int k_last_offset = rHL_iq.last_updated();
00083     if( k_last_offset == IterQuantity::NONE_UPDATED ) {
00084       if( (int)olevel >= (int)PRINT_ALGORITHM_STEPS ) {
00085         out
00086           << "\nNo previous matrix rHL was found!\n"
00087           << "\nReading in the matrix rHL_k from the file \""<<reduced_hessian_input_file_name()<<"\" ...\n";
00088       }
00089       MatrixSymOp &rHL_k = rHL_iq.set_k(0);
00090       Serializable &rHL_serializable = dyn_cast<Serializable>(rHL_k);
00091       std::ifstream reduced_hessian_input_file(reduced_hessian_input_file_name().c_str());
00092       TEST_FOR_EXCEPTION(
00093         !reduced_hessian_input_file, std::logic_error
00094         ,"ReducedHessianSerialization_Step::do_step(...): Error, the file \""<<reduced_hessian_input_file_name()<<"\""
00095         " could not be opened or contains no input!"
00096         );
00097       rHL_serializable.unserialize(reduced_hessian_input_file);
00098       if( (int)ns_olevel >= (int)PRINT_ALGORITHM_STEPS ) {
00099         out << "\nrHL_k.rows() = " << rHL_k.rows() << std::endl;
00100         out << "\nrHL_k.cols() = " << rHL_k.cols() << std::endl;
00101         if(algo.algo_cntr().calc_matrix_norms())
00102           out << "\n||rHL_k||inf    = " << rHL_k.calc_norm(MatrixOp::MAT_NORM_INF).value << std::endl;
00103         if(algo.algo_cntr().calc_conditioning()) {
00104           const MatrixSymOpNonsing
00105             *rHL_ns_k = dynamic_cast<const MatrixSymOpNonsing*>(&rHL_k);
00106           if(rHL_ns_k)
00107             out << "\ncond_inf(rHL_k) = " << rHL_ns_k->calc_cond_num(MatrixOp::MAT_NORM_INF).value << std::endl;
00108         }
00109       }
00110       if( (int)ns_olevel >= (int)PRINT_ITERATION_QUANTITIES )
00111         out << "\nrHL_k = \n" << rHL_k;
00112       // Validate the space
00113       const MatrixOp &Z_k = s.Z().get_k(0);
00114       const VectorSpace &null_space = Z_k.space_rows();
00115       TEST_FOR_EXCEPTION(
00116         !null_space.is_compatible(rHL_k.space_cols()) || !null_space.is_compatible(rHL_k.space_rows())
00117         ,std::runtime_error
00118         ,"ReducedHessianSerialization_Step::do_step(...): Error, the read-in reduced Hessian of dimension "
00119         << rHL_k.rows() << " x " << rHL_k.cols() << " is not compatible with the null space of dimension "
00120         "Z_k.cols() = " << Z_k.cols() << "!"
00121         );
00122     }
00123   }
00124   
00125   return true;
00126   
00127 }
00128 
00129 void ReducedHessianSerialization_Step::finalize_step(
00130   Algorithm& _algo, poss_type step_poss, IterationPack::EDoStepType type
00131   ,poss_type assoc_step_poss
00132   )
00133 {
00134   using Teuchos::dyn_cast;
00135   using SerializationPack::Serializable;
00136 
00137   const NLPAlgo       &algo = rsqp_algo(_algo);
00138   const NLPAlgoState  &s    = algo.rsqp_state();
00139 
00140   EJournalOutputLevel olevel = algo.algo_cntr().journal_output_level();
00141   std::ostream& out = algo.track().journal_out();
00142 
00143   const IterQuantityAccess<MatrixSymOp>  &rHL_iq = s.rHL();
00144   
00145   int k_last_offset = rHL_iq.last_updated();
00146 
00147   if( k_last_offset != IterQuantity::NONE_UPDATED && reduced_hessian_output_file_name().length() ) {
00148     if( (int)olevel >= (int)PRINT_ALGORITHM_STEPS ) {
00149       out
00150         << "\nSerializing the matrix rHL_k("<<k_last_offset<<") to the file "
00151         << "\""<<reduced_hessian_output_file_name()<<"\" ...\n";
00152     }
00153     const MatrixSymOp &rHL_k = rHL_iq.get_k(k_last_offset);
00154     const Serializable &rHL_serializable = dyn_cast<const Serializable>(rHL_k);
00155     std::ofstream reduced_hessian_output_file(reduced_hessian_output_file_name().c_str());
00156     TEST_FOR_EXCEPTION(
00157       !reduced_hessian_output_file, std::logic_error
00158       ,"ReducedHessianSerialization_Step::finalize_step(...): Error, the file \""<<reduced_hessian_output_file_name()<<"\""
00159       " could not be opened!"
00160       );
00161     rHL_serializable.serialize(reduced_hessian_output_file);
00162   }
00163 }
00164 
00165 void ReducedHessianSerialization_Step::print_step(
00166   const Algorithm& algo, poss_type step_poss, IterationPack::EDoStepType type, poss_type assoc_step_poss
00167   ,std::ostream& out, const std::string& L
00168   ) const
00169 {
00170   out
00171     << L << "*** Read in the reduced Hessian of the Lagrangian rHL from a file.\n"
00172     << L << "if (rHL_k is not updated and reduced_hessian_input_file_name != \"\") then\n"
00173     << L << "  k_last_offset = last iteration rHL was updated for\n"
00174     << L << "  if k_last_offset==NONE_UPDATED then\n"
00175     << L << "    rHL_serializable = dyn_cast<Serializable>(rHL_k) *** Throws exception if fails!\n"
00176     << L << "    Unserialize into rHL_serializable from the file \""<<reduced_hessian_input_file_name()<<"\"\n"
00177     << L << "  else\n"
00178     << L << "    *** There is some reduced Hessian that exists for some past iteration so\n"
00179     << L << "    *** we will let some other step object initialize itQ\n"
00180     << L << "  end\n"
00181     << L << "end\n"
00182     << L << "*** Note: On finalization, this step object will serialize rHL_k to the file:\n"
00183     << L << "***   \""<<reduced_hessian_output_file_name()<<"\""
00184     ;
00185 }
00186 
00187 } // end namespace MoochoPack 

Generated on Wed May 12 21:32:13 2010 for MoochoPack : Framework for Large-Scale Optimization Algorithms by  doxygen 1.4.7