ConstrainedOptPack: C++ Tools for Constrained (and Unconstrained) Optimization  Version of the Day
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Groups Pages
ConstrainedOptPack_MatrixSymAddDelBunchKaufman.cpp
1 // @HEADER
2 // ***********************************************************************
3 //
4 // Moocho: Multi-functional Object-Oriented arCHitecture for Optimization
5 // Copyright (2003) Sandia Corporation
6 //
7 // Under terms of Contract DE-AC04-94AL85000, there is a non-exclusive
8 // license for use of this work by or on behalf of the U.S. Government.
9 //
10 // Redistribution and use in source and binary forms, with or without
11 // modification, are permitted provided that the following conditions are
12 // met:
13 //
14 // 1. Redistributions of source code must retain the above copyright
15 // notice, this list of conditions and the following disclaimer.
16 //
17 // 2. Redistributions in binary form must reproduce the above copyright
18 // notice, this list of conditions and the following disclaimer in the
19 // documentation and/or other materials provided with the distribution.
20 //
21 // 3. Neither the name of the Corporation nor the names of the
22 // contributors may be used to endorse or promote products derived from
23 // this software without specific prior written permission.
24 //
25 // THIS SOFTWARE IS PROVIDED BY SANDIA CORPORATION "AS IS" AND ANY
26 // EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
27 // IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR
28 // PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL SANDIA CORPORATION OR THE
29 // CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL,
30 // EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO,
31 // PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR
32 // PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF
33 // LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING
34 // NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
35 // SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
36 //
37 // Questions? Contact Roscoe A. Bartlett (rabartl@sandia.gov)
38 //
39 // ***********************************************************************
40 // @HEADER
41 //
42 
43 #include <assert.h>
44 
45 #include <limits>
46 
47 #include "ConstrainedOptPack_MatrixSymAddDelBunchKaufman.hpp"
48 #include "DenseLinAlgLAPack.hpp"
49 #include "DenseLinAlgPack_DMatrixOut.hpp"
50 #include "DenseLinAlgPack_DMatrixOp.hpp"
51 #include "DenseLinAlgPack_AssertOp.hpp"
52 #include "DenseLinAlgPack_delete_row_col.hpp"
53 
54 namespace ConstrainedOptPack {
55 
57  : S_size_(0), S_indef_(false), fact_updated_(false), fact_in1_(true), inertia_(0,0,0)
58 {}
59 
60 void MatrixSymAddDelBunchKaufman::pivot_tols( PivotTolerances pivot_tols )
61 {
62  S_chol_.pivot_tols(pivot_tols);
63 }
64 
65 MatrixSymAddDelUpdateable::PivotTolerances MatrixSymAddDelBunchKaufman::pivot_tols() const
66 {
67  return S_chol_.pivot_tols();
68 }
69 
70 // Overridden from MatrixSymAddDelUpdateableWithOpNonsingular
71 
72 const MatrixSymOpNonsing& MatrixSymAddDelBunchKaufman::op_interface() const
73 {
74  return *this;
75 }
76 
78 {
79  return *this;
80 }
81 
82 const MatrixSymAddDelUpdateable& MatrixSymAddDelBunchKaufman::update_interface() const
83 {
84  return *this;
85 }
86 
87 // Overridden from MatrixSymAddDelUpdateable
88 
90  value_type alpha
91  ,size_type max_size
92  )
93 {
94  try {
95  // Resize the storage if we have to
96  if( S_store1_.rows() < max_size+1 && S_store1_.cols() < max_size+1 )
97  S_store1_.resize(max_size+1,max_size+1);
98  fact_in1_ = true;
99  // Start out with a p.d. or n.d. matrix and maintain the original
100  S_chol_.init_setup(&S_store1_(),Teuchos::null,0,true,true,true,false,0.0);
101  S_chol_.initialize(alpha,max_size);
102  // Set the state variables:
103  S_size_ = 1;
104  S_indef_ = false; // fact_updated_, fact_in1_ and inertia are meaningless!
105  }
106  catch(...) {
107  S_size_ = 0;
108  throw;
109  }
110 }
111 
113  const DMatrixSliceSym &A
114  ,size_type max_size
115  ,bool force_factorization
116  ,Inertia expected_inertia
117  ,PivotTolerances pivot_tols
118  )
119 {
120  using BLAS_Cpp::upper;
121  using BLAS_Cpp::lower;
122  using DenseLinAlgPack::tri_ele;
123  using DenseLinAlgPack::nonconst_tri_ele;
124  typedef MatrixSymAddDelUpdateable MSADU;
125  typedef MSADU::Inertia Inertia;
126 
127  bool throw_exception = false; // If true then throw exception
128  std::ostringstream omsg; // Will be set if an exception has to be thrown.
129  value_type gamma; // ...
130 
131  const size_type
132  n = A.rows();
133 
134  // Validate proper usage of inertia parameter
135  TEUCHOS_TEST_FOR_EXCEPT( ! ( expected_inertia.zero_eigens == Inertia::UNKNOWN
136  || expected_inertia.zero_eigens == 0 ) );
137 
138  try {
139  // Resize the storage if we have to
140  if( S_store1_.rows() < max_size+1 && S_store1_.cols() < max_size+1 )
141  S_store1_.resize(max_size+1,max_size+1);
142  fact_in1_ = true;
143  // See if the client claims that the matrix is p.d. or n.d.
144  const bool not_indefinite
145  = ( ( expected_inertia.neg_eigens == 0 && expected_inertia.pos_eigens == n )
146  || ( expected_inertia.neg_eigens == n && expected_inertia.pos_eigens == 0 ) );
147  // Initialize the matrix
148  if( not_indefinite ) {
149  // The client says that the matrix is p.d. or n.d. so
150  // we will take their word for it.
151  S_chol_.init_setup(&S_store1_(),Teuchos::null,0,true,true,true,false,0.0);
152  try {
153  S_chol_.initialize(A,max_size,force_factorization,expected_inertia,pivot_tols);
154  }
155  catch(const MSADU::WarnNearSingularUpdateException& except) {
156  throw_exception = true; // Throw this same exception at the end!
157  omsg << except.what();
158  gamma = except.gamma;
159  }
160  // Set the state variables:
161  S_size_ = n;
162  S_indef_ = false; // fact_updated_, fact_in1_ and inertia are meaningless!
163  }
164  else {
165  //
166  // The client did not say that the matrix was p.d. or n.d. so we
167  // must assume that it is indefinite.
168  //
169  bool fact_in1 = !fact_in1_;
170  // Set the new matrix in the unused factorization location
171  DenseLinAlgPack::assign( &DU(n,fact_in1), tri_ele( A.gms(), A.uplo() ) );
172  // Update the factorization in place
173  try {
174  factor_matrix( n, fact_in1 );
175  }
176  catch( const DenseLinAlgLAPack::FactorizationException& excpt ) {
177  omsg
178  << "MatrixSymAddDelBunchKaufman::initialize(...): "
179  << "Error, the matrix A is singular:\n"
180  << excpt.what();
181  throw SingularUpdateException( omsg.str(), -1.0 );
182  }
183  // Compute the inertia and validate that it is correct.
184  Inertia inertia;
185  throw_exception = compute_assert_inertia(
186  n,fact_in1,expected_inertia,"initialize",pivot_tols
187  ,&inertia,&omsg,&gamma);
188  // If the client did not know the inertia of the
189  // matrix but it turns out to be p.d. or n.d. then modify the
190  // DU factor appropriatly and switch to cholesky factorization.
191  if( inertia.pos_eigens == n || inertia.neg_eigens == n ) {
192  TEUCHOS_TEST_FOR_EXCEPT(true); // ToDo Implememnt this!
193  }
194  else {
195  // Indead the new matrix is indefinite
196  // Set the original matrix now
197  DenseLinAlgPack::assign( &DenseLinAlgPack::nonconst_tri_ele( S(n).gms(), BLAS_Cpp::lower)
198  , tri_ele( A.gms(), A.uplo() ) );
199  // Update the state variables:
200  S_size_ = n;
201  S_indef_ = true;
202  fact_updated_ = true;
203  fact_in1_ = fact_in1;
204  inertia_ = inertia;
205  }
206  }
207  }
208  catch(...) {
209  S_size_ = 0;
210  throw;
211  }
212  if( throw_exception )
213  throw WarnNearSingularUpdateException(omsg.str(),gamma);
214 }
215 
217 {
218  return S_store1_.rows() -1; // The center diagonal is used for workspace
219 }
220 
221 MatrixSymAddDelUpdateable::Inertia
223 {
224  return S_indef_ ? inertia_ : S_chol_.inertia();
225 }
226 
228 {
229  S_size_ = 0;
230 }
231 
233  const DVectorSlice *t
234  ,value_type alpha
235  ,bool force_refactorization
236  ,EEigenValType add_eigen_val
237  ,PivotTolerances pivot_tols
238  )
239 {
241  using BLAS_Cpp::no_trans;
242  using DenseLinAlgPack::norm_inf;
244  typedef MatrixSymAddDelUpdateable MSADU;
245 
246  assert_initialized();
247 
248  // Validate the input
249  if (implicit_ref_cast<const MatrixBase>(*this).rows() >= max_size()) {
250  throw MaxSizeExceededException(
251  "MatrixSymAddDelBunchKaufman::augment_update(...) : "
252  "Error, the maximum size would be exceeded." );
253  }
254  if (t && t->dim() != S_size_) {
255  throw std::length_error(
256  "MatrixSymAddDelBunchKaufman::augment_update(...): "
257  "Error, t.dim() must be equal to this->rows()." );
258  }
259  if (!(add_eigen_val == MSADU::EIGEN_VAL_UNKNOWN
260  || add_eigen_val != MSADU::EIGEN_VAL_ZERO))
261  {
262  throw SingularUpdateException(
263  "MatrixSymAddDelBunchKaufman::augment_update(...): "
264  "Error, the client has specified a singular update in add_eigen_val.", -1.0 );
265  }
266 
267  // Try to perform the update
268  bool throw_exception = false; // If true then throw exception
269  std::ostringstream omsg; // Will be set if an exception has to be thrown.
270  value_type gamma; // ...
271  if( !S_indef_ ) {
272  // The current matrix is positive definite or negative definite. If the
273  // new matrix is still p.d. or n.d. when we are good to go. We must first
274  // check if the client has specified that the new matrix will be indefinite
275  // and if so then we will take his/her word for it and do the indefinite
276  // updating.
277  MSADU::Inertia inertia = S_chol_.inertia();
278  const bool
279  new_S_not_indef
280  = ( ( inertia.pos_eigens > 0
281  && ( add_eigen_val == MSADU::EIGEN_VAL_POS || add_eigen_val == MSADU::EIGEN_VAL_UNKNOWN ) )
282  || ( inertia.neg_eigens > 0
283  && ( add_eigen_val == MSADU::EIGEN_VAL_NEG || add_eigen_val == MSADU::EIGEN_VAL_UNKNOWN ) )
284  );
285  if( !new_S_not_indef ) {
286  // We must resort to an indefinite factorization
287  }
288  else {
289  // The new matrix could/should still be p.d. or n.d. so let's try it!
290  bool update_successful = false;
291  try {
292  S_chol_.augment_update(t,alpha,force_refactorization,add_eigen_val,pivot_tols);
293  update_successful = true;
294  }
295  catch(const MSADU::WrongInertiaUpdateException&) {
296  if( add_eigen_val != MSADU::EIGEN_VAL_UNKNOWN )
297  throw; // The client specified the new inertia and it was wrong so throw execepiton.
298  // If the client did not know that inertia then we can't fault them so we will
299  // proceed unto the indefinite factorization.
300  }
301  catch(const MSADU::WarnNearSingularUpdateException& except) {
302  throw_exception = true; // Throw this same exception at the end!
303  update_successful = true;
304  omsg << except.what();
305  gamma = except.gamma;
306  }
307  if( update_successful ) {
308  ++S_size_; // Now we can resize the matrix
309  if(throw_exception)
310  throw MSADU::WarnNearSingularUpdateException(omsg.str(),gamma);
311  else
312  return;
313  }
314  }
315  }
316  //
317  // If we get here then we have no choice but the perform an indefinite factorization.
318  //
319  // ToDo: (2/2/01): We could also try some more fancy updating
320  // procedures and try to get away with some potentially unstable
321  // methods in the interest of time savings.
322  //
323  const size_type n = S_size_;
324  DMatrixSlice S_bar = this->S(n+1).gms();
325  //
326  // Validate that the new matrix will be nonsingular.
327  //
328  // Given any nonsingular matrix S (even unsymmetric) it is easy to show that
329  // beta = alpha - t'*inv(S)*t != 0.0 if [ S, t; t', alpha ] is nonsingular.
330  // This is a cheap O(n^2) way to check that the update is nonsingular before
331  // we go through the O(n^3) refactorization.
332  // In fact, the sign of beta even tells us the sign of the new eigen value
333  // of the updated matrix even before we perform the refactorization.
334  // If the current matrix is not factored then we will just skip this
335  // test and let the full factorization take place to find this out.
336  // We could even cheat a little and actually perform the update with this
337  // diagonal beta and then compute the update to the factors U our selves
338  // in O(n^2) time.
339  //
340  if( force_refactorization && fact_updated_ ) {
341  const value_type
342  beta = alpha - ( t ? transVtInvMtV(*t,*this,no_trans,*t) : 0.0 ),
343  abs_beta = std::fabs(beta),
344  nrm_D_diag = norm_inf(DU(n,fact_in1_).gms().diag()); // ToDo: Consider 2x2 blocks also!
345  gamma = abs_beta / nrm_D_diag;
346  // Check gamma
347  const bool
348  correct_inertia = (
349  add_eigen_val == MSADU::EIGEN_VAL_UNKNOWN
350  || beta > 0.0 && add_eigen_val == MSADU::EIGEN_VAL_POS
351  || beta < 0.0 && add_eigen_val == MSADU::EIGEN_VAL_NEG
352  );
353  PivotTolerances use_pivot_tols = S_chol_.pivot_tols();
354  if( pivot_tols.warning_tol != PivotTolerances::UNKNOWN )
355  use_pivot_tols.warning_tol = pivot_tols.warning_tol;
356  if( pivot_tols.singular_tol != PivotTolerances::UNKNOWN )
357  use_pivot_tols.singular_tol = pivot_tols.singular_tol;
358  if( pivot_tols.wrong_inertia_tol != PivotTolerances::UNKNOWN )
359  use_pivot_tols.wrong_inertia_tol = pivot_tols.wrong_inertia_tol;
360  throw_exception = (
361  gamma == 0.0
362  || correct_inertia && gamma <= use_pivot_tols.singular_tol
363  || !correct_inertia
364  );
365  // Create message and throw exceptions
366  std::ostringstream onum_msg;
367  if(throw_exception) {
368  onum_msg
369  << "gamma = |alpha-t'*inv(S)*t)|/||diag(D)||inf = |" << beta << "|/" << nrm_D_diag
370  << " = " << gamma;
371  omsg
372  << "MatrixSymAddDelBunchKaufman::augment_update(...): ";
373  if( correct_inertia && gamma <= use_pivot_tols.singular_tol ) {
374  omsg
375  << "Singular update!\n" << onum_msg.str() << " <= singular_tol = "
376  << use_pivot_tols.singular_tol;
377  throw SingularUpdateException( omsg.str(), gamma );
378  }
379  else if( !correct_inertia && gamma <= use_pivot_tols.singular_tol ) {
380  omsg
381  << "Singular update!\n" << onum_msg.str() << " <= wrong_inertia_tol = "
382  << use_pivot_tols.wrong_inertia_tol;
383  throw SingularUpdateException( omsg.str(), gamma );
384  }
385 
386  else if( !correct_inertia ) {
387  omsg
388  << "Indefinite update!\n" << onum_msg.str() << " >= wrong_inertia_tol = "
389  << use_pivot_tols.wrong_inertia_tol;
390  throw WrongInertiaUpdateException( omsg.str(), gamma );
391  }
392  else {
393  TEUCHOS_TEST_FOR_EXCEPT(true); // Only local programming error?
394  }
395  }
396  }
397 
398  // Add new row to the lower part of the original symmetric matrix.
399  if(t)
400  S_bar.row(n+1)(1,n) = *t;
401  else
402  S_bar.row(n+1)(1,n) = 0.0;
403  S_bar(n+1,n+1) = alpha;
404 
405  //
406  // Update the factorization
407  //
408  if( force_refactorization ) {
409  // Determine where to copy the original matrix to
410  bool fact_in1 = false;
411  if( S_indef_ ) {
412  // S is already indefinite so let's copy the original into the storage
413  // location other than the current one in case the newly nonsingular matrix
414  // is singular or has the wrong inertia.
415  fact_in1 = !fact_in1_;
416  }
417  else {
418  // S is currently p.d. or n.d. so let's copy the new matrix
419  // into the second storage location so as not to overwrite
420  // the current cholesky factor in case the new matrix is singular
421  // or has the wrong inertia.
422  fact_in1 = false;
423  }
424  // Copy and factor the new matrix
425  try {
426  copy_and_factor_matrix(n+1,fact_in1);
427  }
428  catch( const DenseLinAlgLAPack::FactorizationException& excpt ) {
429  std::ostringstream omsg;
430  omsg
431  << "MatrixSymAddDelBunchKaufman::augment_update(...): "
432  << "Error, singular update but the original matrix was be maintianed:\n"
433  << excpt.what();
434  throw SingularUpdateException( omsg.str(), -1.0 );
435  }
436  // Compute the expected inertia
437  Inertia expected_inertia = this->inertia();
438  if( expected_inertia.zero_eigens == Inertia::UNKNOWN || add_eigen_val == MSADU::EIGEN_VAL_UNKNOWN )
439  expected_inertia = Inertia(); // Unknow inertia!
440  else if( add_eigen_val == MSADU::EIGEN_VAL_NEG )
441  ++expected_inertia.neg_eigens;
442  else if( add_eigen_val == MSADU::EIGEN_VAL_POS )
443  ++expected_inertia.pos_eigens;
444  else
445  TEUCHOS_TEST_FOR_EXCEPT(true); // Should not happen!
446  // Compute the actually inertia and validate that it is what is expected
447  Inertia inertia;
448  throw_exception = compute_assert_inertia(
449  n+1,fact_in1,expected_inertia,"augment_update",pivot_tols
450  ,&inertia,&omsg,&gamma);
451  // Unset S_chol so that there is no chance of accedental modification.
452  if(!S_indef_)
453  S_chol_.init_setup(NULL);
454  // Update the state variables
455  ++S_size_;
456  S_indef_ = true;
457  fact_updated_ = true;
458  fact_in1_ = fact_in1;
459  inertia_ = inertia;
460  }
461  else {
462  //
463  // Don't update the factorization yet
464  //
465  if(!S_indef_) // We must set the inertia if it was definite!
466  inertia_ = S_chol_.inertia();
467  ++S_size_;
468  S_indef_ = true;
469  fact_updated_ = false; // fact_in1_ is meaningless now
470  // We need to keep track of the expected inertia!
471  if( inertia_.zero_eigens == Inertia::UNKNOWN || add_eigen_val == MSADU::EIGEN_VAL_UNKNOWN )
472  inertia_ = Inertia(); // Unknow inertia!
473  else if( add_eigen_val == MSADU::EIGEN_VAL_NEG )
474  ++inertia_.neg_eigens;
475  else if( add_eigen_val == MSADU::EIGEN_VAL_POS )
476  ++inertia_.pos_eigens;
477  else
478  TEUCHOS_TEST_FOR_EXCEPT(true); // Should not happen!
479  }
480  if( throw_exception )
481  throw WarnNearSingularUpdateException(omsg.str(),gamma);
482 }
483 
485  size_type jd
486  ,bool force_refactorization
487  ,EEigenValType drop_eigen_val
488  ,PivotTolerances pivot_tols
489  )
490 {
491  using BLAS_Cpp::upper;
492  using BLAS_Cpp::lower;
493  using DenseLinAlgPack::tri_ele;
494  using DenseLinAlgPack::nonconst_tri_ele;
495  typedef MatrixSymAddDelUpdateable MSADU;
496 
497  assert_initialized();
498 
499  if( jd < 1 || S_size_ < jd )
500  throw std::out_of_range(
501  "MatrixSymAddDelBunchKaufman::delete_update(jd,...): "
502  "Error, the indice jd must be 1 <= jd <= rows()" );
503 
504  bool throw_exception = false; // If true then throw exception
505  std::ostringstream omsg; // Will be set if an exception has to be thrown.
506  value_type gamma; // ...
507 
508  if( !S_indef_ ) {
509  //
510  // The current matrix is p.d. or n.d. and so should the new matrix.
511  //
512  S_chol_.delete_update(jd,force_refactorization,drop_eigen_val,pivot_tols);
513  --S_size_;
514  }
515  else {
516  //
517  // The current matrix is indefinite but the new matrix
518  // could be p.d. or n.d. in which case we could switch
519  // to the cholesky factorization. If the user says the
520  // new matrix will be p.d. or n.d. then we will take
521  // his/her word for it and try the cholesky factorization
522  // otherwise just recompute the indefinite factorization.
523  //
524  // ToDo: (2/2/01): We could also try some more fancy updating
525  // procedures and try to get away with some potentially unstable
526  // methods and resort to the full indefinite factorization if needed.
527  // The sign of the dropped eigen value might tell us what we
528  // might get away with?
529  //
530  // Update the factorization
531  //
532  Inertia inertia = S_chol_.inertia();
533  if( (drop_eigen_val == MSADU::EIGEN_VAL_POS && inertia.pos_eigens == 1 )
534  || (drop_eigen_val == MSADU::EIGEN_VAL_NEG && inertia.neg_eigens == 1 )
535  || S_size_ == 2 )
536  {
537  // Lets take the users word for it and switch to a cholesky factorization.
538  // To do this we will just let S_chol_ do it for us. If the new matrix
539  // turns out not to be what the user says, then we will resize the matrix
540  // to zero and thrown an exception.
541  try {
542  // Delete row and column jd from M
543  DMatrixSlice S = this->S(S_size_).gms();
544  DenseLinAlgPack::delete_row_col( jd, &DenseLinAlgPack::nonconst_tri_ele(S,BLAS_Cpp::lower) );
545  // Setup S_chol and factor the thing
546  S_chol_.init_setup(&S_store1_(),Teuchos::null,0,true,true,true,false,0.0);
547  S_chol_.initialize( this->S(S_size_-1), S_store1_.rows()-1
548  , force_refactorization, Inertia(), PivotTolerances() );
549  }
550  catch( const SingularUpdateException& excpt ) {
551  S_size_ = 0;
552  throw SingularUpdateException(
553  std::string(
554  "MatrixSymAddDelBunchKaufman::delete_update(...) : "
555  "Error, the client incorrectly specified that the new "
556  "matrix would be nonsingular and not indefinite:\n" )
557  + excpt.what()
558  , excpt.gamma
559  );
560  }
561  catch( const WrongInertiaUpdateException& excpt ) {
562  S_size_ = 0;
563  throw WrongInertiaUpdateException(
564  std::string(
565  "MatrixSymAddDelBunchKaufman::delete_update(...) : "
566  "Error, the client incorrectly specified that the new "
567  "matrix would not be indefinite:\n" )
568  + excpt.what()
569  , excpt.gamma
570  );
571  }
572  // If we get here the update succeeded and the new matrix is p.d. or n.d.
573  --S_size_;
574  S_indef_ = false;
575  }
576  else {
577  //
578  // We have been given no indication that the new matrix is p.d. or n.d.
579  // so we will assume it is indefinite and go from there.
580  //
581  DMatrixSlice S = this->S(S_size_).gms();
582  if( force_refactorization ) {
583  //
584  // Perform the refactorization carefully
585  //
586  const bool fact_in1 = !fact_in1_;
587  // Copy the original into the unused storage location
588  DMatrixSliceTriEle DU = this->DU(S_size_,fact_in1);
589  DenseLinAlgPack::assign( &DU, tri_ele(S,lower) );
590  // Delete row and column jd from the storage location for DU
591  DenseLinAlgPack::delete_row_col( jd, &DU );
592  // Now factor the matrix inplace
593  try {
594  factor_matrix(S_size_-1,fact_in1);
595  }
596  catch( const DenseLinAlgLAPack::FactorizationException& excpt ) {
597  omsg
598  << "MatrixSymAddDelBunchKaufman::delete_update(...): "
599  << "Error, singular update but the original matrix was maintianed:\n"
600  << excpt.what();
601  throw SingularUpdateException( omsg.str(), -1.0 );
602  }
603  // Compute the expected inertia
604  Inertia expected_inertia = this->inertia();
605  if( expected_inertia.zero_eigens == Inertia::UNKNOWN || drop_eigen_val == MSADU::EIGEN_VAL_UNKNOWN )
606  expected_inertia = Inertia(); // Unknow inertia!
607  else if( drop_eigen_val == MSADU::EIGEN_VAL_NEG )
608  --expected_inertia.neg_eigens;
609  else if( drop_eigen_val == MSADU::EIGEN_VAL_POS )
610  --expected_inertia.pos_eigens;
611  else
612  TEUCHOS_TEST_FOR_EXCEPT(true); // Should not happen!
613  // Compute the exacted inertia and validate that it is what is expected
614  Inertia inertia;
615  throw_exception = compute_assert_inertia(
616  S_size_-1,fact_in1,expected_inertia,"delete_update",pivot_tols
617  ,&inertia,&omsg,&gamma);
618  // If we get here the factorization worked out.
619  --S_size_;
620  S_indef_ = true;
621  fact_updated_ = true;
622  fact_in1_ = fact_in1;
623  inertia_ = inertia;
624  }
625  // Delete the row and column jd from the original
626  DenseLinAlgPack::delete_row_col( jd, &nonconst_tri_ele(S,lower) );
627  if( !force_refactorization ) {
628  //
629  // The refactorization was not forced
630  //
631  --S_size_;
632  S_indef_ = true;
633  fact_updated_ = false; // fact_in1_ is meaningless now
634  // We need to keep track of the expected inertia!
635  if( inertia_.zero_eigens == Inertia::UNKNOWN || drop_eigen_val == MSADU::EIGEN_VAL_UNKNOWN )
636  inertia_ = Inertia(); // Unknow inertia!
637  else if( drop_eigen_val == MSADU::EIGEN_VAL_NEG )
638  --inertia_.neg_eigens;
639  else if( drop_eigen_val == MSADU::EIGEN_VAL_POS )
640  --inertia_.pos_eigens;
641  else
642  TEUCHOS_TEST_FOR_EXCEPT(true); // Should not happen!
643  }
644  }
645  }
646  if( throw_exception )
647  throw WarnNearSingularUpdateException(omsg.str(),gamma);
648 }
649 
650 // Overridden from MatrixSymOpNonsingSerial
651 
653 {
654  return S_size_;
655 }
656 
657 std::ostream& MatrixSymAddDelBunchKaufman::output(std::ostream& out) const
658 {
659  if( S_size_ ) {
660  out << "Unfactored symmetric matrix stored as lower triangle (ignore upper nonzeros):\n"
661  << S(S_size_).gms();
662  if( S_indef_ ) {
663  out << "Upper symmetric indefinite factor DU returned from sytrf(...) (ignore lower nonzeros):\n"
664  << DU(S_size_,fact_in1_).gms();
665  out << "Permutation array IPIV returned from sytrf(...):\n";
666  for( size_type i = 0; i < S_size_; ++i )
667  out << " " << IPIV_[i];
668  out << std::endl;
669  }
670  else {
671  out << "Upper cholesky factor (ignore lower nonzeros):\n" << DU(S_size_,true).gms();
672  }
673  }
674  else {
675  out << "0 0\n";
676  }
677  return out;
678 }
679 
681  DVectorSlice* y, value_type a, BLAS_Cpp::Transp M_trans
682  ,const DVectorSlice& x, value_type b
683  ) const
684 {
685  DenseLinAlgPack::Vp_MtV_assert_sizes( y->dim(), S_size_, S_size_, M_trans, x.dim() );
686  // Use the unfactored matrix since this is more accurate!
687  DenseLinAlgPack::Vp_StMtV( y, a, S(S_size_), BLAS_Cpp::no_trans, x, b );
688 }
689 
691  DVectorSlice* y, BLAS_Cpp::Transp M_trans
692  ,const DVectorSlice& x
693  ) const
694 {
695  const size_type n = S_size_;
696  DenseLinAlgPack::Vp_MtV_assert_sizes( y->dim(), n, n, M_trans, x.dim() );
697  if( S_indef_ ) {
698  // Update the factorzation if needed
699  if(!fact_updated_) {
700  const bool fact_in1 = true;
702  *nc_this = const_cast<MatrixSymAddDelBunchKaufman*>(this);
703  nc_this->copy_and_factor_matrix(S_size_,fact_in1); // May throw exceptions!
704  nc_this->fact_updated_ = true;
705  nc_this->fact_in1_ = fact_in1;
706  }
707  *y = x;
708  DenseLinAlgLAPack::sytrs(
709  DU(S_size_,fact_in1_), &const_cast<IPIV_t&>(IPIV_)[0]
710  , &DMatrixSlice(y->raw_ptr(),n,n,n,1), &WORK_() );
711  }
712  else {
713  AbstractLinAlgPack::V_InvMtV( y, S_chol_, M_trans, x );
714  }
715 }
716 
717 // Private
718 
719 void MatrixSymAddDelBunchKaufman::assert_initialized() const
720 {
721  if(!S_size_)
722  throw std::logic_error(
723  "MatrixSymAddDelBunchKaufman::assert_initialized() : "
724  "Error, this matrix is not initialized yet" );
725 }
726 
727 void MatrixSymAddDelBunchKaufman::resize_DU_store( bool in_store1 )
728 {
729  if(!in_store1 && S_store2_.rows() < S_store1_.rows() )
730  S_store2_.resize( S_store1_.rows(), S_store1_.cols() );
731 }
732 
733 void MatrixSymAddDelBunchKaufman::copy_and_factor_matrix(
734  size_type S_size, bool fact_in1 )
735 {
737  &DU(S_size,fact_in1)
738  , DenseLinAlgPack::tri_ele(S(S_size).gms(),BLAS_Cpp::lower) );
739  factor_matrix(S_size,fact_in1);
740 }
741 
742 void MatrixSymAddDelBunchKaufman::factor_matrix( size_type S_size, bool fact_in1 )
743 {
744  // Resize the workspace first
745  const size_type opt_block_size = 64; // This is an estimate (see sytrf(...))
746  if( WORK_.dim() < S_store1_.rows() * opt_block_size )
747  WORK_.resize(S_store1_.rows()*opt_block_size);
748  if( IPIV_.size() < S_store1_.rows() )
749  IPIV_.resize(S_store1_.rows());
750  // Factor the matrix (will throw FactorizationException if singular)
751  DenseLinAlgLAPack::sytrf( &DU(S_size,fact_in1), &IPIV_[0], &WORK_() );
752 }
753 
754 bool MatrixSymAddDelBunchKaufman::compute_assert_inertia(
755  size_type S_size, bool fact_in1, const Inertia& exp_inertia, const char func_name[]
756  ,PivotTolerances pivot_tols, Inertia* comp_inertia, std::ostringstream* omsg, value_type* gamma )
757 {
758  bool throw_exception = false;
759  *gamma = 0.0;
760  // Here we will compute the inertia given IPIV[] and D[] as described in the documentation
761  // for dsytrf(...) (see the source code).
762  const DMatrixSlice DU = this->DU(S_size,fact_in1).gms();
763  const size_type n = DU.rows();
764  Inertia inertia(0,0,0);
765  value_type max_diag = 0.0;
766  value_type min_diag = std::numeric_limits<value_type>::max();
767  for( size_type k = 1; k <= n; ) {
768  const FortranTypes::f_int k_p = IPIV_[k-1];
769  if( k_p > 0 ) {
770  // D(k,k) is a 1x1 block.
771  // Lets get the eigen value from the sign of D(k,k)
772  const value_type D_k_k = DU(k,k), abs_D_k_k = std::fabs(D_k_k);
773  if( D_k_k > 0.0 )
774  ++inertia.pos_eigens;
775  else
776  ++inertia.neg_eigens;
777  if(abs_D_k_k > max_diag) max_diag = abs_D_k_k;
778  if(abs_D_k_k < min_diag) min_diag = abs_D_k_k;
779  k++;
780  }
781  else {
782  // D(k:k+1,k:k+1) is a 2x2 block.
783  // This represents one positive eigen value and
784  // on negative eigen value
785  TEUCHOS_TEST_FOR_EXCEPT( !( IPIV_[k] == k_p ) ); // This is what the documentation for xSYTRF(...) says!
786  ++inertia.pos_eigens;
787  ++inertia.neg_eigens;
788  // To find the largest and smallest diagonals of U for L*U we must perform Gaussian
789  // elimination on this 2x2 block:
790  const value_type // [ a b ] = D(k:k+1,k:k+1)
791  a = DU(k,k), b = DU(k,k+1), c = DU(k+1,k+1), // [ b c ]
792  abs_a = std::fabs(a), abs_b = std::fabs(b);
793  value_type pivot_1, pivot_2;
794  if( abs_a > abs_b ) { // Pivot on a = D(k,k)
795  pivot_1 = abs_a; // [ 1 ] * [ a b ] = [ a b ]
796  pivot_2 = std::fabs(c - b*b/a); // [ -b/a 1 ] [ b c ] [ 0 c - b*b/a ]
797  }
798  else { // Pivot on b = D(k+1,k) = D(k,k+1)
799  pivot_1 = abs_b; // [ 1 ] * [ b c ] = [ b c ]
800  pivot_2 = std::fabs(b - a*c/b); // [ -a/b 1 ] [ a b ] [ 0 b - a*c/b ]
801  }
802  if(pivot_1 > max_diag) max_diag = pivot_1;
803  if(pivot_1 < min_diag) min_diag = pivot_1;
804  if(pivot_2 > max_diag) max_diag = pivot_2;
805  if(pivot_2 < min_diag) min_diag = pivot_2;
806  k+=2;
807  }
808  }
809  // gamma = min(|diag(i)|)/max(|diag(i)|)
810  *gamma = min_diag / max_diag;
811  // Now validate that the inertia is what is expected
812  const bool
813  wrong_inertia =
814  ( exp_inertia.neg_eigens != Inertia::UNKNOWN
815  && exp_inertia.neg_eigens != inertia.neg_eigens )
816  || ( exp_inertia.pos_eigens != Inertia::UNKNOWN
817  && exp_inertia.pos_eigens != inertia.pos_eigens ) ;
818  PivotTolerances use_pivot_tols = S_chol_.pivot_tols();
819  if( pivot_tols.warning_tol != PivotTolerances::UNKNOWN )
820  use_pivot_tols.warning_tol = pivot_tols.warning_tol;
821  if( pivot_tols.singular_tol != PivotTolerances::UNKNOWN )
822  use_pivot_tols.singular_tol = pivot_tols.singular_tol;
823  if( pivot_tols.wrong_inertia_tol != PivotTolerances::UNKNOWN )
824  use_pivot_tols.wrong_inertia_tol = pivot_tols.wrong_inertia_tol;
825  throw_exception = (
826  *gamma == 0.0
827  || !wrong_inertia && *gamma <= use_pivot_tols.warning_tol
828  || !wrong_inertia && *gamma <= use_pivot_tols.singular_tol
829  || wrong_inertia
830  );
831  // Create message and throw exceptions
832  std::ostringstream onum_msg;
833  if(throw_exception) {
834  if(wrong_inertia) {
835  onum_msg
836  << "inertia = ("
837  << inertia.neg_eigens << "," << inertia.zero_eigens << "," << inertia.pos_eigens
838  << ") != expected_inertia = ("
839  << exp_inertia.neg_eigens << "," << exp_inertia.zero_eigens << "," << exp_inertia.pos_eigens << ")"
840  << std::endl;
841  }
842  onum_msg
843  << "gamma = min(|diag(D)(k)|)/max(|diag(D)(k)|) = " << min_diag << "/" << max_diag
844  << " = " << *gamma;
845  *omsg
846  << "MatrixSymAddDelBunchKaufman::"<<func_name<<"(...): ";
847  if( *gamma == 0.0 || (!wrong_inertia && *gamma <= use_pivot_tols.singular_tol) ) {
848  *omsg
849  << "Singular update!\n" << onum_msg.str() << " <= singular_tol = "
850  << use_pivot_tols.singular_tol;
851  throw SingularUpdateException( omsg->str(), *gamma );
852  }
853  else if( wrong_inertia && *gamma <= use_pivot_tols.singular_tol ) {
854  *omsg
855  << "Singular update!\n" << onum_msg.str() << " <= wrong_inertia_tol = "
856  << use_pivot_tols.wrong_inertia_tol;
857  throw SingularUpdateException( omsg->str(), *gamma );
858  }
859  else if( wrong_inertia ) {
860  *omsg
861  << "Indefinite update!\n" << onum_msg.str() << " >= wrong_inertia_tol = "
862  << use_pivot_tols.wrong_inertia_tol;
863  throw WrongInertiaUpdateException( omsg->str(), *gamma );
864  }
865  else if( !wrong_inertia && *gamma <= use_pivot_tols.warning_tol ) {
866  *omsg
867  << "Warning, near singular update!\nsingular_tol = " << use_pivot_tols.singular_tol
868  << " < " << onum_msg.str() << " <= warning_tol = "
869  << use_pivot_tols.warning_tol;
870  // Don't throw the exception till the very end!
871  }
872  else {
873  TEUCHOS_TEST_FOR_EXCEPT(true); // Only local programming error?
874  }
875  }
876  // Return
877  *comp_inertia = inertia;
878  return throw_exception;
879 }
880 
881 } // end namespace ConstrainedOptPack
void Vp_StMtV(DVectorSlice *vs_lhs, value_type alpha, BLAS_Cpp::Transp trans_rhs1, const DVectorSlice &vs_rhs2, value_type beta) const
value_type transVtInvMtV(const Vector &v_rhs1, const MatrixNonsing &M_rhs2, BLAS_Cpp::Transp trans_rhs2, const Vector &v_rhs3)
TypeTo & implicit_ref_cast(TypeFrom &t)
This class maintains the factorization of symmetric indefinite matrix using a Bunch & Kaufman factori...
void assign(MatrixOp *M_lhs, const MatrixOp &M_rhs, BLAS_Cpp::Transp trans_rhs)
size_t size_type
void augment_update(const DVectorSlice *t, value_type alpha, bool force_refactorization, EEigenValType add_eigen_val, PivotTolerances pivot_tols)
void V_InvMtV(VectorMutable *v_lhs, const MatrixNonsing &M_rhs1, BLAS_Cpp::Transp trans_rhs1, const Vector &v_rhs2)
void Vp_StMtV(VectorMutable *v_lhs, value_type alpha, const MatrixOp &M_rhs1, BLAS_Cpp::Transp trans_rhs1, const Vector &v_rhs2, value_type beta=1.0)
MatrixSymAddDelBunchKaufman()
Initializes with 0x0 and pivot_tols == (0.0,0.0,0.0).
Transp
void delete_update(size_type jd, bool force_refactorization, EEigenValType drop_eigen_val, PivotTolerances pivot_tols)
void V_InvMtV(DVectorSlice *vs_lhs, BLAS_Cpp::Transp trans_rhs1, const DVectorSlice &vs_rhs2) const
#define TEUCHOS_TEST_FOR_EXCEPT(throw_exception_test)