MOOCHO (Single Doxygen Collection)  Version of the Day
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
AbstractLinAlgPack_MatrixSymPosDefCholFactor.cpp
Go to the documentation of this file.
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 #include <assert.h>
43 
44 #include <limits>
45 #include <vector>
46 
62 #include "DenseLinAlgLAPack.hpp"
70 #include "Teuchos_Assert.hpp"
71 #include "Teuchos_FancyOStream.hpp"
72 
73 #ifdef HAVE_MOOCHO_FORTRAN
74 # define ALAP_DCHUD_DECL FORTRAN_FUNC_DECL_UL( void, DCHUD, dchud )
75 # define ALAP_DCHUD_CALL FORTRAN_FUNC_CALL_UL( DCHUD, dchud )
76 #else
77 # define ALAP_DCHUD_DECL void dchud_c
78 # define ALAP_DCHUD_CALL dchud_c
79 #endif
80 
81 // Helper functions
82 extern "C" {
84  , const FortranTypes::f_int& LDR, const FortranTypes::f_int& P
86  , const FortranTypes::f_int& LDZ, const FortranTypes::f_int& NZ
89 } // end extern "C"
90 
91 namespace {
92 //
93 template< class T >
94 inline
95 T my_max( const T& v1, const T& v2 ) { return v1 > v2 ? v1 : v2; }
96 //
97 template< class T >
98 inline
99 T my_min( const T& v1, const T& v2 ) { return v1 < v2 ? v1 : v2; }
100 } // end namespace
101 
102 namespace AbstractLinAlgPack {
103 
104 // Constructors/initalizers
105 
106 MatrixSymPosDefCholFactor::MatrixSymPosDefCholFactor()
107  : maintain_original_(true), maintain_factor_(false)
108  , factor_is_updated_(false), allocates_storage_(true)
109  , max_size_(0), M_size_(0), M_l_r_(1), M_l_c_(1)
110  , U_l_r_(1), U_l_c_(1), scale_(0.0), is_diagonal_(false)
111  , pivot_tols_(0.0,0.0,0.0) // This is what DPOTRF(...) uses!
112 {}
113 
114 MatrixSymPosDefCholFactor::MatrixSymPosDefCholFactor(
115  DMatrixSlice *MU_store
116  ,const release_resource_ptr_t& release_resource_ptr
117  ,size_type max_size
118  ,bool maintain_original
119  ,bool maintain_factor
120  ,bool allow_factor
121  ,bool set_full_view
122  ,value_type scale
123  )
124  : pivot_tols_(0.0,0.0,0.0) // This is what DPOTRF(...) uses!
125 {
126  init_setup(MU_store,release_resource_ptr,max_size,maintain_original
127  ,maintain_factor,allow_factor,set_full_view,scale);
128 }
129 
130 void MatrixSymPosDefCholFactor::init_setup(
131  DMatrixSlice *MU_store
132  ,const release_resource_ptr_t& release_resource_ptr
133  ,size_type max_size
134  ,bool maintain_original
135  ,bool maintain_factor
136  ,bool allow_factor
137  ,bool set_full_view
138  ,value_type scale
139  )
140 {
141  TEUCHOS_TEST_FOR_EXCEPT( !( maintain_original || maintain_factor ) );
142  if( MU_store == NULL ) {
143  maintain_original_ = maintain_original;
144  maintain_factor_ = maintain_factor;
145  factor_is_updated_ = maintain_factor;
146  allocates_storage_ = true; // We will be able to allocate our own storage!
147  release_resource_ptr_ = Teuchos::null; // Free any bound resource
148  MU_store_.bind( DMatrixSlice(NULL,0,0,0,0) ); // Unbind this!
149  max_size_ = max_size;
150  M_size_ = 0;
151  M_l_r_ = M_l_c_ = 1;
152  if( !maintain_factor && !allow_factor )
153  U_l_r_ = 0; // Do not allow the factor to be computed
154  else
155  U_l_r_ = 1; // Allow the factor to be computed
156  scale_ = +1.0;
157  is_diagonal_ = false;
158  }
159  else {
160  TEUCHOS_TEST_FOR_EXCEPT( !( MU_store->rows() ) );
161  allocates_storage_ = false; // The client allocated the storage!
162  MU_store_.bind(*MU_store);
163  release_resource_ptr_ = release_resource_ptr;
164  max_size_ = my_min( MU_store->rows(), MU_store->cols() ) - 1;
165  if( set_full_view ) {
166  TEUCHOS_TEST_FOR_EXCEPT( !( scale != 0.0 ) );
167  this->set_view(
168  max_size_
169  ,scale,maintain_original,1,1
170  ,maintain_factor, allow_factor ? 1: 0, allow_factor ? 1 : 0
171  );
172  }
173  else {
174  this->set_view(
175  0
176  ,0.0,maintain_original,1,1
177  ,maintain_factor, allow_factor ? 1: 0, allow_factor ? 1 : 0
178  );
179  }
180  }
181 }
182 
183 void MatrixSymPosDefCholFactor::set_view(
184  size_t M_size
185  ,value_type scale
186  ,bool maintain_original
187  ,size_t M_l_r
188  ,size_t M_l_c
189  ,bool maintain_factor
190  ,size_t U_l_r
191  ,size_t U_l_c
192  )
193 {
194  TEUCHOS_TEST_FOR_EXCEPT( !( maintain_original || maintain_factor ) );
195  if( max_size_ )
196  allocate_storage(max_size_);
197  else
198  allocate_storage( my_max( M_l_r + M_size, M_l_c + M_size ) - 1 );
199  // Check the preconditions
200  if( maintain_original ) {
201  TEUCHOS_TEST_FOR_EXCEPT( !( 1 <= M_l_r && M_l_r <= M_l_c ) );
202  TEUCHOS_TEST_FOR_EXCEPT( !( M_l_r+M_size <= MU_store_.rows() ) );
203  }
204  if( maintain_factor ) {
205  TEUCHOS_TEST_FOR_EXCEPT( !( 1 <= U_l_r && U_l_r >= U_l_c ) );
206  TEUCHOS_TEST_FOR_EXCEPT( !( U_l_c+M_size <= MU_store_.cols() ) );
207  }
208  // Set the members
209  maintain_original_ = maintain_original;
210  maintain_factor_ = maintain_factor;
211  is_diagonal_ = false;
212  if( M_size ) {
213  max_size_ = my_min( MU_store_.rows() - U_l_r, MU_store_.cols() - U_l_c );
214  M_size_ = M_size;
215  scale_ = scale;
216  M_l_r_ = M_l_r;
217  M_l_c_ = M_l_c;
218  U_l_r_ = U_l_r;
219  U_l_c_ = U_l_c;
220  factor_is_updated_ = maintain_factor;
221  }
222  else {
223  max_size_ = 0;
224  M_size_ = 0;
225  scale_ = 0.0;
226  M_l_r_ = 0;
227  M_l_c_ = 0;
228  U_l_r_ = U_l_r;
229  U_l_c_ = U_l_r;
230  factor_is_updated_ = maintain_factor;
231  }
232 }
233 
234 void MatrixSymPosDefCholFactor::pivot_tols( PivotTolerances pivot_tols )
235 {
236  pivot_tols_ = pivot_tols;
237 }
238 
239 MatrixSymAddDelUpdateable::PivotTolerances MatrixSymPosDefCholFactor::pivot_tols() const
240 {
241  return pivot_tols_;
242 }
243 
244 // Overridden from MatrixBase
245 
247 {
248  return M_size_;
249 }
250 
251 // Overridden from MatrixOp
252 
253 void MatrixSymPosDefCholFactor::zero_out()
254 {
255  this->init_identity(this->space_cols(),0.0);
256 }
257 
258 std::ostream& MatrixSymPosDefCholFactor::output(std::ostream& out_arg) const
259 {
260  Teuchos::RCP<Teuchos::FancyOStream> out = Teuchos::getFancyOStream(Teuchos::rcp(&out_arg,false));
261  Teuchos::OSTab tab(out);
262  if( M_size_ ) {
263  if( maintain_original_ ) {
264  *out
265  << "Unfactored symmetric matrix stored as lower triangle (ignore upper nonzeros):\n"
266  << M().gms();
267  }
268  if( factor_is_updated_ ) {
269  *out
270  << "Matrix scaling M = scale*U'*U, scale = " << scale_ << std::endl
271  << "Upper cholesky factor U (ignore lower nonzeros):\n"
272  << U().gms();
273  }
274  }
275  else {
276  *out << "0 0\n";
277  }
278  return out_arg;
279 }
280 
282  MatrixOp* m_lhs, value_type alpha
283  ,BLAS_Cpp::Transp trans_rhs
284  ) const
285 {
286  MatrixSymOpGetGMSSymMutable
287  *symwo_gms_lhs = dynamic_cast<MatrixSymOpGetGMSSymMutable*>(m_lhs);
288  if(!symwo_gms_lhs)
289  return false;
290  MatrixDenseSymMutableEncap sym_lhs(symwo_gms_lhs);
291  const DMatrixSliceSym M = this->M();
293  &DMatrixSliceTriEle(sym_lhs().gms(),sym_lhs().uplo())
294  ,alpha
295  ,DMatrixSliceTriEle(M.gms(),M.uplo())
296  );
297 
298  return true;
299 }
300 
302  value_type alpha,const MatrixOp& M_rhs, BLAS_Cpp::Transp trans_rhs
303  )
304 {
306  !maintain_original_, std::logic_error
307  ,"MatrixSymPosDefCholFactor::Mp_StM(alpha,M_rhs,trans_rhs): Error, Current implementation "
308  "can not perform this operation unless the original matrix is being maintained." );
309  // Perform the operation
310  bool did_op = false;
311  bool diag_op = false;
312  if(const MatrixSymOpGetGMSSym *symwo_gms_rhs = dynamic_cast<const MatrixSymOpGetGMSSym*>(&M_rhs)) {
313  DMatrixSliceSym M = this->M();
314  MatrixDenseSymEncap sym_rhs(*symwo_gms_rhs);
316  &DMatrixSliceTriEle(M.gms(),M.uplo())
317  ,alpha
318  ,DMatrixSliceTriEle(sym_rhs().gms(),sym_rhs().uplo())
319  );
320  did_op = true;
321  diag_op = false;
322  }
323  else if(const MatrixSymDiag *symwo_diag_rhs = dynamic_cast<const MatrixSymDiag*>(&M_rhs)) {
324  DMatrixSliceSym M = this->M();
325  VectorDenseEncap sym_rhs_diag(symwo_diag_rhs->diag());
326  LinAlgOpPack::Vp_StV( &M.gms().diag(), alpha, sym_rhs_diag() );
327  did_op = true;
328  diag_op = true;
329  }
330  else if(const MatrixSymOp *symwo_rhs = dynamic_cast<const MatrixSymOp*>(&M_rhs)) {
331  // ToDo: Implement this!
332  }
333  // Properly update the state of *this.
334  // If only the original is updated
335  if(did_op) {
336  if( diag_op && is_diagonal_ )
337  this->init_diagonal(VectorMutableDense(this->M().gms().diag(),Teuchos::null));
338  else
339  this->initialize(this->M());
340  return true;
341  }
342  return false;
343 }
344 
346  const MatrixOp &mwo_rhs
347  ,BLAS_Cpp::Transp M_trans
348  ,value_type alpha
349  ,value_type beta
350  )
351 {
352  MatrixDenseSymMutableEncap sym_gms_lhs(this);
353  const MatrixOpGetGMS *mwo_rhs_gms = dynamic_cast<const MatrixOpGetGMS*>(&mwo_rhs);
354  if(mwo_rhs_gms) {
355  TEUCHOS_TEST_FOR_EXCEPT(true); // ToDo: Implement
356  return true;
357  }
358  else {
359  // Here we will give up on symmetry and just compute the whole product S = op(mwo_rhs)*op(mwo_rhs')
360  DenseLinAlgPack::DMatrixSliceTriEle tri_ele_gms_lhs = DenseLinAlgPack::tri_ele(sym_gms_lhs().gms(),sym_gms_lhs().uplo());
361  if(beta==0.0) DenseLinAlgPack::assign( &tri_ele_gms_lhs, 0.0 );
362  else if(beta!=1.0) DenseLinAlgPack::Mt_S( &tri_ele_gms_lhs, beta );
363  const VectorSpace &spc = mwo_rhs.space_rows();
364  const index_type m = spc.dim();
365  VectorSpace::multi_vec_mut_ptr_t S = spc.create_members(m);
366  S->zero_out();
367  LinAlgOpPack::M_MtM( S.get(), mwo_rhs, M_trans, mwo_rhs, BLAS_Cpp::trans_not(M_trans) );
368  // Copy S into sym_ghs_lhs
369  if( sym_gms_lhs().uplo() == BLAS_Cpp::lower ) {
370  for( index_type j = 1; j <= m; ++j ) {
371  LinAlgOpPack::Vp_StV( &sym_gms_lhs().gms().col(j)(j,m), alpha, VectorDenseEncap(*S->col(j)->sub_view(j,m))() );
372  }
373  }
374  else {
375  for( index_type j = 1; j <= m; ++j ) {
376  LinAlgOpPack::Vp_StV( &sym_gms_lhs().gms().col(j)(1,j), alpha, VectorDenseEncap(*S->col(j)->sub_view(1,j))() );
377  }
378  }
379  return true;
380  }
381  return false;
382 }
383 
384 // Overridden from MatrixOpSerial
385 
388  ,const DVectorSlice& x, value_type b
389  ) const
390 {
392  using BLAS_Cpp::no_trans;
393  using BLAS_Cpp::trans;
394 #ifdef PROFILE_HACK_ENABLED
395  ProfileHackPack::ProfileTiming profile_timing( "MatrixSymPosDefCholFactor::Vp_StMtV(...DVectorSlice...)" );
396 #endif
397  assert_initialized();
398 
400  implicit_ref_cast<const MatrixBase>(*this).rows(), cols(), no_trans, x.dim() );
401 
402  if( maintain_original_ ) {
403  //
404  // M = symmetric
405  //
406  // y = b*y + a*M*x
407  //
408  DenseLinAlgPack::Vp_StMtV( y, a, M(), no_trans, x, b );
409  }
410  else {
411  //
412  // M = scale*U'*U
413  //
414  // y = b*y + a*op(M)*x
415  // = b*y = scale*a*U'*U*x
416  //
417  const DMatrixSliceTri
418  U = this->U();
419 
420  if( b == 0.0 ) {
421  // No temporary needed
422  //
423  // y = U*x
424  DenseLinAlgPack::V_MtV( y, U, no_trans, x );
425  // y = U'*y
426  DenseLinAlgPack::V_MtV( y, U, trans, *y );
427  // y *= scale*a
428  if( a != 1.0 || scale_ != 1.0 )
429  DenseLinAlgPack::Vt_S( y, scale_*a );
430  }
431  else {
432  // We need a temporary
433  DVector t;
434  // t = U*x
435  DenseLinAlgPack::V_MtV( &t, U, no_trans, x );
436  // t = U'*t
437  DenseLinAlgPack::V_MtV( &t(), U, trans, t() );
438  // y *= b
439  if(b != 1.0)
440  DenseLinAlgPack::Vt_S( y, b );
441  // y += scale*a*t
442  DenseLinAlgPack::Vp_StV( y, scale_*a, t() );
443  }
444  }
445 }
446 
449  ,const SpVectorSlice& x, value_type b
450  ) const
451 {
452  using BLAS_Cpp::no_trans;
453  using BLAS_Cpp::trans;
454 #ifdef PROFILE_HACK_ENABLED
455  ProfileHackPack::ProfileTiming profile_timing( "MatrixSymPosDefCholFactor::Vp_StMtV(...SpVectorSlice...)" );
456 #endif
457  assert_initialized();
458  if( maintain_original_ ) {
459  const DMatrixSlice M = this->M().gms(); // This is lower triangular!
460  const size_type n = M.rows();
461  DenseLinAlgPack::Vp_MtV_assert_sizes( y->dim(), n, n, no_trans, x.dim() );
462  DenseLinAlgPack::Vt_S(y,b); // y = b*y
463  //
464  // Compute product column by column corresponding to x_itr->index() + x.offset()
465  //
466  // y += a * M * e(i) * x(i)
467  //
468  // [ y(1:i-1) ] += a * x(i) * [ ... M(1:i-1,i) ... ] stored as M(i,1:i-1)
469  // [ y(i:n) ] [ ... M(i:n,i) ... ] stored as M(i:n,i)
470  //
471  for( SpVectorSlice::const_iterator x_itr = x.begin(); x_itr != x.end(); ++x_itr ) {
472  const size_type i = x_itr->index() + x.offset();
473  if( i > 1 )
474  DenseLinAlgPack::Vp_StV( &(*y)(1,i-1), a * x_itr->value(), M.row(i)(1,i-1) );
475  DenseLinAlgPack::Vp_StV( &(*y)(i,n), a * x_itr->value(), M.col(i)(i,n) );
476  }
477  }
478  else {
479  MatrixOpSerial::Vp_StMtV(y,a,M_trans,x,b); // ToDo: Specialize when needed!
480  }
481 }
482 
484  DVectorSlice* y, value_type a, const GenPermMatrixSlice& P, BLAS_Cpp::Transp P_trans
485  , BLAS_Cpp::Transp H_trans, const DVectorSlice& x, value_type b
486  ) const
487 {
488 #ifdef PROFILE_HACK_ENABLED
489  ProfileHackPack::ProfileTiming profile_timing( "MatrixSymPosDefCholFactor::Vp_StPtMtV(...DVectorSlice...)" );
490 #endif
491  assert_initialized();
492  MatrixOpSerial::Vp_StPtMtV(y,a,P,P_trans,H_trans,x,b); // ToDo: Specialize when needed!
493 }
494 
496  DVectorSlice* y, value_type a, const GenPermMatrixSlice& P, BLAS_Cpp::Transp P_trans
497  , BLAS_Cpp::Transp H_trans, const SpVectorSlice& x, value_type b
498  ) const
499 {
500 #ifdef PROFILE_HACK_ENABLED
501  ProfileHackPack::ProfileTiming profile_timing( "MatrixSymPosDefCholFactor::Vp_StPtMtV(...SpVectorSlice...)" );
502 #endif
503  assert_initialized();
504  if( maintain_original_ ) {
505  DenseLinAlgPack::Vt_S(y,b); // y = b*y
506  const DMatrixSlice M = this->M().gms(); // This is lower triangular!
507  // Compute product column by corresponding to x_itr->index() + x.offset()
508  /*
509  if( P.is_identity() ) {
510  TEUCHOS_TEST_FOR_EXCEPT(true); // ToDo: Implement
511  }
512  else {
513  for( SpVectorSlice::const_iterator x_itr = x.begin(); x_itr != x.end(); ++x_itr ) {
514  const size_type i = x_itr->index() + x.offset();
515 
516 
517  }
518  }
519  */
520  MatrixOpSerial::Vp_StPtMtV(y,a,P,P_trans,H_trans,x,b); // ToDo: Specialize when needed!
521  }
522  else {
523  MatrixOpSerial::Vp_StPtMtV(y,a,P,P_trans,H_trans,x,b); // ToDo: Specialize when needed!
524  }
525 }
526 
527 // Overridden from MatrixSymOpSerial
528 
531  , EMatRhsPlaceHolder dummy_place_holder
532  , const GenPermMatrixSlice& P, BLAS_Cpp::Transp P_trans
533  , value_type b ) const
534 {
535  using BLAS_Cpp::no_trans;
536  using BLAS_Cpp::trans;
537 #ifdef PROFILE_HACK_ENABLED
538  ProfileHackPack::ProfileTiming profile_timing( "MatrixSymPosDefCholFactor::Mp_StPtMtP(...)" );
539 #endif
540  assert_initialized();
541  if( !maintain_original_ ) {
542  MatrixSymOpSerial::Mp_StPtMtP(S,a,dummy_place_holder,P,P_trans,b);
543  }
544  else {
545  MatrixSymOpSerial::Mp_StPtMtP(S,a,dummy_place_holder,P,P_trans,b); // ToDo: Override when needed!
546  }
547 }
548 
549 // Overridden from MatrixNonsingSerial
550 
552  DVectorSlice* y, BLAS_Cpp::Transp M_trans, const DVectorSlice& x
553  ) const
554 {
555  using BLAS_Cpp::no_trans;
556  using BLAS_Cpp::trans;
557 #ifdef PROFILE_HACK_ENABLED
558  ProfileHackPack::ProfileTiming profile_timing( "MatrixSymPosDefCholFactor::V_InvMtV(...DVectorSlice...)" );
559 #endif
560  assert_initialized();
561 
562  //
563  // M = scale*U'*U
564  //
565  // y = inv(op(M))*x
566  // =>
567  // op(M)*y = x
568  // =>
569  // scale*U'*U*y = x
570  // =>
571  // y = (1/scale)*inv(U)*inv(U')*x
572  //
573  update_factorization();
574  const DMatrixSliceTri
575  U = this->U();
576  DenseLinAlgPack::Vp_MtV_assert_sizes( y->dim(), U.rows(), U.cols(), no_trans, x.dim() );
577  // y = inv(U')*x
578  DenseLinAlgPack::V_InvMtV( y, U, trans, x );
579  // y = inv(U)*y
580  DenseLinAlgPack::V_InvMtV( y, U, no_trans, *y );
581  // y *= (1/scale)
582  if( scale_ != 1.0 )
583  DenseLinAlgPack::Vt_S( y, 1.0/scale_ );
584 }
585 
587  DVectorSlice* y, BLAS_Cpp::Transp M_trans, const SpVectorSlice& x
588  ) const
589 {
590 #ifdef PROFILE_HACK_ENABLED
591  ProfileHackPack::ProfileTiming profile_timing( "MatrixSymPosDefCholFactor::V_InvMtV(...SpVectorSlice...)" );
592 #endif
593  assert_initialized();
594  MatrixNonsingSerial::V_InvMtV(y,M_trans,x);
595 }
596 
597 // Overridden from MatrixSymNonsingSerial
598 
600  DMatrixSliceSym* S, value_type a, const MatrixOpSerial& B
601  , BLAS_Cpp::Transp B_trans, EMatrixDummyArg dummy_arg
602  ) const
603 {
604 #ifdef PROFILE_HACK_ENABLED
605  ProfileHackPack::ProfileTiming profile_timing( "MatrixSymPosDefCholFactor::M_StMtInvMtM(...)" );
606 #endif
607 
608 // // Uncomment to use the defalut implementation (for debugging)
609 // MatrixSymFactorized::M_StMtInvMtM(S,a,B,B_trans,dummy_arg); return;
610 
611  using BLAS_Cpp::trans;
612  using BLAS_Cpp::no_trans;
613  using BLAS_Cpp::trans_not;
614  using BLAS_Cpp::upper;
615  using BLAS_Cpp::nonunit;
616  using DenseLinAlgPack::tri;
617  using DenseLinAlgPack::syrk;
619  using LinAlgOpPack::assign;
620 
621  assert_initialized();
622  update_factorization();
623 
624  //
625  // S = a * op(B) * inv(M) * op(B)'
626  //
627  // M = scale*U'*U
628  // =>
629  // inv(M) = scale*inv(U'*U) = scale*inv(U)*inv(U')
630  // =>
631  // S = scale*a * op(B) * inv(U) * inv(U') * op(B)'
632  //
633  // T = op(B)'
634  //
635  // T = inv(U') * T (inplace with BLAS)
636  //
637  // S = scale*a * T' * T
638  //
640  rows(), cols(), no_trans
641  ,B.rows(), B.cols(), trans_not(B_trans)
642  );
644  S->rows(), S->cols(), no_trans
645  ,B.rows(), B.cols(), B_trans
646  ,B.rows(), B.cols(), trans_not(B_trans)
647  );
648  // T = op(B)'
649  DMatrix T;
650  assign( &T, B, trans_not(B_trans) );
651  // T = inv(U') * T (inplace with BLAS)
652  M_StInvMtM( &T(), 1.0, this->U(), trans, T(), no_trans );
653  // S = scale*a * T' * T
654  syrk( trans, scale_*a, T(), 0.0, S );
655 }
656 
657 // Overridden from MatrixSymDenseInitialize
658 
659 void MatrixSymPosDefCholFactor::initialize( const DMatrixSliceSym& M )
660 {
661  // Initialize without knowing the inertia but is must be p.d.
662  this->initialize(
663  M, M.rows(), maintain_factor_
664  , MatrixSymAddDelUpdateable::Inertia()
665  , MatrixSymAddDelUpdateable::PivotTolerances()
666  );
667 }
668 
669 // Overridden from MatrixSymOpGetGMSSym
670 
671 const DenseLinAlgPack::DMatrixSliceSym MatrixSymPosDefCholFactor::get_sym_gms_view() const
672 {
674  !maintain_original_, std::logic_error
675  ,"MatrixSymPosDefCholFactor::get_sym_gms_view(): Error, maintain_original must be "
676  "true in order to call this method!" );
677  return this->M();
678 }
679 
680 void MatrixSymPosDefCholFactor::free_sym_gms_view(const DenseLinAlgPack::DMatrixSliceSym* sym_gms_view) const
681 {
683  !maintain_original_, std::logic_error
684  ,"MatrixSymPosDefCholFactor::free_sym_gms_view(...): Error, maintain_original must be "
685  "true in order to call this method!" );
686  // Nothing todo
687 }
688 
689 // Overridden from MatrixSymOpGetGMSSymMutable
690 
691 DenseLinAlgPack::DMatrixSliceSym MatrixSymPosDefCholFactor::get_sym_gms_view()
692 {
694  !maintain_original_, std::logic_error
695  ,"MatrixSymPosDefCholFactor::get_sym_gms_view(): Error, maintain_original must be "
696  "true in order to call this method!" );
697  return this->M();
698 }
699 
700 void MatrixSymPosDefCholFactor::commit_sym_gms_view(DenseLinAlgPack::DMatrixSliceSym* sym_gms_view)
701 {
703  !maintain_original_, std::logic_error
704  ,"MatrixSymPosDefCholFactor::commit_sym_gms_view(...): Error, maintain_original must be "
705  "true in order to call this method!" );
706  this->initialize(*sym_gms_view);
707 }
708 
709 // Overridden from MatrixExtractInvCholFactor
710 
711 void MatrixSymPosDefCholFactor::extract_inv_chol( DMatrixSliceTriEle* InvChol ) const
712 {
713  assert_initialized();
714  update_factorization();
715  //
716  // The matrix is represented as the upper cholesky factor:
717  // M = scale * U' * U
718  //
719  // inv(M) = inv(scale*U*U') = (1/sqrt(scale))*inv(U)*(1/sqrt(scale))*inv(U')
720  // = UInv * UInv'
721  // =>
722  // UInv = (1/sqrt(scale))*inv(U)
723  //
724  // Here scale > 0 or an exception will be thrown.
725  //
726  // Compute the inverse cholesky factor as:
727  //
728  // Upper cholesky:
729  // sqrt(scale) * U * UInv = I => InvChol = UInv = (1/sqrt(scale))*inv(U) * I
730  //
731  // Lower cholesky:
732  // sqrt(scale) * L * LInv = I => InvChol = LInv = (1/sqrt(scale))*inv(U) * inv(U') * I
733  //
735  scale_ < 0.0, std::logic_error
736  ,"MatrixSymPosDefCholFactor::extract_inv_chol(...) : "
737  "Error, we can not compute the inverse cholesky factor "
738  "af a negative definite matrix." );
739  DenseLinAlgPack::assign( &InvChol->gms(), 0.0 ); // Set InvChol to identity first.
740  InvChol->gms().diag() = 1.0;
741  DenseLinAlgPack::M_StInvMtM( // Comput InvChol using Level-3 BLAS
742  &InvChol->gms(), 1.0 / std::sqrt(scale_), U()
743  , InvChol->uplo() == BLAS_Cpp::upper ? BLAS_Cpp::no_trans : BLAS_Cpp::trans
744  , InvChol->gms(), BLAS_Cpp::no_trans );
745 }
746 
747 // Overridden from MatrixSymSecantUpdateble
748 
749 void MatrixSymPosDefCholFactor::init_identity( const VectorSpace& space_diag, value_type alpha )
750 {
751  const size_type n = space_diag.dim();
752  allocate_storage( max_size_ ? max_size_ : n );
753  //
754  // B = alpha*I = = alpha*I*I = scale*U'*U
755  // =>
756  // U = sqrt(|alpha|) * I
757  //
758  // Here we will set scale = sign(alpha)
759  //
760  const value_type scale = alpha > 0.0 ? +1.0: -1.0; // Explicitly set the scale
761  resize_and_zero_off_diagonal(n,scale);
762  if( maintain_original_ ) {
763  M().gms().diag()(1,n) = alpha;
764  }
765  if( maintain_factor_ ) {
766  U().gms().diag()(1,n) = std::sqrt(std::fabs(alpha));
767  factor_is_updated_ = true;
768  }
769  is_diagonal_ = true;
770 }
771 
772 void MatrixSymPosDefCholFactor::init_diagonal( const Vector& diag_in )
773 {
774  VectorDenseEncap diag_encap(diag_in);
775  const DVectorSlice diag = diag_encap(); // When diag_encap is destroyed, bye-bye view!
776 
777  allocate_storage( max_size_ ? max_size_ : diag.dim() );
778  //
779  // M = scale * U' * U = scale * (1/scale)*diag^(1/2) * (1/scale)*diag^(1/2)
780  //
781  // Here we will set scale = sign(diag(1)) and validate the rest
782  //
783  if( diag.dim() == 0 ) {
784  M_size_ = 0;
785  return; // We are unsizing this thing
786  }
787  const value_type scale = diag(1) > 0.0 ? +1.0: -1.0;
788  resize_and_zero_off_diagonal(diag.dim(),scale);
789  if( maintain_original_ ) {
790  M().gms().diag() = diag;
791  // ToDo: validate that scale*diag > 0
792  }
793  if( maintain_factor_ ) {
794  DVectorSlice U_diag = U().gms().diag();
795  U_diag = diag;
796  if( scale_ != 1.0 )
797  DenseLinAlgPack::Vt_S( &U_diag, 1.0/scale_ );
798  DenseLinAlgPack::sqrt( &U_diag, U_diag );
799  DenseLinAlgPack::assert_print_nan_inf( U_diag, "(1/scale)*diag", true, NULL );
800  factor_is_updated_ = true;
801  }
802  is_diagonal_ = true;
803 }
804 
805 void MatrixSymPosDefCholFactor::secant_update(
806  VectorMutable* s_in, VectorMutable* y_in, VectorMutable* Bs_in
807  )
808 {
809  using BLAS_Cpp::no_trans;
810  using BLAS_Cpp::trans;
811  using DenseLinAlgPack::dot;
814  namespace rcp = MemMngPack;
815 
816  assert_initialized();
817 
818  // Validate the input
819  TEUCHOS_TEST_FOR_EXCEPT( !( s_in && y_in ) );
820  DenseLinAlgPack::Vp_MtV_assert_sizes( y_in->dim(), M_size_, M_size_, no_trans, s_in->dim() );
821 
822  // Get the serial vectors
823  VectorDenseMutableEncap s_encap(*s_in);
824  VectorDenseMutableEncap y_encap(*y_in);
825  VectorDenseMutableEncap Bs_encap( Bs_in ? *Bs_in : *y_in); // Must pass something on
827  *s = &s_encap(), // When s_encap, y_encap and Bs_encap are destroyed
828  *y = &y_encap(), // these views go bye-bye!
829  *Bs = ( Bs_in ? &Bs_encap() : NULL );
830 
831  // Check skipping the BFGS update
832  const value_type
833  sTy = dot(*s,*y),
834  sTy_scale = sTy*scale_;
835  std::ostringstream omsg;
836  if( !BFGS_sTy_suff_p_d(
837  VectorMutableDense((*s)(),Teuchos::null)
838  ,VectorMutableDense((*y)(),Teuchos::null)
839  ,&sTy_scale,&omsg,"\nMatrixSymPosDefCholFactor::secant_update(...)"
840  )
841  )
842  {
843  throw UpdateSkippedException( omsg.str() );
844  }
845  // Compute Bs if it was not passed in
846  DVector Bs_store;
847  DVectorSlice Bs_view;
848  if( !Bs ) {
849  LinAlgOpPack::V_MtV( &Bs_store, *this, no_trans, *s );
850  Bs_view.bind( Bs_store() );
851  Bs = &Bs_view;
852  }
853  // Check that s'*Bs is positive and if not then throw exception
854  const value_type sTBs = dot(*s,*Bs);
856  scale_*sTBs <= 0.0 && scale_ > 0.0, UpdateFailedException
857  ,"MatrixSymPosDefCholFactor::secant_update(...) : "
858  "Error, B can't be positive definite if s'*Bs <= 0.0" );
860  scale_*sTBs <= 0.0 && scale_ <= 0.0, UpdateFailedException
861  ,"MatrixSymPosDefCholFactor::secant_update(...) : "
862  "Error, B can't be negative definite if s'*Bs >= 0.0" );
863  if( maintain_original_ ) {
864  //
865  // Compute the BFGS update of the original, nonfactored matrix
866  //
867  // Just preform two symmetric updates.
868  //
869  // B = B + (-1/s'*Bs) * Bs*Bs' + (1/s'*y) * y*y'
870  //
871  DMatrixSliceSym M = this->M();
872  DenseLinAlgPack::syr( -1.0/sTBs, *Bs, &M );
873  DenseLinAlgPack::syr( 1.0/sTy, *y, &M );
874  }
875  if( maintain_factor_ ) {
876  //
877  // Compute the BFGS update for the cholesky factor
878  //
879  // If this implementation is based on the one in Section 9.2, page 198-201 of:
880  //
881  // Dennis, J.E., R.B. Schnabel
882  // "Numerical Methods for Unconstrained Optimization"
883  //
884  // Given that we have B = scale*U'*U and the BFGS update:
885  //
886  // B_new = B + y*y'/(y'*s) - (B*s)*(B*s)'/(s'*B*s)
887  //
888  // We can rewrite it in the following form:
889  //
890  // B_new = scale*(U' + a*u*v')*(U + a*v*u') = scale*U_new'*U_new
891  //
892  // where:
893  // v = sqrt(y'*s/(s'*B*s))*U*s
894  // u = y - U'*v
895  // a = 1/(v'*v)
896  //
897  DMatrixSliceTri U = this->U();
898  // v = sqrt(y'*s/(s'*B*s))*U*s
899  DVectorSlice v = *s; // Reuse s as storage for v
900  DenseLinAlgPack::V_MtV( &v, U, no_trans, v ); // Direct call to xSYMV(...)
901  DenseLinAlgPack::Vt_S( &v, std::sqrt( sTy / sTBs ) );
902  // u = (y - U'*v)
903  DVectorSlice u = *y; // Reuse y as storage for u
904  DenseLinAlgPack::Vp_StMtV( &u, -1.0, U, trans, v );
905  // a = 1/(v'*v)
906  const value_type a = 1.0/dot(v,v);
907  // Perform Givens rotations to make Q*(U' + a*u*v') -> U_new upper triangular:
908  //
909  // B_new = scale*(U' + a*u*v')*Q'*Q*(U + a*v*u') = scale*U_new'*U_new
911  a, &v, u, v.dim() > 1 ? &U.gms().diag(-1) : NULL
912  , &DMatrixSliceTriEle(U.gms(),BLAS_Cpp::upper), no_trans );
913  }
914  else {
915  factor_is_updated_ = false;
916  }
917  is_diagonal_ = false;
918 }
919 
920 // Overridden from MatrixSymAddDelUpdateble
921 
922 void MatrixSymPosDefCholFactor::initialize(
923  value_type alpha
924  ,size_type max_size
925  )
926 {
927 #ifdef PROFILE_HACK_ENABLED
928  ProfileHackPack::ProfileTiming profile_timing( "MatrixSymPosDefCholFactor::initialize(alpha,max_size)" );
929 #endif
930  allocate_storage(max_size);
931 
932  if( alpha == 0.0 )
933  throw SingularUpdateException(
934  "MatrixSymPosDefCholFactor::initialize(...): "
935  "Error, alpha == 0.0, matrix is singular.", 0.0 );
936 
937  // Resize the view
938  if( maintain_original_ ) {
939  M_l_r_ = 1;
940  M_l_c_ = 1;
941  }
942  if( U_l_r_ ) {
943  U_l_r_ = 1;
944  U_l_c_ = 1;
945  }
946  M_size_ = 1;
947  max_size_ = my_min( MU_store_.rows(), MU_store_.cols() ) - 1;
948  scale_ = alpha > 0.0 ? +1.0 : -1.0;
949  // Update the matrix
950  if( maintain_original_ ) {
951  M().gms()(1,1) = alpha;
952  }
953  if( U_l_r_ ) {
954  U().gms()(1,1) = std::sqrt( scale_ * alpha );
955  factor_is_updated_ = true;
956  }
957  is_diagonal_ = false;
958 }
959 
960 void MatrixSymPosDefCholFactor::initialize(
961  const DMatrixSliceSym &A
962  ,size_type max_size
963  ,bool force_factorization
964  ,Inertia inertia
965  ,PivotTolerances pivot_tols
966  )
967 {
968 #ifdef PROFILE_HACK_ENABLED
969  ProfileHackPack::ProfileTiming profile_timing( "MatrixSymPosDefCholFactor::initialize(A,max_size...)" );
970 #endif
971  typedef MatrixSymAddDelUpdateable::Inertia Inertia;
972 
973  allocate_storage(max_size);
974 
975  const size_type
976  n = A.rows();
977 
978  // Validate proper usage of inertia parameter
979  TEUCHOS_TEST_FOR_EXCEPT( !( inertia.zero_eigens == Inertia::UNKNOWN || inertia.zero_eigens == 0 ) );
980  TEUCHOS_TEST_FOR_EXCEPT( !( (inertia.neg_eigens == Inertia::UNKNOWN && inertia.pos_eigens == Inertia::UNKNOWN ) )
981  || ( inertia.neg_eigens == n && inertia.pos_eigens == 0 )
982  || ( inertia.neg_eigens == 0 && inertia.pos_eigens == n )
983  );
984 
985  // We can infer if the matrix is p.d. or n.d. by the sign of the diagonal
986  // elements. If a matrix is s.p.d. (s.n.d) then A(i,i) > 0 (A(i,i) < 0)
987  // for all i = 1...n so we will just check the first i.
988  const value_type
989  a_11 = A.gms()(1,1);
990  const int sign_a_11 = ( a_11 == 0.0 ? 0 : ( a_11 > 0.0 ? +1 : -1 ) );
991  if( sign_a_11 == 0.0 )
992  std::logic_error(
993  "MatrixSymPosDefCholFactor::initialize(...) : "
994  "Error, A can not be positive definite or negative definete "
995  "if A(1,1) == 0.0" );
996  if( inertia.pos_eigens == n && sign_a_11 < 0 )
997  std::logic_error(
998  "MatrixSymPosDefCholFactor::initialize(...) : "
999  "Error, A can not be positive definite "
1000  "if A(1,1) < 0.0" );
1001  if( inertia.neg_eigens == n && sign_a_11 > 0 )
1002  std::logic_error(
1003  "MatrixSymPosDefCholFactor::initialize(...) : "
1004  "Error, A can not be negative definite "
1005  "if A(1,1) > 0.0" );
1006  // Now we have got the sign
1007  const value_type
1008  scale = (value_type)sign_a_11;
1009  // Setup the view
1010  set_view(
1011  n,scale,maintain_original_,1,1
1012  ,maintain_factor_, U_l_r_ ? 1 : 0, U_l_r_ ? 1 : 0
1013  );
1014  // Now set the matrix and update the factors
1015  if( maintain_original_ ) {
1016  // Set M = S
1018  &DMatrixSliceTriEle( M().gms(), BLAS_Cpp::lower )
1019  ,DMatrixSliceTriEle( A.gms(), A.uplo() )
1020  );
1021  }
1022  if( maintain_factor_ || force_factorization ) {
1023  // Copy S into U for an inplace factorization.
1025  DenseLinAlgPack::assign( &U_ele, DMatrixSliceTriEle( A.gms(), A.uplo() ) );
1026  if( sign_a_11 < 0 )
1027  DenseLinAlgPack::Mt_S( &U_ele, -1.0 );
1028  try {
1029  DenseLinAlgLAPack::potrf( &U_ele );
1030  factor_is_updated_ = true;
1031  }
1032  catch( const DenseLinAlgLAPack::FactorizationException& excpt ) {
1033  M_size_ = 0; // set unsized
1034  throw SingularUpdateException( excpt.what(), 0.0 );
1035  }
1036  catch(...) {
1037  M_size_ = 0;
1038  throw;
1039  }
1040  // Validate that the tolerances are met and throw appropriate
1041  // exceptions. We already know that the matrix is technically
1042  // p.d. or n.d. Now we must determine gamma = (min(|diag|)/max(|diag|))^2
1043  value_type
1045  max_diag = 0.0;
1046  DVectorSlice::iterator
1047  U_itr = U_ele.gms().diag().begin(),
1048  U_end = U_ele.gms().diag().end();
1049  while( U_itr != U_end ) {
1050  const value_type U_abs = std::abs(*U_itr++);
1051  if(U_abs < min_diag) min_diag = U_abs;
1052  if(U_abs > max_diag) max_diag = U_abs;
1053  }
1054  const value_type gamma = (min_diag*min_diag)/(max_diag*max_diag);
1055  // Validate gamma
1056  PivotTolerances use_pivot_tols = pivot_tols_;
1057  if( pivot_tols.warning_tol != PivotTolerances::UNKNOWN )
1058  use_pivot_tols.warning_tol = pivot_tols.warning_tol;
1059  if( pivot_tols.singular_tol != PivotTolerances::UNKNOWN )
1060  use_pivot_tols.singular_tol = pivot_tols.singular_tol;
1061  if( pivot_tols.wrong_inertia_tol != PivotTolerances::UNKNOWN )
1062  use_pivot_tols.wrong_inertia_tol = pivot_tols.wrong_inertia_tol;
1063  const bool throw_exception = (
1064  gamma == 0.0
1065  || gamma <= use_pivot_tols.warning_tol
1066  || gamma <= use_pivot_tols.singular_tol
1067  );
1068  // Create message and throw exceptions
1069  std::ostringstream onum_msg, omsg;
1070  if(throw_exception) {
1071  onum_msg
1072  << "gamma = (min(|diag(U)(k)|)/|max(|diag(U)(k)|))^2 = (" << min_diag <<"/"
1073  << max_diag << ")^2 = " << gamma;
1074  omsg
1075  << "MatrixSymPosDefCholFactor::initialize(...): ";
1076  if( gamma <= use_pivot_tols.singular_tol ) {
1077  M_size_ = 0; // The initialization failed!
1078  omsg
1079  << "Singular update!\n" << onum_msg.str() << " <= singular_tol = "
1080  << use_pivot_tols.singular_tol;
1081  throw SingularUpdateException( omsg.str(), gamma );
1082  }
1083  else if( gamma <= use_pivot_tols.warning_tol ) {
1084  omsg
1085  << "Warning, near singular update!\nsingular_tol = " << use_pivot_tols.singular_tol
1086  << " < " << onum_msg.str() << " <= warning_tol = "
1087  << use_pivot_tols.warning_tol;
1088  throw WarnNearSingularUpdateException( omsg.str(), gamma ); // The initialization still succeeded through
1089  }
1090  else {
1091  TEUCHOS_TEST_FOR_EXCEPT(true); // Only local programming error?
1092  }
1093  }
1094  }
1095  else {
1096  factor_is_updated_ = false; // The factor is not updated!
1097  }
1098 }
1099 
1100 size_type MatrixSymPosDefCholFactor::max_size() const
1101 {
1102  return max_size_;
1103 }
1104 
1105 MatrixSymAddDelUpdateable::Inertia
1106 MatrixSymPosDefCholFactor::inertia() const
1107 {
1108  typedef MatrixSymAddDelUpdateable MSADU;
1109  typedef MSADU::Inertia Inertia;
1110  return ( M_size_
1111  ? ( scale_ > 0.0
1112  ? Inertia(0,0,M_size_)
1113  : Inertia(M_size_,0,0) )
1114  : Inertia(0,0,0) );
1115 }
1116 
1117 void MatrixSymPosDefCholFactor::set_uninitialized()
1118 {
1119  M_size_ = 0;
1120 }
1121 
1122 void MatrixSymPosDefCholFactor::augment_update(
1123  const DVectorSlice *t
1124  ,value_type alpha
1125  ,bool force_refactorization
1126  ,EEigenValType add_eigen_val
1127  ,PivotTolerances pivot_tols
1128  )
1129 {
1130 #ifdef PROFILE_HACK_ENABLED
1131  ProfileHackPack::ProfileTiming profile_timing( "MatrixSymPosDefCholFactor::augment_udpate(...)" );
1132 #endif
1134  using DenseLinAlgPack::dot;
1136  typedef MatrixSymAddDelUpdateable MSADU;
1137 
1138  assert_initialized();
1139 
1140  // Validate the input
1142  implicit_ref_cast<const MatrixBase>(*this).rows() >= max_size(),
1143  MaxSizeExceededException,
1144  "MatrixSymPosDefCholFactor::augment_update(...) : "
1145  "Error, the maximum size would be exceeded." );
1147  t && t->dim() != M_size_, std::length_error
1148  ,"MatrixSymPosDefCholFactor::augment_update(...): "
1149  "Error, t.dim() must be equal to this->rows()." );
1150  if( !(add_eigen_val == MSADU::EIGEN_VAL_UNKNOWN || add_eigen_val != MSADU::EIGEN_VAL_ZERO ) )
1151  throw SingularUpdateException(
1152  "MatrixSymPosDefCholFactor::augment_update(...): "
1153  "Error, the client has specified a singular update in add_eigen_val.", -1.0 );
1154  if( alpha == 0.0 )
1155  throw SingularUpdateException(
1156  "MatrixSymPosDefCholFactor::augment_update(...): "
1157  "Error, alpha == 0.0 so the matrix is not positive definite "
1158  "or negative definite.", -1.0 );
1159  if( scale_ > 0.0 && alpha < 0.0 )
1160  throw WrongInertiaUpdateException(
1161  "MatrixSymPosDefCholFactor::augment_update(...): "
1162  "Error, alpha < 0.0 so the matrix is not postivie definite ", -1.0 );
1163  if( scale_ < 0.0 && alpha > 0.0 )
1164  throw WrongInertiaUpdateException(
1165  "MatrixSymPosDefCholFactor::augment_update(...): "
1166  "Error, alpha > 0.0 so the matrix is not negative definite ", -1.0 );
1167  if( scale_ > 0.0 && !( add_eigen_val == MSADU::EIGEN_VAL_POS || add_eigen_val == MSADU::EIGEN_VAL_UNKNOWN ) )
1168  throw WrongInertiaUpdateException(
1169  "MatrixSymPosDefCholFactor::augment_update(...): "
1170  "Error, alpha > 0.0 but the user specified a non-positive definite new matrix.", -1.0 );
1171  if( scale_ < 0.0 && !( add_eigen_val == MSADU::EIGEN_VAL_NEG || add_eigen_val == MSADU::EIGEN_VAL_UNKNOWN ) )
1172  throw WrongInertiaUpdateException(
1173  "MatrixSymPosDefCholFactor::augment_update(...): "
1174  "Error, alpha > 0.0 but the user specified a non-positive definite new matrix.", -1.0 );
1175  // First try to augment the factor to verify that the matrix is still p.d. or n.d.
1176  bool throw_exception = false; // If true then throw exception
1177  std::ostringstream omsg; // Will be set if an exception has to be thrown.
1178  value_type gamma; // ...
1179  if( maintain_factor_ ) {
1180  //
1181  // The update is:
1182  //
1183  // B_new = [ B t ]
1184  // [ t' alpha ]
1185  //
1186  // = scale * [ U'*U (1/scale)*t ]
1187  // [ (1/scale)*t' (1/scale)*alpha ]
1188  //
1189  // We seek to find a new cholesky factor of the form:
1190  //
1191  // U_new = [ U11 u12 ]
1192  // [ 0 u22 ]
1193  //
1194  // B_new = scale*U_new'*U_new
1195  //
1196  // = scale * [ U11' 0 ] * [ U11 u12 ]
1197  // [ u12' u22 ] [ 0 u22 ]
1198  //
1199  // = scale * [ U11'*U11 U11'*u12 ]
1200  // [ u12'*U11 u12'*u12 + u22^2 ]
1201  //
1202  // From the above we see that:
1203  // =>
1204  // U11 = U
1205  // u12 = inv(U') * (1/scale) * t
1206  // u22 = sqrt( (1/scale)*alpha - u12'*u12 );
1207  //
1208  // We must compute gamma relative to the LU factorization
1209  // so we must square ||U11.diag()||inf.
1210 
1211  // Get references to the storage for the to-be-updated parts for the new factor.
1212  DVectorSlice u12 = MU_store_.col(U_l_c_+M_size_+1)(U_l_r_,U_l_r_+M_size_-1);
1213  value_type &u22 = MU_store_(U_l_r_+M_size_,U_l_c_+M_size_+1);
1214  // u12 = inv(U') * (1/scale) * t
1215  if(t) {
1216  DenseLinAlgPack::V_InvMtV( &u12, U(), BLAS_Cpp::trans, *t );
1217  if( scale_ != 1.0 ) DenseLinAlgPack::Vt_S( &u12, 1.0/scale_ );
1218  }
1219  else {
1220  u12 = 0.0;
1221  }
1222  // u22^2 = (1/scale)*alpha - u12'*u12;
1223  const value_type
1224  u22sqr = (1/scale_) * alpha - ( t ? dot( u12, u12 ) : 0.0 ),
1225  u22sqrabs = std::abs(u22sqr),
1226  nrm_U_diag = norm_inf(U().gms().diag()),
1227  sqr_nrm_U_diag = nrm_U_diag * nrm_U_diag;
1228  // Calculate gamma in proper context
1229  gamma = u22sqrabs / sqr_nrm_U_diag;
1230  // Check gamma
1231  const bool
1232  correct_inertia = u22sqr > 0.0;
1233  PivotTolerances use_pivot_tols = pivot_tols_;
1234  if( pivot_tols.warning_tol != PivotTolerances::UNKNOWN )
1235  use_pivot_tols.warning_tol = pivot_tols.warning_tol;
1236  if( pivot_tols.singular_tol != PivotTolerances::UNKNOWN )
1237  use_pivot_tols.singular_tol = pivot_tols.singular_tol;
1238  if( pivot_tols.wrong_inertia_tol != PivotTolerances::UNKNOWN )
1239  use_pivot_tols.wrong_inertia_tol = pivot_tols.wrong_inertia_tol;
1240  throw_exception = (
1241  gamma == 0.0
1242  || correct_inertia && gamma <= use_pivot_tols.warning_tol
1243  || correct_inertia && gamma <= use_pivot_tols.singular_tol
1244  || !correct_inertia
1245  );
1246  // Create message and throw exceptions
1247  std::ostringstream onum_msg;
1248  if(throw_exception) {
1249  onum_msg
1250  << "gamma = |u22^2|/(||diag(U11)||inf)^2 = |" << u22sqr <<"|/("
1251  << nrm_U_diag << ")^2 = " << gamma;
1252  omsg
1253  << "MatrixSymPosDefCholFactor::augment_update(...): ";
1254  if( correct_inertia && gamma <= use_pivot_tols.singular_tol ) {
1255  omsg
1256  << "Singular update!\n" << onum_msg.str() << " <= singular_tol = "
1257  << use_pivot_tols.singular_tol;
1258  throw SingularUpdateException( omsg.str(), gamma );
1259  }
1260  else if( !correct_inertia && gamma <= use_pivot_tols.singular_tol ) {
1261  omsg
1262  << "Singular update!\n" << onum_msg.str() << " <= wrong_inertia_tol = "
1263  << use_pivot_tols.wrong_inertia_tol;
1264  throw SingularUpdateException( omsg.str(), gamma );
1265  }
1266 
1267  else if( !correct_inertia ) {
1268  omsg
1269  << "Indefinite update!\n" << onum_msg.str() << " >= wrong_inertia_tol = "
1270  << use_pivot_tols.wrong_inertia_tol;
1271  throw WrongInertiaUpdateException( omsg.str(), gamma );
1272  }
1273  else if( correct_inertia && gamma <= use_pivot_tols.warning_tol ) {
1274  omsg
1275  << "Warning, near singular update!\nsingular_tol = " << use_pivot_tols.singular_tol
1276  << " < " << onum_msg.str() << " <= warning_tol = "
1277  << use_pivot_tols.warning_tol;
1278  // Don't throw the exception till the very end!
1279  }
1280  else {
1281  TEUCHOS_TEST_FOR_EXCEPT(true); // Only local programming error?
1282  }
1283  }
1284  // u22 = sqrt(u22^2)
1285  u22 = std::sqrt(u22sqrabs);
1286  }
1287  else {
1288  factor_is_updated_ = false;
1289  }
1290  // Now augment the original
1291  if( maintain_original_ ) {
1292  //
1293  // M_new = [ M t ]
1294  // [ t' alpha ]
1295  //
1296  DVectorSlice M12 = MU_store_.row(M_l_r_+M_size_+1)(M_l_c_,M_l_c_+M_size_-1);
1297  if(t)
1298  M12 = *t;
1299  else
1300  M12 = 0.0;
1301  MU_store_(M_l_r_+M_size_+1,M_l_c_+M_size_) = alpha;
1302  }
1303  ++M_size_; // Enlarge the matrix by one
1304  is_diagonal_ = false;
1305  if( throw_exception )
1306  throw WarnNearSingularUpdateException(omsg.str(),gamma);
1307 }
1308 
1309 void MatrixSymPosDefCholFactor::delete_update(
1310  size_type jd
1311  ,bool force_refactorization
1312  ,EEigenValType drop_eigen_val
1313  ,PivotTolerances pivot_tols
1314  )
1315 {
1316 #ifdef PROFILE_HACK_ENABLED
1317  ProfileHackPack::ProfileTiming profile_timing( "MatrixSymPosDefCholFactor::delete_udpate(...)" );
1318 #endif
1319  typedef MatrixSymAddDelUpdateable MSADU;
1320 
1322  jd < 1 || M_size_ < jd, std::out_of_range
1323  ,"MatrixSymPosDefCholFactor::delete_update(jd,...): "
1324  "Error, the indice jd must be 1 <= jd <= rows()" );
1325 
1326  TEUCHOS_TEST_FOR_EXCEPT( !( drop_eigen_val == MSADU::EIGEN_VAL_UNKNOWN
1327  || (scale_ > 0.0 && drop_eigen_val == MSADU::EIGEN_VAL_POS)
1328  || (scale_ < 0.0 && drop_eigen_val == MSADU::EIGEN_VAL_NEG)
1329  ) );
1330 
1331  if( maintain_original_ ) {
1332  //
1333  // Here we have the lower portion of M partitioned as:
1334  //
1335  // 1 |\
1336  // | \
1337  // | \
1338  // | M11 \
1339  // |________\ _
1340  // jd |_________|_|
1341  // | | |\
1342  // | | | \
1343  // | | | \
1344  // | M31 | | M33 \
1345  // n | | | \
1346  // ----------------------
1347  // 1 jd n
1348  //
1349  // We want to move M31 up one row and M33 up and to the left.
1350  //
1352  }
1353  if( maintain_factor_ ) {
1354  //
1355  // Here we have U partitioned as:
1356  //
1357  // 1 jd n
1358  // ------------------------- 1
1359  // \ | | |
1360  // \ U11 | | U13 |
1361  // \ | | |
1362  // \ u12->| | u23' |
1363  // \ | | | |
1364  // \ | | \./ |
1365  // \ |_|_______|
1366  // u22 -> |_|_______| jd
1367  // \ |
1368  // \ U33 |
1369  // \ |
1370  // \ | n
1371  //
1372  // To preform the update we need to first update U33 as:
1373  // U33'*U33 + u23'*u23 ==> U33'*U33
1374  // Then we need to move U12 (one column at a time) to the
1375  // left one column and overwrite u12.
1376  // Then we will move the updated U33 (one column at a time)
1377  // up and to the left one position to overwrite u22 and
1378  // the left part of u23'. We then decrease the dimension
1379  // of U by one and we are finished updating the factorization.
1380  //
1381  // See RAB notes from 1/21/99 and 1/26/99 for details on this update.
1382  //
1383  const size_type n = M_size_;
1384  // Resize workspace if it has not been done yet.
1385  size_type work_size = 3 * max_size_;
1386  if(work_.dim() < work_size)
1387  work_.resize(work_size);
1388  // Update the factors
1389  {
1390  DMatrixSlice U = this->U().gms();
1391  // Update U33 where it sits.
1392  if(jd < n) {
1393  size_type _size = n-jd; // Set storage for u23, c and s
1394  Range1D rng(1,_size);
1395  DVectorSlice
1396  u23 = work_(rng),
1397  c = work_(rng+_size),
1398  s = work_(rng+2*_size);
1399  Range1D U_rng(jd+1,n); // Set U33 and u23
1400  DMatrixSlice U33 = U(U_rng,U_rng);
1401  u23 = U.row(jd)(U_rng);
1402  // Update U33
1403  value_type dummy;
1404  ALAP_DCHUD_CALL (
1405  U33.col_ptr(1), U33.max_rows()
1406  ,U_rng.size(), u23.start_ptr(), &dummy, 1, 0, &dummy
1407  ,&dummy, c.start_ptr(), s.start_ptr() );
1408  }
1409  // Move U13 and U33 to delete row and column jd
1411  }
1412  }
1413  else {
1414  factor_is_updated_ = false;
1415  }
1416  // Strink the size of M and U
1417  --M_size_;
1418 }
1419 
1420 // Overridden from Serializable
1421 
1422 // ToDo: Refactor this code and create an external utility matrix
1423 // serialization class that will convert from matrix type to matrix
1424 // type.
1425 
1426 void MatrixSymPosDefCholFactor::serialize( std::ostream &out ) const
1427 {
1428  // Write key words on top line
1429  out << build_serialization_string() << std::endl;
1430  // Set the precision (very important!)
1431  out.precision(std::numeric_limits<value_type>::digits10+4);
1432  // Write the dimmension
1433  out << M_size_ << std::endl;
1434  if(M_size_) {
1435  // Write the matrix values
1436  if( maintain_original_ ) {
1437  const DMatrixSliceSym M = this->M();
1438  write_matrix( M.gms(), M.uplo(), out );
1439  }
1440  else {
1441  const DMatrixSliceTri U = this->U();
1442  write_matrix( U.gms(), U.uplo(), out );
1443  }
1444  }
1445  // ToDo: You need to write both M and U if both are computed!
1446 }
1447 
1448 void MatrixSymPosDefCholFactor::unserialize( std::istream &in )
1449 {
1450  // Get the keywords for the matrix type
1451  std::string keywords;
1452  std::getline( in, keywords, '\n' );
1453  // For now make sure the types are exactly the same!
1454  const std::string this_keywords = build_serialization_string();
1456  this_keywords != keywords, std::logic_error
1457  ,"MatrixSymPosDefCholFactor::unserialize(...): Error, the matrix type being read in from file of "
1458  "\'"<<keywords<<"\' does not equal the type expected of \'"<<this_keywords<<"\'!"
1459  );
1460  // Read in the dimension of the matrix
1461  in >> M_size_;
1463  M_size_ < 0, std::logic_error
1464  ,"MatrixSymPosDefCholFactor::unserialize(...): Error, read in a size of M_size = "<<M_size_<<" < 0!"
1465  );
1466  allocate_storage(M_size_);
1467  // Read in the matrix into storage
1468  if(maintain_original_) {
1469  M_l_r_ = M_l_c_ = 1;
1470  DMatrixSliceSym M = this->M();
1471  read_matrix( in, M.uplo(), &M.gms() );
1472  }
1473  else {
1474  U_l_r_ = U_l_c_ = 1;
1475  DMatrixSliceTri U = this->U();
1476  read_matrix( in, U.uplo(), &U.gms() );
1477  }
1478 }
1479 
1480 // Private member functions
1481 
1482 void MatrixSymPosDefCholFactor::assert_storage() const
1483 {
1484  TEUCHOS_TEST_FOR_EXCEPT( !( MU_store_.rows() ) );
1485 }
1486 
1487 void MatrixSymPosDefCholFactor::allocate_storage(size_type max_size) const
1488 {
1489  namespace rcp = MemMngPack;
1490  if( allocates_storage_ && MU_store_.rows() < max_size + 1 ) {
1491  // We have the right to allocate storage so lets just do it.
1493  MU_store = Teuchos::rcp(new DMatrix( max_size + 1, max_size + 1 ));
1495  const_cast<MatrixSymPosDefCholFactor*>(this)->release_resource_ptr_ = Teuchos::rcp(new ptr_t(MU_store));
1496  const_cast<MatrixSymPosDefCholFactor*>(this)->MU_store_.bind( (*MU_store)() );
1497  const_cast<MatrixSymPosDefCholFactor*>(this)->max_size_ = max_size;
1498  }
1499  else {
1500  TEUCHOS_TEST_FOR_EXCEPT( !( MU_store_.rows() >= max_size + 1 ) );
1501  }
1502 }
1503 
1504 void MatrixSymPosDefCholFactor::assert_initialized() const
1505 {
1506  TEUCHOS_TEST_FOR_EXCEPT( !( M_size_ ) );
1507 }
1508 
1509 void MatrixSymPosDefCholFactor::resize_and_zero_off_diagonal(size_type n, value_type scale)
1510 {
1513 
1514  TEUCHOS_TEST_FOR_EXCEPT( !( n <= my_min( MU_store_.rows(), MU_store_.cols() ) - 1 ) );
1515 
1516  // Resize the views
1517  set_view(
1518  n, scale, maintain_original_, 1, 1, maintain_factor_
1519  ,U_l_r_ ? 1 : 0, U_l_r_ ? 1 : 0 );
1520  if( maintain_original_ ) {
1521  if(!is_diagonal_ && n > 1)
1522  assign( &nonconst_tri_ele( M().gms()(2,n,1,n-1), BLAS_Cpp::lower ), 0.0 );
1523  }
1524  if( U_l_r_ ) {
1525  if(!is_diagonal_ && n > 1)
1526  assign( &nonconst_tri_ele( U().gms()(1,n-1,2,n), BLAS_Cpp::upper ), 0.0 );
1527  }
1528 }
1529 
1530 void MatrixSymPosDefCholFactor::update_factorization() const
1531 {
1532  if( factor_is_updated_ ) return; // The factor should already be updated.
1534  U_l_r_ == 0, std::logic_error
1535  ,"MatrixSymPosDefCholFactor::update_factorization() : "
1536  "Error, U_l_r == 0 was set in MatrixSymPosDefCholFactor::set_view(...) "
1537  "and therefore we can not update the factorization in the provided storage space." );
1538  // Update the factorization from scratch!
1539 #ifdef PROFILE_HACK_ENABLED
1540  ProfileHackPack::ProfileTiming profile_timing( "MatrixSymPosDefCholFactor::update_factorization(...) ... update" );
1541 #endif
1542  MatrixSymPosDefCholFactor
1543  *nc_this = const_cast<MatrixSymPosDefCholFactor*>(this);
1545  DenseLinAlgPack::assign( &U, DenseLinAlgPack::tri_ele( M().gms(), BLAS_Cpp::lower ) ); // Copy in the original
1546  {
1547 #ifdef PROFILE_HACK_ENABLED
1548  ProfileHackPack::ProfileTiming profile_timing( "MatrixSymPosDefCholFactor::update_factorization(...) ... potrf" );
1549 #endif
1551  }
1552  nc_this->factor_is_updated_ = true;
1553 }
1554 
1555 std::string MatrixSymPosDefCholFactor::build_serialization_string() const
1556 {
1557  std::string str = "SYMMETRIC POS_DEF";
1558  if( !maintain_original_ )
1559  str.append(" CHOL_FACTOR");
1560  if( maintain_original_ )
1561  str.append(" LOWER");
1562  else
1563  str.append(" UPPER");
1564  return str;
1565 }
1566 
1567 void MatrixSymPosDefCholFactor::write_matrix( const DMatrixSlice &Q, BLAS_Cpp::Uplo Q_uplo, std::ostream &out )
1568 {
1569  const int Q_dim = Q.rows();
1570  if( Q_uplo == BLAS_Cpp::lower ) {
1571  for( int i = 1; i <= Q_dim; ++i ) {
1572  for( int j = 1; j <= i; ++j ) {
1573  out << " " << Q(i,j);
1574  }
1575  out << std::endl;
1576  }
1577  }
1578  else {
1579  for( int i = 1; i <= Q_dim; ++i ) {
1580  for( int j = i; j <= Q_dim; ++j ) {
1581  out << " " << Q(i,j);
1582  }
1583  out << std::endl;
1584  }
1585  }
1586 }
1587 
1588 void MatrixSymPosDefCholFactor::read_matrix( std::istream &in, BLAS_Cpp::Uplo Q_uplo, DMatrixSlice *Q_out )
1589 {
1590  DMatrixSlice &Q = *Q_out;
1591  const int Q_dim = Q.rows();
1592  if( Q_uplo == BLAS_Cpp::lower ) {
1593  for( int i = 1; i <= Q_dim; ++i ) {
1594  for( int j = 1; j <= i; ++j ) {
1595 #ifdef TEUCHOS_DEBUG
1597  in.eof(), std::logic_error
1598  ,"MatrixSymPosDefCholFactor::read_matrix(in,lower,Q_out): Error, not finished reading in matrix yet (i="<<i<<",j="<<j<<")!"
1599  );
1600 #endif
1601  in >> Q(i,j);
1602  }
1603  }
1604  }
1605  else {
1606  for( int i = 1; i <= Q_dim; ++i ) {
1607  for( int j = i; j <= Q_dim; ++j ) {
1608 #ifdef TEUCHOS_DEBUG
1610  in.eof(), std::logic_error
1611  ,"MatrixSymPosDefCholFactor::read_matrix(in,upper,Q_out): Error, not finished reading in matrix yet (i="<<i<<",j="<<j<<")!"
1612  );
1613 #endif
1614  in >> Q(i,j);
1615  }
1616  }
1617  }
1618 }
1619 
1620 } // end namespace AbstractLinAlgPack
f_dbl_prec const f_int f_dbl_prec * Y
void V_MtV(DVector *v_lhs, const DMatrixSliceTri &tri_rhs1, BLAS_Cpp::Transp trans_rhs1, const DVectorSlice &vs_rhs2)
AbstractLinAlgPack::size_type size_type
void Vp_StV(DVectorSlice *vs_lhs, value_type alpha, const DVectorSlice &vs_rhs)
vs_lhs += alpha * vs_rhs (BLAS xAXPY)
FortranTypes::f_int f_int
void assign(DMatrix *gm_lhs, value_type alpha)
gm_lhs = alpha (elementwise)
void MtM_assert_sizes(size_type m_rhs1_rows, size_type m_rhs1_cols, BLAS_Cpp::Transp trans_rhs1, size_type m_rhs2_rows, size_type m_rhs2_cols, BLAS_Cpp::Transp trans_rhs2)
op(m_lhs) += op(m_rhs1)
const f_int const f_int const f_dbl_prec const f_dbl_prec * X
void Vp_StPtMtV(VectorMutable *v_lhs, value_type alpha, const GenPermMatrixSlice &P_rhs1, BLAS_Cpp::Transp P_rhs1_trans, const MatrixOp &M_rhs2, BLAS_Cpp::Transp M_rhs2_trans, const Vector &v_rhs3, value_type beta=1.0)
v_lhs = alpha * op(P_rhs1) * op(M_rhs2) * v_rhs3 + beta * v_rhs
const DMatrixSliceTriEle tri_ele(const DMatrixSlice &gms, BLAS_Cpp::Uplo uplo)
#define TEUCHOS_TEST_FOR_EXCEPTION(throw_exception_test, Exception, msg)
void Vp_StV(VectorMutable *v_lhs, const value_type &alpha, const Vector &v_rhs)
v_lhs = alpha * v_rhs + v_lhs
void abs(DVectorSlice *vs_lhs, const DVectorSlice &vs_rhs)
vs_lhs = abs(vs_rhs)
size_type rows(size_type rows, size_type cols, BLAS_Cpp::Transp _trans)
Return rows of a possible transposed matrix.
void assign(VectorMutable *v_lhs, const V &V_rhs)
v_lhs = V_rhs.
std::ostream & output(std::ostream &o, const COOMatrix &coom)
Output stream function for COOMatrix.
void Mp_MtM_assert_sizes(size_type m_lhs_rows, size_type m_lhs_cols, BLAS_Cpp::Transp trans_lhs, size_type m_rhs1_rows, size_type m_rhs1_cols, BLAS_Cpp::Transp trans_rhs1, size_type m_rhs2_rows, size_type m_rhs2_cols, BLAS_Cpp::Transp trans_rhs2)
op(m_lhs) += op(m_rhs1) * op(m_rhs2)
Helper class that takes care of timing.
void M_MtM(MatrixOp *M_lhs, const MatrixOp &M_rhs1, BLAS_Cpp::Transp trans_rhs1, const MatrixOp &M_rhs2, BLAS_Cpp::Transp trans_rhs2)
M_lhs = op(M_rhs1) * op(M_rhs2).
bool assert_print_nan_inf(const value_type &val, const std::string &name, bool throw_excpt, std::ostream *out)
TypeTo & implicit_ref_cast(TypeFrom &t)
void syrk(BLAS_Cpp::Transp trans, value_type alpha, const DMatrixSlice &gms_rhs, value_type beta, DMatrixSliceSym *sym_lhs)
Transposed.
const DMatrixSliceTri tri(const DMatrixSlice &gms, BLAS_Cpp::Uplo uplo, BLAS_Cpp::Diag diag)
void V_InvMtV(DVector *v_lhs, const DMatrixSliceTri &tri_rhs1, BLAS_Cpp::Transp trans_rhs1, const DVectorSlice &vs_rhs2)
void Mp_StPtMtP(MatrixOp *mwo_lhs, value_type alpha, const GenPermMatrixSlice &P_rhs1, BLAS_Cpp::Transp P_rhs1_trans, const MatrixOp &M_rhs, BLAS_Cpp::Transp trans_rhs, const GenPermMatrixSlice &P_rhs2, BLAS_Cpp::Transp P_rhs2_trans)
mwo_lhs += alpha * op(P_rhs1) * op(M_rhs) * op(P_rhs2).
Not transposed.
value_type norm_2(const DVectorSlice &vs_rhs)
result = ||vs_rhs||2 (BLAS xNRM2)
TEUCHOS_DEPRECATED RCP< T > rcp(T *p, Dealloc_T dealloc, bool owns_mem)
void syr(value_type alpha, const DVectorSlice &vs_rhs, DMatrixSliceSym *sym_lhs)
void syrk(const MatrixOp &mwo_rhs, BLAS_Cpp::Transp M_trans, value_type alpha, value_type beta, MatrixSymOp *sym_lhs)
Perform a rank-k update of a symmetric matrix of the form:
void Mt_S(DMatrixSlice *gms_lhs, value_type alpha)
gms_lhs *= alpha (BLAS xSCAL)
std::ostream * out
DenseLinAlgPack::DMatrixSliceTriEle DMatrixSliceTriEle
void rank_2_chol_update(const value_type a, DVectorSlice *u, const DVectorSlice &v, DVectorSlice *w, DMatrixSliceTriEle *R, BLAS_Cpp::Transp R_trans)
Update the cholesky factor of symmetric positive definite matrix for a rank-2 update.
const LAPACK_C_Decl::f_int & M
void V_InvMtV(VectorMutable *v_lhs, const MatrixNonsing &M_rhs1, BLAS_Cpp::Transp trans_rhs1, const Vector &v_rhs2)
v_lhs = inv(op(M_rhs1)) * v_rhs2
f_dbl_prec f_dbl_prec f_dbl_prec * S
void Mp_StM(MatrixOp *mwo_lhs, value_type alpha, const MatrixOp &M_rhs, BLAS_Cpp::Transp trans_rhs)
void potrf(DMatrixSliceTriEle *A)
Calls xPOTRF to compute the cholesky factorization of a symmetric positive definte matrix...
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)
v_lhs = alpha * op(M_rhs1) * v_rhs2 + beta * v_lhs (BLAS xGEMV)
DVectorSlice col(DMatrixSlice &gms, BLAS_Cpp::Transp trans, size_type j)
value_type norm_inf(const SparseVectorSlice< T_Ele > &sv_rhs)
result = ||sv_rhs||inf (BLAS IxAMAX)
bool BFGS_sTy_suff_p_d(const Vector &s, const Vector &y, const value_type *sTy=NULL, std::ostream *out=NULL, const char func_name[]=NULL)
Check that s'*y is sufficiently positive and print the result if it is not.
DMatrixSliceTriEle nonconst_tri_ele(DMatrixSlice gms, BLAS_Cpp::Uplo uplo)
Return a triangular element-wise matrix.
Template class that implements ReleaseResource interface for a RCP<T> object.
value_type dot(const Vector &v_rhs1, const Vector &v_rhs2)
result = v_rhs1' * v_rhs2
void sqrt(DVectorSlice *vs_lhs, const DVectorSlice &vs_rhs)
vs_lhs = sqrt(vs_rhs)
SparseVectorSlice< SparseElement< index_type, value_type > > SpVectorSlice
Transp trans_not(Transp _trans)
Return the opposite of the transpose argument.
DenseLinAlgPack::VectorSliceTmpl< value_type > DVectorSlice
size_type rows() const
Return the number of rows.
void Vt_S(DVectorSlice *vs_lhs, value_type alpha)
vs_lhs *= alpha (BLAS xSCAL) (*** Note that alpha == 0.0 is handeled as vs_lhs = 0.0)
int read_matrix(const std::string &filename, Epetra_Comm &Comm, const Epetra_Map *rowmap, Epetra_Map *colmap, const Epetra_Map *rangemap, const Epetra_Map *domainmap, Epetra_CrsMatrix *&mat)
void V_MtV(VectorMutable *v_lhs, const MatrixOp &M_rhs1, BLAS_Cpp::Transp trans_rhs1, const V &V_rhs2)
v_lhs = op(M_rhs1) * V_rhs2.
AbstractLinAlgPack::value_type value_type
DenseLinAlgPack::DMatrixSliceTri DMatrixSliceTri
void Mp_StM(DMatrixSliceTriEle *tri_lhs, value_type alpha, const DMatrixSliceTriEle &tri_rhs)
tri_lhs += alpha * tri_rhs (BLAS xAXPY)
void Vp_MtV_assert_sizes(size_type v_lhs_size, size_type m_rhs1_rows, size_type m_rhs1_cols, BLAS_Cpp::Transp trans_rhs1, size_type v_rhs2_size)
v_lhs += op(m_rhs1) * v_rhs2
void Vp_StMtV(DVectorSlice *vs_lhs, value_type alpha, const DMatrixSlice &gms_rhs1, BLAS_Cpp::Transp trans_rhs1, const DVectorSlice &vs_rhs2, value_type beta=1.0)
vs_lhs = alpha * op(gms_rhs1) * vs_rhs2 + beta * vs_lhs (BLAS xGEMV)
void M_StInvMtM(MatrixOp *m_lhs, value_type alpha, const MatrixNonsing &M_rhs1, BLAS_Cpp::Transp trans_rhs1, const MatrixOp &mwo_rhs2, BLAS_Cpp::Transp trans_rhs2)
m_lhs = alpha * inv(op(mwo_rhs1)) * op(mwo_rhs2) (right)
value_type norm_inf(const DVectorSlice &vs_rhs)
result = ||vs_rhs||infinity (BLAS IxAMAX)
void delete_row_col(size_type kd, DMatrixSliceTriEle *M)
DenseLinAlgPack::DMatrixSlice DMatrixSlice
void M_StMtInvMtM(MatrixSymOp *sym_gms_lhs, value_type alpha, const MatrixOp &mwo, BLAS_Cpp::Transp mwo_trans, const MatrixSymNonsing &mswof, MatrixSymNonsing::EMatrixDummyArg mwo_rhs)
sym_gms_lhs = alpha * op(mwo) * inv(mswof) * op(mwo)'
DenseLinAlgPack::DMatrixSliceSym DMatrixSliceSym
Transp
TRANS.
FortranTypes::f_dbl_prec f_dbl_prec
DenseLinAlgPack::DMatrix DMatrix
int n
size_type cols(size_type rows, size_type cols, BLAS_Cpp::Transp _trans)
Return columns of a possible transposed matrix.
const f_int f_dbl_prec const f_int f_int const f_int f_int const f_dbl_prec & u
#define TEUCHOS_TEST_FOR_EXCEPT(throw_exception_test)
void M_StInvMtM(DMatrix *gm_lhs, value_type alpha, const DMatrixSliceTri &tri_rhs1, BLAS_Cpp::Transp trans_rhs1, const DMatrixSlice &gms_rhs2, BLAS_Cpp::Transp trans_rhs2)
value_type dot(const DVectorSlice &vs_rhs1, const DVectorSlice &vs_rhs2)
result = vs_rhs1' * vs_rhs2 (BLAS xDOT)