MoochoPack_UpdateReducedSigma_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 <ostream>
00030 #include <typeinfo>
00031 #include <iostream>
00032 
00033 #include "AbstractLinAlgPack_MatrixSymDiagStd.hpp"
00034 #include "AbstractLinAlgPack_MatrixSymOpNonsing.hpp"
00035 #include "AbstractLinAlgPack_MatrixOpOut.hpp"
00036 #include "AbstractLinAlgPack_MultiVectorMutable.hpp"
00037 #include "AbstractLinAlgPack_VectorStdOps.hpp"
00038 #include "AbstractLinAlgPack_VectorOut.hpp"
00039 #include "AbstractLinAlgPack_VectorAuxiliaryOps.hpp"
00040 #include "AbstractLinAlgPack_LinAlgOpPack.hpp"
00041 #include "AbstractLinAlgPack_assert_print_nan_inf.hpp"
00042 #include "ConstrainedOptPack_MatrixIdentConcat.hpp"
00043 #include "IterationPack_print_algorithm_step.hpp"
00044 #include "MoochoPack_IpState.hpp"
00045 #include "MoochoPack_UpdateReducedSigma_Step.hpp"
00046 #include "MoochoPack_moocho_algo_conversion.hpp"
00047 
00048 #include "OptionsFromStreamPack_StringToIntMap.hpp"
00049 
00050 #include "Teuchos_dyn_cast.hpp"
00051 
00052 namespace MoochoPack {
00053 
00054 UpdateReducedSigma_Step::UpdateReducedSigma_Step(
00055   const e_update_methods update_method
00056   )
00057   :
00058   update_method_(update_method)
00059   {}
00060 
00061 bool UpdateReducedSigma_Step::do_step(
00062   Algorithm& _algo, poss_type step_poss, IterationPack::EDoStepType type
00063   ,poss_type assoc_step_poss
00064   )
00065   {
00066   using Teuchos::dyn_cast;
00067   using IterationPack::print_algorithm_step;
00068 
00069   NLPAlgo            &algo   = dyn_cast<NLPAlgo>(_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     {
00078     using IterationPack::print_algorithm_step;
00079     print_algorithm_step( _algo, step_poss, type, assoc_step_poss, out );
00080     }
00081 
00082   switch (update_method_)
00083     {
00084     case ALWAYS_EXPLICIT:
00085       {
00086       if( (int)olevel >= (int)PRINT_ALGORITHM_STEPS ) 
00087         {
00088         out << "\nupdate_method is always_explicit, forming Reduced Sigma Explicitly ...\n";
00089         }
00090       FormReducedSigmaExplicitly(algo,s,olevel,out);
00091       }
00092       break;
00093     case BFGS_PRIMAL:
00094     case BFGS_DUAL_NO_CORRECTION:
00095     case BFGS_DUAL_EXPLICIT_CORRECTION:
00096     case BFGS_DUAL_SCALING_CORRECTION:
00097       {
00098       TEST_FOR_EXCEPTION(true, std::logic_error, "Option BFGS_primal not handled yet in UpdateReducedSigma\n");
00099       }
00100       break;
00101     default:
00102       TEST_FOR_EXCEPT(true); // local error ?
00103     };
00104 
00105   if( (int)olevel >= (int)PRINT_ITERATION_QUANTITIES ) 
00106     {
00107     out << "\nrHB_k =\n" << s.rHB().get_k(0);
00108     }
00109 
00110   return true;
00111   }
00112 
00113 
00114 void UpdateReducedSigma_Step::print_step(
00115   const Algorithm& _algo, poss_type step_poss, IterationPack::EDoStepType type
00116   ,poss_type assoc_step_poss, std::ostream& out, const std::string& L
00117   ) const
00118   {
00119   out << L << "*** Update Z'*Sigma*Z\n"
00120     << L << "if (update_method == always_explicit) then\n"
00121     << L << "  Sigma_k = invXl*Vl-invXu*Vu\n"
00122     << L << "  Sigma_I = Sigma_k.sub_view(Z.I_rng)\n"
00123     << L << "  Sigma_D_sqrt = (Sigma_k.sub_view(Z.D_rng))^1/2\n"
00124     << L << "  J = Sigma_D_sqrt*Z\n"
00125     << L << "  rHB_k = J'*J + Sigma_I\n"
00126     << L << "elsif ( update_method == BFGS_??? ) then\n"
00127     << L << "  NOT IMPLEMENTED YET!\n"
00128     << L << "end\n";
00129   }
00130 
00131 void UpdateReducedSigma_Step::FormReducedSigmaExplicitly(
00132   NLPAlgo& algo, IpState& s, EJournalOutputLevel olevel,  std::ostream& out
00133   )
00134   {
00135   namespace mmp = MemMngPack;
00136   using Teuchos::dyn_cast;
00137   using AbstractLinAlgPack::ele_wise_prod;
00138   using AbstractLinAlgPack::ele_wise_sqrt;
00139    using LinAlgOpPack::Mp_M;
00140    using LinAlgOpPack::Mp_MtM;
00141    using LinAlgOpPack::M_MtM;
00142    using LinAlgOpPack::M_StM;
00143   using LinAlgOpPack::V_MtV;
00144   using LinAlgOpPack::assign;
00145 
00146   // Calculate Reduced Sigma directly from
00147   // Sigma = invXl*Vl + invXu*Vu
00148   // Z_kT*Sigma*Z_k
00149 
00150   // Get the iteration quantities
00151   const MatrixIdentConcat     &Z     = dyn_cast<MatrixIdentConcat>(s.Z().get_k(0));
00152   const MatrixSymDiagStd  &invXl = s.invXl().get_k(0);
00153   const MatrixSymDiagStd  &invXu = s.invXu().get_k(0);
00154   const MatrixSymDiagStd  &Vl    = s.Vl().get_k(0);
00155   const MatrixSymDiagStd  &Vu    = s.Vu().get_k(0);
00156   
00157   MatrixSymDiagStd  &Sigma = s.Sigma().set_k(0);
00158 
00159   MatrixSymOpNonsing& rHB = dyn_cast<MatrixSymOpNonsing>(s.rHB().set_k(0));
00160   if (rHB.cols() != Z.cols())
00161     {
00162     // Initialize space in rHB
00163     dyn_cast<MatrixSymInitDiag>(rHB).init_identity(Z.space_rows(), 0.0);
00164     }
00165   
00166   // Calculate Sigma = invXl*Vl + invXu*Vu
00167 
00168   ele_wise_prod(1.0, invXl.diag(), Vl.diag(), &(Sigma.diag() = 0.0));
00169   
00170   if( (int)olevel >= (int)PRINT_ALGORITHM_STEPS )
00171     out << "\n||Sigma_l||inf = " << Sigma.diag().norm_inf() << std::endl;
00172   if( (int)olevel >= (int)PRINT_VECTORS )
00173     out << "\nSigma_l =\n" << Sigma.diag();
00174   
00175   ele_wise_prod(1.0, invXu.diag(), Vu.diag(), &Sigma.diag() );
00176 
00177   if( (int)olevel >= (int)PRINT_ALGORITHM_STEPS )
00178     out << "\n||Sigma_k||inf = ||Sigma_l + Sigma_u||inf = " << Sigma.diag().norm_inf() << std::endl;
00179   if( (int)olevel >= (int)PRINT_VECTORS ) 
00180     out << "\nSigma_k = Sigma_l + Sigma_u =\n" << Sigma.diag();
00181 
00182   // Calculate the cross term (Z'*Sigma*Ypy) first
00183   VectorSpace::vec_mut_ptr_t temp = Z.space_cols().create_member(0.0);
00184   ele_wise_prod(1.0, s.Ypy().get_k(0), Sigma.diag(), temp.get());
00185   VectorMutable& w_sigma = s.w_sigma().set_k(0);
00186   V_MtV(&w_sigma, Z, BLAS_Cpp::trans, *temp);
00187 
00188   if( (int)olevel >= (int)PRINT_ALGORITHM_STEPS )
00189     out << "\n||w_sigma_k||inf = " << w_sigma.norm_inf() << std::endl;
00190   if( (int)olevel >= (int)PRINT_VECTORS ) 
00191     out << "\nw_sigma_k = \n" << w_sigma;
00192   
00193   // Calculate Reduced Sigma
00194   // Try sigma^1/2 making use of dependent and independent variables
00195 
00196   Teuchos::RCP<const VectorMutable>
00197     Sigma_D_diag = Sigma.diag().sub_view(Z.D_rng()),
00198     Sigma_I_diag = Sigma.diag().sub_view(Z.I_rng());
00199   const size_type
00200     Sigma_D_nz = Sigma_D_diag->nz(),
00201     Sigma_I_nz = Sigma_I_diag->nz();
00202 
00203   if( (int)olevel >= (int)PRINT_ALGORITHM_STEPS )
00204     {
00205     out << "\nSigma_D.diag().nz() = " << Sigma_D_nz;
00206     out << "\nSigma_I.diag().nz() = " << Sigma_I_nz << std::endl;
00207     }
00208 
00209   if( Sigma_D_nz || Sigma_I_nz )
00210     {
00211     if( (int)olevel >= (int)PRINT_ALGORITHM_STEPS )
00212       {
00213       out << "\nForming explicit, nonzero rHB_k = Z_k'*Sigma_k*Z_k ...\n";
00214       }
00215     if( Sigma_D_nz )
00216       {
00217 
00218       MatrixSymDiagStd Sigma_D_sqrt(Sigma_D_diag->clone());
00219 
00220       ele_wise_sqrt(&Sigma_D_sqrt.diag());
00221 
00222       if( static_cast<int>(olevel) >= static_cast<int>(PRINT_VECTORS) ) 
00223         {
00224         out << "\nSigma_D_sqrt =\n" << Sigma_D_sqrt;
00225         }
00226   
00227       Teuchos::RCP<MultiVectorMutable>
00228         J = Sigma_D_sqrt.space_cols().create_members(Z.cols());
00229       M_MtM(
00230         static_cast<MatrixOp*>(J.get())
00231         ,Sigma_D_sqrt, BLAS_Cpp::no_trans, Z.D(), BLAS_Cpp::no_trans);
00232 
00233       LinAlgOpPack::syrk( *J, BLAS_Cpp::trans, 1.0, 0.0, &rHB );
00234 
00235       if( static_cast<int>(olevel) >= static_cast<int>(PRINT_ITERATION_QUANTITIES) ) 
00236         {
00237         out << "\nJ =\n" << *J;
00238         out << "\nJ'*J =\n" << rHB;
00239         }
00240 
00241       }
00242 
00243     if( Sigma_I_nz )
00244       {
00245 
00246       const MatrixSymDiagStd Sigma_I(
00247         Teuchos::rcp_const_cast<VectorMutable>(Sigma_I_diag)
00248         );
00249 
00250       if(Sigma_D_nz)
00251         {
00252         Mp_M( &rHB, Sigma_I, BLAS_Cpp::no_trans );
00253         }
00254       else
00255         {
00256         assign( &rHB, Sigma_I, BLAS_Cpp::no_trans );
00257         }
00258 
00259       }
00260     }
00261   else
00262     {
00263     if( (int)olevel >= (int)PRINT_ALGORITHM_STEPS )
00264       {
00265       out << "\nSigma_k is zero so setting rHB_k = 0.0 ...\n";
00266       }
00267     rHB.zero_out();
00268     }
00269 
00270   /*
00271    // Try the full unspecialised calculation, this is expensive, but will
00272   // serve as a debug for the more efficient calculations.
00273 
00274   VectorSpace::multi_vec_mut_ptr_t Sigma_Z = Z.space_cols().create_members(Z.cols());
00275   M_MtM((MatrixOp*)Sigma_Z.get(), Sigma, BLAS_Cpp::no_trans, Z, BLAS_Cpp::no_trans);
00276 
00277   //std::cout << "Sigma_Z\n";
00278   //Sigma_Z->output(std::cout);
00279 
00280   VectorSpace::multi_vec_mut_ptr_t ZT_Sigma_Z = Z.space_rows().create_members(Z.cols());
00281   M_MtM((MatrixOp*)ZT_Sigma_Z.get(), (MatrixOp&)Z, BLAS_Cpp::trans, (MatrixOp&)*Sigma_Z, BLAS_Cpp::no_trans);
00282 
00283   std::cout << "ZT_Sigma_Z=\n";
00284   ZT_Sigma_Z->output(std::cout);
00285   */
00286   }
00287 
00288 
00289 namespace {
00290 
00291 const int local_num_options = 1;
00292 
00293 enum local_EOptions 
00294   {
00295   UPDATE_METHOD
00296   };
00297 
00298 const char* local_SOptions[local_num_options] = 
00299   {
00300   "update_method",
00301   };
00302 
00303 const int num_update_methods = 5;
00304 
00305 const char* s_update_methods[num_update_methods] = 
00306   {
00307   "always_explicit",
00308   "BFGS_primal",
00309   "BFGS_dual_no_correction",
00310   "BFGS_dual_explicit_correction",
00311   "BFGS_dual_scaling_correction"
00312   };
00313 
00314 }
00315 
00316  
00317 UpdateReducedSigma_StepSetOptions::UpdateReducedSigma_StepSetOptions(
00318   UpdateReducedSigma_Step* target
00319   , const char opt_grp_name[] )
00320   :
00321   OptionsFromStreamPack::SetOptionsFromStreamNode(
00322     opt_grp_name, local_num_options, local_SOptions ),
00323   OptionsFromStreamPack::SetOptionsToTargetBase< UpdateReducedSigma_Step >( target )
00324   {
00325   }
00326 
00327 void UpdateReducedSigma_StepSetOptions::setOption( 
00328   int option_num, const std::string& option_value )
00329   {
00330   using OptionsFromStreamPack::StringToIntMap;
00331 
00332   typedef UpdateReducedSigma_Step target_t;
00333   switch( (local_EOptions)option_num ) 
00334     {
00335     case UPDATE_METHOD:
00336       {
00337       StringToIntMap  config_map( UpdateReducedSigma_opt_grp_name, num_update_methods, s_update_methods );
00338       target().update_method( (UpdateReducedSigma_Step::e_update_methods) config_map( option_value.c_str() ) );
00339       }
00340       break;
00341     default:
00342       TEST_FOR_EXCEPT(true);  // Local error only?
00343     }
00344   }
00345 
00346 } // end namespace MoochoPack

Generated on Tue Jul 13 09:29:32 2010 for MoochoPack : Framework for Large-Scale Optimization Algorithms by  doxygen 1.4.7