MOOCHO (Single Doxygen Collection) Version of the Day
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 // Redistribution and use in source and binary forms, with or without
00013 // modification, are permitted provided that the following conditions are
00014 // met:
00015 //
00016 // 1. Redistributions of source code must retain the above copyright
00017 // notice, this list of conditions and the following disclaimer.
00018 //
00019 // 2. Redistributions in binary form must reproduce the above copyright
00020 // notice, this list of conditions and the following disclaimer in the
00021 // documentation and/or other materials provided with the distribution.
00022 //
00023 // 3. Neither the name of the Corporation nor the names of the
00024 // contributors may be used to endorse or promote products derived from
00025 // this software without specific prior written permission.
00026 //
00027 // THIS SOFTWARE IS PROVIDED BY SANDIA CORPORATION "AS IS" AND ANY
00028 // EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
00029 // IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR
00030 // PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL SANDIA CORPORATION OR THE
00031 // CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL,
00032 // EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO,
00033 // PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR
00034 // PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF
00035 // LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING
00036 // NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
00037 // SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
00038 //
00039 // Questions? Contact Roscoe A. Bartlett (rabartl@sandia.gov) 
00040 // 
00041 // ***********************************************************************
00042 // @HEADER
00043 
00044 #include "ConstrainedOptPack_MatrixSymPosDefInvCholFactor.hpp"
00045 #include "SymInvCholMatrixOp.hpp"
00046 #include "AbstractLinAlgPack_SpVectorOp.hpp"
00047 #include "DenseLinAlgPack_DVectorOp.hpp"
00048 #include "DenseLinAlgPack_LinAlgOpPack.hpp"
00049 #include "DenseLinAlgPack_DMatrixOp.hpp"
00050 #include "DenseLinAlgPack_DMatrixOut.hpp"
00051 
00052 namespace LinAlgOpPack {
00053 
00054 using AbstractLinAlgPack::Vp_StV;
00055 using AbstractLinAlgPack::Vp_StMtV;
00056 using AbstractLinAlgPack::Mp_StM;
00057 using ConstrainedOptPack::Vp_StMtV;
00058 
00059 } // end namespace LinAlgOpPack
00060 
00061 namespace ConstrainedOptPack {
00062 
00063 // Overridden from Matrix 
00064 
00065 size_type MatrixSymPosDefInvCholFactor::cols() const
00066 {
00067   return rows();
00068 }
00069 
00070 // Overridden from MatrixOp
00071 
00072 MatrixOp& MatrixSymPosDefInvCholFactor::operator=(const MatrixOp& m)
00073 {
00074   return MatrixWithOpConcreteEncap<SymInvCholMatrix>::operator=(m);
00075 }
00076 
00077 std::ostream& MatrixSymPosDefInvCholFactor::output(std::ostream& out) const
00078 { return out << "\n*** Inverse Cholesky factor:\n" << m().UInv(); }
00079 
00080 // Level-2 BLAS
00081 
00082 void MatrixSymPosDefInvCholFactor::Vp_StMtV(DVectorSlice* vs_lhs, value_type alpha, BLAS_Cpp::Transp trans_rhs1
00083   , const DVectorSlice& vs_rhs2, value_type beta) const
00084 {
00085   ConstrainedOptPack::Vp_StMtV(vs_lhs,alpha,m(),trans_rhs1,vs_rhs2,beta);
00086 }
00087 
00088 void MatrixSymPosDefInvCholFactor::Vp_StMtV(DVectorSlice* vs_lhs, value_type alpha, BLAS_Cpp::Transp trans_rhs1
00089   , const SpVectorSlice& sv_rhs2, value_type beta) const
00090 {
00091   using LinAlgOpPack::assign;
00092   DVector vs_rhs2;
00093   assign(&vs_rhs2,sv_rhs2);
00094   ConstrainedOptPack::Vp_StMtV(vs_lhs,alpha,m(),trans_rhs1,vs_rhs2,beta);
00095 }
00096 
00097 value_type MatrixSymPosDefInvCholFactor::transVtMtV(const DVectorSlice& vs_rhs1, BLAS_Cpp::Transp trans_rhs2
00098   , const DVectorSlice& vs_rhs3) const
00099 {
00100   return ConstrainedOptPack::transVtMtV(vs_rhs1,m(),vs_rhs3);
00101 }
00102 
00103 value_type MatrixSymPosDefInvCholFactor::transVtMtV(const SpVectorSlice& sv_rhs1, BLAS_Cpp::Transp trans_rhs2
00104   , const SpVectorSlice& sv_rhs3) const
00105 {
00106   using LinAlgOpPack::assign;
00107   DVector vs_rhs1, vs_rhs3;
00108   assign(&vs_rhs1,sv_rhs1);
00109   assign(&vs_rhs3,sv_rhs3);
00110   return ConstrainedOptPack::transVtMtV(vs_rhs1,m(),vs_rhs3);
00111 }
00112 
00113 // Overridden from MatrixFactorized
00114 
00115 void MatrixSymPosDefInvCholFactor::V_InvMtV(DVector* v_lhs, BLAS_Cpp::Transp trans_rhs1
00116   , const DVectorSlice& vs_rhs2) const
00117 {
00118   ConstrainedOptPack::V_InvMtV(v_lhs,m(),vs_rhs2);
00119 }
00120 
00121 void MatrixSymPosDefInvCholFactor::V_InvMtV(DVectorSlice* vs_lhs, BLAS_Cpp::Transp trans_rhs1
00122   , const DVectorSlice& vs_rhs2) const
00123 {
00124   ConstrainedOptPack::V_InvMtV(vs_lhs,m(),vs_rhs2);
00125 }
00126 
00127 void MatrixSymPosDefInvCholFactor::V_InvMtV(DVector* v_lhs, BLAS_Cpp::Transp trans_rhs1
00128   , const SpVectorSlice& sv_rhs2) const
00129 {
00130   ConstrainedOptPack::V_InvMtV(v_lhs,m(),sv_rhs2);
00131 }
00132 
00133 void MatrixSymPosDefInvCholFactor::V_InvMtV(DVectorSlice* vs_lhs, BLAS_Cpp::Transp trans_rhs1
00134   , const SpVectorSlice& sv_rhs2) const
00135 {
00136   ConstrainedOptPack::V_InvMtV(vs_lhs,m(),sv_rhs2);
00137 }
00138 
00139 value_type MatrixSymPosDefInvCholFactor::transVtInvMtV(const DVectorSlice& vs_rhs1
00140   , BLAS_Cpp::Transp trans_rhs2, const DVectorSlice& vs_rhs3) const
00141 {
00142   return ConstrainedOptPack::transVtInvMtV(vs_rhs1,m(),vs_rhs3);
00143 }
00144 
00145 value_type MatrixSymPosDefInvCholFactor::transVtInvMtV(const SpVectorSlice& sv_rhs1
00146   , BLAS_Cpp::Transp trans_rhs2, const SpVectorSlice& sv_rhs3) const
00147 {
00148   return ConstrainedOptPack::transVtInvMtV(sv_rhs1,m(),sv_rhs3);}
00149 
00150 // Overridden from MatrixSymFactorized
00151 
00152 void MatrixSymPosDefInvCholFactor::M_StMtInvMtM(
00153     DMatrixSliceSym* S, value_type a, const MatrixOp& B
00154   , BLAS_Cpp::Transp B_trans, EMatrixDummyArg dummy_arg ) const
00155 {
00156 //  // Uncomment to use the defalut implementation (for debugging)
00157 //  MatrixSymFactorized::M_StMtInvMtM(S,a,B,B_trans,dummy_arg); return;
00158 
00159   using BLAS_Cpp::trans;
00160   using BLAS_Cpp::no_trans;
00161   using BLAS_Cpp::trans_not;
00162   using BLAS_Cpp::upper;
00163   using BLAS_Cpp::nonunit;
00164   using AbstractLinAlgPack::M_StInvMtM;
00165   using DenseLinAlgPack::tri;
00166   using DenseLinAlgPack::syrk;
00167   using DenseLinAlgPack::M_StInvMtM;
00168   using LinAlgOpPack::M_StMtM;
00169   using LinAlgOpPack::assign;
00170 
00171   DenseLinAlgPack::MtM_assert_sizes( rows(), cols(), no_trans
00172     , B.rows(), B.cols(), trans_not(B_trans) );
00173   DenseLinAlgPack::Mp_MtM_assert_sizes( S->rows(), S->cols(), no_trans
00174     , B.rows(), B.cols(), B_trans
00175     , B.rows(), B.cols(), trans_not(B_trans) );
00176   //
00177   // S = a * op(B) * inv(M) * op(B)'
00178   // 
00179   // M = L * L'
00180   // inv(M) = inv(L * L') = inv(L') * inv(L) = UInv * UInv'
00181   // 
00182   // S = a * op(B) * UInv * UInv' * op(B)'
00183   // 
00184   // T = op(B)'
00185   // 
00186   // T = UInv' * T (inplace with BLAS)
00187   // 
00188   // S = a * T' * T
00189   // 
00190 
00191   // T = op(B)'
00192   DMatrix T;
00193   assign( &T, B, trans_not(B_trans) );
00194   // T = UInv' * T (inplace with BLAS)
00195   M_StMtM( &T(), 1.0, tri(m().UInv(),upper,nonunit), trans, T(), no_trans );
00196   // S = a * T' * T
00197   syrk( trans, a, T(), 0.0, S );
00198 }
00199 
00200 // Overridden from MatrixSymSecant
00201 
00202 void MatrixSymPosDefInvCholFactor::init_identity(size_type n, value_type alpha)
00203 {
00204   if( alpha <= 0.0 ) {
00205     std::ostringstream omsg;
00206     omsg  << "MatrixSymPosDefInvCholFactor::init_identity(...) : Error, alpha = " << alpha
00207         << " <= 0.0 and therefore this is not a positive definite matrix.";
00208     throw UpdateSkippedException( omsg.str() ); 
00209   }
00210   m().resize(n);
00211   m().UInv() = 0.0;
00212   m().UInv().diag() = 1.0 / ::sqrt( alpha );
00213 }
00214 
00215 void MatrixSymPosDefInvCholFactor::init_diagonal( const DVectorSlice& diag )
00216 {
00217   DVectorSlice::const_iterator
00218     min_ele_ptr = std::min_element( diag.begin(), diag.end() );
00219   if( *min_ele_ptr <= 0.0 ) {
00220     std::ostringstream omsg;
00221     omsg  << "MatrixSymPosDefInvCholFactor::init_diagonal(...) : Error, "
00222         << "diag(" << min_ele_ptr - diag.begin() + 1 << " ) = "
00223         << (*min_ele_ptr) << " <= 0.0.\n"
00224         << "Therefore this is not a positive definite matrix.";
00225     throw UpdateSkippedException( omsg.str() ); 
00226   }
00227   const size_type n = diag.size();
00228   m().resize(n);
00229   m().UInv() = 0.0;
00230 
00231   DVectorSlice::const_iterator
00232     diag_itr = diag.begin();
00233   DVectorSlice::iterator
00234     inv_fact_diag_itr = m().UInv().diag().begin();
00235 
00236   while( diag_itr != diag.end() )
00237     *inv_fact_diag_itr++ = 1.0 / ::sqrt( *diag_itr++ );
00238 }
00239 
00240 void MatrixSymPosDefInvCholFactor::secant_update(DVectorSlice* s, DVectorSlice* y, DVectorSlice* _Bs)
00241 {
00242   using LinAlgOpPack::V_MtV;
00243   try {
00244     if(!_Bs) {
00245       DVector Bs;
00246       V_MtV( &Bs, *this, BLAS_Cpp::no_trans, *s );
00247       ConstrainedOptPack::BFGS_update(&m(),s,y,&Bs());
00248     }
00249     else {
00250       ConstrainedOptPack::BFGS_update(&m(),s,y,_Bs);
00251     }
00252   }
00253   catch(const BFGSUpdateSkippedException& excpt) {
00254     throw UpdateSkippedException( excpt.what() );
00255   }
00256 }
00257 
00258 // Overridden from MatrixExtractInvCholFactor
00259 
00260 void MatrixSymPosDefInvCholFactor::extract_inv_chol( DMatrixSliceTriEle* InvChol ) const
00261 {
00262   DenseLinAlgPack::assign( InvChol, DenseLinAlgPack::tri_ele( m().UInv(), BLAS_Cpp::upper ) );
00263 }
00264 
00265 // Overridden from Serializable
00266 
00267 void MatrixSymPosDefInvCholFactor::serialize( std::ostream &out ) const
00268 {
00269   TEUCHOS_TEST_FOR_EXCEPT(true);
00270 }
00271 
00272 void MatrixSymPosDefInvCholFactor::unserialize( std::istream &in )
00273 {
00274   TEUCHOS_TEST_FOR_EXCEPT(true);
00275 }
00276 
00277 } // end namespace ConstrainedOptPack
00278 
00279 #endif // 0
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Defines