MoochoPack : Framework for Large-Scale Optimization Algorithms Version of the Day
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 // Redistribution and use in source and binary forms, with or without
00011 // modification, are permitted provided that the following conditions are
00012 // met:
00013 //
00014 // 1. Redistributions of source code must retain the above copyright
00015 // notice, this list of conditions and the following disclaimer.
00016 //
00017 // 2. Redistributions in binary form must reproduce the above copyright
00018 // notice, this list of conditions and the following disclaimer in the
00019 // documentation and/or other materials provided with the distribution.
00020 //
00021 // 3. Neither the name of the Corporation nor the names of the
00022 // contributors may be used to endorse or promote products derived from
00023 // this software without specific prior written permission.
00024 //
00025 // THIS SOFTWARE IS PROVIDED BY SANDIA CORPORATION "AS IS" AND ANY
00026 // EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
00027 // IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR
00028 // PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL SANDIA CORPORATION OR THE
00029 // CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL,
00030 // EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO,
00031 // PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR
00032 // PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF
00033 // LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING
00034 // NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
00035 // SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
00036 //
00037 // Questions? Contact Roscoe A. Bartlett (rabartl@sandia.gov) 
00038 // 
00039 // ***********************************************************************
00040 // @HEADER
00041 
00042 #include <ostream>
00043 #include <typeinfo>
00044 #include <iostream>
00045 
00046 #include "AbstractLinAlgPack_MatrixSymDiagStd.hpp"
00047 #include "AbstractLinAlgPack_MatrixSymOpNonsing.hpp"
00048 #include "AbstractLinAlgPack_MatrixOpOut.hpp"
00049 #include "AbstractLinAlgPack_MultiVectorMutable.hpp"
00050 #include "AbstractLinAlgPack_VectorStdOps.hpp"
00051 #include "AbstractLinAlgPack_VectorOut.hpp"
00052 #include "AbstractLinAlgPack_VectorAuxiliaryOps.hpp"
00053 #include "AbstractLinAlgPack_LinAlgOpPack.hpp"
00054 #include "AbstractLinAlgPack_assert_print_nan_inf.hpp"
00055 #include "ConstrainedOptPack_MatrixIdentConcat.hpp"
00056 #include "IterationPack_print_algorithm_step.hpp"
00057 #include "MoochoPack_IpState.hpp"
00058 #include "MoochoPack_UpdateReducedSigma_Step.hpp"
00059 #include "MoochoPack_moocho_algo_conversion.hpp"
00060 
00061 #include "OptionsFromStreamPack_StringToIntMap.hpp"
00062 
00063 #include "Teuchos_dyn_cast.hpp"
00064 
00065 namespace MoochoPack {
00066 
00067 UpdateReducedSigma_Step::UpdateReducedSigma_Step(
00068   const e_update_methods update_method
00069   )
00070   :
00071   update_method_(update_method)
00072   {}
00073 
00074 bool UpdateReducedSigma_Step::do_step(
00075   Algorithm& _algo, poss_type step_poss, IterationPack::EDoStepType type
00076   ,poss_type assoc_step_poss
00077   )
00078   {
00079   using Teuchos::dyn_cast;
00080   using IterationPack::print_algorithm_step;
00081 
00082   NLPAlgo            &algo   = dyn_cast<NLPAlgo>(_algo);
00083   IpState             &s      = dyn_cast<IpState>(_algo.state());
00084   
00085   EJournalOutputLevel olevel  = algo.algo_cntr().journal_output_level();
00086   std::ostream        &out    = algo.track().journal_out();
00087   
00088   // print step header.
00089   if( static_cast<int>(olevel) >= static_cast<int>(PRINT_ALGORITHM_STEPS) ) 
00090     {
00091     using IterationPack::print_algorithm_step;
00092     print_algorithm_step( _algo, step_poss, type, assoc_step_poss, out );
00093     }
00094 
00095   switch (update_method_)
00096     {
00097     case ALWAYS_EXPLICIT:
00098       {
00099       if( (int)olevel >= (int)PRINT_ALGORITHM_STEPS ) 
00100         {
00101         out << "\nupdate_method is always_explicit, forming Reduced Sigma Explicitly ...\n";
00102         }
00103       FormReducedSigmaExplicitly(algo,s,olevel,out);
00104       }
00105       break;
00106     case BFGS_PRIMAL:
00107     case BFGS_DUAL_NO_CORRECTION:
00108     case BFGS_DUAL_EXPLICIT_CORRECTION:
00109     case BFGS_DUAL_SCALING_CORRECTION:
00110       {
00111       TEUCHOS_TEST_FOR_EXCEPTION(true, std::logic_error, "Option BFGS_primal not handled yet in UpdateReducedSigma\n");
00112       }
00113       break;
00114     default:
00115       TEUCHOS_TEST_FOR_EXCEPT(true); // local error ?
00116     };
00117 
00118   if( (int)olevel >= (int)PRINT_ITERATION_QUANTITIES ) 
00119     {
00120     out << "\nrHB_k =\n" << s.rHB().get_k(0);
00121     }
00122 
00123   return true;
00124   }
00125 
00126 
00127 void UpdateReducedSigma_Step::print_step(
00128   const Algorithm& _algo, poss_type step_poss, IterationPack::EDoStepType type
00129   ,poss_type assoc_step_poss, std::ostream& out, const std::string& L
00130   ) const
00131   {
00132   out << L << "*** Update Z'*Sigma*Z\n"
00133     << L << "if (update_method == always_explicit) then\n"
00134     << L << "  Sigma_k = invXl*Vl-invXu*Vu\n"
00135     << L << "  Sigma_I = Sigma_k.sub_view(Z.I_rng)\n"
00136     << L << "  Sigma_D_sqrt = (Sigma_k.sub_view(Z.D_rng))^1/2\n"
00137     << L << "  J = Sigma_D_sqrt*Z\n"
00138     << L << "  rHB_k = J'*J + Sigma_I\n"
00139     << L << "elsif ( update_method == BFGS_??? ) then\n"
00140     << L << "  NOT IMPLEMENTED YET!\n"
00141     << L << "end\n";
00142   }
00143 
00144 void UpdateReducedSigma_Step::FormReducedSigmaExplicitly(
00145   NLPAlgo& algo, IpState& s, EJournalOutputLevel olevel,  std::ostream& out
00146   )
00147   {
00148   namespace mmp = MemMngPack;
00149   using Teuchos::dyn_cast;
00150   using AbstractLinAlgPack::ele_wise_prod;
00151   using AbstractLinAlgPack::ele_wise_sqrt;
00152    using LinAlgOpPack::Mp_M;
00153    using LinAlgOpPack::Mp_MtM;
00154    using LinAlgOpPack::M_MtM;
00155    using LinAlgOpPack::M_StM;
00156   using LinAlgOpPack::V_MtV;
00157   using LinAlgOpPack::assign;
00158 
00159   // Calculate Reduced Sigma directly from
00160   // Sigma = invXl*Vl + invXu*Vu
00161   // Z_kT*Sigma*Z_k
00162 
00163   // Get the iteration quantities
00164   const MatrixIdentConcat     &Z     = dyn_cast<MatrixIdentConcat>(s.Z().get_k(0));
00165   const MatrixSymDiagStd  &invXl = s.invXl().get_k(0);
00166   const MatrixSymDiagStd  &invXu = s.invXu().get_k(0);
00167   const MatrixSymDiagStd  &Vl    = s.Vl().get_k(0);
00168   const MatrixSymDiagStd  &Vu    = s.Vu().get_k(0);
00169   
00170   MatrixSymDiagStd  &Sigma = s.Sigma().set_k(0);
00171 
00172   MatrixSymOpNonsing& rHB = dyn_cast<MatrixSymOpNonsing>(s.rHB().set_k(0));
00173   if (rHB.cols() != Z.cols())
00174     {
00175     // Initialize space in rHB
00176     dyn_cast<MatrixSymInitDiag>(rHB).init_identity(Z.space_rows(), 0.0);
00177     }
00178   
00179   // Calculate Sigma = invXl*Vl + invXu*Vu
00180 
00181   ele_wise_prod(1.0, invXl.diag(), Vl.diag(), &(Sigma.diag() = 0.0));
00182   
00183   if( (int)olevel >= (int)PRINT_ALGORITHM_STEPS )
00184     out << "\n||Sigma_l||inf = " << Sigma.diag().norm_inf() << std::endl;
00185   if( (int)olevel >= (int)PRINT_VECTORS )
00186     out << "\nSigma_l =\n" << Sigma.diag();
00187   
00188   ele_wise_prod(1.0, invXu.diag(), Vu.diag(), &Sigma.diag() );
00189 
00190   if( (int)olevel >= (int)PRINT_ALGORITHM_STEPS )
00191     out << "\n||Sigma_k||inf = ||Sigma_l + Sigma_u||inf = " << Sigma.diag().norm_inf() << std::endl;
00192   if( (int)olevel >= (int)PRINT_VECTORS ) 
00193     out << "\nSigma_k = Sigma_l + Sigma_u =\n" << Sigma.diag();
00194 
00195   // Calculate the cross term (Z'*Sigma*Ypy) first
00196   VectorSpace::vec_mut_ptr_t temp = Z.space_cols().create_member(0.0);
00197   ele_wise_prod(1.0, s.Ypy().get_k(0), Sigma.diag(), temp.get());
00198   VectorMutable& w_sigma = s.w_sigma().set_k(0);
00199   V_MtV(&w_sigma, Z, BLAS_Cpp::trans, *temp);
00200 
00201   if( (int)olevel >= (int)PRINT_ALGORITHM_STEPS )
00202     out << "\n||w_sigma_k||inf = " << w_sigma.norm_inf() << std::endl;
00203   if( (int)olevel >= (int)PRINT_VECTORS ) 
00204     out << "\nw_sigma_k = \n" << w_sigma;
00205   
00206   // Calculate Reduced Sigma
00207   // Try sigma^1/2 making use of dependent and independent variables
00208 
00209   Teuchos::RCP<const VectorMutable>
00210     Sigma_D_diag = Sigma.diag().sub_view(Z.D_rng()),
00211     Sigma_I_diag = Sigma.diag().sub_view(Z.I_rng());
00212   const size_type
00213     Sigma_D_nz = Sigma_D_diag->nz(),
00214     Sigma_I_nz = Sigma_I_diag->nz();
00215 
00216   if( (int)olevel >= (int)PRINT_ALGORITHM_STEPS )
00217     {
00218     out << "\nSigma_D.diag().nz() = " << Sigma_D_nz;
00219     out << "\nSigma_I.diag().nz() = " << Sigma_I_nz << std::endl;
00220     }
00221 
00222   if( Sigma_D_nz || Sigma_I_nz )
00223     {
00224     if( (int)olevel >= (int)PRINT_ALGORITHM_STEPS )
00225       {
00226       out << "\nForming explicit, nonzero rHB_k = Z_k'*Sigma_k*Z_k ...\n";
00227       }
00228     if( Sigma_D_nz )
00229       {
00230 
00231       MatrixSymDiagStd Sigma_D_sqrt(Sigma_D_diag->clone());
00232 
00233       ele_wise_sqrt(&Sigma_D_sqrt.diag());
00234 
00235       if( static_cast<int>(olevel) >= static_cast<int>(PRINT_VECTORS) ) 
00236         {
00237         out << "\nSigma_D_sqrt =\n" << Sigma_D_sqrt;
00238         }
00239   
00240       Teuchos::RCP<MultiVectorMutable>
00241         J = Sigma_D_sqrt.space_cols().create_members(Z.cols());
00242       M_MtM(
00243         static_cast<MatrixOp*>(J.get())
00244         ,Sigma_D_sqrt, BLAS_Cpp::no_trans, Z.D(), BLAS_Cpp::no_trans);
00245 
00246       LinAlgOpPack::syrk( *J, BLAS_Cpp::trans, 1.0, 0.0, &rHB );
00247 
00248       if( static_cast<int>(olevel) >= static_cast<int>(PRINT_ITERATION_QUANTITIES) ) 
00249         {
00250         out << "\nJ =\n" << *J;
00251         out << "\nJ'*J =\n" << rHB;
00252         }
00253 
00254       }
00255 
00256     if( Sigma_I_nz )
00257       {
00258 
00259       const MatrixSymDiagStd Sigma_I(
00260         Teuchos::rcp_const_cast<VectorMutable>(Sigma_I_diag)
00261         );
00262 
00263       if(Sigma_D_nz)
00264         {
00265         Mp_M( &rHB, Sigma_I, BLAS_Cpp::no_trans );
00266         }
00267       else
00268         {
00269         assign( &rHB, Sigma_I, BLAS_Cpp::no_trans );
00270         }
00271 
00272       }
00273     }
00274   else
00275     {
00276     if( (int)olevel >= (int)PRINT_ALGORITHM_STEPS )
00277       {
00278       out << "\nSigma_k is zero so setting rHB_k = 0.0 ...\n";
00279       }
00280     rHB.zero_out();
00281     }
00282 
00283   /*
00284    // Try the full unspecialised calculation, this is expensive, but will
00285   // serve as a debug for the more efficient calculations.
00286 
00287   VectorSpace::multi_vec_mut_ptr_t Sigma_Z = Z.space_cols().create_members(Z.cols());
00288   M_MtM((MatrixOp*)Sigma_Z.get(), Sigma, BLAS_Cpp::no_trans, Z, BLAS_Cpp::no_trans);
00289 
00290   //std::cout << "Sigma_Z\n";
00291   //Sigma_Z->output(std::cout);
00292 
00293   VectorSpace::multi_vec_mut_ptr_t ZT_Sigma_Z = Z.space_rows().create_members(Z.cols());
00294   M_MtM((MatrixOp*)ZT_Sigma_Z.get(), (MatrixOp&)Z, BLAS_Cpp::trans, (MatrixOp&)*Sigma_Z, BLAS_Cpp::no_trans);
00295 
00296   std::cout << "ZT_Sigma_Z=\n";
00297   ZT_Sigma_Z->output(std::cout);
00298   */
00299   }
00300 
00301 
00302 namespace {
00303 
00304 const int local_num_options = 1;
00305 
00306 enum local_EOptions 
00307   {
00308   UPDATE_METHOD
00309   };
00310 
00311 const char* local_SOptions[local_num_options] = 
00312   {
00313   "update_method",
00314   };
00315 
00316 const int num_update_methods = 5;
00317 
00318 const char* s_update_methods[num_update_methods] = 
00319   {
00320   "always_explicit",
00321   "BFGS_primal",
00322   "BFGS_dual_no_correction",
00323   "BFGS_dual_explicit_correction",
00324   "BFGS_dual_scaling_correction"
00325   };
00326 
00327 }
00328 
00329  
00330 UpdateReducedSigma_StepSetOptions::UpdateReducedSigma_StepSetOptions(
00331   UpdateReducedSigma_Step* target
00332   , const char opt_grp_name[] )
00333   :
00334   OptionsFromStreamPack::SetOptionsFromStreamNode(
00335     opt_grp_name, local_num_options, local_SOptions ),
00336   OptionsFromStreamPack::SetOptionsToTargetBase< UpdateReducedSigma_Step >( target )
00337   {
00338   }
00339 
00340 void UpdateReducedSigma_StepSetOptions::setOption( 
00341   int option_num, const std::string& option_value )
00342   {
00343   using OptionsFromStreamPack::StringToIntMap;
00344 
00345   typedef UpdateReducedSigma_Step target_t;
00346   switch( (local_EOptions)option_num ) 
00347     {
00348     case UPDATE_METHOD:
00349       {
00350       StringToIntMap  config_map( UpdateReducedSigma_opt_grp_name, num_update_methods, s_update_methods );
00351       target().update_method( (UpdateReducedSigma_Step::e_update_methods) config_map( option_value.c_str() ) );
00352       }
00353       break;
00354     default:
00355       TEUCHOS_TEST_FOR_EXCEPT(true);  // Local error only?
00356     }
00357   }
00358 
00359 } // end namespace MoochoPack
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends