ConstrainedOptPack: C++ Tools for Constrained (and Unconstrained) Optimization Version of the Day
ConstrainedOptPack_MatrixSymAddDelBunchKaufman.cpp
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 
00043 #include <assert.h>
00044 
00045 #include <limits>
00046 
00047 #include "ConstrainedOptPack_MatrixSymAddDelBunchKaufman.hpp"
00048 #include "DenseLinAlgLAPack.hpp"
00049 #include "DenseLinAlgPack_DMatrixOut.hpp"
00050 #include "DenseLinAlgPack_DMatrixOp.hpp"
00051 #include "DenseLinAlgPack_AssertOp.hpp"
00052 #include "DenseLinAlgPack_delete_row_col.hpp"
00053 
00054 namespace ConstrainedOptPack {
00055 
00056 MatrixSymAddDelBunchKaufman::MatrixSymAddDelBunchKaufman()
00057   : S_size_(0), S_indef_(false), fact_updated_(false), fact_in1_(true), inertia_(0,0,0)
00058 {}
00059 
00060 void MatrixSymAddDelBunchKaufman::pivot_tols( PivotTolerances pivot_tols )
00061 {
00062   S_chol_.pivot_tols(pivot_tols);
00063 }
00064 
00065 MatrixSymAddDelUpdateable::PivotTolerances MatrixSymAddDelBunchKaufman::pivot_tols() const
00066 {
00067   return S_chol_.pivot_tols();
00068 }
00069 
00070 // Overridden from MatrixSymAddDelUpdateableWithOpNonsingular
00071 
00072 const MatrixSymOpNonsing& MatrixSymAddDelBunchKaufman::op_interface() const
00073 {
00074   return *this;
00075 }
00076 
00077 MatrixSymAddDelUpdateable& MatrixSymAddDelBunchKaufman::update_interface()
00078 {
00079   return *this;
00080 }
00081 
00082 const MatrixSymAddDelUpdateable& MatrixSymAddDelBunchKaufman::update_interface() const
00083 {
00084   return *this;
00085 }
00086 
00087 // Overridden from MatrixSymAddDelUpdateable
00088 
00089 void MatrixSymAddDelBunchKaufman::initialize(
00090   value_type         alpha
00091   ,size_type         max_size
00092   )
00093 {
00094   try {
00095     // Resize the storage if we have to
00096     if( S_store1_.rows() < max_size+1 && S_store1_.cols() < max_size+1 )
00097       S_store1_.resize(max_size+1,max_size+1);
00098     fact_in1_ = true;
00099     // Start out with a p.d. or n.d. matrix and maintain the original
00100     S_chol_.init_setup(&S_store1_(),Teuchos::null,0,true,true,true,false,0.0);
00101     S_chol_.initialize(alpha,max_size);
00102     // Set the state variables:
00103     S_size_  = 1;
00104     S_indef_ = false; // fact_updated_, fact_in1_ and inertia are meaningless!
00105   }
00106   catch(...) {
00107     S_size_ = 0;
00108     throw;
00109   }
00110 }
00111 
00112 void MatrixSymAddDelBunchKaufman::initialize(
00113   const DMatrixSliceSym      &A
00114   ,size_type         max_size
00115   ,bool              force_factorization
00116   ,Inertia           expected_inertia
00117   ,PivotTolerances   pivot_tols
00118   )
00119 {
00120   using BLAS_Cpp::upper;
00121   using BLAS_Cpp::lower;
00122   using DenseLinAlgPack::tri_ele;
00123   using DenseLinAlgPack::nonconst_tri_ele;
00124   typedef MatrixSymAddDelUpdateable  MSADU;
00125   typedef MSADU::Inertia   Inertia;
00126 
00127   bool                throw_exception = false; // If true then throw exception
00128   std::ostringstream  omsg;                    // Will be set if an exception has to be thrown.
00129   value_type          gamma;                   // ...
00130 
00131   const size_type
00132     n = A.rows();
00133 
00134   // Validate proper usage of inertia parameter
00135   TEUCHOS_TEST_FOR_EXCEPT( ! ( expected_inertia.zero_eigens == Inertia::UNKNOWN
00136     || expected_inertia.zero_eigens == 0 ) );
00137   
00138   try {
00139     // Resize the storage if we have to
00140     if( S_store1_.rows() < max_size+1 && S_store1_.cols() < max_size+1 )
00141       S_store1_.resize(max_size+1,max_size+1);
00142     fact_in1_ = true;
00143     // See if the client claims that the matrix is p.d. or n.d.
00144     const bool not_indefinite
00145       = ( ( expected_inertia.neg_eigens == 0 && expected_inertia.pos_eigens == n )
00146         || ( expected_inertia.neg_eigens == n && expected_inertia.pos_eigens == 0 ) );
00147     // Initialize the matrix
00148     if( not_indefinite ) {
00149       // The client says that the matrix is p.d. or n.d. so
00150       // we will take their word for it.
00151       S_chol_.init_setup(&S_store1_(),Teuchos::null,0,true,true,true,false,0.0);
00152       try {
00153         S_chol_.initialize(A,max_size,force_factorization,expected_inertia,pivot_tols);
00154       }
00155       catch(const MSADU::WarnNearSingularUpdateException& except) {
00156         throw_exception = true; // Throw this same exception at the end!
00157         omsg << except.what();
00158         gamma = except.gamma;
00159       }
00160       // Set the state variables:
00161       S_size_       = n;
00162       S_indef_      = false; // fact_updated_, fact_in1_ and inertia are meaningless!
00163     }
00164     else {
00165       //
00166       // The client did not say that the matrix was p.d. or n.d. so we
00167       // must assume that it is indefinite.
00168       //
00169       bool fact_in1 = !fact_in1_;
00170       // Set the new matrix in the unused factorization location
00171       DenseLinAlgPack::assign(  &DU(n,fact_in1), tri_ele( A.gms(), A.uplo() ) );
00172       // Update the factorization in place
00173       try {
00174         factor_matrix( n, fact_in1 );
00175       }
00176       catch( const DenseLinAlgLAPack::FactorizationException& excpt ) {
00177         omsg
00178           << "MatrixSymAddDelBunchKaufman::initialize(...): "
00179           << "Error, the matrix A is singular:\n"
00180           << excpt.what();
00181         throw SingularUpdateException( omsg.str(), -1.0 );
00182       }
00183       // Compute the inertia and validate that it is correct.
00184       Inertia inertia;
00185       throw_exception = compute_assert_inertia(
00186         n,fact_in1,expected_inertia,"initialize",pivot_tols
00187         ,&inertia,&omsg,&gamma);
00188       // If the client did not know the inertia of the
00189       // matrix but it turns out to be p.d. or n.d. then modify the
00190       // DU factor appropriatly and switch to cholesky factorization.
00191       if( inertia.pos_eigens == n || inertia.neg_eigens == n ) {
00192         TEUCHOS_TEST_FOR_EXCEPT(true); // ToDo Implememnt this!
00193       }
00194       else {
00195         // Indead the new matrix is indefinite
00196         // Set the original matrix now
00197         DenseLinAlgPack::assign(  &DenseLinAlgPack::nonconst_tri_ele( S(n).gms(), BLAS_Cpp::lower)
00198                   , tri_ele( A.gms(), A.uplo() ) );
00199         // Update the state variables:
00200         S_size_       = n;
00201         S_indef_      = true;
00202         fact_updated_ = true;
00203         fact_in1_     = fact_in1;
00204         inertia_      = inertia;
00205       }
00206     }
00207   }
00208   catch(...) {
00209     S_size_ = 0;
00210     throw;
00211   }
00212   if( throw_exception )
00213     throw WarnNearSingularUpdateException(omsg.str(),gamma);
00214 }
00215 
00216 size_type MatrixSymAddDelBunchKaufman::max_size() const
00217 {
00218   return S_store1_.rows() -1;  // The center diagonal is used for workspace
00219 }
00220 
00221 MatrixSymAddDelUpdateable::Inertia
00222 MatrixSymAddDelBunchKaufman::inertia() const
00223 {
00224   return S_indef_ ? inertia_ : S_chol_.inertia();
00225 }
00226 
00227 void MatrixSymAddDelBunchKaufman::set_uninitialized()
00228 {
00229   S_size_ = 0;
00230 }
00231 
00232 void MatrixSymAddDelBunchKaufman::augment_update(
00233   const DVectorSlice  *t
00234   ,value_type        alpha
00235   ,bool              force_refactorization
00236   ,EEigenValType     add_eigen_val
00237   ,PivotTolerances   pivot_tols
00238   )
00239 {
00240   using BLAS_Cpp::no_trans;
00241   using DenseLinAlgPack::norm_inf;
00242   using AbstractLinAlgPack::transVtInvMtV;
00243   typedef MatrixSymAddDelUpdateable  MSADU;
00244 
00245   assert_initialized();
00246 
00247   // Validate the input
00248   if( rows() >= max_size() )
00249     throw MaxSizeExceededException(
00250       "MatrixSymAddDelBunchKaufman::augment_update(...) : "
00251       "Error, the maximum size would be exceeded." );
00252   if( t && t->dim() != S_size_ )
00253     throw std::length_error(
00254       "MatrixSymAddDelBunchKaufman::augment_update(...): "
00255       "Error, t.dim() must be equal to this->rows()." );
00256   if( !(add_eigen_val == MSADU::EIGEN_VAL_UNKNOWN || add_eigen_val != MSADU::EIGEN_VAL_ZERO) )
00257     throw SingularUpdateException(
00258       "MatrixSymAddDelBunchKaufman::augment_update(...): "
00259       "Error, the client has specified a singular update in add_eigen_val.", -1.0 );
00260 
00261   // Try to perform the update
00262   bool throw_exception = false; // If true then throw exception
00263   std::ostringstream omsg;      // Will be set if an exception has to be thrown.
00264   value_type gamma;             // ...
00265   if( !S_indef_ ) {
00266     // The current matrix is positive definite or negative definite.  If the
00267     // new matrix is still p.d. or n.d. when we are good to go.  We must first
00268     // check if the client has specified that the new matrix will be indefinite
00269     // and if so then we will take his/her word for it and do the indefinite
00270     // updating.
00271     MSADU::Inertia inertia = S_chol_.inertia();
00272     const bool
00273       new_S_not_indef
00274       = ( ( inertia.pos_eigens > 0
00275         && ( add_eigen_val == MSADU::EIGEN_VAL_POS || add_eigen_val == MSADU::EIGEN_VAL_UNKNOWN ) )
00276         || ( inertia.neg_eigens > 0
00277         && ( add_eigen_val == MSADU::EIGEN_VAL_NEG || add_eigen_val == MSADU::EIGEN_VAL_UNKNOWN ) )
00278         );
00279     if( !new_S_not_indef ) {
00280       // We must resort to an indefinite factorization
00281     }
00282     else {
00283       // The new matrix could/should still be p.d. or n.d. so let's try it!
00284       bool update_successful = false;
00285       try {
00286         S_chol_.augment_update(t,alpha,force_refactorization,add_eigen_val,pivot_tols);
00287         update_successful = true;
00288       }
00289       catch(const MSADU::WrongInertiaUpdateException&) {
00290         if( add_eigen_val != MSADU::EIGEN_VAL_UNKNOWN )
00291           throw; // The client specified the new inertia and it was wrong so throw execepiton.
00292         // If the client did not know that inertia then we can't fault them so we will
00293         // proceed unto the indefinite factorization.
00294       }
00295       catch(const MSADU::WarnNearSingularUpdateException& except) {
00296         throw_exception = true; // Throw this same exception at the end!
00297         update_successful = true;
00298         omsg << except.what();
00299         gamma = except.gamma;
00300       }
00301       if( update_successful ) {
00302         ++S_size_;     // Now we can resize the matrix
00303         if(throw_exception)
00304           throw MSADU::WarnNearSingularUpdateException(omsg.str(),gamma);
00305         else
00306           return;
00307       }
00308     }
00309   }
00310   //
00311   // If we get here then we have no choice but the perform an indefinite factorization.
00312   //
00313   // ToDo: (2/2/01): We could also try some more fancy updating
00314   // procedures and try to get away with some potentially unstable
00315   // methods in the interest of time savings.
00316   //
00317   const size_type n     = S_size_;
00318   DMatrixSlice  S_bar = this->S(n+1).gms();
00319   //
00320   // Validate that the new matrix will be nonsingular.
00321   //
00322   // Given any nonsingular matrix S (even unsymmetric) it is easy to show that
00323   // beta = alpha - t'*inv(S)*t != 0.0 if [ S, t; t', alpha ] is nonsingular.
00324   // This is a cheap O(n^2) way to check that the update is nonsingular before
00325   // we go through the O(n^3) refactorization.
00326   // In fact, the sign of beta even tells us the sign of the new eigen value
00327   // of the updated matrix even before we perform the refactorization.
00328   // If the current matrix is not factored then we will just skip this
00329   // test and let the full factorization take place to find this out.
00330   // We could even cheat a little and actually perform the update with this
00331   // diagonal beta and then compute the update to the factors U our selves
00332   // in O(n^2) time.
00333   //
00334   if( force_refactorization && fact_updated_ ) {
00335     const value_type
00336       beta        = alpha - ( t ? transVtInvMtV(*t,*this,no_trans,*t) : 0.0 ),
00337       abs_beta    = std::fabs(beta),
00338       nrm_D_diag  = norm_inf(DU(n,fact_in1_).gms().diag()); // ToDo: Consider 2x2 blocks also!
00339     gamma = abs_beta / nrm_D_diag;
00340     // Check gamma
00341     const bool
00342       correct_inertia = (
00343         add_eigen_val == MSADU::EIGEN_VAL_UNKNOWN
00344         || beta > 0.0 && add_eigen_val == MSADU::EIGEN_VAL_POS
00345         || beta < 0.0 && add_eigen_val == MSADU::EIGEN_VAL_NEG
00346         );
00347     PivotTolerances use_pivot_tols = S_chol_.pivot_tols();
00348     if( pivot_tols.warning_tol != PivotTolerances::UNKNOWN )
00349       use_pivot_tols.warning_tol = pivot_tols.warning_tol;
00350     if( pivot_tols.singular_tol != PivotTolerances::UNKNOWN )
00351       use_pivot_tols.singular_tol = pivot_tols.singular_tol;
00352     if( pivot_tols.wrong_inertia_tol != PivotTolerances::UNKNOWN )
00353       use_pivot_tols.wrong_inertia_tol = pivot_tols.wrong_inertia_tol;
00354     throw_exception = (
00355       gamma == 0.0
00356       || correct_inertia && gamma <= use_pivot_tols.singular_tol
00357       || !correct_inertia
00358       );
00359     // Create message and throw exceptions
00360     std::ostringstream onum_msg;
00361     if(throw_exception) {
00362       onum_msg
00363         << "gamma = |alpha-t'*inv(S)*t)|/||diag(D)||inf = |" <<  beta << "|/" << nrm_D_diag
00364         << " = " << gamma;
00365       omsg
00366         << "MatrixSymAddDelBunchKaufman::augment_update(...): ";
00367       if( correct_inertia && gamma <= use_pivot_tols.singular_tol ) {
00368         omsg
00369           << "Singular update!\n" << onum_msg.str() << " <= singular_tol = "
00370           << use_pivot_tols.singular_tol;
00371         throw SingularUpdateException( omsg.str(), gamma );
00372       }
00373       else if( !correct_inertia && gamma <= use_pivot_tols.singular_tol ) {
00374         omsg
00375           << "Singular update!\n" << onum_msg.str() << " <= wrong_inertia_tol = "
00376           << use_pivot_tols.wrong_inertia_tol;
00377         throw SingularUpdateException( omsg.str(), gamma );
00378       }
00379       
00380       else if( !correct_inertia ) {
00381         omsg
00382           << "Indefinite update!\n" << onum_msg.str() << " >= wrong_inertia_tol = "
00383           << use_pivot_tols.wrong_inertia_tol;
00384         throw WrongInertiaUpdateException( omsg.str(), gamma );
00385       }
00386       else {
00387         TEUCHOS_TEST_FOR_EXCEPT(true); // Only local programming error?
00388       }
00389     }
00390   }
00391 
00392   // Add new row to the lower part of the original symmetric matrix.
00393   if(t)
00394     S_bar.row(n+1)(1,n) = *t;
00395   else
00396     S_bar.row(n+1)(1,n) = 0.0;
00397   S_bar(n+1,n+1) = alpha;
00398 
00399   //
00400   // Update the factorization
00401   //
00402   if( force_refactorization ) {
00403     // Determine where to copy the original matrix to
00404     bool fact_in1 = false;
00405     if( S_indef_ ) {
00406       // S is already indefinite so let's copy the original into the storage
00407       // location other than the current one in case the newly nonsingular matrix
00408       // is singular or has the wrong inertia.
00409       fact_in1 = !fact_in1_;
00410     }
00411     else {
00412       // S is currently p.d. or n.d. so let's copy the new matrix
00413       // into the second storage location so as not to overwrite
00414       // the current cholesky factor in case the new matrix is singular
00415       // or has the wrong inertia.
00416       fact_in1 = false;
00417     }
00418     // Copy and factor the new matrix
00419     try {
00420       copy_and_factor_matrix(n+1,fact_in1);
00421     }
00422     catch( const DenseLinAlgLAPack::FactorizationException& excpt ) {
00423       std::ostringstream omsg;
00424       omsg
00425         << "MatrixSymAddDelBunchKaufman::augment_update(...): "
00426         << "Error, singular update but the original matrix was be maintianed:\n"
00427         << excpt.what();
00428       throw SingularUpdateException( omsg.str(), -1.0 );
00429     }
00430     // Compute the expected inertia
00431     Inertia expected_inertia = this->inertia();
00432     if( expected_inertia.zero_eigens == Inertia::UNKNOWN || add_eigen_val == MSADU::EIGEN_VAL_UNKNOWN )
00433       expected_inertia = Inertia(); // Unknow inertia!
00434     else if( add_eigen_val == MSADU::EIGEN_VAL_NEG )
00435       ++expected_inertia.neg_eigens;
00436     else if( add_eigen_val == MSADU::EIGEN_VAL_POS )
00437       ++expected_inertia.pos_eigens;
00438     else
00439       TEUCHOS_TEST_FOR_EXCEPT(true); // Should not happen!
00440     // Compute the actually inertia and validate that it is what is expected
00441     Inertia inertia;
00442     throw_exception = compute_assert_inertia(
00443       n+1,fact_in1,expected_inertia,"augment_update",pivot_tols
00444       ,&inertia,&omsg,&gamma);
00445     // Unset S_chol so that there is no chance of accedental modification.
00446     if(!S_indef_)
00447       S_chol_.init_setup(NULL);
00448     // Update the state variables
00449     ++S_size_;
00450     S_indef_      = true;
00451     fact_updated_ = true;
00452     fact_in1_     = fact_in1;
00453     inertia_      = inertia;
00454   }
00455   else {
00456     //
00457     // Don't update the factorization yet
00458     //
00459     if(!S_indef_) // We must set the inertia if it was definite!
00460       inertia_ = S_chol_.inertia();
00461     ++S_size_;
00462     S_indef_      = true;
00463     fact_updated_ = false;  // fact_in1_ is meaningless now
00464     // We need to keep track of the expected inertia!
00465     if( inertia_.zero_eigens == Inertia::UNKNOWN || add_eigen_val == MSADU::EIGEN_VAL_UNKNOWN )
00466       inertia_ = Inertia(); // Unknow inertia!
00467     else if( add_eigen_val == MSADU::EIGEN_VAL_NEG )
00468       ++inertia_.neg_eigens;
00469     else if( add_eigen_val == MSADU::EIGEN_VAL_POS )
00470       ++inertia_.pos_eigens;
00471     else
00472       TEUCHOS_TEST_FOR_EXCEPT(true); // Should not happen!
00473   }
00474   if( throw_exception )
00475     throw WarnNearSingularUpdateException(omsg.str(),gamma);
00476 }
00477 
00478 void MatrixSymAddDelBunchKaufman::delete_update(
00479   size_type          jd
00480   ,bool              force_refactorization
00481   ,EEigenValType     drop_eigen_val
00482   ,PivotTolerances   pivot_tols
00483   )
00484 {
00485   using BLAS_Cpp::upper;
00486   using BLAS_Cpp::lower;
00487   using DenseLinAlgPack::tri_ele;
00488   using DenseLinAlgPack::nonconst_tri_ele;
00489   typedef MatrixSymAddDelUpdateable  MSADU;
00490 
00491   assert_initialized();
00492 
00493   if( jd < 1 || S_size_ < jd )
00494     throw std::out_of_range(
00495       "MatrixSymAddDelBunchKaufman::delete_update(jd,...): "
00496       "Error, the indice jd must be 1 <= jd <= rows()" );
00497 
00498   bool                throw_exception = false; // If true then throw exception
00499   std::ostringstream  omsg;                    // Will be set if an exception has to be thrown.
00500   value_type          gamma;                   // ...
00501 
00502   if( !S_indef_ ) {
00503     //
00504     // The current matrix is p.d. or n.d. and so should the new matrix.
00505     //
00506     S_chol_.delete_update(jd,force_refactorization,drop_eigen_val,pivot_tols);
00507     --S_size_;
00508   }
00509   else {
00510     //
00511     // The current matrix is indefinite but the new matrix
00512     // could be p.d. or n.d. in which case we could switch
00513     // to the cholesky factorization.  If the user says the
00514     // new matrix will be p.d. or n.d. then we will take
00515     // his/her word for it and try the cholesky factorization
00516     // otherwise just recompute the indefinite factorization.
00517     //
00518     // ToDo: (2/2/01): We could also try some more fancy updating
00519     // procedures and try to get away with some potentially unstable
00520     // methods and resort to the full indefinite factorization if needed.
00521     // The sign of the dropped eigen value might tell us what we
00522     // might get away with?
00523     //
00524     // Update the factorization
00525     //
00526     Inertia inertia = S_chol_.inertia();
00527     if( (drop_eigen_val == MSADU::EIGEN_VAL_POS && inertia.pos_eigens == 1 )
00528       || (drop_eigen_val == MSADU::EIGEN_VAL_NEG && inertia.neg_eigens == 1 )
00529       || S_size_ == 2 )
00530     {
00531       // Lets take the users word for it and switch to a cholesky factorization.
00532       // To do this we will just let S_chol_ do it for us.  If the new matrix
00533       // turns out not to be what the user says, then we will resize the matrix
00534       // to zero and thrown an exception.
00535       try {
00536         // Delete row and column jd from M
00537         DMatrixSlice S = this->S(S_size_).gms();
00538         DenseLinAlgPack::delete_row_col( jd, &DenseLinAlgPack::nonconst_tri_ele(S,BLAS_Cpp::lower) ); 
00539         // Setup S_chol and factor the thing
00540         S_chol_.init_setup(&S_store1_(),Teuchos::null,0,true,true,true,false,0.0);
00541         S_chol_.initialize( this->S(S_size_-1), S_store1_.rows()-1
00542                   , force_refactorization, Inertia(), PivotTolerances() );
00543       }
00544       catch( const SingularUpdateException& excpt ) {
00545         S_size_ = 0;
00546         throw SingularUpdateException(
00547           std::string(
00548             "MatrixSymAddDelBunchKaufman::delete_update(...) : "
00549             "Error, the client incorrectly specified that the new "
00550             "matrix would be nonsingular and not indefinite:\n" )
00551           + excpt.what()
00552           , excpt.gamma
00553           );
00554       }
00555       catch( const WrongInertiaUpdateException& excpt ) {
00556         S_size_ = 0;
00557         throw WrongInertiaUpdateException(
00558           std::string(
00559             "MatrixSymAddDelBunchKaufman::delete_update(...) : "
00560             "Error, the client incorrectly specified that the new "
00561             "matrix would not be indefinite:\n" )
00562           + excpt.what()
00563           , excpt.gamma
00564           );
00565       }
00566       // If we get here the update succeeded and the new matrix is p.d. or n.d.
00567       --S_size_;
00568       S_indef_ = false;
00569     }
00570     else {
00571       // 
00572       // We have been given no indication that the new matrix is p.d. or n.d.
00573       // so we will assume it is indefinite and go from there.
00574       //
00575       DMatrixSlice  S  = this->S(S_size_).gms();
00576       if( force_refactorization ) {
00577         // 
00578         // Perform the refactorization carefully
00579         //
00580         const bool fact_in1 = !fact_in1_;
00581         // Copy the original into the unused storage location
00582         DMatrixSliceTriEle  DU = this->DU(S_size_,fact_in1);
00583         DenseLinAlgPack::assign(  &DU, tri_ele(S,lower) );
00584         // Delete row and column jd from the storage location for DU
00585         DenseLinAlgPack::delete_row_col( jd, &DU );
00586         // Now factor the matrix inplace
00587         try {
00588           factor_matrix(S_size_-1,fact_in1);
00589         }
00590         catch( const DenseLinAlgLAPack::FactorizationException& excpt ) {
00591           omsg
00592             << "MatrixSymAddDelBunchKaufman::delete_update(...): "
00593             << "Error, singular update but the original matrix was maintianed:\n"
00594             << excpt.what();
00595           throw SingularUpdateException( omsg.str(), -1.0 );
00596         }
00597         // Compute the expected inertia
00598         Inertia expected_inertia = this->inertia();
00599         if( expected_inertia.zero_eigens == Inertia::UNKNOWN || drop_eigen_val == MSADU::EIGEN_VAL_UNKNOWN )
00600           expected_inertia = Inertia(); // Unknow inertia!
00601         else if( drop_eigen_val == MSADU::EIGEN_VAL_NEG )
00602           --expected_inertia.neg_eigens;
00603         else if( drop_eigen_val == MSADU::EIGEN_VAL_POS )
00604           --expected_inertia.pos_eigens;
00605         else
00606           TEUCHOS_TEST_FOR_EXCEPT(true); // Should not happen!
00607         // Compute the exacted inertia and validate that it is what is expected
00608         Inertia inertia;
00609         throw_exception = compute_assert_inertia(
00610           S_size_-1,fact_in1,expected_inertia,"delete_update",pivot_tols
00611           ,&inertia,&omsg,&gamma);
00612         // If we get here the factorization worked out.
00613         --S_size_;
00614         S_indef_      = true;
00615         fact_updated_ = true;
00616         fact_in1_     = fact_in1;
00617         inertia_      = inertia;
00618       }
00619       // Delete the row and column jd from the original
00620       DenseLinAlgPack::delete_row_col( jd, &nonconst_tri_ele(S,lower) );
00621       if( !force_refactorization ) {
00622         //
00623         // The refactorization was not forced
00624         //
00625         --S_size_;
00626         S_indef_      = true;
00627         fact_updated_ = false;  // fact_in1_ is meaningless now
00628         // We need to keep track of the expected inertia!
00629         if( inertia_.zero_eigens == Inertia::UNKNOWN || drop_eigen_val == MSADU::EIGEN_VAL_UNKNOWN )
00630           inertia_ = Inertia(); // Unknow inertia!
00631         else if( drop_eigen_val == MSADU::EIGEN_VAL_NEG )
00632           --inertia_.neg_eigens;
00633         else if( drop_eigen_val == MSADU::EIGEN_VAL_POS )
00634           --inertia_.pos_eigens;
00635         else
00636           TEUCHOS_TEST_FOR_EXCEPT(true); // Should not happen!
00637       }
00638     }
00639   }
00640   if( throw_exception )
00641     throw WarnNearSingularUpdateException(omsg.str(),gamma);
00642 }
00643 
00644 // Overridden from MatrixSymOpNonsingSerial
00645 
00646 size_type MatrixSymAddDelBunchKaufman::rows() const
00647 {
00648   return S_size_;
00649 }
00650 
00651 std::ostream& MatrixSymAddDelBunchKaufman::output(std::ostream& out) const
00652 {
00653   if( S_size_ ) {
00654     out << "Unfactored symmetric matrix stored as lower triangle (ignore upper nonzeros):\n"
00655       << S(S_size_).gms();
00656     if( S_indef_ ) {
00657       out << "Upper symmetric indefinite factor DU returned from sytrf(...) (ignore lower nonzeros):\n"
00658         << DU(S_size_,fact_in1_).gms();
00659       out << "Permutation array IPIV returned from sytrf(...):\n";
00660       for( size_type i = 0; i < S_size_; ++i )
00661         out << "  " << IPIV_[i];
00662       out << std::endl;
00663     }
00664     else {
00665       out << "Upper cholesky factor (ignore lower nonzeros):\n" << DU(S_size_,true).gms();
00666     }
00667   }
00668   else {
00669     out << "0 0\n";
00670   }
00671   return out;
00672 }
00673 
00674 void MatrixSymAddDelBunchKaufman::Vp_StMtV(
00675   DVectorSlice* y, value_type a, BLAS_Cpp::Transp M_trans
00676   ,const DVectorSlice& x, value_type b
00677   ) const
00678 {
00679   DenseLinAlgPack::Vp_MtV_assert_sizes( y->dim(), S_size_, S_size_, M_trans, x.dim() );
00680   // Use the unfactored matrix since this is more accurate!
00681   DenseLinAlgPack::Vp_StMtV( y, a, S(S_size_), BLAS_Cpp::no_trans, x, b );
00682 }
00683 
00684 void MatrixSymAddDelBunchKaufman::V_InvMtV(
00685   DVectorSlice* y, BLAS_Cpp::Transp M_trans
00686   ,const DVectorSlice& x
00687   ) const
00688 {
00689   const size_type n = S_size_;
00690   DenseLinAlgPack::Vp_MtV_assert_sizes( y->dim(), n, n, M_trans, x.dim() );
00691   if( S_indef_ ) {
00692     // Update the factorzation if needed
00693     if(!fact_updated_) {
00694       const bool fact_in1 = true;
00695       MatrixSymAddDelBunchKaufman
00696         *nc_this = const_cast<MatrixSymAddDelBunchKaufman*>(this);
00697       nc_this->copy_and_factor_matrix(S_size_,fact_in1); // May throw exceptions!
00698       nc_this->fact_updated_ = true;
00699       nc_this->fact_in1_     = fact_in1;
00700     }
00701     *y = x;
00702     DenseLinAlgLAPack::sytrs(
00703       DU(S_size_,fact_in1_), &const_cast<IPIV_t&>(IPIV_)[0]
00704       , &DMatrixSlice(y->raw_ptr(),n,n,n,1), &WORK_() );
00705   }
00706   else {
00707     AbstractLinAlgPack::V_InvMtV( y, S_chol_, M_trans, x );
00708   }
00709 }
00710 
00711 // Private
00712 
00713 void MatrixSymAddDelBunchKaufman::assert_initialized() const
00714 {
00715   if(!S_size_)
00716     throw std::logic_error(
00717       "MatrixSymAddDelBunchKaufman::assert_initialized() : "
00718       "Error, this matrix is not initialized yet" );
00719 }
00720 
00721 void MatrixSymAddDelBunchKaufman::resize_DU_store( bool in_store1 )
00722 {
00723   if(!in_store1 && S_store2_.rows() < S_store1_.rows() )
00724     S_store2_.resize( S_store1_.rows(), S_store1_.cols() );
00725 }
00726 
00727 void MatrixSymAddDelBunchKaufman::copy_and_factor_matrix(
00728   size_type S_size, bool fact_in1 )
00729 {
00730   DenseLinAlgPack::assign(
00731     &DU(S_size,fact_in1)
00732     , DenseLinAlgPack::tri_ele(S(S_size).gms(),BLAS_Cpp::lower) );
00733   factor_matrix(S_size,fact_in1);
00734 }
00735 
00736 void MatrixSymAddDelBunchKaufman::factor_matrix( size_type S_size, bool fact_in1 )
00737 {
00738   // Resize the workspace first
00739   const size_type opt_block_size = 64; // This is an estimate (see sytrf(...))
00740   if( WORK_.dim() < S_store1_.rows() * opt_block_size )
00741     WORK_.resize(S_store1_.rows()*opt_block_size);
00742   if( IPIV_.size() < S_store1_.rows() )
00743     IPIV_.resize(S_store1_.rows());
00744   // Factor the matrix (will throw FactorizationException if singular)
00745   DenseLinAlgLAPack::sytrf( &DU(S_size,fact_in1), &IPIV_[0], &WORK_() );
00746 }
00747 
00748 bool MatrixSymAddDelBunchKaufman::compute_assert_inertia(
00749   size_type S_size, bool fact_in1, const Inertia& exp_inertia, const char func_name[]
00750   ,PivotTolerances pivot_tols, Inertia* comp_inertia, std::ostringstream* omsg, value_type* gamma )
00751 {
00752   bool throw_exception = false;
00753   *gamma = 0.0;
00754   // Here we will compute the inertia given IPIV[] and D[] as described in the documentation
00755   // for dsytrf(...) (see the source code).
00756   const DMatrixSlice DU = this->DU(S_size,fact_in1).gms();
00757   const size_type      n = DU.rows();
00758   Inertia inertia(0,0,0);
00759   value_type max_diag = 0.0;
00760   value_type min_diag = std::numeric_limits<value_type>::max();
00761   for( size_type k = 1; k <= n; ) {
00762     const FortranTypes::f_int k_p = IPIV_[k-1];
00763     if( k_p > 0 ) {
00764       // D(k,k) is a 1x1 block.
00765       // Lets get the eigen value from the sign of D(k,k)
00766       const value_type D_k_k = DU(k,k), abs_D_k_k = std::fabs(D_k_k);
00767       if( D_k_k > 0.0 )
00768         ++inertia.pos_eigens;
00769       else
00770         ++inertia.neg_eigens;
00771       if(abs_D_k_k > max_diag) max_diag = abs_D_k_k;
00772       if(abs_D_k_k < min_diag) min_diag = abs_D_k_k;
00773       k++;
00774     }
00775     else {
00776       // D(k:k+1,k:k+1) is a 2x2 block.
00777       // This represents one positive eigen value and
00778       // on negative eigen value
00779       TEUCHOS_TEST_FOR_EXCEPT( !(  IPIV_[k] == k_p ) ); // This is what the documentation for xSYTRF(...) says!
00780       ++inertia.pos_eigens;
00781       ++inertia.neg_eigens;
00782       // To find the largest and smallest diagonals of U for L*U we must perform Gaussian
00783       // elimination on this 2x2 block:
00784       const value_type                                   // [ a   b ] = D(k:k+1,k:k+1) 
00785         a = DU(k,k), b = DU(k,k+1), c = DU(k+1,k+1),   // [ b   c ]
00786         abs_a = std::fabs(a), abs_b = std::fabs(b);  
00787       value_type pivot_1, pivot_2;
00788       if( abs_a > abs_b ) { // Pivot on a = D(k,k)
00789         pivot_1 = abs_a;              // [   1      ] * [ a   b ] = [ a      b     ]
00790         pivot_2 = std::fabs(c - b*b/a);  // [ -b/a  1  ]   [ b   c ]   [ 0  c - b*b/a ]
00791       }
00792       else {                // Pivot on b = D(k+1,k) = D(k,k+1)
00793         pivot_1 = abs_b;              // [   1      ] * [ b   c ] = [ b      c     ]
00794         pivot_2 = std::fabs(b - a*c/b);  // [ -a/b  1  ]   [ a   b ]   [ 0  b - a*c/b ]
00795       }
00796       if(pivot_1 > max_diag) max_diag = pivot_1;
00797       if(pivot_1 < min_diag) min_diag = pivot_1;
00798       if(pivot_2 > max_diag) max_diag = pivot_2;
00799       if(pivot_2 < min_diag) min_diag = pivot_2;
00800       k+=2;
00801     }
00802   }
00803   // gamma = min(|diag(i)|)/max(|diag(i)|)
00804   *gamma = min_diag / max_diag;
00805   // Now validate that the inertia is what is expected
00806     const bool
00807     wrong_inertia =
00808     ( exp_inertia.neg_eigens != Inertia::UNKNOWN
00809       && exp_inertia.neg_eigens != inertia.neg_eigens )
00810     || ( exp_inertia.pos_eigens != Inertia::UNKNOWN
00811        && exp_inertia.pos_eigens != inertia.pos_eigens ) ;
00812   PivotTolerances use_pivot_tols = S_chol_.pivot_tols();
00813   if( pivot_tols.warning_tol != PivotTolerances::UNKNOWN )
00814     use_pivot_tols.warning_tol = pivot_tols.warning_tol;
00815   if( pivot_tols.singular_tol != PivotTolerances::UNKNOWN )
00816     use_pivot_tols.singular_tol = pivot_tols.singular_tol;
00817   if( pivot_tols.wrong_inertia_tol != PivotTolerances::UNKNOWN )
00818     use_pivot_tols.wrong_inertia_tol = pivot_tols.wrong_inertia_tol;
00819   throw_exception = (
00820     *gamma == 0.0
00821     || !wrong_inertia && *gamma <= use_pivot_tols.warning_tol
00822     || !wrong_inertia && *gamma <= use_pivot_tols.singular_tol
00823     || wrong_inertia
00824     );
00825   // Create message and throw exceptions
00826   std::ostringstream onum_msg;
00827   if(throw_exception) {
00828     if(wrong_inertia) {
00829       onum_msg
00830         << "inertia = ("
00831         << inertia.neg_eigens << "," << inertia.zero_eigens << "," << inertia.pos_eigens
00832         << ") != expected_inertia = ("
00833         << exp_inertia.neg_eigens << "," << exp_inertia.zero_eigens << "," << exp_inertia.pos_eigens << ")"
00834         << std::endl;
00835     }
00836     onum_msg
00837       << "gamma = min(|diag(D)(k)|)/max(|diag(D)(k)|) = " <<  min_diag << "/" << max_diag
00838       << " = " << *gamma;
00839     *omsg
00840       << "MatrixSymAddDelBunchKaufman::"<<func_name<<"(...): ";
00841     if( *gamma == 0.0 || (!wrong_inertia && *gamma <= use_pivot_tols.singular_tol) ) {
00842       *omsg
00843         << "Singular update!\n" << onum_msg.str() << " <= singular_tol = "
00844         << use_pivot_tols.singular_tol;
00845       throw SingularUpdateException( omsg->str(), *gamma );
00846     }
00847     else if( wrong_inertia && *gamma <= use_pivot_tols.singular_tol ) {
00848       *omsg
00849         << "Singular update!\n" << onum_msg.str() << " <= wrong_inertia_tol = "
00850         << use_pivot_tols.wrong_inertia_tol;
00851       throw SingularUpdateException( omsg->str(), *gamma );
00852     }
00853     else if( wrong_inertia ) {
00854       *omsg
00855         << "Indefinite update!\n" << onum_msg.str() << " >= wrong_inertia_tol = "
00856         << use_pivot_tols.wrong_inertia_tol;
00857       throw WrongInertiaUpdateException( omsg->str(), *gamma );
00858     }
00859     else if( !wrong_inertia && *gamma <= use_pivot_tols.warning_tol ) {
00860       *omsg
00861         << "Warning, near singular update!\nsingular_tol = " << use_pivot_tols.singular_tol
00862         << " < " << onum_msg.str() << " <= warning_tol = "
00863         << use_pivot_tols.warning_tol;
00864       // Don't throw the exception till the very end!
00865     }
00866     else {
00867       TEUCHOS_TEST_FOR_EXCEPT(true); // Only local programming error?
00868     }
00869   }
00870   // Return
00871   *comp_inertia = inertia;
00872   return throw_exception;
00873 }
00874 
00875 } // end namespace ConstrainedOptPack
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends