MOOCHO (Single Doxygen Collection) Version of the Day
ConstrainedOptPack_MatrixSymHessianRelaxNonSing.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 #include <assert.h>
00043 
00044 #include "ConstrainedOptPack_MatrixSymHessianRelaxNonSing.hpp"
00045 #include "AbstractLinAlgPack_VectorDenseEncap.hpp"
00046 #include "AbstractLinAlgPack_LinAlgOpPackHack.hpp"
00047 #include "AbstractLinAlgPack_SpVectorClass.hpp"
00048 #include "AbstractLinAlgPack_GenPermMatrixSlice.hpp"
00049 #include "AbstractLinAlgPack_VectorSpaceBlocked.hpp"
00050 #include "AbstractLinAlgPack_LinAlgOpPack.hpp"
00051 #include "DenseLinAlgPack_AssertOp.hpp"
00052 #include "DenseLinAlgPack_LinAlgOpPack.hpp"
00053 #include "ProfileHackPack_profile_hack.hpp"
00054 #include "Teuchos_Assert.hpp"
00055 
00056 namespace {
00057 
00058 //
00059 template<class V>
00060 void Vp_StPtMtV_imp( 
00061   DenseLinAlgPack::DVectorSlice* y, DenseLinAlgPack::value_type a
00062   ,const AbstractLinAlgPack::GenPermMatrixSlice& P, BLAS_Cpp::Transp P_trans
00063   ,const ConstrainedOptPack::MatrixSymHessianRelaxNonSing& H, BLAS_Cpp::Transp H_trans
00064   ,const V& x, DenseLinAlgPack::value_type b
00065   )
00066 {
00067   using BLAS_Cpp::no_trans;
00068   using BLAS_Cpp::trans;
00069   using BLAS_Cpp::trans_not;
00070   using AbstractLinAlgPack::Vp_StMtV;
00071   using AbstractLinAlgPack::Vp_StPtMtV;
00072   using AbstractLinAlgPack::GenPermMatrixSlice;
00073   using AbstractLinAlgPack::MatrixOp;
00074   namespace GPMSIP = AbstractLinAlgPack::GenPermMatrixSliceIteratorPack;
00075 
00076 #ifdef TEUCHOS_DEBUG
00077   TEUCHOS_TEST_FOR_EXCEPT(y==NULL);
00078 #endif
00079 
00080   const DenseLinAlgPack::size_type
00081     no = H.G().rows(),  // number of original variables
00082     nr = H.M().rows(),  // number of relaxation variables
00083     nd = no + nr,       // total number of variables
00084     ydim = y->dim();    // y->dim() == number of rows in op(P)
00085 
00086   DenseLinAlgPack::Vp_MtV_assert_sizes(
00087     y->dim(),P.rows(),P.cols(),P_trans
00088     ,BLAS_Cpp::rows( nd, nd, H_trans) );
00089   DenseLinAlgPack::Vp_MtV_assert_sizes(
00090     BLAS_Cpp::cols( P.rows(), P.cols(), P_trans)
00091     ,nd, nd, H_trans, x.dim() );
00092 
00093   //
00094   // y = b*y + a * op(P) * H * x
00095   //
00096   // y = b*y + a * [op(P1)  op(P2) ] * [ G  0 ] * [ x1 ]
00097   //                                   [ 0  M ]   [ x2 ]
00098   //
00099   // =>
00100   //
00101   // y = b*y + a*op(P1)*G*x1 + a*op(P2)*H*x2
00102   //
00103   // For this to work op(P) must be sorted by column.
00104   //
00105   if(  ( P.ordered_by() == GPMSIP::BY_ROW && P_trans == BLAS_Cpp::no_trans )
00106        || ( P.ordered_by() == GPMSIP::BY_COL && P_trans == BLAS_Cpp::trans )
00107        ||  ( P.ordered_by() == GPMSIP::UNORDERED ) )
00108   {
00109     // Call the default implementation
00110     //H.MatrixOp::Vp_StPtMtV(y,a,P,P_trans,H_trans,x,b);
00111     TEUCHOS_TEST_FOR_EXCEPT(true); // ToDo: Implement!
00112     return;
00113   }
00114   const DenseLinAlgPack::Range1D
00115     o_rng(1,no),
00116     r_rng(no+1,no+nr);
00117   const AbstractLinAlgPack::GenPermMatrixSlice
00118     P1 = ( P.is_identity() 
00119          ? GenPermMatrixSlice(
00120            P_trans == no_trans ? ydim : no 
00121            ,P_trans == no_trans ? no : ydim
00122            ,GenPermMatrixSlice::IDENTITY_MATRIX )
00123          : P.create_submatrix(o_rng,P_trans==trans?GPMSIP::BY_ROW:GPMSIP::BY_COL)
00124       ),
00125     P2 = ( P.is_identity()
00126          ? GenPermMatrixSlice(
00127            P_trans == no_trans ? ydim : nr
00128            ,P_trans == no_trans ? nr : ydim
00129            ,GenPermMatrixSlice::ZERO_MATRIX )
00130          : P.create_submatrix(r_rng,P_trans==trans?GPMSIP::BY_ROW:GPMSIP::BY_COL)
00131       );
00132   const V
00133     x1 = x(o_rng),
00134     x2 = x(r_rng);
00135   // y = b*y
00136   LinAlgOpPack::Vt_S(y,b);
00137   // y += a*op(P1)*G*x1
00138   if( P1.nz() )
00139     LinAlgOpPack::Vp_StPtMtV( y, a, P1, P_trans, H.G(), H_trans, x1, b );
00140   // y += a*op(P2)*M*x2
00141   if( P2.nz() )
00142     LinAlgOpPack::Vp_StPtMtV( y, a, P2, P_trans, H.M(), H_trans, x2, 1.0 );
00143 }
00144 
00145 } // end namespace
00146 
00147 namespace ConstrainedOptPack {
00148 
00149 MatrixSymHessianRelaxNonSing::MatrixSymHessianRelaxNonSing()
00150   : vec_space_(Teuchos::null)
00151 {}
00152 
00153 MatrixSymHessianRelaxNonSing::MatrixSymHessianRelaxNonSing(
00154   const G_ptr_t         &G_ptr
00155   ,const vec_mut_ptr_t  &M_diag_ptr
00156   ,const space_ptr_t    &space
00157   )
00158   : vec_space_(Teuchos::null)
00159 {
00160   initialize(G_ptr,M_diag_ptr,space);
00161 }
00162 
00163 void MatrixSymHessianRelaxNonSing::initialize(
00164   const G_ptr_t         &G_ptr
00165   ,const vec_mut_ptr_t  &M_diag_ptr
00166   ,const space_ptr_t    &space
00167   )
00168 {
00169   namespace mmp = MemMngPack;
00170 #ifdef TEUCHOS_DEBUG
00171   const char err_msg_head[] = "MatrixSymHessianRelaxNonSing::initialize(...) : Error!";
00172   TEUCHOS_TEST_FOR_EXCEPTION(G_ptr.get()==NULL, std::invalid_argument, err_msg_head);
00173   TEUCHOS_TEST_FOR_EXCEPTION(M_diag_ptr.get()==NULL, std::invalid_argument, err_msg_head);
00174   const size_type G_rows = G_ptr->rows(), M_diag_dim = M_diag_ptr->dim();
00175   TEUCHOS_TEST_FOR_EXCEPTION(G_rows==0, std::invalid_argument, err_msg_head);
00176   TEUCHOS_TEST_FOR_EXCEPTION(M_diag_dim==0, std::invalid_argument, err_msg_head);
00177 #endif
00178   if( space.get() ) {
00179 #ifdef TEUCHOS_DEBUG
00180     const size_type space_dim = space->dim();
00181     TEUCHOS_TEST_FOR_EXCEPTION(space_dim != G_rows + M_diag_dim, std::invalid_argument, err_msg_head);
00182 #endif
00183     vec_space_ = space;
00184   }
00185   else {
00186     VectorSpace::space_ptr_t spaces[]
00187       = { Teuchos::rcp(&G_ptr->space_cols(),false), Teuchos::rcp(&M_diag_ptr->space(),false) };
00188     vec_space_ = Teuchos::rcp(new VectorSpaceBlocked( spaces, 2 ) );
00189   }
00190   G_ptr_ = G_ptr;
00191   M_.initialize(M_diag_ptr);
00192 }
00193   
00194 // Overridden from MatrixOp
00195 
00196 const VectorSpace& MatrixSymHessianRelaxNonSing::space_cols() const
00197 {
00198   assert_initialized();
00199   return *vec_space_;
00200 }
00201 
00202 bool MatrixSymHessianRelaxNonSing::Mp_StM(
00203   MatrixOp* C, value_type a, BLAS_Cpp::Transp H_trans
00204   ) const
00205 {
00206 #ifdef PROFILE_HACK_ENABLED
00207   ProfileHackPack::ProfileTiming profile_timing( "MatrixSymHessianRelaxNonSing::Mp_StM(...)" );
00208 #endif
00209   assert_initialized();
00210   return MatrixOp::Mp_StM(C,a,H_trans); // ToDo: Update below code!
00211 /*
00212   const size_type
00213     nG = G_ptr_->rows(),
00214     nM = M_.rows();
00215   AbstractLinAlgPack::Mp_StM( &(*C)(1,nG,1,nG), a, *G_ptr_, H_trans);
00216   AbstractLinAlgPack::Mp_StM( &(*C)(nG+1,nG+nM,nG+1,nG+nM), a, M_, H_trans);
00217 */
00218 }
00219 
00220 void MatrixSymHessianRelaxNonSing::Vp_StMtV(
00221   VectorMutable* y, value_type a, BLAS_Cpp::Transp H_trans
00222   ,const Vector& x, value_type b
00223   ) const
00224 {
00225 #ifdef PROFILE_HACK_ENABLED
00226   ProfileHackPack::ProfileTiming profile_timing( "MatrixSymHessianRelaxNonSing::Vp_StMtV(...Vector...)" );
00227 #endif
00228   assert_initialized();
00229   const size_type
00230     nG = G_ptr_->rows(),
00231     nM = M_.rows();
00232   AbstractLinAlgPack::Vt_S(y,b);
00233   AbstractLinAlgPack::Vp_StMtV( y->sub_view(1,nG).get(), a, *G_ptr_, H_trans, *x.sub_view(1,nG) );
00234   AbstractLinAlgPack::Vp_StMtV( y->sub_view(nG+1,nG+nM).get(), a, M_, H_trans, *x.sub_view(nG+1,nG+nM) );
00235 }
00236 
00237 void MatrixSymHessianRelaxNonSing::Vp_StMtV(
00238   VectorMutable* y, value_type a, BLAS_Cpp::Transp H_trans
00239   ,const SpVectorSlice& x, value_type b
00240   ) const
00241 {
00242 #ifdef PROFILE_HACK_ENABLED
00243   ProfileHackPack::ProfileTiming profile_timing( "MatrixSymHessianRelaxNonSing::Vp_StMtV(...SpVectorSlice...)" );
00244 #endif
00245   assert_initialized();
00246   const size_type
00247     nG = G_ptr_->rows(),
00248     nM = M_.rows();
00249   AbstractLinAlgPack::Vt_S(y,b); // Takes care of b == 0.0 and y uninitialized
00250   AbstractLinAlgPack::Vp_StMtV( y->sub_view(1,nG).get(), a, *G_ptr_, H_trans, x(1,nG) );
00251   AbstractLinAlgPack::Vp_StMtV( y->sub_view(nG+1,nG+nM).get(), a, M_, H_trans, x(nG+1,nG+nM) );
00252 }
00253 
00254 void MatrixSymHessianRelaxNonSing::Vp_StPtMtV(
00255   VectorMutable* y, value_type a, const GenPermMatrixSlice& P, BLAS_Cpp::Transp P_trans
00256   ,BLAS_Cpp::Transp H_trans, const Vector& x, value_type b
00257   ) const
00258 {
00259 #ifdef PROFILE_HACK_ENABLED
00260   ProfileHackPack::ProfileTiming profile_timing( "MatrixSymHessianRelaxNonSing::Vp_StPtMtV(...Vector...)" );
00261 #endif
00262   assert_initialized();
00263   //MatrixOp::Vp_StPtMtV(y,a,P,P_trans,H_trans,x,b); // Uncomment for this default implementation
00264   AbstractLinAlgPack::VectorDenseMutableEncap y_d(*y);
00265   AbstractLinAlgPack::VectorDenseEncap x_d(x);
00266   Vp_StPtMtV_imp(&y_d(),a,P,P_trans,*this,H_trans,x_d(),b);
00267 }
00268 
00269 void MatrixSymHessianRelaxNonSing::Vp_StPtMtV(
00270   VectorMutable* y, value_type a, const GenPermMatrixSlice& P, BLAS_Cpp::Transp P_trans
00271   ,BLAS_Cpp::Transp H_trans, const SpVectorSlice& x, value_type b
00272   ) const
00273 {
00274 #ifdef PROFILE_HACK_ENABLED
00275   ProfileHackPack::ProfileTiming profile_timing( "MatrixSymHessianRelaxNonSing::Vp_StPtMtV(...SpVectorSlice...)" );
00276 #endif
00277   assert_initialized();
00278   //MatrixOp::Vp_StPtMtV(y,a,P,P_trans,H_trans,x,b); // Uncomment for this default implementation
00279   AbstractLinAlgPack::VectorDenseMutableEncap y_d(*y);
00280   Vp_StPtMtV_imp(&y_d(),a,P,P_trans,*this,H_trans,x,b);
00281 }
00282 
00283 // Overridden form MatrixSymOp
00284 
00285 void MatrixSymHessianRelaxNonSing::Mp_StPtMtP(
00286   MatrixSymOp* S, value_type a
00287   ,EMatRhsPlaceHolder dummy_place_holder
00288   ,const GenPermMatrixSlice& P, BLAS_Cpp::Transp P_trans
00289   ,value_type b
00290   ) const
00291 {
00292   using BLAS_Cpp::no_trans;
00293   using BLAS_Cpp::trans;
00294   using BLAS_Cpp::trans_not;
00295   namespace GPMSIP = AbstractLinAlgPack::GenPermMatrixSliceIteratorPack;
00296 #ifdef PROFILE_HACK_ENABLED
00297   ProfileHackPack::ProfileTiming profile_timing( "MatrixSymHessianRelaxNonSing::Mp_StPtMtP(...)" );
00298 #endif
00299   assert_initialized();
00300 
00301   MatrixSymOp::Mp_StPtMtP(S,a,dummy_place_holder,P,P_trans,b); // ToDo: Override when needed!
00302   return;
00303 /* ToDo: Update below code!
00304   const DenseLinAlgPack::size_type
00305     no = G().rows(),     // number of original variables
00306     nr = M().rows(),     // number of relaxation variables
00307     nd = no + nr;        // total number of variables
00308 
00309   DenseLinAlgPack::Mp_MtM_assert_sizes( S->rows(), S->cols(), no_trans
00310                    , P.rows(), P.cols(), trans_not(P_trans)
00311                    , P.rows(), P.cols(), P_trans );
00312   DenseLinAlgPack::Vp_V_assert_sizes( BLAS_Cpp::rows( P.rows(), P.cols(), P_trans), nd );
00313 
00314   //
00315   // S = b*S + a * op(P)' * H * op(P)
00316   //
00317   // S = b*S + a * [op(P1)'  op(P2)' ] * [ G  0 ] * [ op(P1) ]
00318   //                                     [ 0  M ]   [ op(P2) ]
00319   //
00320   // =>
00321   //
00322   // S = b*S
00323   // S1 += op(P1)' * G * op(P1)
00324   // S2 += op(P2)' * M * op(P2)
00325   //
00326   // For this to work op(P) must be sorted by row.
00327   //
00328   if(   ( P.ordered_by() == GPMSIP::BY_ROW && P_trans == BLAS_Cpp::trans )
00329       ||  ( P.ordered_by() == GPMSIP::BY_COL && P_trans == BLAS_Cpp::no_trans )
00330     ||  ( P.ordered_by() == GPMSIP::UNORDERED ) )
00331   {
00332     // Call the default implementation
00333     MatrixSymOp::Mp_StPtMtP(S,a,dummy_place_holder,P,P_trans,b);
00334     return;
00335   }
00336   const DenseLinAlgPack::Range1D
00337     o_rng(1,no),
00338     r_rng(no+1,no+nr);
00339   const AbstractLinAlgPack::GenPermMatrixSlice
00340     P1 = ( P.is_identity() 
00341          ? GenPermMatrixSlice(
00342            P_trans == no_trans ? nd : no 
00343            ,P_trans == no_trans ? no : nd
00344            ,GenPermMatrixSlice::IDENTITY_MATRIX )
00345          : P.create_submatrix(o_rng,P_trans==no_trans?GPMSIP::BY_ROW:GPMSIP::BY_COL)
00346       ),
00347     P2 = ( P.is_identity()
00348          ? GenPermMatrixSlice(
00349            P_trans == no_trans ? nd : nr
00350            ,P_trans == no_trans ? nr : nd
00351            ,GenPermMatrixSlice::ZERO_MATRIX )
00352          : P.create_submatrix(r_rng,P_trans==no_trans?GPMSIP::BY_ROW:GPMSIP::BY_COL)
00353       );
00354   // S = b*S
00355   DenseLinAlgPack::Mt_S( &DMatrixSliceTriEle(S->gms(),S->uplo()),b); // Handles b == 0.0 properly!
00356 
00357   // S1 += a*op(P1)'*G*op(P1)
00358   if( P1.nz() )
00359     AbstractLinAlgPack::Mp_StPtMtP(
00360       &DMatrixSliceSym( S->gms()(1,no,1,no), S->uplo() )
00361       , a, dummy_place_holder, G(), P1, P_trans );
00362   // S2 += a*op(P2)'*M*op(P2)
00363   if( P2.nz() )
00364     AbstractLinAlgPack::Mp_StPtMtP(
00365       &DMatrixSliceSym( S->gms()(no+1,nd,no+1,nd), S->uplo() )
00366       , a, dummy_place_holder, M(), P2, P_trans );
00367 */
00368 }
00369 
00370 // Overridden from MatrixOpNonsing
00371 
00372 void MatrixSymHessianRelaxNonSing::V_InvMtV(
00373   VectorMutable* y, BLAS_Cpp::Transp H_trans, const Vector& x
00374   ) const
00375 {
00376 #ifdef PROFILE_HACK_ENABLED
00377   ProfileHackPack::ProfileTiming profile_timing( "MatrixSymHessianRelaxNonSing::V_InvMtV(...Vector...)" );
00378 #endif
00379   assert_initialized();
00380   const size_type
00381     nG = G_ptr_->rows(),
00382     nM = M_.rows();
00383   AbstractLinAlgPack::V_InvMtV( y->sub_view(1,nG).get(), *G_ptr_, H_trans, *x.sub_view(1,nG) );
00384   AbstractLinAlgPack::V_InvMtV( y->sub_view(nG+1,nG+nM).get(), M_, H_trans, *x.sub_view(nG+1,nG+nM) );
00385 }
00386 
00387 void MatrixSymHessianRelaxNonSing::V_InvMtV(
00388   VectorMutable* y, BLAS_Cpp::Transp H_trans, const SpVectorSlice& x
00389   ) const
00390 {
00391 #ifdef PROFILE_HACK_ENABLED
00392   ProfileHackPack::ProfileTiming profile_timing( "MatrixSymHessianRelaxNonSing::V_InvMtV(...SpVectorSlice...)" );
00393 #endif
00394   assert_initialized();
00395   const size_type
00396     nG = G_ptr_->rows(),
00397     nM = M_.rows();
00398   AbstractLinAlgPack::V_InvMtV( y->sub_view(1,nG).get(), *G_ptr_, H_trans, x(1,nG) );
00399   AbstractLinAlgPack::V_InvMtV( y->sub_view(nG+1,nG+nM).get(), M_, H_trans, x(nG+1,nG+nM) );
00400 }
00401 
00402 // private
00403 
00404 void MatrixSymHessianRelaxNonSing::assert_initialized() const
00405 {
00406   TEUCHOS_TEST_FOR_EXCEPTION(
00407     G_ptr_.get() == NULL, std::logic_error
00408     ,"MatrixSymHessianRelaxNonSing::assert_initialized(): Error, Not initalized yet!" );
00409 }
00410 
00411 } // end namespace ConstrainedOptPack
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Defines