MOOCHO (Single Doxygen Collection) Version of the Day
ConstrainedOptPack_MatrixSymPosDefLBFGS.cpp
Go to the documentation of this file.
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 // Let's define a compact representation for the matrix B^{k} and
00043 // its inverse H^{k} = inv(B^{k}).
00044 //
00045 // Bk = (1/gk)*I - [ (1/gk)*S  Y ] * inv( [ (1/gk)*S'S   L ]     [ (1/gk)*S' ]
00046 //                                        [    L'       -D ] ) * [    Y'     ]
00047 //                                        \________________/
00048 //                                                Q
00049 //
00050 // Hk = gk*I + [ S  gk*Y ] * [ inv(R')*(D+gk*Y'Y)*inv(R)     -inv(R') ] * [   S'  ]
00051 //                           [            -inv(R)                0    ]   [ gk*Y' ]
00052 //
00053 // where:
00054 //
00055 // gk = gamma_k <: R
00056 //
00057 //      [ s^{1}'*y^{1}  s^{1}'*y^{2}  ...   s^{1}'*y^{m}  ]
00058 //  S'Y = [ s^{2}'*y^{1}  s^{2}'*y^{2}  ...   s^{2}'*y^{m}  ] <: R^(m x m)
00059 //      [ .       .           .       ]
00060 //      [ s^{m}'*y^{1}  s^{m}'*y^{2}  ...   s^{m}'*y^{m}  ]
00061 // 
00062 //      [ s^{1}'*y^{1}  0       ...   0       ]
00063 //  D =   [ 0       s^{2}'*y^{2}  ...   0       ] <: R^(m x m)
00064 //      [ .       .           .       ]
00065 //      [ 0       0       ...   s^{m}'*y^{m}  ]
00066 //
00067 //  R = upper triangular part of S'Y
00068 // 
00069 //  L = lower tirangular part of S'Y with zeros on the diagonal
00070 //
00071 
00072 #include <assert.h>
00073 
00074 #include <typeinfo>
00075 
00076 #include "ConstrainedOptPack_MatrixSymPosDefLBFGS.hpp"
00077 #include "AbstractLinAlgPack_BFGS_helpers.hpp"
00078 #include "AbstractLinAlgPack_VectorDenseEncap.hpp"
00079 #include "AbstractLinAlgPack_MatrixOpOut.hpp"
00080 #include "AbstractLinAlgPack_VectorStdOps.hpp"
00081 #include "AbstractLinAlgPack_LinAlgOpPack.hpp"
00082 #include "DenseLinAlgPack_LinAlgOpPack.hpp"
00083 #include "DenseLinAlgPack_DMatrixOut.hpp"
00084 #include "DenseLinAlgLAPack.hpp"
00085 #include "Teuchos_Workspace.hpp"
00086 #include "Teuchos_Assert.hpp"
00087 
00088 namespace {
00089 
00090   using DenseLinAlgPack::DVectorSlice;
00091   using DenseLinAlgPack::DMatrixSlice;
00092 
00100   void comp_Cb( const DMatrixSlice& Lb, const DVectorSlice& Db_diag
00101     , DMatrixSlice* Cb );
00102 
00103 } // end namespace
00104 
00105 namespace ConstrainedOptPack {
00106 
00107 // /////////////////////////////////
00108 // Inline private member functions
00109 
00110 inline
00111 const DMatrixSliceTri MatrixSymPosDefLBFGS::R() const
00112 {
00113   return DenseLinAlgPack::tri( STY_(1,m_bar_,1,m_bar_), BLAS_Cpp::upper, BLAS_Cpp::nonunit );
00114 }
00115 
00116 inline
00117 const DMatrixSliceTri MatrixSymPosDefLBFGS::Lb() const
00118 {
00119   return DenseLinAlgPack::tri( STY_(2,m_bar_,1,m_bar_-1), BLAS_Cpp::lower, BLAS_Cpp::nonunit );
00120 }
00121 
00122 inline
00123 DMatrixSlice MatrixSymPosDefLBFGS::STY()
00124 {
00125   return STY_(1,m_bar_,1,m_bar_);
00126 }
00127 
00128 inline
00129 const DMatrixSlice MatrixSymPosDefLBFGS::STY() const
00130 {
00131   return STY_(1,m_bar_,1,m_bar_);
00132 }
00133 
00134 inline
00135 DMatrixSliceSym MatrixSymPosDefLBFGS::STS()
00136 {
00137   return DenseLinAlgPack::nonconst_sym( STSYTY_(2,m_bar_+1,1,m_bar_),BLAS_Cpp::lower );
00138 }
00139 
00140 inline
00141 const DMatrixSliceSym MatrixSymPosDefLBFGS::STS() const
00142 {
00143   return DenseLinAlgPack::sym( STSYTY_(2,m_bar_+1,1,m_bar_),BLAS_Cpp::lower );
00144 }
00145 
00146 inline
00147 DMatrixSliceSym MatrixSymPosDefLBFGS::YTY()
00148 {
00149   return DenseLinAlgPack::nonconst_sym( STSYTY_(1,m_bar_,2,m_bar_+1),BLAS_Cpp::upper );
00150 }
00151 
00152 inline
00153 const DMatrixSliceSym MatrixSymPosDefLBFGS::YTY() const
00154 {
00155   return DenseLinAlgPack::sym( STSYTY_(1,m_bar_,2,m_bar_+1),BLAS_Cpp::upper );
00156 }
00157 
00158 // ///////////////////////
00159 // Nonlinined functions
00160 
00161 MatrixSymPosDefLBFGS::MatrixSymPosDefLBFGS(
00162     size_type   m
00163   ,bool       maintain_original
00164   ,bool       maintain_inverse
00165   ,bool       auto_rescaling
00166   )
00167 {
00168   initial_setup(m,maintain_original,maintain_inverse,auto_rescaling);
00169 }
00170 
00171 void MatrixSymPosDefLBFGS::initial_setup(
00172   size_type m,
00173   bool maintain_original,
00174   bool maintain_inverse,
00175   bool auto_rescaling
00176   )
00177 {
00178   // Validate input
00179   TEUCHOS_TEST_FOR_EXCEPTION(
00180     !maintain_original && !maintain_inverse, std::invalid_argument
00181     ,"MatrixSymPosDefLBFGS::initial_setup(...) : "
00182     "Error, both maintain_original and maintain_inverse can not both be false!" );
00183   TEUCHOS_TEST_FOR_EXCEPTION(
00184     m < 1, std::invalid_argument
00185     ,"MatrixSymPosDefLBFGS::set_num_updates_stored(m) : "
00186     "Error, the number of storage locations must be > 0" );
00187   vec_spc_           = Teuchos::null;
00188   maintain_original_ = maintain_original;
00189   maintain_inverse_  = maintain_inverse;
00190   m_                 = m;
00191   n_                 = 0; // make uninitialized
00192   num_secant_updates_= 0;
00193   auto_rescaling_    = auto_rescaling;
00194 }
00195 
00196 // Overridden from MatrixOp
00197 
00198 const VectorSpace& MatrixSymPosDefLBFGS::space_cols() const
00199 {
00200   return *vec_spc_;
00201 }
00202 
00203 std::ostream& MatrixSymPosDefLBFGS::output(std::ostream& out) const
00204 {
00205   assert_initialized();
00206   out << "\n*** Limited Memory BFGS matrix\n";
00207   out << "\nConversion to dense =\n";
00208   MatrixOp::output(out);
00209   out << "\nStored quantities\n"
00210     << "\nn       = " << n_
00211     << "\nm       = " << m_
00212     << "\nm_bar   = " << m_bar_
00213     << "\ngamma_k = " << gamma_k_ << std::endl;
00214   if( m_bar_ ) {
00215     out << "\nS =\n" << *S()
00216       << "\nY =\n" << *Y()
00217       << "\nS'Y =\n" << STY_(1,m_bar_,1,m_bar_)
00218       << "\nlower(S'S) \\ zero diagonal \\ upper(Y'Y) =\n"
00219         << STSYTY_(1,m_bar_+1,1,m_bar_+1)
00220       << "\nQ updated? = " << Q_updated_ << std::endl;
00221     if(Q_updated_)
00222       out << "\nCholesky of schur complement of Q, QJ =\n" << QJ_(1,m_bar_,1,m_bar_);
00223   }
00224   return out;
00225 }
00226 
00227 MatrixOp& MatrixSymPosDefLBFGS::operator=(const MatrixOp& mwo)
00228 { 
00229   const MatrixSymPosDefLBFGS *p_m = dynamic_cast<const MatrixSymPosDefLBFGS*>(&mwo);
00230   if(p_m) {
00231     if( p_m == this ) return *this; // assignment to self
00232     // Important: Assign all members here.
00233     auto_rescaling_      = p_m->auto_rescaling_;
00234     maintain_original_   = p_m->maintain_original_;
00235     original_is_updated_ = p_m->original_is_updated_;
00236     maintain_inverse_    = p_m->maintain_inverse_;
00237     inverse_is_updated_  = p_m->inverse_is_updated_;
00238     vec_spc_             = p_m->vec_spc_.get() ? p_m->vec_spc_->clone() : Teuchos::null;
00239     n_               = p_m->n_;
00240     m_               = p_m->m_;
00241     m_bar_             = p_m->m_bar_;
00242     num_secant_updates_  = p_m->num_secant_updates_;
00243     gamma_k_           = p_m->gamma_k_;
00244     S_               = p_m->S_;
00245     Y_               = p_m->Y_;
00246     STY_             = p_m->STY_;
00247     STSYTY_            = p_m->STSYTY_;
00248     Q_updated_           = p_m->Q_updated_;
00249     QJ_              = p_m->QJ_;
00250   }
00251   else {
00252     TEUCHOS_TEST_FOR_EXCEPTION(
00253       true,std::invalid_argument
00254       ,"MatrixSymPosDefLBFGS::operator=(const MatrixOp& mwo) : Error, "
00255       "The concrete type of mwo \'" << typeName(mwo) << "\' is not "
00256       "as subclass of MatrixSymPosDefLBFGS as required" );
00257   }
00258   return *this;
00259 }
00260 
00261 // Level-2 BLAS
00262 
00263 void MatrixSymPosDefLBFGS::Vp_StMtV(
00264     VectorMutable* y, value_type alpha, BLAS_Cpp::Transp trans_rhs1
00265   , const Vector& x, value_type beta
00266   ) const
00267 {
00268   using AbstractLinAlgPack::Vt_S;
00269   using AbstractLinAlgPack::Vp_StV;
00270   using AbstractLinAlgPack::Vp_StMtV;
00271   using LinAlgOpPack::V_StMtV;
00272   using LinAlgOpPack::V_MtV;
00273   typedef VectorDenseEncap         vde;
00274   typedef VectorDenseMutableEncap  vdme;
00275   using Teuchos::Workspace;
00276   Teuchos::WorkspaceStore* wss = Teuchos::get_default_workspace_store().get();
00277 
00278   assert_initialized();
00279 
00280   TEUCHOS_TEST_FOR_EXCEPT( !(  original_is_updated_  ) ); // For now just always update
00281 
00282   // y = b*y + Bk * x
00283   //
00284   // y = b*y + (1/gk)*x - [ (1/gk)*S  Y ] * inv(Q) * [ (1/gk)*S' ] * x
00285   //                                                 [     Y'    ]
00286   // Perform the following operations (in order):
00287   //
00288   // y = b*y
00289   //
00290   // y += (1/gk)*x
00291   //
00292   // t1 = [ (1/gk)*S'*x ]   <: R^(2*m)
00293   //    [      Y'*x   ]
00294   //
00295   // t2 = inv(Q) * t1     <: R^(2*m)
00296   //
00297   // y += -(1/gk) * S * t2(1:m)
00298   //
00299   // y += -1.0 * Y * t2(m+1,2m)
00300 
00301   const value_type
00302     invgk = 1.0 / gamma_k_;
00303 
00304   // y = b*y
00305   Vt_S( y, beta );
00306 
00307   // y += (1/gk)*x
00308   Vp_StV( y, invgk, x );
00309 
00310   if( !m_bar_ )
00311     return; // No updates have been added yet.
00312 
00313   const multi_vec_ptr_t
00314     S = this->S(),
00315     Y = this->Y();
00316 
00317   // Get workspace
00318 
00319   const size_type
00320     mb = m_bar_;
00321 
00322   Workspace<value_type>  t1_ws(wss,2*mb);
00323   DVectorSlice                 t1(&t1_ws[0],t1_ws.size());
00324   Workspace<value_type>  t2_ws(wss,2*mb);
00325   DVectorSlice                 t2(&t2_ws[0],t2_ws.size());
00326 
00327   VectorSpace::vec_mut_ptr_t
00328     t = S->space_rows().create_member();
00329 
00330   // t1 = [ (1/gk)*S'*x ]
00331   //    [      Y'*x   ]
00332 
00333   V_StMtV( t.get(), invgk, *S, BLAS_Cpp::trans, x );
00334   t1(1,mb) = vde(*t)();
00335   V_MtV( t.get(), *Y, BLAS_Cpp::trans, x );
00336   t1(mb+1,2*mb) = vde(*t)();
00337 
00338   // t2 = inv(Q) * t1
00339   V_invQtV( &t2, t1 );
00340 
00341   // y += -(1/gk) * S * t2(1:m)
00342   (vdme(*t)() = t2(1,mb));
00343   Vp_StMtV( y, -invgk, *S, BLAS_Cpp::no_trans, *t );
00344 
00345   // y += -1.0 * Y * t2(m+1,2m
00346   (vdme(*t)() = t2(mb+1,2*mb));
00347   Vp_StMtV( y, -1.0, *Y, BLAS_Cpp::no_trans, *t );
00348 
00349 }
00350 
00351 // Overridden from MatrixOpNonsing
00352 
00353 void MatrixSymPosDefLBFGS::V_InvMtV(
00354   VectorMutable* y, BLAS_Cpp::Transp trans_rhs1
00355   , const Vector& x
00356   ) const
00357 {
00358   using AbstractLinAlgPack::Vp_StMtV;
00359   using DenseLinAlgPack::V_InvMtV;
00360   using LinAlgOpPack::V_mV;
00361   using LinAlgOpPack::V_StV;
00362   using LinAlgOpPack::Vp_V;
00363   using LinAlgOpPack::V_MtV;
00364   using LinAlgOpPack::V_StMtV;
00365   using LinAlgOpPack::Vp_MtV;
00366   using DenseLinAlgPack::Vp_StMtV;
00367   typedef VectorDenseEncap         vde;
00368   typedef VectorDenseMutableEncap  vdme;
00369   using Teuchos::Workspace;
00370   Teuchos::WorkspaceStore* wss = Teuchos::get_default_workspace_store().get();
00371 
00372   assert_initialized();
00373 
00374   TEUCHOS_TEST_FOR_EXCEPT( !(  inverse_is_updated_  ) ); // For now just always update
00375 
00376   // y = inv(Bk) * x = Hk * x
00377   //
00378   // = gk*x + [S gk*Y] * [ inv(R')*(D+gk*Y'Y)*inv(R)     -inv(R') ] * [   S'  ] * x
00379   //                     [            -inv(R)                0    ]   [ gk*Y' ]
00380   //
00381   // Perform in the following (in order):
00382   //
00383   // y = gk*x
00384   //
00385   // t1 = [   S'*x  ]         <: R^(2*m)
00386   //      [ gk*Y'*x ]
00387   //
00388   // t2 = inv(R) * t1(1:m)      <: R^(m)
00389   //
00390   // t3 = - inv(R') * t1(m+1,2*m)   <: R^(m)
00391   //
00392   // t4 = gk * Y'Y * t2       <: R^(m)
00393   //
00394   // t4 += D*t2
00395   //
00396   // t5 = inv(R') * t4        <: R^(m)
00397   //
00398   // t5 += t3
00399   //
00400   // y += S*t5
00401   //
00402   // y += -gk*Y*t2
00403 
00404   // y = gk*x
00405   V_StV( y, gamma_k_, x );
00406 
00407   const size_type
00408     mb = m_bar_;  
00409   
00410   if( !mb )
00411     return; // No updates have been performed.
00412 
00413   const multi_vec_ptr_t
00414     S = this->S(),
00415     Y = this->Y();
00416 
00417   // Get workspace
00418 
00419   Workspace<value_type>    t1_ws(wss,2*mb);
00420   DVectorSlice                   t1(&t1_ws[0],t1_ws.size());
00421   Workspace<value_type>    t2_ws(wss,mb);
00422   DVectorSlice                   t2(&t2_ws[0],t2_ws.size());
00423   Workspace<value_type>    t3_ws(wss,mb);
00424   DVectorSlice                   t3(&t3_ws[0],t3_ws.size());
00425   Workspace<value_type>    t4_ws(wss,mb);
00426   DVectorSlice                   t4(&t4_ws[0],t4_ws.size());
00427   Workspace<value_type>    t5_ws(wss,mb);
00428   DVectorSlice                   t5(&t5_ws[0],t5_ws.size());
00429 
00430   VectorSpace::vec_mut_ptr_t
00431     t = S->space_rows().create_member();
00432 
00433   const DMatrixSliceTri
00434     &R = this->R();
00435 
00436   const DMatrixSliceSym
00437     &YTY = this->YTY();
00438 
00439   // t1 = [   S'*x  ]
00440   //      [ gk*Y'*x ]
00441   V_MtV( t.get(), *S, BLAS_Cpp::trans, x );
00442   t1(1,mb) = vde(*t)();
00443   V_StMtV( t.get(), gamma_k_, *Y, BLAS_Cpp::trans, x );
00444   t1(mb+1,2*mb) = vde(*t)();
00445 
00446   // t2 = inv(R) * t1(1:m)
00447   V_InvMtV( &t2, R, BLAS_Cpp::no_trans, t1(1,mb) );
00448 
00449   // t3 = - inv(R') * t1(m+1,2*m)
00450   V_mV( &t3, t1(mb+1,2*mb) );
00451   V_InvMtV( &t3, R, BLAS_Cpp::trans, t3 );
00452 
00453   // t4 = gk * Y'Y * t2
00454   V_StMtV( &t4, gamma_k_, YTY, BLAS_Cpp::no_trans, t2 );
00455 
00456   // t4 += D*t2
00457   Vp_DtV( &t4, t2 );
00458 
00459   // t5 = inv(R') * t4
00460   V_InvMtV( &t5, R, BLAS_Cpp::trans, t4 );
00461 
00462   // t5 += t3
00463   Vp_V( &t5, t3 );
00464 
00465   // y += S*t5
00466   (vdme(*t)() = t5);
00467   Vp_MtV( y, *S, BLAS_Cpp::no_trans, *t );
00468 
00469   // y += -gk*Y*t2
00470   (vdme(*t)() = t2);
00471   Vp_StMtV( y, -gamma_k_, *Y, BLAS_Cpp::no_trans, *t );
00472 
00473 }
00474 
00475 // Overridden from MatrixSymSecant
00476 
00477 void MatrixSymPosDefLBFGS::init_identity( const VectorSpace& space_diag, value_type alpha )
00478 {
00479   // Validate input
00480   TEUCHOS_TEST_FOR_EXCEPTION(
00481     alpha <= 0.0, std::invalid_argument
00482     ,"MatrixSymPosDefLBFGS::init_identity(n,alpha) : Error, "
00483     "alpha = " << alpha << " <= 0 is not allowed!" );
00484   
00485   // Set the vector space
00486   vec_spc_ = space_diag.clone();
00487   vec_spc_.get();
00488 
00489   // Set storage
00490   S_ = vec_spc_->create_members(m_);
00491   Y_ = vec_spc_->create_members(m_);
00492   TEUCHOS_TEST_FOR_EXCEPT( !( S_.get() ) );
00493   TEUCHOS_TEST_FOR_EXCEPT( !( Y_.get() ) );
00494   STY_.resize( m_, m_ );
00495   STSYTY_.resize( m_+1, m_+1 );
00496   STSYTY_.diag(0) = 0.0;
00497 
00498   gamma_k_ = 1.0/alpha;
00499 
00500   // Initialize counters
00501   m_bar_  = 0;
00502 
00503   n_ = vec_spc_->dim();        // initialized;
00504   original_is_updated_ = true; // This will never change for now
00505   inverse_is_updated_  = true; // This will never change for now
00506   num_secant_updates_  = 0;    // reset this to zero
00507 }
00508 
00509 void MatrixSymPosDefLBFGS::init_diagonal(const Vector& diag)
00510 {
00511   init_identity( diag.space(), diag.norm_inf() );
00512 }
00513 
00514 void MatrixSymPosDefLBFGS::secant_update(
00515   VectorMutable* s, VectorMutable* y, VectorMutable* Bs
00516   )
00517 {
00518   using AbstractLinAlgPack::BFGS_sTy_suff_p_d;
00519   using AbstractLinAlgPack::dot;
00520   using LinAlgOpPack::V_MtV;
00521   using Teuchos::Workspace;
00522   Teuchos::WorkspaceStore* wss = Teuchos::get_default_workspace_store().get();
00523 
00524   assert_initialized();
00525 
00526   // Check skipping the BFGS update
00527   const value_type
00528     sTy       = dot(*s,*y);
00529   std::ostringstream omsg;
00530   if( !BFGS_sTy_suff_p_d(*s,*y,&sTy,&omsg,"MatrixSymPosDefLBFGS::secant_update(...)" ) ) {
00531     throw UpdateSkippedException( omsg.str() ); 
00532   }
00533 
00534   try {
00535 
00536   // Update counters
00537   if( m_bar_ == m_ ) {
00538     // We are at the end of the storage so remove the oldest stored update
00539     // and move updates to make room for the new update.  This has to be done for the
00540     // the matrix to behave properly
00541     {for( size_type k = 1; k <= m_-1; ++k ) {
00542       S_->col(k) = S_->col(k+1);                            // Shift S.col() to the left
00543       Y_->col(k) = Y_->col(k+1);                            // Shift Y.col() to the left
00544       STY_.col(k)(1,m_-1) = STY_.col(k+1)(2,m_);            // Move submatrix STY(2,m-1,2,m-1) up and left
00545       STSYTY_.col(k)(k+1,m_) = STSYTY_.col(k+1)(k+2,m_+1);  // Move triangular submatrix STS(2,m-1,2,m-1) up and left
00546       STSYTY_.col(k+1)(1,k) = STSYTY_.col(k+2)(2,k+1);      // Move triangular submatrix YTY(2,m-1,2,m-1) up and left
00547     }}
00548     // ToDo: Create an abstract interface, call it MultiVectorShiftVecs, to rearrange S and Y all at once.
00549     // This will be important for parallel performance.
00550   }
00551   else {
00552     m_bar_++;
00553   }
00554   // Set the update vectors
00555   *S_->col(m_bar_) = *s;
00556   *Y_->col(m_bar_) = *y;
00557 
00558   // /////////////////////////////////////////////////////////////////////////////////////
00559   // Update S'Y
00560   //
00561   // Update the row and column m_bar
00562   //
00563   //  S'Y = 
00564   //
00565   //  [ s(1)'*y(1)    ...   s(1)'*y(m_bar)    ...   s(1)'*y(m_bar)    ]
00566   //  [ .           .             .         ] row
00567   //  [ s(m_bar)'*y(1)  ...   s(m_bar)'*y(m_bar)  ...   s(m_bar)'*y(m_bar)  ] m_bar
00568   //  [ .           .             .         ]
00569   //  [ s(m_bar)'*y(1)  ...   s(m_bar)'*y(m_bar)  ...   s(m_bar)'*y(m_bar)  ]
00570   //
00571   //                col m_bar
00572   //
00573   // Therefore we set:
00574   //  (S'Y)(:,m_bar) =  S'*y(m_bar)
00575   //  (S'Y)(m_bar,:) =  s(m_bar)'*Y
00576 
00577   const multi_vec_ptr_t
00578     S = this->S(),
00579     Y = this->Y();
00580 
00581   VectorSpace::vec_mut_ptr_t
00582     t = S->space_rows().create_member();  // temporary workspace
00583 
00584   //  (S'Y)(:,m_bar) =  S'*y(m_bar)
00585   V_MtV( t.get(), *S, BLAS_Cpp::trans, *y );
00586   STY_.col(m_bar_)(1,m_bar_) = VectorDenseEncap(*t)();
00587 
00588   //  (S'Y)(m_bar,:)' =  Y'*s(m_bar)
00589   V_MtV( t.get(), *Y, BLAS_Cpp::trans, *s );
00590   STY_.row(m_bar_)(1,m_bar_) = VectorDenseEncap(*t)();
00591 
00592   // /////////////////////////////////////////////////////////////////
00593   // Update S'S
00594   //
00595   //  S'S = 
00596   //
00597   //  [ s(1)'*s(1)    ...   symmetric         symmetric     ]
00598   //  [ .           .             .         ] row
00599   //  [ s(m_bar)'*s(1)  ...   s(m_bar)'*s(m_bar)  ...   symmetric     ] m_bar
00600   //  [ .           .             .         ]
00601   //  [ s(m_bar)'*s(1)  ...   s(m_bar)'*s(m_bar)  ...   s(m_bar)'*s(m_bar)  ]
00602   //
00603   //                col m_bar
00604   //
00605   // Here we will update the lower triangular part of S'S.  To do this we
00606   // only need to compute:
00607   //    t = S'*s(m_bar) = { s(m_bar)' * [ s(1),..,s(m_bar),..,s(m_bar) ]  }'
00608   // then set the appropriate rows and columns of S'S.
00609 
00610   Workspace<value_type>   work_ws(wss,m_bar_);
00611   DVectorSlice                  work(&work_ws[0],work_ws.size());
00612 
00613   // work = S'*s(m_bar)
00614   V_MtV( t.get(), *S, BLAS_Cpp::trans, *s );
00615   work = VectorDenseEncap(*t)();
00616 
00617   // Set row elements
00618   STSYTY_.row(m_bar_+1)(1,m_bar_) = work;
00619   // Set column elements
00620   STSYTY_.col(m_bar_)(m_bar_+1,m_bar_+1) = work(m_bar_,m_bar_);
00621 
00622   // /////////////////////////////////////////////////////////////////////////////////////
00623   // Update Y'Y
00624   //
00625   // Update the row and column m_bar
00626   //
00627   //  Y'Y = 
00628   //
00629   //  [ y(1)'*y(1)    ...   y(1)'*y(m_bar)    ...   y(1)'*y(m_bar)    ]
00630   //  [ .           .             .         ] row
00631   //  [ symmetric   ...   y(m_bar)'*y(m_bar)  ...   y(m_bar)'*y(m_bar)  ] m_bar
00632   //  [ .           .             .         ]
00633   //  [ symmetric   ...   symmetric     ...   y(m_bar)'*y(m_bar)  ]
00634   //
00635   //                col m_bar
00636   //
00637   // Here we will update the upper triangular part of Y'Y.  To do this we
00638   // only need to compute:
00639   //    t = Y'*y(m_bar) = { y(m_bar)' * [ y(1),..,y(m_bar),..,y(m_bar) ]  }'
00640   // then set the appropriate rows and columns of Y'Y.
00641 
00642   // work = Y'*y(m_bar)
00643   V_MtV( t.get(), *Y, BLAS_Cpp::trans, *y );
00644   work = VectorDenseEncap(*t)();
00645 
00646   // Set row elements
00647   STSYTY_.col(m_bar_+1)(1,m_bar_) = work;
00648   // Set column elements
00649   STSYTY_.row(m_bar_)(m_bar_+1,m_bar_+1) = work(m_bar_,m_bar_);
00650 
00651   // /////////////////////////////
00652   // Update gamma_k
00653 
00654   // gamma_k = s'*y / y'*y
00655   if(auto_rescaling_)
00656     gamma_k_ = STY_(m_bar_,m_bar_) / STSYTY_(m_bar_,m_bar_+1);
00657 
00658   // We do not initially update Q unless we have to form a matrix-vector
00659   // product later.
00660   
00661   Q_updated_ = false;
00662   num_secant_updates_++;
00663 
00664   } //  end try
00665   catch(...) {
00666     // If we throw any exception the we should make the matrix uninitialized
00667     // so that we do not leave this object in an inconsistant state.
00668     n_ = 0;
00669     throw;
00670   }
00671 
00672 }
00673 
00674 // Private member functions
00675 
00676 void MatrixSymPosDefLBFGS::Vp_DtV( DVectorSlice* y, const DVectorSlice& x ) const
00677 {
00678   DenseLinAlgPack::Vp_MtV_assert_sizes(
00679     y->dim(), m_bar_, m_bar_, BLAS_Cpp::no_trans, x.dim() );
00680 
00681   DVectorSlice::const_iterator
00682     d_itr = STY_.diag(0).begin(),
00683     x_itr = x.begin();
00684   DVectorSlice::iterator
00685     y_itr = y->begin();
00686 
00687   while( y_itr != y->end() )
00688     *y_itr++ += (*d_itr++) * (*x_itr++);    
00689 }
00690 
00691 //
00692 // We need to update the factorizations to solve for:
00693 //
00694 // x = inv(Q) * y   =>   Q * x = y
00695 //
00696 //  [ (1/gk)*S'S   L  ] * [ x1 ] = [ y1 ]
00697 //  [      L'   -D  ]   [ x2 ]   [ y2 ]
00698 //
00699 // We will solve the above system using a schur complement:
00700 //
00701 // C = (1/gk)*S'S + L*inv(D)*L'
00702 //
00703 // According to the referenced paper, C is p.d. so:
00704 //
00705 // C = J*J'
00706 //
00707 // We then compute the solution as:
00708 //
00709 // x1 = inv(C) * ( y1 + L*inv(D)*y2 )
00710 // x2 = - inv(D) * ( y2 - L'*x1 )
00711 //
00712 // Therefore we will just update the factorization C = J*J'
00713 // where the factor J is stored in QJ_.
00714 //
00715 
00716 void MatrixSymPosDefLBFGS::update_Q() const
00717 {
00718   using DenseLinAlgPack::tri;
00719   using DenseLinAlgPack::tri_ele;
00720   using DenseLinAlgPack::Mp_StM;
00721 
00722   //
00723   // We need update the factorizations to solve for:
00724   //
00725   // x = inv(Q) * y
00726   //
00727   //  [ y1 ]  = [ (1/gk)*S'S   L  ] * [ x1 ]
00728   //  [ y2 ]    [      L'   -D  ]   [ x2 ]
00729   //
00730   // We will solve the above system using the schur complement:
00731   //
00732   // C = (1/gk)*S'S + L*inv(D)*L'
00733   //
00734   // According to the referenced paper, C is p.d. so:
00735   //
00736   // C = J*J'
00737   //
00738   // We then compute the solution as:
00739   //
00740   // x1 = inv(C) * ( y1 + L*inv(D)*y2 )
00741   // x2 = - inv(D) * ( y2 - L'*x1 )
00742   //
00743   // Therefore we will just update the factorization C = J*J'
00744   //
00745 
00746   // Form the upper triangular part of C which will become J
00747   // which we are using storage of QJ
00748 
00749   if( QJ_.rows() < m_ )
00750     QJ_.resize( m_, m_ );
00751 
00752   const size_type
00753     mb = m_bar_;
00754 
00755   DMatrixSlice
00756     C = QJ_(1,mb,1,mb);
00757 
00758   // C = L * inv(D) * L'
00759   //
00760   // Here L is a strictly lower triangular (zero diagonal) matrix where:
00761   //
00762   // L = [ 0  0 ]
00763   //     [ Lb 0 ]
00764   //
00765   // Lb is lower triangular (nonzero diagonal)
00766   //
00767   // Therefore we compute L*inv(D)*L' as:
00768   //
00769   // C = [ 0  0 ] * [ Db  0  ] * [ 0  Lb' ]
00770   //     [ Lb 0 ]   [ 0   db ]   [ 0   0  ]
00771   //
00772   //   = [ 0  0  ] = [ 0      0     ]
00773   //     [ 0  Cb ]   [ 0  Lb*Db*Lb' ]
00774   //
00775   // We need to compute the upper triangular part of Cb.
00776 
00777   C.row(1) = 0.0;
00778   if( mb > 1 )
00779     comp_Cb( STY_(2,mb,1,mb-1), STY_.diag(0)(1,mb-1), &C(2,mb,2,mb) );
00780 
00781   // C += (1/gk)*S'S
00782 
00783   const DMatrixSliceSym &STS = this->STS();
00784   Mp_StM( &C, (1/gamma_k_), tri( STS.gms(), STS.uplo(), BLAS_Cpp::nonunit )
00785     , BLAS_Cpp::trans );
00786 
00787   // Now perform a cholesky factorization of C
00788   // After this factorization the upper triangular part of QJ
00789   // (through C) will contain the cholesky factor.
00790 
00791   DMatrixSliceTriEle C_upper = tri_ele( C, BLAS_Cpp::upper );
00792   try {
00793     DenseLinAlgLAPack::potrf( &C_upper );
00794   }
00795   catch( const DenseLinAlgLAPack::FactorizationException &fe ) {
00796     TEUCHOS_TEST_FOR_EXCEPTION(
00797       true, UpdateFailedException
00798       ,"Error, the factorization of Q which should be s.p.d. failed with"
00799       " the error message: {" << fe.what() << "}";
00800       );
00801   }
00802 
00803   Q_updated_ = true;
00804 }
00805 
00806 void MatrixSymPosDefLBFGS::V_invQtV( DVectorSlice* x, const DVectorSlice& y ) const
00807 {
00808   using DenseLinAlgPack::sym;
00809   using DenseLinAlgPack::tri;
00810   using DenseLinAlgPack::Vp_StV;
00811   using DenseLinAlgPack::V_InvMtV;
00812 
00813   using LinAlgOpPack::Vp_V;
00814   using LinAlgOpPack::V_MtV;
00815 
00816 
00817   // Solve the system 
00818   //
00819   // Q * x = y
00820   //
00821   // Using the schur complement factorization as described above.
00822 
00823   const size_type
00824     mb = m_bar_;
00825 
00826   if(!Q_updated_) {
00827     update_Q();
00828   }
00829 
00830   DVectorSlice
00831     x1 = (*x)(1,mb),
00832     x2 = (*x)(mb+1,2*mb);
00833 
00834   const DVectorSlice
00835     y1 = y(1,mb),
00836     y2 = y(mb+1,2*mb);
00837 
00838   // //////////////////////////////////////
00839   // x1 = inv(C) * ( y1 + L*inv(D)*y2 )
00840   //     = inv(J'*J) * r
00841   //     = inv(J) * inv(J') * r
00842 
00843   { // x1 = inv(D) * y2
00844     DVectorSlice::const_iterator
00845       d_itr = STY_.diag(0).begin(),
00846       y2_itr = y2.begin();
00847     DVectorSlice::iterator
00848       x1_itr = x1.begin();
00849     while( x1_itr != x1.end() )
00850       *x1_itr++ = *y2_itr++ / *d_itr++;
00851   }
00852 
00853   // x1 = L * x1
00854   //
00855   //    = [ 0  0 ] * [ x1(1:mb-1) ]
00856   //      [ Lb 0 ]   [ x1(mb)     ]
00857   //
00858   //    = [ 0             ]
00859   //      [ Lb*x1(1:mb-1) ]
00860   //
00861   if( mb > 1 ) {
00862     // x1(2,mb) = x1(1,mb-1) ( copy from mb-1 to mb then mb-2 to mb-1
00863     // etc. so that we don't overwrite the elements we need to copy ).
00864     DVectorSlice
00865       x1a = x1(1,mb-1),
00866       x1b = x1(2,mb);
00867     std::copy( x1a.rbegin(), x1a.rend(), x1b.rbegin() );
00868     V_MtV( &x1b, Lb(), BLAS_Cpp::no_trans, x1b );
00869   }
00870   x1(1) = 0.0;
00871 
00872   // r = x1 += y1
00873   Vp_V( &x1, y1 );
00874 
00875   // x1 = inv(J') * r
00876   const DMatrixSliceTri J = tri( QJ_(1,mb,1,mb), BLAS_Cpp::upper, BLAS_Cpp::nonunit );
00877   V_InvMtV( &x1, J, BLAS_Cpp::trans, x1 );
00878 
00879   // x1 = inv(J) * x1
00880   V_InvMtV( &x1, J, BLAS_Cpp::no_trans, x1 );
00881 
00882   // /////////////////////////////////////
00883   // x2 = inv(D) * ( - y2 + L'*x1 )
00884 
00885   // x2 = L'*x1
00886   //
00887   //    = [ 0  Lb' ] * [ x1(1)    ]
00888   //      [ 0  0   ]   [ x1(2,mb) ]
00889   //
00890   //    = [ Lb'*x1(2,mb) ]
00891   //      [      0       ]
00892   //
00893   if( mb > 1 ) {
00894     V_MtV( &x2(1,mb-1), Lb(), BLAS_Cpp::trans, x1(2,mb) );
00895   }
00896   x2(mb) = 0.0;
00897 
00898   // x2 += -y2
00899   Vp_StV( &x2, -1.0, y2 );
00900 
00901   // x2 = inv(D) * x2
00902   {
00903     DVectorSlice::const_iterator
00904       d_itr = STY_.diag(0).begin();
00905     DVectorSlice::iterator
00906       x2_itr = x2.begin();
00907     while( x2_itr != x2.end() )
00908       *x2_itr++ /= *d_itr++;
00909   }
00910 }
00911 
00912 void MatrixSymPosDefLBFGS::assert_initialized() const
00913 {
00914   if(!n_)
00915     throw std::logic_error( "MatrixSymPosDefLBFGS::assert_initialized() : "
00916       "Error, matrix not initialized" );
00917 }
00918 
00919 } // end namespace ConstrainedOptPack 
00920 
00921 namespace {
00922 
00923 void comp_Cb(
00924   const DMatrixSlice& Lb, const DVectorSlice& Db_diag
00925   ,DMatrixSlice* Cb
00926   )
00927 {
00928   // Lb * inv(Db) * Lb =
00929   //
00930   // [ l(1,1)           ]   [ dd(1)         ]   [ l(1,1)  l(2,1)  ... l(p,1)  ]
00931   // [ l(2,1) l(2,2)        ]   [   dd(2)     ]   [     l(2,2)  ... l(p,2)  ]
00932   // [ .    .        .      ] * [     .     ] * [         . .   ]
00933   // [ l(p,1) l(p,2)  ... l(p,p)  ]   [       dd(p) ]   [           l(p,p)  ]
00934   //
00935   //
00936   // [ l(1,1)*dd(1)*l(1,1)  l(1,1)*dd(1)*l(2,1)             ...   l(1,1)*dd(1)*l(p,1)             ]
00937   // [ symmetric        l(2,1)*dd(1)*l(2,1) + l(2,2)*dd(2)*l(2,2) ...   l(2,1)*dd(1)*l(p,1) + l(2,2)*dd(2)*l(p,2) ]   
00938   // [ .            .                     ...
00939   // [ symmetric        symmetric                 ...   sum( l(p,i)*dd(i)*l(p,i), i=1,..,p )    ]
00940   //
00941   // Therefore we can express the upper triangular elemetns of Cb as:
00942   //
00943   // Cb(i,j) = sum( l(i,k)*dd(k)*l(j,k), k = 1,..,i )
00944 
00945   typedef DenseLinAlgPack::size_type size_type;
00946   typedef DenseLinAlgPack::value_type value_type;
00947 
00948   TEUCHOS_TEST_FOR_EXCEPT( !(  Lb.rows() == Cb->rows() && Cb->rows() == Db_diag.dim()  ) ); // only a local error!
00949 
00950   const size_type p = Db_diag.dim();
00951 
00952   for( size_type i = 1; i <= p; ++i ) {
00953     for( size_type j = i; j <= p; ++j ) {
00954       value_type &c = (*Cb)(i,j) = 0.0;
00955       for( size_type k = 1; k <= i; ++k ) {
00956         c += Lb(i,k) * Lb(j,k) / Db_diag(k);
00957       }
00958     }
00959   }
00960 
00961   // ToDo: Make the above operation more efficent if needed! (i.e. write
00962   // it in fortran or something?).
00963 }
00964 
00965 } // end namespace
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Defines