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