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