AbstractLinAlgPack: C++ Interfaces For Vectors, Matrices And Related Linear Algebra Objects  Version of the Day
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Groups Pages
AbstractLinAlgPack_MatrixSymPosDefCholFactor.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_CHOL_FACTOR_H
43 #define MATRIX_SYM_POS_DEF_CHOL_FACTOR_H
44 
45 #include "AbstractLinAlgPack_Types.hpp"
46 #include "AbstractLinAlgPack_MatrixExtractInvCholFactor.hpp"
47 #include "AbstractLinAlgPack_MatrixSymAddDelUpdateable.hpp"
48 #include "AbstractLinAlgPack_MatrixSymOpNonsingSerial.hpp"
49 #include "AbstractLinAlgPack_MatrixSymDenseInitialize.hpp"
50 #include "AbstractLinAlgPack_MatrixSymOpGetGMSSymMutable.hpp"
51 #include "AbstractLinAlgPack_MatrixSymSecant.hpp"
52 #include "DenseLinAlgPack_DMatrixClass.hpp"
53 #include "DenseLinAlgPack_DMatrixAsTriSym.hpp"
54 #include "SerializationPack_Serializable.hpp"
55 #include "Teuchos_RCP.hpp"
56 #include "ReleaseResource.hpp"
57 
58 namespace AbstractLinAlgPack {
119 class MatrixSymPosDefCholFactor
120  : virtual public AbstractLinAlgPack::MatrixSymOpNonsingSerial // doxygen needs full name
121  , virtual public AbstractLinAlgPack::MatrixSymDenseInitialize // ""
123  , virtual public MatrixExtractInvCholFactor
124  , virtual public MatrixSymSecant
125  , virtual public MatrixSymAddDelUpdateable
126  , virtual public SerializationPack::Serializable
127 {
128 public:
129 
132 
134  typedef Teuchos::RCP<
135  MemMngPack::ReleaseResource> release_resource_ptr_t;
136 
139  class PostMod {
140  public:
142  PostMod(
143  bool maintain_original = true
144  ,bool maintain_factor = false
145  ,bool allow_factor = true
146  )
147  :maintain_original_(maintain_original)
148  ,maintain_factor_(maintain_factor)
149  ,allow_factor_(allow_factor)
150  {}
152  void initialize(MatrixSymPosDefCholFactor* p) const
153  {
154  p->init_setup(
155  NULL,Teuchos::null,0 // Allocate own storage
156  ,maintain_original_,maintain_factor_,allow_factor_
157  );
158  }
159  private:
160  bool maintain_original_;
161  bool maintain_factor_;
162  bool allow_factor_;
163  }; // end PostMod
164 
166 
169 
176  MatrixSymPosDefCholFactor();
177 
182  MatrixSymPosDefCholFactor(
183  DMatrixSlice *MU_store
184  ,const release_resource_ptr_t& release_resource_ptr = Teuchos::null
185  ,size_type max_size = 0
186  ,bool maintain_original = true
187  ,bool maintain_factor = false
188  ,bool allow_factor = true
189  ,bool set_full_view = true
190  ,value_type scale = 1.0
191  );
192 
246  void init_setup(
247  DMatrixSlice *MU_store
248  ,const release_resource_ptr_t& release_resource_ptr = Teuchos::null
249  ,size_type max_size = 0
250  ,bool maintain_original = true
251  ,bool maintain_factor = false
252  ,bool allow_factor = true
253  ,bool set_full_view = true
254  ,value_type scale = 1.0
255  );
299  void set_view(
300  size_t M_size
301  ,value_type scale
302  ,bool maintain_original
303  ,size_t M_l_r
304  ,size_t M_l_c
305  ,bool maintain_factor
306  ,size_t U_l_r
307  ,size_t U_l_c
308  );
311  void pivot_tols( PivotTolerances pivot_tols );
313  PivotTolerances pivot_tols() const;
314 
316 
319 
322  void get_view_setup(
323  size_t *M_size
324  ,value_type *scale
325  ,bool *maintain_original
326  ,size_t *M_l_r
327  ,size_t *M_l_c
328  ,bool *maintain_factor
329  ,size_t *U_l_r
330  ,size_t *U_l_c
331  ) const;
334  bool allocates_storage() const;
335 
338  DMatrixSlice& MU_store();
340  const DMatrixSlice& MU_store() const;
343  DMatrixSliceTri U();
345  const DMatrixSliceTri U() const;
348  DMatrixSliceSym M();
350  const DMatrixSliceSym M() const;
351 
353 
356 
358  size_type rows() const;
359 
361 
364 
366  void zero_out();
368  std::ostream& output(std::ostream& out) const;
370  bool Mp_StM(
371  MatrixOp* m_lhs, value_type alpha
372  ,BLAS_Cpp::Transp trans_rhs
373  ) const;
375  bool Mp_StM(
376  value_type alpha,const MatrixOp& M_rhs, BLAS_Cpp::Transp trans_rhs
377  );
379  bool syrk(
380  const MatrixOp &mwo_rhs
381  ,BLAS_Cpp::Transp M_trans
382  ,value_type alpha
383  ,value_type beta
384  );
385 
387 
390 
392  void Vp_StMtV(DVectorSlice* vs_lhs, value_type alpha, BLAS_Cpp::Transp trans_rhs1
393  , const DVectorSlice& vs_rhs2, value_type beta) const;
395  void Vp_StMtV(DVectorSlice* vs_lhs, value_type alpha, BLAS_Cpp::Transp trans_rhs1
396  , const SpVectorSlice& vs_rhs2, value_type beta) const;
398  void Vp_StPtMtV(DVectorSlice* vs_lhs, value_type alpha
399  , const GenPermMatrixSlice& P_rhs1, BLAS_Cpp::Transp P_rhs1_trans
400  , BLAS_Cpp::Transp M_rhs2_trans
401  , const DVectorSlice& vs_rhs3, value_type beta) const;
403  void Vp_StPtMtV(DVectorSlice* vs_lhs, value_type alpha
404  , const GenPermMatrixSlice& P_rhs1, BLAS_Cpp::Transp P_rhs1_trans
405  , BLAS_Cpp::Transp M_rhs2_trans
406  , const SpVectorSlice& sv_rhs3, value_type beta) const;
407 
409 
412 
413  void Mp_StPtMtP( DMatrixSliceSym* sym_lhs, value_type alpha
414  , EMatRhsPlaceHolder dummy_place_holder
415  , const GenPermMatrixSlice& gpms_rhs, BLAS_Cpp::Transp gpms_rhs_trans
416  , value_type beta ) const;
417 
419 
422 
424  void V_InvMtV(DVectorSlice* vs_lhs, BLAS_Cpp::Transp trans_rhs1
425  , const DVectorSlice& vs_rhs2) const;
427  void V_InvMtV(DVectorSlice* vs_lhs, BLAS_Cpp::Transp trans_rhs1
428  , const SpVectorSlice& sv_rhs2) const;
429 
431 
434 
436  void M_StMtInvMtM(
437  DMatrixSliceSym* sym_gms_lhs, value_type alpha
438  , const MatrixOpSerial& mwo, BLAS_Cpp::Transp mwo_trans, EMatrixDummyArg
439  ) const;
440 
442 
445 
447  void initialize( const DMatrixSliceSym& M );
448 
450 
453 
455  const DenseLinAlgPack::DMatrixSliceSym get_sym_gms_view() const;
457  void free_sym_gms_view(const DenseLinAlgPack::DMatrixSliceSym* sym_gms_view) const;
458 
460 
463 
465  DenseLinAlgPack::DMatrixSliceSym get_sym_gms_view();
467  void commit_sym_gms_view(DenseLinAlgPack::DMatrixSliceSym* sym_gms_view);
468 
470 
473 
475  void extract_inv_chol( DMatrixSliceTriEle* InvChol ) const;
476 
478 
481 
483  void init_identity( const VectorSpace& space_diag, value_type alpha );
485  void init_diagonal( const Vector& diag );
487  void secant_update(VectorMutable* s, VectorMutable* y, VectorMutable* Bs);
488 
490 
493 
495  void initialize(
496  value_type alpha
497  ,size_type max_size
498  );
500  void initialize(
501  const DMatrixSliceSym &A
502  ,size_type max_size
503  ,bool force_factorization
504  ,Inertia inertia
505  ,PivotTolerances pivot_tols
506  );
508  size_type max_size() const;
510  Inertia inertia() const;
512  void set_uninitialized();
514  void augment_update(
515  const DVectorSlice *t
516  ,value_type alpha
517  ,bool force_refactorization
518  ,EEigenValType add_eigen_val
519  ,PivotTolerances pivot_tols
520  );
522  void delete_update(
523  size_type jd
524  ,bool force_refactorization
525  ,EEigenValType drop_eigen_val
526  ,PivotTolerances pivot_tols
527  );
528 
530 
531 
534 
536  void serialize( std::ostream &out ) const;
538  void unserialize( std::istream &in );
539 
541 
542 private:
543 
544  // /////////////////////////////
545  // Private data members
546 
547  bool maintain_original_;
548  bool maintain_factor_;
549  bool factor_is_updated_; // Is set to true if maintain_factor_=true
550  bool allocates_storage_; // If true then this object allocates the storage
551  release_resource_ptr_t release_resource_ptr_;
552  DMatrixSlice MU_store_;
553  size_t max_size_;
554  size_t M_size_, // M_size == 0 is flag that we are completely uninitialized
555  M_l_r_,
556  M_l_c_,
557  U_l_r_,
558  U_l_c_;
559  value_type scale_;
560  bool is_diagonal_;
561  PivotTolerances pivot_tols_;
562  DVector work_; // workspace.
563 
564  // /////////////////////////////
565  // Private member functions
566 
567  void assert_storage() const;
568  void allocate_storage(size_type max_size) const;
569  void assert_initialized() const;
570  void resize_and_zero_off_diagonal(size_type n, value_type scale);
571  void update_factorization() const;
572  std::string build_serialization_string() const;
573  static void write_matrix( const DMatrixSlice &Q, BLAS_Cpp::Uplo Q_uplo, std::ostream &out );
574  static void read_matrix( std::istream &in, BLAS_Cpp::Uplo Q_uplo, DMatrixSlice *Q );
575 
576 }; // end class MatrixSymPosDefCholFactor
577 
578 // ///////////////////////////////////////////////////////
579 // Inline members for MatrixSymPosDefCholFactor
580 
581 inline
582 bool MatrixSymPosDefCholFactor::allocates_storage() const
583 {
584  return allocates_storage_;
585 }
586 
587 inline
588 DMatrixSlice& MatrixSymPosDefCholFactor::MU_store()
589 {
590  return MU_store_;
591 }
592 
593 inline
594 const DMatrixSlice& MatrixSymPosDefCholFactor::MU_store() const
595 {
596  return MU_store_;
597 }
598 
599 inline
600 void MatrixSymPosDefCholFactor::get_view_setup(
601  size_t *M_size
602  ,value_type *scale
603  ,bool *maintain_original
604  ,size_t *M_l_r
605  ,size_t *M_l_c
606  ,bool *maintain_factor
607  ,size_t *U_l_r
608  ,size_t *U_l_c
609  ) const
610 {
611  *M_size = M_size_;
612  *scale = scale_;
613  *maintain_original = maintain_original_;
614  *M_l_r = maintain_original_ ? M_l_r_ : 0;
615  *M_l_c = maintain_original_ ? M_l_c_ : 0;
616  *maintain_factor = maintain_factor_;
617  *U_l_r = maintain_factor_ ? U_l_r_ : 0;
618  *U_l_c = maintain_factor_ ? U_l_c_ : 0;
619 }
620 
621 inline
622 DMatrixSliceTri MatrixSymPosDefCholFactor::U()
623 {
624  return DenseLinAlgPack::nonconst_tri(
625  MU_store_(U_l_r_,U_l_r_+M_size_-1,U_l_c_+1,U_l_c_+M_size_)
626  , BLAS_Cpp::upper, BLAS_Cpp::nonunit
627  );
628 }
629 
630 inline
631 const DMatrixSliceTri MatrixSymPosDefCholFactor::U() const
632 {
633  return DenseLinAlgPack::tri(
634  MU_store_(U_l_r_,U_l_r_+M_size_-1,U_l_c_+1,U_l_c_+M_size_)
635  , BLAS_Cpp::upper, BLAS_Cpp::nonunit
636  );
637 }
638 
639 inline
640 DMatrixSliceSym MatrixSymPosDefCholFactor::M()
641 {
642  return DenseLinAlgPack::nonconst_sym(
643  MU_store_(M_l_r_+1,M_l_r_+M_size_,M_l_c_,M_l_c_+M_size_-1)
644  , BLAS_Cpp::lower
645  );
646 }
647 
648 inline
649 const DMatrixSliceSym MatrixSymPosDefCholFactor::M() const
650 {
651  return DenseLinAlgPack::sym(
652  MU_store_(M_l_r_+1,M_l_r_+M_size_,M_l_c_,M_l_c_+M_size_-1)
653  , BLAS_Cpp::lower
654  );
655 }
656 
657 } // end namespace AbstractLinAlgPack
658 
659 #endif // MATRIX_SYM_POS_DEF_CHOL_FACTOR_H
Uplo
void Vp_StPtMtV(VectorMutable *v_lhs, value_type alpha, const GenPermMatrixSlice &P_rhs1, BLAS_Cpp::Transp P_rhs1_trans, const MatrixOp &M_rhs2, BLAS_Cpp::Transp M_rhs2_trans, const Vector &v_rhs3, value_type beta=1.0)
v_lhs = alpha * op(P_rhs1) * op(M_rhs2) * v_rhs3 + beta * v_rhs
size_type rows(size_type rows, size_type cols, BLAS_Cpp::Transp _trans)
Abstract interface that allows the extraction of a non-const DenseLinAlgPack::DMatrixSliceSym view of...
Abstract base class for all serial polymorphic symmetric nonsingular matrices that can be used to com...
void Mp_StPtMtP(MatrixOp *mwo_lhs, value_type alpha, const GenPermMatrixSlice &P_rhs1, BLAS_Cpp::Transp P_rhs1_trans, const MatrixOp &M_rhs, BLAS_Cpp::Transp trans_rhs, const GenPermMatrixSlice &P_rhs2, BLAS_Cpp::Transp P_rhs2_trans)
mwo_lhs += alpha * op(P_rhs1) * op(M_rhs) * op(P_rhs2).
size_t size_type
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:
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
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 M_StMtInvMtM(MatrixSymOp *sym_gms_lhs, value_type alpha, const MatrixOp &mwo, BLAS_Cpp::Transp mwo_trans, const MatrixSymNonsing &mswof, MatrixSymNonsing::EMatrixDummyArg mwo_rhs)
sym_gms_lhs = alpha * op(mwo) * inv(mswof) * op(mwo)'
Transp
Mix-in Interface for initializing a matrix with a dense symmetric matrix.