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

Generated on Wed May 12 21:52:28 2010 for MOOCHO (Single Doxygen Collection) by  doxygen 1.4.7