MOOCHO (Single Doxygen Collection)  Version of the Day
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
ConstrainedOptPack_MatrixSymPosDefInvCholFactor.cpp
Go to the documentation of this file.
1 #if 0
2 
3 // @HEADER
4 // ***********************************************************************
5 //
6 // Moocho: Multi-functional Object-Oriented arCHitecture for Optimization
7 // Copyright (2003) Sandia Corporation
8 //
9 // Under terms of Contract DE-AC04-94AL85000, there is a non-exclusive
10 // license for use of this work by or on behalf of the U.S. Government.
11 //
12 // Redistribution and use in source and binary forms, with or without
13 // modification, are permitted provided that the following conditions are
14 // met:
15 //
16 // 1. Redistributions of source code must retain the above copyright
17 // notice, this list of conditions and the following disclaimer.
18 //
19 // 2. Redistributions in binary form must reproduce the above copyright
20 // notice, this list of conditions and the following disclaimer in the
21 // documentation and/or other materials provided with the distribution.
22 //
23 // 3. Neither the name of the Corporation nor the names of the
24 // contributors may be used to endorse or promote products derived from
25 // this software without specific prior written permission.
26 //
27 // THIS SOFTWARE IS PROVIDED BY SANDIA CORPORATION "AS IS" AND ANY
28 // EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
29 // IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR
30 // PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL SANDIA CORPORATION OR THE
31 // CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL,
32 // EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO,
33 // PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR
34 // PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF
35 // LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING
36 // NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
37 // SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
38 //
39 // Questions? Contact Roscoe A. Bartlett (rabartl@sandia.gov)
40 //
41 // ***********************************************************************
42 // @HEADER
43 
45 #include "SymInvCholMatrixOp.hpp"
51 
52 namespace LinAlgOpPack {
53 
58 
59 } // end namespace LinAlgOpPack
60 
61 namespace ConstrainedOptPack {
62 
63 // Overridden from Matrix
64 
66 {
67  return rows();
68 }
69 
70 // Overridden from MatrixOp
71 
72 MatrixOp& MatrixSymPosDefInvCholFactor::operator=(const MatrixOp& m)
73 {
74  return MatrixWithOpConcreteEncap<SymInvCholMatrix>::operator=(m);
75 }
76 
77 std::ostream& MatrixSymPosDefInvCholFactor::output(std::ostream& out) const
78 { return out << "\n*** Inverse Cholesky factor:\n" << m().UInv(); }
79 
80 // Level-2 BLAS
81 
83  , const DVectorSlice& vs_rhs2, value_type beta) const
84 {
85  ConstrainedOptPack::Vp_StMtV(vs_lhs,alpha,m(),trans_rhs1,vs_rhs2,beta);
86 }
87 
89  , const SpVectorSlice& sv_rhs2, value_type beta) const
90 {
92  DVector vs_rhs2;
93  assign(&vs_rhs2,sv_rhs2);
94  ConstrainedOptPack::Vp_StMtV(vs_lhs,alpha,m(),trans_rhs1,vs_rhs2,beta);
95 }
96 
98  , const DVectorSlice& vs_rhs3) const
99 {
100  return ConstrainedOptPack::transVtMtV(vs_rhs1,m(),vs_rhs3);
101 }
102 
104  , const SpVectorSlice& sv_rhs3) const
105 {
106  using LinAlgOpPack::assign;
107  DVector vs_rhs1, vs_rhs3;
108  assign(&vs_rhs1,sv_rhs1);
109  assign(&vs_rhs3,sv_rhs3);
110  return ConstrainedOptPack::transVtMtV(vs_rhs1,m(),vs_rhs3);
111 }
112 
113 // Overridden from MatrixFactorized
114 
116  , const DVectorSlice& vs_rhs2) const
117 {
118  ConstrainedOptPack::V_InvMtV(v_lhs,m(),vs_rhs2);
119 }
120 
122  , const DVectorSlice& vs_rhs2) const
123 {
124  ConstrainedOptPack::V_InvMtV(vs_lhs,m(),vs_rhs2);
125 }
126 
128  , const SpVectorSlice& sv_rhs2) const
129 {
130  ConstrainedOptPack::V_InvMtV(v_lhs,m(),sv_rhs2);
131 }
132 
134  , const SpVectorSlice& sv_rhs2) const
135 {
136  ConstrainedOptPack::V_InvMtV(vs_lhs,m(),sv_rhs2);
137 }
138 
140  , BLAS_Cpp::Transp trans_rhs2, const DVectorSlice& vs_rhs3) const
141 {
142  return ConstrainedOptPack::transVtInvMtV(vs_rhs1,m(),vs_rhs3);
143 }
144 
146  , BLAS_Cpp::Transp trans_rhs2, const SpVectorSlice& sv_rhs3) const
147 {
148  return ConstrainedOptPack::transVtInvMtV(sv_rhs1,m(),sv_rhs3);}
149 
150 // Overridden from MatrixSymFactorized
151 
153  DMatrixSliceSym* S, value_type a, const MatrixOp& B
154  , BLAS_Cpp::Transp B_trans, EMatrixDummyArg dummy_arg ) const
155 {
156 // // Uncomment to use the defalut implementation (for debugging)
157 // MatrixSymFactorized::M_StMtInvMtM(S,a,B,B_trans,dummy_arg); return;
158 
159  using BLAS_Cpp::trans;
160  using BLAS_Cpp::no_trans;
161  using BLAS_Cpp::trans_not;
162  using BLAS_Cpp::upper;
163  using BLAS_Cpp::nonunit;
165  using DenseLinAlgPack::tri;
166  using DenseLinAlgPack::syrk;
168  using LinAlgOpPack::M_StMtM;
169  using LinAlgOpPack::assign;
170 
172  , B.rows(), B.cols(), trans_not(B_trans) );
173  DenseLinAlgPack::Mp_MtM_assert_sizes( S->rows(), S->cols(), no_trans
174  , B.rows(), B.cols(), B_trans
175  , B.rows(), B.cols(), trans_not(B_trans) );
176  //
177  // S = a * op(B) * inv(M) * op(B)'
178  //
179  // M = L * L'
180  // inv(M) = inv(L * L') = inv(L') * inv(L) = UInv * UInv'
181  //
182  // S = a * op(B) * UInv * UInv' * op(B)'
183  //
184  // T = op(B)'
185  //
186  // T = UInv' * T (inplace with BLAS)
187  //
188  // S = a * T' * T
189  //
190 
191  // T = op(B)'
192  DMatrix T;
193  assign( &T, B, trans_not(B_trans) );
194  // T = UInv' * T (inplace with BLAS)
195  M_StMtM( &T(), 1.0, tri(m().UInv(),upper,nonunit), trans, T(), no_trans );
196  // S = a * T' * T
197  syrk( trans, a, T(), 0.0, S );
198 }
199 
200 // Overridden from MatrixSymSecant
201 
203 {
204  if( alpha <= 0.0 ) {
205  std::ostringstream omsg;
206  omsg << "MatrixSymPosDefInvCholFactor::init_identity(...) : Error, alpha = " << alpha
207  << " <= 0.0 and therefore this is not a positive definite matrix.";
208  throw UpdateSkippedException( omsg.str() );
209  }
210  m().resize(n);
211  m().UInv() = 0.0;
212  m().UInv().diag() = 1.0 / ::sqrt( alpha );
213 }
214 
216 {
217  DVectorSlice::const_iterator
218  min_ele_ptr = std::min_element( diag.begin(), diag.end() );
219  if( *min_ele_ptr <= 0.0 ) {
220  std::ostringstream omsg;
221  omsg << "MatrixSymPosDefInvCholFactor::init_diagonal(...) : Error, "
222  << "diag(" << min_ele_ptr - diag.begin() + 1 << " ) = "
223  << (*min_ele_ptr) << " <= 0.0.\n"
224  << "Therefore this is not a positive definite matrix.";
225  throw UpdateSkippedException( omsg.str() );
226  }
227  const size_type n = diag.size();
228  m().resize(n);
229  m().UInv() = 0.0;
230 
231  DVectorSlice::const_iterator
232  diag_itr = diag.begin();
233  DVectorSlice::iterator
234  inv_fact_diag_itr = m().UInv().diag().begin();
235 
236  while( diag_itr != diag.end() )
237  *inv_fact_diag_itr++ = 1.0 / ::sqrt( *diag_itr++ );
238 }
239 
241 {
242  using LinAlgOpPack::V_MtV;
243  try {
244  if(!_Bs) {
245  DVector Bs;
246  V_MtV( &Bs, *this, BLAS_Cpp::no_trans, *s );
247  ConstrainedOptPack::BFGS_update(&m(),s,y,&Bs());
248  }
249  else {
250  ConstrainedOptPack::BFGS_update(&m(),s,y,_Bs);
251  }
252  }
253  catch(const BFGSUpdateSkippedException& excpt) {
254  throw UpdateSkippedException( excpt.what() );
255  }
256 }
257 
258 // Overridden from MatrixExtractInvCholFactor
259 
261 {
263 }
264 
265 // Overridden from Serializable
266 
267 void MatrixSymPosDefInvCholFactor::serialize( std::ostream &out ) const
268 {
270 }
271 
272 void MatrixSymPosDefInvCholFactor::unserialize( std::istream &in )
273 {
275 }
276 
277 } // end namespace ConstrainedOptPack
278 
279 #endif // 0
void M_StMtM(MatrixOp *M_lhs, value_type alpha, const MatrixOp &M_rhs1, BLAS_Cpp::Transp trans_rhs1, const MatrixOp &M_rhs2, BLAS_Cpp::Transp trans_rhs2)
M_lhs = alpha * op(M_rhs1) * op(M_rhs2).
AbstractLinAlgPack::size_type size_type
value_type transVtInvMtV(const Vector &v_rhs1, const MatrixNonsing &M_rhs2, BLAS_Cpp::Transp trans_rhs2, const Vector &v_rhs3)
result = v_rhs1' * inv(op(M_rhs2)) * v_rhs3
void assign(DMatrix *gm_lhs, value_type alpha)
gm_lhs = alpha (elementwise)
void MtM_assert_sizes(size_type m_rhs1_rows, size_type m_rhs1_cols, BLAS_Cpp::Transp trans_rhs1, size_type m_rhs2_rows, size_type m_rhs2_cols, BLAS_Cpp::Transp trans_rhs2)
op(m_lhs) += op(m_rhs1)
void M_StMtInvMtM(DMatrixSliceSym *sym_gms_lhs, value_type alpha, const MatrixOp &mwo, BLAS_Cpp::Transp mwo_trans, EMatrixDummyArg) const
const DMatrixSliceTriEle tri_ele(const DMatrixSlice &gms, BLAS_Cpp::Uplo uplo)
void init_diagonal(const DVectorSlice &diag)
void Vp_StV(VectorMutable *v_lhs, const value_type &alpha, const Vector &v_rhs)
v_lhs = alpha * v_rhs + v_lhs
size_type rows(size_type rows, size_type cols, BLAS_Cpp::Transp _trans)
Return rows of a possible transposed matrix.
void assign(VectorMutable *v_lhs, const V &V_rhs)
v_lhs = V_rhs.
void serialize(std::ostream &out) const
void Mp_MtM_assert_sizes(size_type m_lhs_rows, size_type m_lhs_cols, BLAS_Cpp::Transp trans_lhs, size_type m_rhs1_rows, size_type m_rhs1_cols, BLAS_Cpp::Transp trans_rhs1, size_type m_rhs2_rows, size_type m_rhs2_cols, BLAS_Cpp::Transp trans_rhs2)
op(m_lhs) += op(m_rhs1) * op(m_rhs2)
void secant_update(DVectorSlice *s, DVectorSlice *y, DVectorSlice *Bs)
void Vp_StMtV(DVectorSlice *vs_lhs, value_type alpha, BLAS_Cpp::Transp trans_rhs1, const DVectorSlice &vs_rhs2, value_type beta) const
void syrk(BLAS_Cpp::Transp trans, value_type alpha, const DMatrixSlice &gms_rhs, value_type beta, DMatrixSliceSym *sym_lhs)
value_type transVtMtV(const DVectorSlice &vs_rhs1, BLAS_Cpp::Transp trans_rhs2, const DVectorSlice &vs_rhs3) const
Transposed.
const DMatrixSliceTri tri(const DMatrixSlice &gms, BLAS_Cpp::Uplo uplo, BLAS_Cpp::Diag diag)
value_type transVtMtV(const Vector &v_rhs1, const MatrixOp &M_rhs2, BLAS_Cpp::Transp trans_rhs2, const Vector &v_rhs3)
result = v_rhs1' * op(M_rhs2) * v_rhs3
Not transposed.
std::ostream & output(std::ostream &out) const
value_type transVtInvMtV(const DVectorSlice &vs_rhs1, BLAS_Cpp::Transp trans_rhs2, const DVectorSlice &vs_rhs3) const
void syrk(const MatrixOp &mwo_rhs, BLAS_Cpp::Transp M_trans, value_type alpha, value_type beta, MatrixSymOp *sym_lhs)
Perform a rank-k update of a symmetric matrix of the form:
std::ostream * out
DenseLinAlgPack::DMatrixSliceTriEle DMatrixSliceTriEle
void V_InvMtV(VectorMutable *v_lhs, const MatrixNonsing &M_rhs1, BLAS_Cpp::Transp trans_rhs1, const Vector &v_rhs2)
v_lhs = inv(op(M_rhs1)) * v_rhs2
f_dbl_prec f_dbl_prec f_dbl_prec * S
void Mp_StM(MatrixOp *mwo_lhs, value_type alpha, const MatrixOp &M_rhs, BLAS_Cpp::Transp trans_rhs)
void Vp_StMtV(VectorMutable *v_lhs, value_type alpha, const MatrixOp &M_rhs1, BLAS_Cpp::Transp trans_rhs1, const Vector &v_rhs2, value_type beta=1.0)
v_lhs = alpha * op(M_rhs1) * v_rhs2 + beta * v_lhs (BLAS xGEMV)
void extract_inv_chol(DMatrixSliceTriEle *InvChol) const
void init_identity(size_type n, value_type alpha)
void sqrt(DVectorSlice *vs_lhs, const DVectorSlice &vs_rhs)
vs_lhs = sqrt(vs_rhs)
SparseVectorSlice< SparseElement< index_type, value_type > > SpVectorSlice
Transp trans_not(Transp _trans)
Return the opposite of the transpose argument.
DenseLinAlgPack::VectorSliceTmpl< value_type > DVectorSlice
void V_InvMtV(DVector *v_lhs, BLAS_Cpp::Transp trans_rhs1, const DVectorSlice &vs_rhs2) const
void V_MtV(VectorMutable *v_lhs, const MatrixOp &M_rhs1, BLAS_Cpp::Transp trans_rhs1, const V &V_rhs2)
v_lhs = op(M_rhs1) * V_rhs2.
AbstractLinAlgPack::value_type value_type
void M_StInvMtM(MatrixOp *m_lhs, value_type alpha, const MatrixNonsing &M_rhs1, BLAS_Cpp::Transp trans_rhs1, const MatrixOp &mwo_rhs2, BLAS_Cpp::Transp trans_rhs2)
m_lhs = alpha * inv(op(mwo_rhs1)) * op(mwo_rhs2) (right)
DenseLinAlgPack::DMatrixSliceSym DMatrixSliceSym
Transp
TRANS.
#define TEUCHOS_TEST_FOR_EXCEPT(throw_exception_test)
void M_StInvMtM(DMatrix *gm_lhs, value_type alpha, const DMatrixSliceTri &tri_rhs1, BLAS_Cpp::Transp trans_rhs1, const DMatrixSlice &gms_rhs2, BLAS_Cpp::Transp trans_rhs2)