ConstrainedOptPack_MatrixSymPosDefInvCholFactor.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 "ConstrainedOptPack_MatrixSymPosDefInvCholFactor.hpp"
00030 #include "SymInvCholMatrixOp.hpp"
00031 #include "AbstractLinAlgPack_SpVectorOp.hpp"
00032 #include "DenseLinAlgPack_DVectorOp.hpp"
00033 #include "DenseLinAlgPack_LinAlgOpPack.hpp"
00034 #include "DenseLinAlgPack_DMatrixOp.hpp"
00035 #include "DenseLinAlgPack_DMatrixOut.hpp"
00036 
00037 namespace LinAlgOpPack {
00038 
00039 using AbstractLinAlgPack::Vp_StV;
00040 using AbstractLinAlgPack::Vp_StMtV;
00041 using AbstractLinAlgPack::Mp_StM;
00042 using ConstrainedOptPack::Vp_StMtV;
00043 
00044 } // end namespace LinAlgOpPack
00045 
00046 namespace ConstrainedOptPack {
00047 
00048 // Overridden from Matrix 
00049 
00050 size_type MatrixSymPosDefInvCholFactor::cols() const
00051 {
00052   return rows();
00053 }
00054 
00055 // Overridden from MatrixOp
00056 
00057 MatrixOp& MatrixSymPosDefInvCholFactor::operator=(const MatrixOp& m)
00058 {
00059   return MatrixWithOpConcreteEncap<SymInvCholMatrix>::operator=(m);
00060 }
00061 
00062 std::ostream& MatrixSymPosDefInvCholFactor::output(std::ostream& out) const
00063 { return out << "\n*** Inverse Cholesky factor:\n" << m().UInv(); }
00064 
00065 // Level-2 BLAS
00066 
00067 void MatrixSymPosDefInvCholFactor::Vp_StMtV(DVectorSlice* vs_lhs, value_type alpha, BLAS_Cpp::Transp trans_rhs1
00068   , const DVectorSlice& vs_rhs2, value_type beta) const
00069 {
00070   ConstrainedOptPack::Vp_StMtV(vs_lhs,alpha,m(),trans_rhs1,vs_rhs2,beta);
00071 }
00072 
00073 void MatrixSymPosDefInvCholFactor::Vp_StMtV(DVectorSlice* vs_lhs, value_type alpha, BLAS_Cpp::Transp trans_rhs1
00074   , const SpVectorSlice& sv_rhs2, value_type beta) const
00075 {
00076   using LinAlgOpPack::assign;
00077   DVector vs_rhs2;
00078   assign(&vs_rhs2,sv_rhs2);
00079   ConstrainedOptPack::Vp_StMtV(vs_lhs,alpha,m(),trans_rhs1,vs_rhs2,beta);
00080 }
00081 
00082 value_type MatrixSymPosDefInvCholFactor::transVtMtV(const DVectorSlice& vs_rhs1, BLAS_Cpp::Transp trans_rhs2
00083   , const DVectorSlice& vs_rhs3) const
00084 {
00085   return ConstrainedOptPack::transVtMtV(vs_rhs1,m(),vs_rhs3);
00086 }
00087 
00088 value_type MatrixSymPosDefInvCholFactor::transVtMtV(const SpVectorSlice& sv_rhs1, BLAS_Cpp::Transp trans_rhs2
00089   , const SpVectorSlice& sv_rhs3) const
00090 {
00091   using LinAlgOpPack::assign;
00092   DVector vs_rhs1, vs_rhs3;
00093   assign(&vs_rhs1,sv_rhs1);
00094   assign(&vs_rhs3,sv_rhs3);
00095   return ConstrainedOptPack::transVtMtV(vs_rhs1,m(),vs_rhs3);
00096 }
00097 
00098 // Overridden from MatrixFactorized
00099 
00100 void MatrixSymPosDefInvCholFactor::V_InvMtV(DVector* v_lhs, BLAS_Cpp::Transp trans_rhs1
00101   , const DVectorSlice& vs_rhs2) const
00102 {
00103   ConstrainedOptPack::V_InvMtV(v_lhs,m(),vs_rhs2);
00104 }
00105 
00106 void MatrixSymPosDefInvCholFactor::V_InvMtV(DVectorSlice* vs_lhs, BLAS_Cpp::Transp trans_rhs1
00107   , const DVectorSlice& vs_rhs2) const
00108 {
00109   ConstrainedOptPack::V_InvMtV(vs_lhs,m(),vs_rhs2);
00110 }
00111 
00112 void MatrixSymPosDefInvCholFactor::V_InvMtV(DVector* v_lhs, BLAS_Cpp::Transp trans_rhs1
00113   , const SpVectorSlice& sv_rhs2) const
00114 {
00115   ConstrainedOptPack::V_InvMtV(v_lhs,m(),sv_rhs2);
00116 }
00117 
00118 void MatrixSymPosDefInvCholFactor::V_InvMtV(DVectorSlice* vs_lhs, BLAS_Cpp::Transp trans_rhs1
00119   , const SpVectorSlice& sv_rhs2) const
00120 {
00121   ConstrainedOptPack::V_InvMtV(vs_lhs,m(),sv_rhs2);
00122 }
00123 
00124 value_type MatrixSymPosDefInvCholFactor::transVtInvMtV(const DVectorSlice& vs_rhs1
00125   , BLAS_Cpp::Transp trans_rhs2, const DVectorSlice& vs_rhs3) const
00126 {
00127   return ConstrainedOptPack::transVtInvMtV(vs_rhs1,m(),vs_rhs3);
00128 }
00129 
00130 value_type MatrixSymPosDefInvCholFactor::transVtInvMtV(const SpVectorSlice& sv_rhs1
00131   , BLAS_Cpp::Transp trans_rhs2, const SpVectorSlice& sv_rhs3) const
00132 {
00133   return ConstrainedOptPack::transVtInvMtV(sv_rhs1,m(),sv_rhs3);}
00134 
00135 // Overridden from MatrixSymFactorized
00136 
00137 void MatrixSymPosDefInvCholFactor::M_StMtInvMtM(
00138     DMatrixSliceSym* S, value_type a, const MatrixOp& B
00139   , BLAS_Cpp::Transp B_trans, EMatrixDummyArg dummy_arg ) const
00140 {
00141 //  // Uncomment to use the defalut implementation (for debugging)
00142 //  MatrixSymFactorized::M_StMtInvMtM(S,a,B,B_trans,dummy_arg); return;
00143 
00144   using BLAS_Cpp::trans;
00145   using BLAS_Cpp::no_trans;
00146   using BLAS_Cpp::trans_not;
00147   using BLAS_Cpp::upper;
00148   using BLAS_Cpp::nonunit;
00149   using AbstractLinAlgPack::M_StInvMtM;
00150   using DenseLinAlgPack::tri;
00151   using DenseLinAlgPack::syrk;
00152   using DenseLinAlgPack::M_StInvMtM;
00153   using LinAlgOpPack::M_StMtM;
00154   using LinAlgOpPack::assign;
00155 
00156   DenseLinAlgPack::MtM_assert_sizes( rows(), cols(), no_trans
00157     , B.rows(), B.cols(), trans_not(B_trans) );
00158   DenseLinAlgPack::Mp_MtM_assert_sizes( S->rows(), S->cols(), no_trans
00159     , B.rows(), B.cols(), B_trans
00160     , B.rows(), B.cols(), trans_not(B_trans) );
00161   //
00162   // S = a * op(B) * inv(M) * op(B)'
00163   // 
00164   // M = L * L'
00165   // inv(M) = inv(L * L') = inv(L') * inv(L) = UInv * UInv'
00166   // 
00167   // S = a * op(B) * UInv * UInv' * op(B)'
00168   // 
00169   // T = op(B)'
00170   // 
00171   // T = UInv' * T (inplace with BLAS)
00172   // 
00173   // S = a * T' * T
00174   // 
00175 
00176   // T = op(B)'
00177   DMatrix T;
00178   assign( &T, B, trans_not(B_trans) );
00179   // T = UInv' * T (inplace with BLAS)
00180   M_StMtM( &T(), 1.0, tri(m().UInv(),upper,nonunit), trans, T(), no_trans );
00181   // S = a * T' * T
00182   syrk( trans, a, T(), 0.0, S );
00183 }
00184 
00185 // Overridden from MatrixSymSecant
00186 
00187 void MatrixSymPosDefInvCholFactor::init_identity(size_type n, value_type alpha)
00188 {
00189   if( alpha <= 0.0 ) {
00190     std::ostringstream omsg;
00191     omsg  << "MatrixSymPosDefInvCholFactor::init_identity(...) : Error, alpha = " << alpha
00192         << " <= 0.0 and therefore this is not a positive definite matrix.";
00193     throw UpdateSkippedException( omsg.str() ); 
00194   }
00195   m().resize(n);
00196   m().UInv() = 0.0;
00197   m().UInv().diag() = 1.0 / ::sqrt( alpha );
00198 }
00199 
00200 void MatrixSymPosDefInvCholFactor::init_diagonal( const DVectorSlice& diag )
00201 {
00202   DVectorSlice::const_iterator
00203     min_ele_ptr = std::min_element( diag.begin(), diag.end() );
00204   if( *min_ele_ptr <= 0.0 ) {
00205     std::ostringstream omsg;
00206     omsg  << "MatrixSymPosDefInvCholFactor::init_diagonal(...) : Error, "
00207         << "diag(" << min_ele_ptr - diag.begin() + 1 << " ) = "
00208         << (*min_ele_ptr) << " <= 0.0.\n"
00209         << "Therefore this is not a positive definite matrix.";
00210     throw UpdateSkippedException( omsg.str() ); 
00211   }
00212   const size_type n = diag.size();
00213   m().resize(n);
00214   m().UInv() = 0.0;
00215 
00216   DVectorSlice::const_iterator
00217     diag_itr = diag.begin();
00218   DVectorSlice::iterator
00219     inv_fact_diag_itr = m().UInv().diag().begin();
00220 
00221   while( diag_itr != diag.end() )
00222     *inv_fact_diag_itr++ = 1.0 / ::sqrt( *diag_itr++ );
00223 }
00224 
00225 void MatrixSymPosDefInvCholFactor::secant_update(DVectorSlice* s, DVectorSlice* y, DVectorSlice* _Bs)
00226 {
00227   using LinAlgOpPack::V_MtV;
00228   try {
00229     if(!_Bs) {
00230       DVector Bs;
00231       V_MtV( &Bs, *this, BLAS_Cpp::no_trans, *s );
00232       ConstrainedOptPack::BFGS_update(&m(),s,y,&Bs());
00233     }
00234     else {
00235       ConstrainedOptPack::BFGS_update(&m(),s,y,_Bs);
00236     }
00237   }
00238   catch(const BFGSUpdateSkippedException& excpt) {
00239     throw UpdateSkippedException( excpt.what() );
00240   }
00241 }
00242 
00243 // Overridden from MatrixExtractInvCholFactor
00244 
00245 void MatrixSymPosDefInvCholFactor::extract_inv_chol( DMatrixSliceTriEle* InvChol ) const
00246 {
00247   DenseLinAlgPack::assign( InvChol, DenseLinAlgPack::tri_ele( m().UInv(), BLAS_Cpp::upper ) );
00248 }
00249 
00250 // Overridden from Serializable
00251 
00252 void MatrixSymPosDefInvCholFactor::serialize( std::ostream &out ) const
00253 {
00254   TEST_FOR_EXCEPT(true);
00255 }
00256 
00257 void MatrixSymPosDefInvCholFactor::unserialize( std::istream &in )
00258 {
00259   TEST_FOR_EXCEPT(true);
00260 }
00261 
00262 } // end namespace ConstrainedOptPack

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