ConstrainedOptPack: C++ Tools for Constrained (and Unconstrained) Optimization  Version of the Day
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Groups Pages
ConstrainedOptPack_MatrixSymPosDefLBFGS.hpp
1 // @HEADER
2 // ***********************************************************************
3 //
4 // Moocho: Multi-functional Object-Oriented arCHitecture for Optimization
5 // Copyright (2003) Sandia Corporation
6 //
7 // Under terms of Contract DE-AC04-94AL85000, there is a non-exclusive
8 // license for use of this work by or on behalf of the U.S. Government.
9 //
10 // Redistribution and use in source and binary forms, with or without
11 // modification, are permitted provided that the following conditions are
12 // met:
13 //
14 // 1. Redistributions of source code must retain the above copyright
15 // notice, this list of conditions and the following disclaimer.
16 //
17 // 2. Redistributions in binary form must reproduce the above copyright
18 // notice, this list of conditions and the following disclaimer in the
19 // documentation and/or other materials provided with the distribution.
20 //
21 // 3. Neither the name of the Corporation nor the names of the
22 // contributors may be used to endorse or promote products derived from
23 // this software without specific prior written permission.
24 //
25 // THIS SOFTWARE IS PROVIDED BY SANDIA CORPORATION "AS IS" AND ANY
26 // EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
27 // IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR
28 // PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL SANDIA CORPORATION OR THE
29 // CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL,
30 // EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO,
31 // PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR
32 // PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF
33 // LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING
34 // NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
35 // SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
36 //
37 // Questions? Contact Roscoe A. Bartlett (rabartl@sandia.gov)
38 //
39 // ***********************************************************************
40 // @HEADER
41 
42 #ifndef MATRIX_SYM_POS_DEF_LBFGS_H
43 #define MATRIX_SYM_POS_DEF_LBFGS_H
44 
45 #include <vector>
46 
47 #include "ConstrainedOptPack_Types.hpp"
48 #include "AbstractLinAlgPack_MatrixSymSecant.hpp"
49 #include "AbstractLinAlgPack_MatrixSymOpNonsing.hpp"
50 #include "AbstractLinAlgPack_MultiVectorMutable.hpp"
51 #include "AbstractLinAlgPack_VectorSpace.hpp"
52 #include "DenseLinAlgPack_DMatrixAsTriSym.hpp"
54 
55 namespace ConstrainedOptPack {
56 
111  : public AbstractLinAlgPack::MatrixSymOpNonsing // doxygen needs full path
112  , public MatrixSymSecant
113 {
114 public:
115 
118 
121 
124  class PostMod {
125  public:
128  size_type m = 10
129  ,bool maintain_original = true
130  ,bool maintain_inverse = true
131  ,bool auto_rescaling = false
132  )
133  :m_(m)
134  ,maintain_original_(maintain_original)
135  ,maintain_inverse_(maintain_inverse)
136  ,auto_rescaling_(auto_rescaling)
137  {}
140  {
141  p->initial_setup(m_,maintain_original_,maintain_inverse_,auto_rescaling_);
142  }
143  private:
144  size_type m_;
145  bool maintain_original_;
146  bool maintain_inverse_;
147  bool auto_rescaling_;
148  }; // end PostMod
149 
151 
154 
157  size_type m = 10
158  ,bool maintain_original = true
159  ,bool maintain_inverse = true
160  ,bool auto_rescaling = false
161  );
162 
168  STANDARD_MEMBER_COMPOSITION_MEMBERS( bool, auto_rescaling );
169 
203  void initial_setup(
204  size_type m = 10
205  ,bool maintain_original = true
206  ,bool maintain_inverse = true
207  ,bool auto_rescaling = false
208  );
209 
210  // //////////////////////////////////
211  // Representation access
212 
214  size_type m() const;
216  size_type m_bar() const;
218  value_type gamma_k() const;
220  const multi_vec_ptr_t S() const;
222  const multi_vec_ptr_t Y() const;
224  bool maintain_original() const;
226  bool maintain_inverse() const;
231 
233 
236 
238  const VectorSpace& space_cols() const;
240  std::ostream& output(std::ostream& out) const;
242  MatrixOp& operator=(const MatrixOp& mwo);
244  void Vp_StMtV(
245  VectorMutable* v_lhs, value_type alpha, BLAS_Cpp::Transp trans_rhs1
246  , const Vector& v_rhs2, value_type beta ) const;
247 
249 
252 
254  void V_InvMtV(
255  VectorMutable* v_lhs, BLAS_Cpp::Transp trans_rhs1
256  , const Vector& v_rhs2 ) const;
257 
259 
262 
264  void init_identity( const VectorSpace& space_diag, value_type alpha );
271  void init_diagonal( const Vector& diag );
273  void secant_update(
274  VectorMutable *s
275  ,VectorMutable *y
276  ,VectorMutable *Bs
277  );
278 
280 
281 private:
282 
283  // //////////////////////////////////
284  // Private types
285 
286  typedef VectorSpace::multi_vec_mut_ptr_t multi_vec_mut_ptr_t;
287 
288  // //////////////////////////////////
289  // Private data members
290 
291  bool maintain_original_; // If true, qualities needed for Bk will be maintained
292  bool original_is_updated_;// If true, qualities needed for Bk are already updated
293  bool maintain_inverse_; // If true, quantities needed for Hk will be maintained
294  bool inverse_is_updated_; // If true, quantities needed for Hk are already updated
295 
296  VectorSpace::space_ptr_t
297  vec_spc_; // The vector space that everything is based on.
298 
299  size_type n_, // Size of the matrix. If 0 then is uninitialized
300  m_, // Maximum number of update vectors that can be stored.
301  m_bar_, // Current number of update vectors being stored.
302  // 0 <= m_bar <= m
303  num_secant_updates_; // Records the number of secant updates performed
304  value_type gamma_k_;// Scaling factor for Bo = (1/gamma_k) * I.
305 
306  multi_vec_mut_ptr_t
307  S_, // (n x m) Matrix of stored update vectors = [ s1, ..., sm ]
308  // S(:,m_bar) is the most recently stored s update vector
309  Y_; // (n_max x m) Matrix of stored update vectors = [ y1, ..., ym ]
310  // Y(:,k_bar) is the most recently stored y update vector
311  DMatrix STY_, // (m x m) The matrix S'Y
312  STSYTY_;// ((m+1) x (m+1)) The strictly upper triangular part stores the
313  // upper triangular part Y'Y and the strictly lower triangular
314  // part stores the lower triangular part of S'S. The diagonal
315  // can be used for workspace.
316 
317  mutable bool Q_updated_; // True if Q has been updated for the most current update.
318  mutable DMatrix QJ_; // Used to store factorization of the schur complement of Q.
319 
320  // //////////////////////////////////
321  // Private member functions
322 
323  // Access to important matrices.
324 
326  const DMatrixSliceTri R() const;
328  const DMatrixSliceTri Lb() const;
330  DMatrixSlice STY();
332  const DMatrixSlice STY() const;
334  DMatrixSliceSym STS();
336  const DMatrixSliceSym STS() const;
338  DMatrixSliceSym YTY();
340  const DMatrixSliceSym YTY() const;
342  void V_invQtV( DVectorSlice* y, const DVectorSlice& x ) const;
344  void Vp_DtV( DVectorSlice* y, const DVectorSlice& x ) const;
345 
346  // Updates
347 
349  void update_Q() const;
350 
352  void assert_initialized() const;
353 
354 }; // end class MatrixSymPosDefLBFGS
355 
356 // //////////////////////////////////////////////
357 // Inline member functions
358 
359 inline
361 {
362  return m_;
363 }
364 
365 inline
367 {
368  return m_bar_;
369 }
370 
371 inline
373 {
374  return gamma_k_;
375 }
376 
377 inline
380 {
381  return S_->mv_sub_view(1,n_,1,m_bar_);
382 }
383 
384 inline
387 {
388  return Y_->mv_sub_view(1,n_,1,m_bar_);
389 }
390 
391 inline
393 {
394  return maintain_original_;
395 }
396 
397 inline
399 {
400  return maintain_inverse_;
401 }
402 
403 inline
405 {
406  return num_secant_updates_;
407 }
408 
409 } // end namespace ConstrainedOptPack
410 
411 #endif // MATRIX_SYM_POS_DEF_LBFGS_H
void init_identity(const VectorSpace &space_diag, value_type alpha)
void secant_update(VectorMutable *s, VectorMutable *y, VectorMutable *Bs)
void initial_setup(size_type m=10, bool maintain_original=true, bool maintain_inverse=true, bool auto_rescaling=false)
Initial setup for the matrix.
PostMod(size_type m=10, bool maintain_original=true, bool maintain_inverse=true, bool auto_rescaling=false)
void V_InvMtV(VectorMutable *v_lhs, BLAS_Cpp::Transp trans_rhs1, const Vector &v_rhs2) const
MatrixSymPosDefLBFGS(size_type m=10, bool maintain_original=true, bool maintain_inverse=true, bool auto_rescaling=false)
Calls this->initial_setup()
void init_diagonal(const Vector &diag)
Actually this calls init_identity( diag.space(), diag.norm_inf() ).
PostMod class to use with MemMngPack::AbstractFactorStd.
STANDARD_MEMBER_COMPOSITION_MEMBERS(bool, auto_rescaling)
Set whether automatic rescaling is used or not.
size_t size_type
void Vp_StMtV(VectorMutable *v_lhs, value_type alpha, BLAS_Cpp::Transp trans_rhs1, const Vector &v_rhs2, value_type beta) const
Implementation of limited Memory BFGS matrix for arbitrary vector spaces.
size_type num_secant_updates() const
Returns the total number of successful secant updates performed since this->init_identity() was calle...
Transp