ConstrainedOptPack_MatrixHessianRelaxed.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 // 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 <assert.h>
00030 
00031 #include "ConstrainedOptPack_MatrixHessianRelaxed.hpp"
00032 #include "AbstractLinAlgPack/src/AbstractLinAlgPack_MatrixSymOp.hpp"
00033 #include "AbstractLinAlgPack/src/AbstractLinAlgPack_GenPermMatrixSlice.hpp"
00034 #include "AbstractLinAlgPack/src/AbstractLinAlgPack_SpVectorClass.hpp"
00035 #include "AbstractLinAlgPack_SpVectorOp.hpp"
00036 #include "DenseLinAlgPack_LinAlgOpPack.hpp"
00037 
00038 namespace LinAlgOpPack {
00039   using AbstractLinAlgPack::Vp_StV;
00040   using AbstractLinAlgPack::Vp_StMtV;
00041 }
00042 
00043 namespace ConstrainedOptPack {
00044 
00045 MatrixHessianRelaxed::MatrixHessianRelaxed()
00046   :
00047     n_(0)
00048     ,H_(NULL)
00049     ,bigM_(0.0)
00050 {}
00051 
00052 void MatrixHessianRelaxed::initialize(
00053     const MatrixSymOp &H
00054   , value_type      bigM
00055   )
00056 {
00057   n_  = H.rows();
00058   H_  = &H;
00059   bigM_ = bigM;
00060 }
00061 
00062 // Overridden from Matrix
00063 
00064 size_type MatrixHessianRelaxed::rows() const
00065 {
00066   return n_ ? n_ + 1 : 0;
00067 }
00068 
00069 // Overridden from MatrixOp
00070 
00071 void MatrixHessianRelaxed::Vp_StMtV(
00072     DVectorSlice* y, value_type a, BLAS_Cpp::Transp M_trans
00073   , const DVectorSlice& x, value_type b ) const
00074 {
00075   using BLAS_Cpp::no_trans;
00076   using BLAS_Cpp::trans;
00077   using AbstractLinAlgPack::Vp_StMtV;
00078   //
00079   // y = b*y + a * M * x
00080   // 
00081   //   = b*y + a * [ H  0    ] * [ x1 ]
00082   //               [ 0  bigM ]   [ x2 ]
00083   //               
00084   // =>              
00085   //               
00086   // y1 = b*y1 + a*H*x1
00087   // 
00088   // y2 = b*y2 + bigM * x2
00089   //
00090   LinAlgOpPack::Vp_MtV_assert_sizes(y->size(),rows(),cols(),M_trans,x.size());
00091 
00092   DVectorSlice
00093     y1 = (*y)(1,n_);
00094   value_type
00095     &y2 = (*y)(n_+1);
00096   const DVectorSlice
00097     x1 = x(1,n_);
00098   const value_type
00099     x2 = x(n_+1);
00100 
00101   // y1 = b*y1 + a*H*x1
00102   Vp_StMtV( &y1, a, *H_, no_trans, x1, b );
00103 
00104   // y2 = b*y2 + bigM * x2
00105   if( b == 0.0 )
00106     y2 = bigM_ * x2;
00107   else
00108     y2 = b*y2 + bigM_ * x2;
00109   
00110 }
00111 
00112 void MatrixHessianRelaxed::Vp_StMtV(
00113     DVectorSlice* y, value_type a, BLAS_Cpp::Transp M_trans
00114   , const SpVectorSlice& x, value_type b ) const
00115 {
00116   using BLAS_Cpp::no_trans;
00117   using BLAS_Cpp::trans;
00118   using AbstractLinAlgPack::Vp_StMtV;
00119   //
00120   // y = b*y + a * M * x
00121   // 
00122   //   = b*y + a * [ H  0    ] * [ x1 ]
00123   //               [ 0  bigM ]   [ x2 ]
00124   //               
00125   // =>              
00126   //               
00127   // y1 = b*y1 + a*H*x1
00128   // 
00129   // y2 = b*y2 + bigM * x2
00130   //
00131   LinAlgOpPack::Vp_MtV_assert_sizes(y->size(),rows(),cols(),M_trans,x.size());
00132 
00133   DVectorSlice
00134     y1 = (*y)(1,n_);
00135   value_type
00136     &y2 = (*y)(n_+1);
00137   const SpVectorSlice
00138     x1 = x(1,n_);
00139   const SpVectorSlice::element_type
00140     *x2_ele = x.lookup_element(n_+1);
00141   const value_type
00142     x2 = x2_ele ? x2_ele->value() : 0.0;
00143 
00144   // y1 = b*y1 + a*H*x1
00145   Vp_StMtV( &y1, a, *H_, no_trans, x1, b );
00146 
00147   // y2 = b*y2 + bigM * x2
00148   if( b == 0.0 )
00149     y2 = bigM_ * x2;
00150   else
00151     y2 = b*y2 + bigM_ * x2;
00152   
00153 }
00154 
00155 void MatrixHessianRelaxed::Vp_StPtMtV(
00156   DVectorSlice* y, value_type a
00157   , const GenPermMatrixSlice& P, BLAS_Cpp::Transp P_trans
00158   , BLAS_Cpp::Transp M_trans
00159   , const DVectorSlice& x, value_type b ) const
00160 {
00161   using BLAS_Cpp::no_trans;
00162   using BLAS_Cpp::trans;
00163   namespace GPMSIP = AbstractLinAlgPack::GenPermMatrixSliceIteratorPack;
00164   //
00165   // y = b*y + a * op(P) * M * x
00166   // 
00167   //   = b*y + a * [ op(P1)  op(P2) ] *  [ H   0   ] * [ x1 ]
00168   //                                     [ 0  bigM ]   [ x2 ]
00169   //               
00170   // =>              
00171   //               
00172   // y = b*y + a*op(P1)*H*x1 + a*op(P2)*bigM*x2
00173   //
00174   LinAlgOpPack::Vp_MtV_assert_sizes(y->size(),P.rows(),P.cols(),P_trans
00175     , BLAS_Cpp::rows( rows(), cols(), M_trans) );
00176   LinAlgOpPack::Vp_MtV_assert_sizes( BLAS_Cpp::cols( P.rows(), P.cols(), P_trans)
00177     ,rows(),cols(),M_trans,x.size());
00178 
00179   // For this to work (as shown below) we need to have P sorted by
00180   // row if op(P) = P' or sorted by column if op(P) = P.  If
00181   // P is not sorted properly, we will just use the default
00182   // implementation of this operation.
00183   if(   ( P.ordered_by() == GPMSIP::BY_ROW && P_trans == no_trans )
00184       ||  ( P.ordered_by() == GPMSIP::BY_COL && P_trans == trans ) )
00185   {
00186     // Call the default implementation
00187     MatrixOp::Vp_StPtMtV(y,a,P,P_trans,M_trans,x,b);
00188     return;
00189   }
00190 
00191   if( P.is_identity() )
00192     assert( BLAS_Cpp::rows( P.rows(), P.cols(), P_trans ) == n_ );
00193 
00194   const GenPermMatrixSlice
00195     P1 = ( P.is_identity() 
00196          ? GenPermMatrixSlice( n_, n_, GenPermMatrixSlice::IDENTITY_MATRIX )
00197          : P.create_submatrix(Range1D(1,n_),P_trans==trans?GPMSIP::BY_ROW:GPMSIP::BY_COL)
00198       ),
00199     P2 = ( P.is_identity()
00200          ? GenPermMatrixSlice(
00201            P_trans == no_trans ? n_ : 1
00202            , P_trans == no_trans ? 1 : n_
00203            , GenPermMatrixSlice::ZERO_MATRIX )
00204          : P.create_submatrix(Range1D(n_+1,n_+1),P_trans==trans?GPMSIP::BY_ROW:GPMSIP::BY_COL)
00205       );
00206   
00207   const DVectorSlice
00208     x1 = x(1,n_);
00209   const value_type
00210     x2 = x(n_+1);
00211   // y = b*y + a*op(P1)*H*x1
00212   AbstractLinAlgPack::Vp_StPtMtV( y, a, P1, P_trans, *H_, no_trans, x1, b );
00213   // y += a*op(P2)*bigM*x2
00214   if( P2.nz() ){
00215     assert(P2.nz() == 1);
00216     const size_type
00217       i = P_trans == no_trans ? P2.begin()->row_i() : P2.begin()->col_j();
00218     (*y)(i) += a * bigM_ * x2;
00219   }
00220 }
00221 
00222 void MatrixHessianRelaxed::Vp_StPtMtV(
00223   DVectorSlice* y, value_type a
00224   , const GenPermMatrixSlice& P, BLAS_Cpp::Transp P_trans
00225   , BLAS_Cpp::Transp M_trans
00226   , const SpVectorSlice& x, value_type b ) const
00227 {
00228   using BLAS_Cpp::no_trans;
00229   using BLAS_Cpp::trans;
00230   namespace GPMSIP = AbstractLinAlgPack::GenPermMatrixSliceIteratorPack;
00231   //
00232   // y = b*y + a * op(P) * M * x
00233   // 
00234   //   = b*y + a * [ op(P1)  op(P2) ] *  [ H   0   ] * [ x1 ]
00235   //                                     [ 0  bigM ]   [ x2 ]
00236   //               
00237   // =>              
00238   //               
00239   // y = b*y + a*op(P1)*H*x1 + a*op(P2)*bigM*x2
00240   //
00241   LinAlgOpPack::Vp_MtV_assert_sizes(y->size(),P.rows(),P.cols(),P_trans
00242     , BLAS_Cpp::rows( rows(), cols(), M_trans) );
00243   LinAlgOpPack::Vp_MtV_assert_sizes( BLAS_Cpp::cols( P.rows(), P.cols(), P_trans)
00244     ,rows(),cols(),M_trans,x.size());
00245 
00246   // For this to work (as shown below) we need to have P sorted by
00247   // row if op(P) = P' or sorted by column if op(P) = P.  If
00248   // P is not sorted properly, we will just use the default
00249   // implementation of this operation.
00250   if(   ( P.ordered_by() == GPMSIP::BY_ROW && P_trans == no_trans )
00251       ||  ( P.ordered_by() == GPMSIP::BY_COL && P_trans == trans ) )
00252   {
00253     // Call the default implementation
00254     MatrixOp::Vp_StPtMtV(y,a,P,P_trans,M_trans,x,b);
00255     return;
00256   }
00257 
00258   if( P.is_identity() )
00259     assert( BLAS_Cpp::rows( P.rows(), P.cols(), P_trans ) == n_ );
00260 
00261   const GenPermMatrixSlice
00262     P1 = ( P.is_identity() 
00263          ? GenPermMatrixSlice( n_, n_, GenPermMatrixSlice::IDENTITY_MATRIX )
00264          : P.create_submatrix(Range1D(1,n_),P_trans==trans?GPMSIP::BY_ROW:GPMSIP::BY_COL)
00265       ),
00266     P2 = ( P.is_identity()
00267          ? GenPermMatrixSlice(
00268            P_trans == no_trans ? n_ : 1
00269            , P_trans == no_trans ? 1 : n_
00270            , GenPermMatrixSlice::ZERO_MATRIX )
00271          : P.create_submatrix(Range1D(n_+1,n_+1),P_trans==trans?GPMSIP::BY_ROW:GPMSIP::BY_COL)
00272       );
00273 
00274   const SpVectorSlice
00275     x1 = x(1,n_);
00276   const SpVectorSlice::element_type
00277     *x2_ele = x.lookup_element(n_+1);
00278   const value_type
00279     x2 = x2_ele ? x2_ele->value() : 0.0;
00280   // y = b*y + a*op(P1)*H*x1
00281   AbstractLinAlgPack::Vp_StPtMtV( y, a, P1, P_trans, *H_, no_trans, x1, b );
00282   // y += a*op(P2)*bigM*x2
00283   if( P2.nz() ){
00284     assert(P2.nz() == 1);
00285     const size_type
00286       i = P_trans == no_trans ? P2.begin()->row_i() : P2.begin()->col_j();
00287     (*y)(i) += a * bigM_ * x2;
00288   }
00289 }
00290 
00291 value_type MatrixHessianRelaxed::transVtMtV(
00292   const SpVectorSlice& x1, BLAS_Cpp::Transp M_trans
00293   , const SpVectorSlice& x2 ) const
00294 {
00295   using BLAS_Cpp::no_trans;
00296   //
00297   // a = x1' * M * x2
00298   // 
00299   //   = [ x11'  x12' ] * [ H   0   ] * [ x21 ]
00300   //                      [ 0  bigM ]   [ x22 ]
00301   //               
00302   // =>              
00303   //               
00304   // a = x11'*H*x21 + x12'*bigM*x22
00305   //
00306   DenseLinAlgPack::Vp_MtV_assert_sizes(x1.size(),rows(),cols(),M_trans,x2.size());
00307 
00308   if( &x1 == &x2 ) {
00309     // x1 and x2 are the same sparse vector
00310     const SpVectorSlice
00311       x11 = x1(1,n_);
00312     const SpVectorSlice::element_type
00313       *x12_ele = x1.lookup_element(n_+1);
00314     const value_type
00315       x12 = x12_ele ? x12_ele->value() : 0.0;
00316     return AbstractLinAlgPack::transVtMtV( x11, *H_, no_trans, x11) 
00317       + x12 * bigM_ * x12;
00318   }
00319   else {
00320     // x1 and x2 could be different sparse vectors
00321     TEST_FOR_EXCEPT(true); // ToDo: Implement this!
00322   }
00323   return 0.0; // Will never be executed!
00324 }
00325 
00326 } // end namespace ConstrainedOptPack

Generated on Thu Sep 18 12:35:14 2008 for MOOCHO (Single Doxygen Collection) by doxygen 1.3.9.1