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

Generated on Tue Oct 20 12:51:45 2009 for MOOCHO (Single Doxygen Collection) by doxygen 1.4.7