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

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