ConstrainedOptPack_MatrixSymAddDelBunchKaufman.cpp

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

Generated on Tue Jul 13 09:30:51 2010 for MOOCHO (Single Doxygen Collection) by  doxygen 1.4.7