MOOCHO (Single Doxygen Collection) Version of the Day
MoochoPack_CalcLambdaIndepStd_AddedStep.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 // Redistribution and use in source and binary forms, with or without
00013 // modification, are permitted provided that the following conditions are
00014 // met:
00015 //
00016 // 1. Redistributions of source code must retain the above copyright
00017 // notice, this list of conditions and the following disclaimer.
00018 //
00019 // 2. Redistributions in binary form must reproduce the above copyright
00020 // notice, this list of conditions and the following disclaimer in the
00021 // documentation and/or other materials provided with the distribution.
00022 //
00023 // 3. Neither the name of the Corporation nor the names of the
00024 // contributors may be used to endorse or promote products derived from
00025 // this software without specific prior written permission.
00026 //
00027 // THIS SOFTWARE IS PROVIDED BY SANDIA CORPORATION "AS IS" AND ANY
00028 // EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
00029 // IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR
00030 // PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL SANDIA CORPORATION OR THE
00031 // CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL,
00032 // EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO,
00033 // PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR
00034 // PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF
00035 // LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING
00036 // NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
00037 // SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
00038 //
00039 // Questions? Contact Roscoe A. Bartlett (rabartl@sandia.gov) 
00040 // 
00041 // ***********************************************************************
00042 // @HEADER
00043 
00044 #include <ostream>
00045 
00046 #include "MoochoPack_CalcLambdaIndepStd_AddedStep.hpp"
00047 #include "MoochoPack_moocho_algo_conversion.hpp"
00048 #include "IterationPack_print_algorithm_step.hpp"
00049 #include "ConstrainedOptPack_ComputeMinMult.hpp"
00050 #include "ConstrainedOptPack/src/VectorWithNorms.h"
00051 #include "AbstractLinAlgPack_SpVectorOp.hpp"
00052 #include "AbstractLinAlgPack/src/AbstractLinAlgPack_MatrixOp.hpp"
00053 #include "DenseLinAlgPack_LinAlgOpPack.hpp"
00054 #include "DenseLinAlgPack_DVectorOp.hpp"
00055 #include "DenseLinAlgPack_DVectorOut.hpp"
00056 
00057 namespace LinAlgOpPack {
00058   using AbstractLinAlgPack::Vp_StV;
00059   using AbstractLinAlgPack::Vp_StMtV;
00060 }
00061 
00062 bool MoochoPack::CalcLambdaIndepStd_AddedStep::do_step(Algorithm& _algo
00063   , poss_type step_poss, IterationPack::EDoStepType type, poss_type assoc_step_poss)
00064 {
00065 
00066   using BLAS_Cpp::no_trans;
00067   using BLAS_Cpp::trans;
00068 
00069   using DenseLinAlgPack::Vt_S;
00070   using DenseLinAlgPack::norm_inf;
00071 
00072   using AbstractLinAlgPack::Vp_StMtV;
00073   
00074   using ConstrainedOptPack::min_abs;
00075 
00076   using LinAlgOpPack::Vp_V;
00077   using LinAlgOpPack::V_StMtV;
00078   using LinAlgOpPack::Vp_MtV;
00079 
00080   NLPAlgo &algo = rsqp_algo(_algo);
00081   NLPAlgoState  &s    = algo.rsqp_state();
00082   Range1D   decomp  = s.equ_decomp();
00083   
00084   EJournalOutputLevel olevel = algo.algo_cntr().journal_output_level();
00085   std::ostream& out = algo.track().journal_out();
00086 
00087   // print step header.
00088   if( static_cast<int>(olevel) >= static_cast<int>(PRINT_ALGORITHM_STEPS) ) {
00089     using IterationPack::print_algorithm_step;
00090     print_algorithm_step( algo, step_poss, type, assoc_step_poss, out );
00091   }
00092 
00093   // Compute: lambda(decomp) = inv(Gc(decomp)'* Y)' * ( - Y' * (Gf + nu) - U' * lambda(undecomp) )
00094   // where U = Gc(undecomp)' * Y
00095 
00096   // Must resize lambda here explicitly since we will only be updating a region of it.
00097   // If lambda(undecomp) has already been updated then lambda will have been resized
00098   // already but lambda(decomp) will not be initialized yet.
00099   if( !s.lambda().updated_k(0) ) s.lambda().set_k(0).v().resize( algo.nlp().m() );
00100   
00101   DVectorSlice lambda_decomp = s.lambda().get_k(0).v()(decomp);
00102   
00103   // lambda_decomp_tmp1 = - Y' * (Gf + nu)
00104   if( algo.nlp().has_bounds() ) {
00105     // _tmp = Gf + nu
00106     DVector _tmp = s.Gf().get_k(0)();
00107     DVectorSlice _vs_tmp = _tmp;  // only create this DVectorSlice once
00108     Vp_V( &_vs_tmp, s.nu().get_k(0)() );
00109     // lambda_decomp_tmp1 = - Y' * _tmp
00110     V_StMtV( &lambda_decomp, -1.0, s.Y().get_k(0), trans, _vs_tmp );
00111   }
00112   else {
00113     // lambda_decomp__tmp1 = - Y' * Gf
00114     V_StMtV( &lambda_decomp, -1.0, s.Y().get_k(0), trans, s.Gf().get_k(0)() );
00115   }
00116 
00117   // lambda_decomp_tmp2 = lambda_decomp_tmp1 - U' * lambda(undecomp)
00118   if( algo.nlp().r() < algo.nlp().m() ) {
00119     Range1D undecomp = s.equ_undecomp();
00120     Vp_StMtV( &lambda_decomp, -1.0, s.U().get_k(0), trans, s.lambda().get_k(0).v()(undecomp) );
00121   }
00122   // else lambda(decomp)_tmp2 = lambda(decomp)_tmp1
00123 
00124   // lambda_decomp = inv(Gc(decomp)'* Y)' * lambda_decomp_tmp2
00125   s.decomp_sys().solve_transAtY( lambda_decomp, trans, &lambda_decomp );
00126 
00127   if( static_cast<int>(olevel) >= static_cast<int>(PRINT_ALGORITHM_STEPS) ) {
00128     out << "\nmax(|lambda_k(equ_decomp)(i)|) = " << norm_inf(lambda_decomp)
00129       << "\nmin(|lambda_k(equ_decomp)(i)|) = " << min_abs(lambda_decomp)  << std::endl;
00130   }
00131 
00132   if( static_cast<int>(olevel) >= static_cast<int>(PRINT_VECTORS) ) {
00133     out << "\nlambda_k(equ_decomp) = \n" << lambda_decomp;
00134   }
00135 
00136   return true;
00137 }
00138 
00139 void MoochoPack::CalcLambdaIndepStd_AddedStep::print_step( const Algorithm& algo
00140   , poss_type step_poss, IterationPack::EDoStepType type, poss_type assoc_step_poss
00141   , std::ostream& out, const std::string& L ) const
00142 {
00143   out
00144     << L << "*** Calculate the Lagrange multipliers for the decomposed constraints\n"
00145     << L << "lambda_k(equ_decomp) = - inv(Gc_k(:,equ_decomp)'*Y_k)\n"
00146     << L << "                        * (Y_k'*(Gf_k + nu_k) + U_k'*lambda_k(equ_undecomp))\n";
00147 }
00148 
00149 #endif // 0
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Defines