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_MatrixSymAddDelBunchKaufman.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_BUNCH_KAUFMAN_H
43 #define MATRIX_SYM_POS_DEF_BUNCH_KAUFMAN_H
44 
45 #include <vector>
46 
47 #include "ConstrainedOptPack_MatrixSymAddDelUpdateableWithOpNonsingular.hpp"
48 #include "AbstractLinAlgPack_MatrixSymAddDelUpdateable.hpp"
49 #include "AbstractLinAlgPack_MatrixSymPosDefCholFactor.hpp"
50 #include "AbstractLinAlgPack_MatrixSymOpNonsingSerial.hpp"
51 #include "DenseLinAlgPack_DMatrixAsTriSym.hpp"
52 
53 namespace ConstrainedOptPack {
54 
67  :public virtual MatrixSymOpNonsingSerial
68  ,public virtual MatrixSymAddDelUpdateable
70 {
71 public:
72 
75 
77  void pivot_tols( PivotTolerances pivot_tols );
79  PivotTolerances pivot_tols() const;
80 
83 
85  const MatrixSymOpNonsing& op_interface() const;
87  MatrixSymAddDelUpdateable& update_interface();
89  const MatrixSymAddDelUpdateable& update_interface() const;
90 
92 
95 
97  void initialize(
98  value_type alpha
100  );
102  void initialize(
103  const DMatrixSliceSym &A
104  ,size_type max_size
105  ,bool force_factorization
106  ,Inertia inertia
107  ,PivotTolerances pivot_tols
108  );
110  size_type max_size() const;
112  Inertia inertia() const;
114  void set_uninitialized();
116  void augment_update(
117  const DVectorSlice *t
118  ,value_type alpha
119  ,bool force_refactorization
120  ,EEigenValType add_eigen_val
121  ,PivotTolerances pivot_tols
122  );
124  void delete_update(
125  size_type jd
126  ,bool force_refactorization
127  ,EEigenValType drop_eigen_val
128  ,PivotTolerances pivot_tols
129  );
130 
132 
135 
137  size_type rows() const;
139  std::ostream& output(std::ostream& out) const;
141  void Vp_StMtV(
142  DVectorSlice* vs_lhs, value_type alpha, BLAS_Cpp::Transp trans_rhs1
143  ,const DVectorSlice& vs_rhs2, value_type beta
144  ) const;
146  void V_InvMtV(
147  DVectorSlice* vs_lhs, BLAS_Cpp::Transp trans_rhs1
148  ,const DVectorSlice& vs_rhs2
149  )const;
150 
151 private:
152 
153  // /////////////////////////////////////////////////////
154  // Private types
155 
156  typedef std::vector<FortranTypes::f_int> IPIV_t;
157 
158  // /////////////////////////////////////////////////////
159  // Private data members
160 
161  size_type S_size_; // The size of the current symmetric matrix. Size == 0 if flag for uninitialized.
162  bool S_indef_; // True if the current matrix is indefinite.
163  bool fact_updated_;// True if the factorization for the current matrix is updated. Only meaningful
164  // if S_indef_==true.
165  bool fact_in1_; // If true then the current factorization is in S_store1_
166  // otherwise it is in S_store2_. Only meaningful if S_indef_==true and
167  // fact_updated_==true
168  MatrixSymAddDelUpdateable::Inertia
169  inertia_; // Inertial for the indefinite L*D*L' factorization. If fact_updated_ == false
170  // then this will be UNKNOWN. IF S_indef_==false then this is meaningless.
171  DMatrix S_store1_; // Storage for the factored matrix in the
172  // upper triangle as well as the original matrix
173  // in the lower triangle. This uses the same storage scheme as
174  // in MatrixSymPosDefCholFactor.
175  DMatrix S_store2_; // Storage for the factorization also. We need
176  // two storage locations for the L*D*L factorization
177  // in case an update is singular. This will not
178  // be initialized for a p.d. or n.d. matrix.
179  IPIV_t IPIV_; // Stores permutations computed by LAPACK
180  mutable DVector
181  WORK_; // workspace
182  MatrixSymPosDefCholFactor
183  S_chol_; // Used to factorize the matrix
184  // when it is p.d. or n.d.
185 
186  // /////////////////////////////////////////////////////
187  // Private member funcitons.
188 
191  DMatrixSliceTriEle DU(size_type S_size, bool fact_in1);
193  const DMatrixSliceTriEle DU(size_type S_size, bool fact_in1) const;
196  DMatrixSliceSym S(size_type S_size);
198  const DMatrixSliceSym S(size_type S_size) const;
200  void assert_initialized() const;
202  void resize_DU_store( bool in_store1 );
207  void copy_and_factor_matrix( size_type S_size, bool fact_in1 );
212  void factor_matrix( size_type S_size, bool fact_in1 );
219  bool compute_assert_inertia(
220  size_type S_size, bool fact_in1
221  ,const Inertia& expected_inertia, const char func_name[]
222  ,PivotTolerances pivot_tols, Inertia* comp_inertia, std::ostringstream* err_msg, value_type* gamma );
223 
227 
228 }; // end class MatrixSymAddDelBunchKaufman
229 
230 // ///////////////////////////
231 // Inline member functions
232 
233 inline
234 DMatrixSliceTriEle MatrixSymAddDelBunchKaufman::DU(size_type S_size, bool fact_in1)
235 {
236  resize_DU_store(fact_in1);
237  return DenseLinAlgPack::nonconst_tri_ele(
238  ( fact_in1 ? S_store1_ : S_store2_ )(1,S_size,2,S_size+1)
239  ,BLAS_Cpp::upper );
240 }
241 
242 inline
243 const DMatrixSliceTriEle MatrixSymAddDelBunchKaufman::DU(size_type S_size, bool fact_in1) const
244 {
245  return DenseLinAlgPack::tri_ele(
246  ( fact_in1 ? S_store1_ : S_store2_ )(1,S_size,2,S_size+1)
247  ,BLAS_Cpp::upper);
248 }
249 
250 inline
251 DMatrixSliceSym MatrixSymAddDelBunchKaufman::S(size_type S_size)
252 {
253  return DenseLinAlgPack::nonconst_sym(
254  S_store1_(2,S_size+1,1,S_size)
255  , BLAS_Cpp::lower );
256 }
257 
258 inline
259 const DMatrixSliceSym MatrixSymAddDelBunchKaufman::S(size_type S_size) const
260 {
261  return DenseLinAlgPack::sym(
262  S_store1_(2,S_size+1,1,S_size)
263  , BLAS_Cpp::lower );
264 }
265 
266 } // namespace ConstrainedOptPack
267 
268 #endif // MATRIX_SYM_POS_DEF_BUNCH_KAUFMAN_H
void Vp_StMtV(DVectorSlice *vs_lhs, value_type alpha, BLAS_Cpp::Transp trans_rhs1, const DVectorSlice &vs_rhs2, value_type beta) const
This class maintains the factorization of symmetric indefinite matrix using a Bunch & Kaufman factori...
Interface for updating a symmetric matrix and its factorization by adding and deleting rows and colum...
size_t size_type
void augment_update(const DVectorSlice *t, value_type alpha, bool force_refactorization, EEigenValType add_eigen_val, PivotTolerances pivot_tols)
MatrixSymAddDelBunchKaufman()
Initializes with 0x0 and pivot_tols == (0.0,0.0,0.0).
Transp
void delete_update(size_type jd, bool force_refactorization, EEigenValType drop_eigen_val, PivotTolerances pivot_tols)
void V_InvMtV(DVectorSlice *vs_lhs, BLAS_Cpp::Transp trans_rhs1, const DVectorSlice &vs_rhs2) const