ConstrainedOptPack: C++ Tools for Constrained (and Unconstrained) Optimization Version of the Day
ConstrainedOptPack_MatrixSymAddDelBunchKaufman.hpp
00001 // @HEADER
00002 // ***********************************************************************
00003 // 
00004 // Moocho: Multi-functional Object-Oriented arCHitecture for Optimization
00005 //                  Copyright (2003) Sandia Corporation
00006 // 
00007 // Under terms of Contract DE-AC04-94AL85000, there is a non-exclusive
00008 // license for use of this work by or on behalf of the U.S. Government.
00009 // 
00010 // Redistribution and use in source and binary forms, with or without
00011 // modification, are permitted provided that the following conditions are
00012 // met:
00013 //
00014 // 1. Redistributions of source code must retain the above copyright
00015 // notice, this list of conditions and the following disclaimer.
00016 //
00017 // 2. Redistributions in binary form must reproduce the above copyright
00018 // notice, this list of conditions and the following disclaimer in the
00019 // documentation and/or other materials provided with the distribution.
00020 //
00021 // 3. Neither the name of the Corporation nor the names of the
00022 // contributors may be used to endorse or promote products derived from
00023 // this software without specific prior written permission.
00024 //
00025 // THIS SOFTWARE IS PROVIDED BY SANDIA CORPORATION "AS IS" AND ANY
00026 // EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
00027 // IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR
00028 // PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL SANDIA CORPORATION OR THE
00029 // CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL,
00030 // EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO,
00031 // PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR
00032 // PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF
00033 // LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING
00034 // NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
00035 // SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
00036 //
00037 // Questions? Contact Roscoe A. Bartlett (rabartl@sandia.gov) 
00038 // 
00039 // ***********************************************************************
00040 // @HEADER
00041 
00042 #ifndef MATRIX_SYM_POS_DEF_BUNCH_KAUFMAN_H
00043 #define MATRIX_SYM_POS_DEF_BUNCH_KAUFMAN_H
00044 
00045 #include <vector>
00046 
00047 #include "ConstrainedOptPack_MatrixSymAddDelUpdateableWithOpNonsingular.hpp"
00048 #include "AbstractLinAlgPack_MatrixSymAddDelUpdateable.hpp"
00049 #include "AbstractLinAlgPack_MatrixSymPosDefCholFactor.hpp"
00050 #include "AbstractLinAlgPack_MatrixSymOpNonsingSerial.hpp"
00051 #include "DenseLinAlgPack_DMatrixAsTriSym.hpp"
00052 
00053 namespace ConstrainedOptPack {
00054 
00066 class MatrixSymAddDelBunchKaufman
00067   :public virtual MatrixSymOpNonsingSerial
00068   ,public virtual MatrixSymAddDelUpdateable
00069   ,public virtual MatrixSymAddDelUpdateableWithOpNonsingular
00070 {
00071 public:
00072 
00074   MatrixSymAddDelBunchKaufman();
00075 
00077   void pivot_tols( PivotTolerances pivot_tols );
00079   PivotTolerances pivot_tols() const;
00080 
00083 
00085   const MatrixSymOpNonsing& op_interface() const;
00087   MatrixSymAddDelUpdateable& update_interface();
00089   const MatrixSymAddDelUpdateable& update_interface() const;
00090 
00092 
00095 
00097   void initialize(
00098     value_type         alpha
00099     ,size_type         max_size
00100     );
00102   void initialize(
00103     const DMatrixSliceSym      &A
00104     ,size_type         max_size
00105     ,bool              force_factorization
00106     ,Inertia           inertia
00107     ,PivotTolerances   pivot_tols
00108     );
00110   size_type max_size() const;
00112   Inertia inertia() const;
00114   void set_uninitialized();
00116   void augment_update(
00117     const DVectorSlice  *t
00118     ,value_type        alpha
00119     ,bool              force_refactorization
00120     ,EEigenValType     add_eigen_val
00121     ,PivotTolerances   pivot_tols
00122     );
00124   void delete_update(
00125     size_type          jd
00126     ,bool              force_refactorization
00127     ,EEigenValType     drop_eigen_val
00128     ,PivotTolerances   pivot_tols
00129     );
00130 
00132 
00135 
00137   size_type rows() const;
00139   std::ostream& output(std::ostream& out) const;
00141   void Vp_StMtV(
00142     DVectorSlice* vs_lhs, value_type alpha, BLAS_Cpp::Transp trans_rhs1
00143     ,const DVectorSlice& vs_rhs2, value_type beta
00144     ) const;
00146   void V_InvMtV(
00147     DVectorSlice* vs_lhs, BLAS_Cpp::Transp trans_rhs1
00148     ,const DVectorSlice& vs_rhs2
00149     )const;
00150 
00151 private:
00152 
00153   // /////////////////////////////////////////////////////
00154   // Private types
00155 
00156   typedef std::vector<FortranTypes::f_int> IPIV_t;
00157 
00158   // /////////////////////////////////////////////////////
00159   // Private data members
00160 
00161   size_type S_size_;      // The size of the current symmetric matrix.  Size == 0 if flag for uninitialized.
00162   bool      S_indef_;     // True if the current matrix is indefinite.
00163   bool      fact_updated_;// True if the factorization for the current matrix is updated.  Only meaningful
00164                           // if S_indef_==true.
00165   bool      fact_in1_;    // If true then the current factorization is in S_store1_
00166                           // otherwise it is in S_store2_.  Only meaningful if S_indef_==true and
00167                           // fact_updated_==true
00168   MatrixSymAddDelUpdateable::Inertia
00169             inertia_;     // Inertial for the indefinite L*D*L' factorization.  If fact_updated_ == false
00170                           // then this will be UNKNOWN.  IF S_indef_==false then this is meaningless.
00171   DMatrix S_store1_;    // Storage for the factored matrix in the
00172                           // upper triangle as well as the original matrix
00173                           // in the lower triangle.  This uses the same storage scheme as
00174                           // in MatrixSymPosDefCholFactor.
00175   DMatrix S_store2_;    // Storage for the factorization also.  We need
00176                           // two storage locations for the L*D*L factorization
00177                           // in case an update is singular.  This will not
00178                           // be initialized for a p.d. or n.d. matrix.
00179   IPIV_t    IPIV_;        // Stores permutations computed by LAPACK
00180   mutable DVector
00181               WORK_;        // workspace
00182   MatrixSymPosDefCholFactor
00183           S_chol_;      // Used to factorize the matrix
00184                           // when it is p.d. or n.d.
00185 
00186   // /////////////////////////////////////////////////////
00187   // Private member funcitons.
00188 
00191   DMatrixSliceTriEle DU(size_type S_size, bool fact_in1);
00193   const DMatrixSliceTriEle DU(size_type S_size, bool fact_in1) const;
00196   DMatrixSliceSym S(size_type S_size);
00198   const DMatrixSliceSym S(size_type S_size) const;
00200   void assert_initialized() const;
00202   void resize_DU_store( bool in_store1 );
00207   void copy_and_factor_matrix( size_type S_size, bool fact_in1 );
00212   void factor_matrix( size_type S_size, bool fact_in1 );
00219   bool compute_assert_inertia(
00220     size_type S_size, bool fact_in1
00221     ,const Inertia& expected_inertia, const char func_name[]
00222     ,PivotTolerances pivot_tols, Inertia* comp_inertia, std::ostringstream* err_msg, value_type* gamma );
00223 
00225   MatrixSymAddDelBunchKaufman( const MatrixSymAddDelBunchKaufman& );
00226   MatrixSymAddDelBunchKaufman& operator=( const MatrixSymAddDelBunchKaufman& );
00227 
00228 };  // end class MatrixSymAddDelBunchKaufman
00229 
00230 // ///////////////////////////
00231 // Inline member functions
00232 
00233 inline
00234 DMatrixSliceTriEle MatrixSymAddDelBunchKaufman::DU(size_type S_size, bool fact_in1)
00235 {
00236   resize_DU_store(fact_in1);
00237   return DenseLinAlgPack::nonconst_tri_ele(
00238     ( fact_in1 ? S_store1_ : S_store2_ )(1,S_size,2,S_size+1)
00239     ,BLAS_Cpp::upper );
00240 }
00241 
00242 inline
00243 const DMatrixSliceTriEle MatrixSymAddDelBunchKaufman::DU(size_type S_size, bool fact_in1) const
00244 {
00245   return DenseLinAlgPack::tri_ele(
00246     ( fact_in1 ? S_store1_ : S_store2_ )(1,S_size,2,S_size+1)
00247     ,BLAS_Cpp::upper);
00248 }
00249 
00250 inline
00251 DMatrixSliceSym MatrixSymAddDelBunchKaufman::S(size_type S_size)
00252 {
00253   return DenseLinAlgPack::nonconst_sym(
00254     S_store1_(2,S_size+1,1,S_size)
00255     , BLAS_Cpp::lower );
00256 }
00257 
00258 inline
00259 const DMatrixSliceSym MatrixSymAddDelBunchKaufman::S(size_type S_size) const
00260 {
00261   return DenseLinAlgPack::sym(
00262     S_store1_(2,S_size+1,1,S_size)
00263     , BLAS_Cpp::lower );
00264 }
00265 
00266 } // namespace ConstrainedOptPack 
00267 
00268 #endif  // MATRIX_SYM_POS_DEF_BUNCH_KAUFMAN_H
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends