Teuchos - Trilinos Tools Package  Version of the Day
Teuchos_SerialSpdDenseSolver.hpp
Go to the documentation of this file.
1 
2 // @HEADER
3 // ***********************************************************************
4 //
5 // Teuchos: Common Tools Package
6 // Copyright (2004) Sandia Corporation
7 //
8 // Under terms of Contract DE-AC04-94AL85000, there is a non-exclusive
9 // license for use of this work by or on behalf of the U.S. Government.
10 //
11 // Redistribution and use in source and binary forms, with or without
12 // modification, are permitted provided that the following conditions are
13 // met:
14 //
15 // 1. Redistributions of source code must retain the above copyright
16 // notice, this list of conditions and the following disclaimer.
17 //
18 // 2. Redistributions in binary form must reproduce the above copyright
19 // notice, this list of conditions and the following disclaimer in the
20 // documentation and/or other materials provided with the distribution.
21 //
22 // 3. Neither the name of the Corporation nor the names of the
23 // contributors may be used to endorse or promote products derived from
24 // this software without specific prior written permission.
25 //
26 // THIS SOFTWARE IS PROVIDED BY SANDIA CORPORATION "AS IS" AND ANY
27 // EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
28 // IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR
29 // PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL SANDIA CORPORATION OR THE
30 // CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL,
31 // EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO,
32 // PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR
33 // PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF
34 // LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING
35 // NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
36 // SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
37 //
38 // Questions? Contact Michael A. Heroux (maherou@sandia.gov)
39 //
40 // ***********************************************************************
41 // @HEADER
42 
43 #ifndef TEUCHOS_SERIALSPDDENSESOLVER_H
44 #define TEUCHOS_SERIALSPDDENSESOLVER_H
45 
50 #include "Teuchos_CompObject.hpp"
51 #include "Teuchos_BLAS.hpp"
52 #include "Teuchos_LAPACK.hpp"
53 #include "Teuchos_RCP.hpp"
54 #include "Teuchos_ConfigDefs.hpp"
58 
59 #ifdef HAVE_TEUCHOSNUMERICS_EIGEN
60 #include "Eigen/Dense"
61 #endif
62 
130 namespace Teuchos {
131 
132  template<typename OrdinalType, typename ScalarType>
133  class SerialSpdDenseSolver : public CompObject, public Object, public BLAS<OrdinalType, ScalarType>,
134  public LAPACK<OrdinalType, ScalarType>
135  {
136  public:
137 
138  typedef typename Teuchos::ScalarTraits<ScalarType>::magnitudeType MagnitudeType;
139 #ifdef HAVE_TEUCHOSNUMERICS_EIGEN
140  typedef typename Eigen::Matrix<ScalarType,Eigen::Dynamic,Eigen::Dynamic,Eigen::ColMajor> EigenMatrix;
141  typedef typename Eigen::Matrix<ScalarType,Eigen::Dynamic,1> EigenVector;
142  typedef typename Eigen::InnerStride<Eigen::Dynamic> EigenInnerStride;
143  typedef typename Eigen::OuterStride<Eigen::Dynamic> EigenOuterStride;
144  typedef typename Eigen::Map<EigenVector,0,EigenInnerStride > EigenVectorMap;
145  typedef typename Eigen::Map<const EigenVector,0,EigenInnerStride > EigenConstVectorMap;
146  typedef typename Eigen::Map<EigenMatrix,0,EigenOuterStride > EigenMatrixMap;
147  typedef typename Eigen::Map<const EigenMatrix,0,EigenOuterStride > EigenConstMatrixMap;
148  typedef typename Eigen::PermutationMatrix<Eigen::Dynamic,Eigen::Dynamic,OrdinalType> EigenPermutationMatrix;
149  typedef typename Eigen::Array<OrdinalType,Eigen::Dynamic,1> EigenOrdinalArray;
150  typedef typename Eigen::Map<EigenOrdinalArray> EigenOrdinalArrayMap;
151  typedef typename Eigen::Array<ScalarType,Eigen::Dynamic,1> EigenScalarArray;
152  typedef typename Eigen::Map<EigenScalarArray> EigenScalarArrayMap;
153 #endif
154 
156 
159 
160 
162  virtual ~SerialSpdDenseSolver();
164 
166 
167 
170 
172 
178 
180 
181 
183 
185  void factorWithEquilibration(bool flag) {equilibrate_ = flag; return;}
186 
188  void solveToRefinedSolution(bool flag) {refineSolution_ = flag; return;}
189 
191 
194  void estimateSolutionErrors(bool flag);
196 
198 
199 
201 
204  int factor();
205 
207 
210  int solve();
211 
213 
218  int invert();
219 
221 
225 
227 
230  int equilibrateMatrix();
231 
233 
236  int equilibrateRHS();
237 
238 
240 
243  int applyRefinement();
244 
246 
249  int unequilibrateLHS();
250 
252 
258  int reciprocalConditionEstimate(MagnitudeType & Value);
260 
262 
263 
265  bool transpose() {return(transpose_);}
266 
268  bool factored() {return(factored_);}
269 
271  bool equilibratedA() {return(equilibratedA_);}
272 
274  bool equilibratedB() {return(equilibratedB_);}
275 
277  bool shouldEquilibrate() {computeEquilibrateScaling(); return(shouldEquilibrate_);}
278 
280  bool solutionErrorsEstimated() {return(solutionErrorsEstimated_);}
281 
283  bool inverted() {return(inverted_);}
284 
286  bool reciprocalConditionEstimated() {return(reciprocalConditionEstimated_);}
287 
289  bool solved() {return(solved_);}
290 
292  bool solutionRefined() {return(solutionRefined_);}
294 
296 
297 
300 
303 
306 
309 
311  OrdinalType numRows() const {return(numRowCols_);}
312 
314  OrdinalType numCols() const {return(numRowCols_);}
315 
317  MagnitudeType ANORM() const {return(ANORM_);}
318 
320  MagnitudeType RCOND() const {return(RCOND_);}
321 
323 
325  MagnitudeType SCOND() {return(SCOND_);};
326 
328  MagnitudeType AMAX() const {return(AMAX_);}
329 
331  std::vector<MagnitudeType> FERR() const {return(FERR_);}
332 
334  std::vector<MagnitudeType> BERR() const {return(BERR_);}
335 
337  std::vector<MagnitudeType> R() const {return(R_);}
338 
340 
342 
343  void Print(std::ostream& os) const;
346 
347  protected:
348 
349  void allocateWORK() { LWORK_ = 4*numRowCols_; WORK_.resize( LWORK_ ); return;}
350  void allocateIWORK() { IWORK_.resize( numRowCols_ ); return;}
351  void resetMatrix();
352  void resetVectors();
353 
354  bool equilibrate_;
355  bool shouldEquilibrate_;
356  bool equilibratedA_;
357  bool equilibratedB_;
358  bool transpose_;
359  bool factored_;
360  bool estimateSolutionErrors_;
361  bool solutionErrorsEstimated_;
362  bool solved_;
363  bool inverted_;
364  bool reciprocalConditionEstimated_;
365  bool refineSolution_;
366  bool solutionRefined_;
367 
368  OrdinalType numRowCols_;
369  OrdinalType LDA_;
370  OrdinalType LDAF_;
371  OrdinalType INFO_;
372  OrdinalType LWORK_;
373 
374  std::vector<int> IWORK_;
375 
376  MagnitudeType ANORM_;
377  MagnitudeType RCOND_;
378  MagnitudeType SCOND_;
379  MagnitudeType AMAX_;
380 
385 
386  ScalarType * A_;
387  ScalarType * AF_;
388  std::vector<MagnitudeType> FERR_;
389  std::vector<MagnitudeType> BERR_;
390  std::vector<ScalarType> WORK_;
391  std::vector<MagnitudeType> R_;
392 #ifdef HAVE_TEUCHOSNUMERICS_EIGEN
393  Eigen::LLT<EigenMatrix,Eigen::Upper> lltu_;
394  Eigen::LLT<EigenMatrix,Eigen::Lower> lltl_;
395 #endif
396 
397  private:
398  // SerialSpdDenseSolver copy constructor (put here because we don't want user access)
399 
402 
403  };
404 
405  // Helper traits to distinguish work arrays for real and complex-valued datatypes.
406  using namespace details;
407 
408 //=============================================================================
409 
410 template<typename OrdinalType, typename ScalarType>
412  : CompObject(),
413  equilibrate_(false),
414  shouldEquilibrate_(false),
415  equilibratedA_(false),
416  equilibratedB_(false),
417  transpose_(false),
418  factored_(false),
419  estimateSolutionErrors_(false),
420  solutionErrorsEstimated_(false),
421  solved_(false),
422  inverted_(false),
423  reciprocalConditionEstimated_(false),
424  refineSolution_(false),
425  solutionRefined_(false),
426  numRowCols_(0),
427  LDA_(0),
428  LDAF_(0),
429  INFO_(0),
430  LWORK_(0),
431  ANORM_(ScalarTraits<MagnitudeType>::zero()),
432  RCOND_(ScalarTraits<MagnitudeType>::zero()),
433  SCOND_(ScalarTraits<MagnitudeType>::zero()),
434  AMAX_(ScalarTraits<MagnitudeType>::zero()),
435  A_(0),
436  AF_(0)
437 {
438  resetMatrix();
439 }
440 //=============================================================================
441 
442 template<typename OrdinalType, typename ScalarType>
444 {}
445 
446 //=============================================================================
447 
448 template<typename OrdinalType, typename ScalarType>
450 {
451  LHS_ = Teuchos::null;
452  RHS_ = Teuchos::null;
453  reciprocalConditionEstimated_ = false;
454  solutionRefined_ = false;
455  solved_ = false;
456  solutionErrorsEstimated_ = false;
457  equilibratedB_ = false;
458 }
459 //=============================================================================
460 
461 template<typename OrdinalType, typename ScalarType>
463 {
464  resetVectors();
465  equilibratedA_ = false;
466  factored_ = false;
467  inverted_ = false;
468  numRowCols_ = 0;
469  LDA_ = 0;
470  LDAF_ = 0;
475  A_ = 0;
476  AF_ = 0;
477  INFO_ = 0;
478  LWORK_ = 0;
479  R_.resize(0);
480 }
481 //=============================================================================
482 
483 template<typename OrdinalType, typename ScalarType>
485 {
486  resetMatrix();
487  Matrix_ = A;
488  Factor_ = A;
489  numRowCols_ = A->numRows();
490  LDA_ = A->stride();
491  LDAF_ = LDA_;
492  A_ = A->values();
493  AF_ = A->values();
494  return(0);
495 }
496 //=============================================================================
497 
498 template<typename OrdinalType, typename ScalarType>
501 {
502  // Check that these new vectors are consistent.
503  TEUCHOS_TEST_FOR_EXCEPTION(B->numRows()!=X->numRows() || B->numCols() != X->numCols(), std::invalid_argument,
504  "SerialSpdDenseSolver<T>::setVectors: X and B are not the same size!");
505  TEUCHOS_TEST_FOR_EXCEPTION(B->values()==0, std::invalid_argument,
506  "SerialSpdDenseSolver<T>::setVectors: B is an empty SerialDenseMatrix<T>!");
507  TEUCHOS_TEST_FOR_EXCEPTION(X->values()==0, std::invalid_argument,
508  "SerialSpdDenseSolver<T>::setVectors: X is an empty SerialDenseMatrix<T>!");
509  TEUCHOS_TEST_FOR_EXCEPTION(B->stride()<1, std::invalid_argument,
510  "SerialSpdDenseSolver<T>::setVectors: B has an invalid stride!");
511  TEUCHOS_TEST_FOR_EXCEPTION(X->stride()<1, std::invalid_argument,
512  "SerialSpdDenseSolver<T>::setVectors: X has an invalid stride!");
513 
514  resetVectors();
515  LHS_ = X;
516  RHS_ = B;
517  return(0);
518 }
519 //=============================================================================
520 
521 template<typename OrdinalType, typename ScalarType>
523 {
524  estimateSolutionErrors_ = flag;
525 
526  // If the errors are estimated, this implies that the solution must be refined
527  refineSolution_ = refineSolution_ || flag;
528 }
529 //=============================================================================
530 
531 template<typename OrdinalType, typename ScalarType>
533 
534  if (factored()) return(0); // Already factored
535 
536  TEUCHOS_TEST_FOR_EXCEPTION(inverted(), std::logic_error,
537  "SerialSpdDenseSolver<T>::factor: Cannot factor an inverted matrix!");
538 
539  ANORM_ = Matrix_->normOne(); // Compute 1-Norm of A
540 
541 
542  // If we want to refine the solution, then the factor must
543  // be stored separatedly from the original matrix
544 
545  if (A_ == AF_)
546  if (refineSolution_ ) {
547  Factor_ = rcp( new SerialSymDenseMatrix<OrdinalType,ScalarType>(*Matrix_) );
548  AF_ = Factor_->values();
549  LDAF_ = Factor_->stride();
550  }
551 
552  int ierr = 0;
553  if (equilibrate_) ierr = equilibrateMatrix();
554 
555  if (ierr!=0) return(ierr);
556 
557  INFO_ = 0;
558 
559 #ifdef HAVE_TEUCHOSNUMERICS_EIGEN
560  EigenMatrixMap aMap( AF_, numRowCols_, numRowCols_, EigenOuterStride(LDAF_) );
561 
562  if (Matrix_->upper())
563  lltu_.compute(aMap);
564  else
565  lltl_.compute(aMap);
566 
567 #else
568  this->POTRF(Matrix_->UPLO(), numRowCols_, AF_, LDAF_, &INFO_);
569 #endif
570 
571  factored_ = true;
572 
573  return(INFO_);
574 }
575 
576 //=============================================================================
577 
578 template<typename OrdinalType, typename ScalarType>
580 
581  // We will call one of four routines depending on what services the user wants and
582  // whether or not the matrix has been inverted or factored already.
583  //
584  // If the matrix has been inverted, use DGEMM to compute solution.
585  // Otherwise, if the user want the matrix to be equilibrated or wants a refined solution, we will
586  // call the X interface.
587  // Otherwise, if the matrix is already factored we will call the TRS interface.
588  // Otherwise, if the matrix is unfactored we will call the SV interface.
589 
590  int ierr = 0;
591  if (equilibrate_) {
592  ierr = equilibrateRHS();
593  equilibratedB_ = true;
594  }
595  if (ierr != 0) return(ierr); // Can't equilibrate B, so return.
596 
597  TEUCHOS_TEST_FOR_EXCEPTION( (equilibratedA_ && !equilibratedB_) || (!equilibratedA_ && equilibratedB_) ,
598  std::logic_error, "SerialSpdDenseSolver<T>::solve: Matrix and vectors must be similarly scaled!");
599  TEUCHOS_TEST_FOR_EXCEPTION( RHS_==Teuchos::null, std::invalid_argument,
600  "SerialSpdDenseSolver<T>::solve: No right-hand side vector (RHS) has been set for the linear system!");
601  TEUCHOS_TEST_FOR_EXCEPTION( LHS_==Teuchos::null, std::invalid_argument,
602  "SerialSpdDenseSolver<T>::solve: No solution vector (LHS) has been set for the linear system!");
603 
604  if (shouldEquilibrate() && !equilibratedA_)
605  std::cout << "WARNING! SerialSpdDenseSolver<T>::solve: System should be equilibrated!" << std::endl;
606 
607  if (inverted()) {
608 
609  TEUCHOS_TEST_FOR_EXCEPTION( RHS_->values() == LHS_->values(), std::invalid_argument,
610  "SerialSpdDenseSolver<T>::solve: X and B must be different vectors if matrix is inverted.");
611 
612  INFO_ = 0;
613  this->GEMM(Teuchos::NO_TRANS, Teuchos::NO_TRANS, numRowCols_, RHS_->numCols(),
614  numRowCols_, 1.0, AF_, LDAF_, RHS_->values(), RHS_->stride(), 0.0,
615  LHS_->values(), LHS_->stride());
616  if (INFO_!=0) return(INFO_);
617  solved_ = true;
618  }
619  else {
620 
621  if (!factored()) factor(); // Matrix must be factored
622 
623  if (RHS_->values()!=LHS_->values()) {
624  (*LHS_) = (*RHS_); // Copy B to X if needed
625  }
626  INFO_ = 0;
627 
628 #ifdef HAVE_TEUCHOSNUMERICS_EIGEN
629  EigenMatrixMap rhsMap( RHS_->values(), RHS_->numRows(), RHS_->numCols(), EigenOuterStride(RHS_->stride()) );
630  EigenMatrixMap lhsMap( LHS_->values(), LHS_->numRows(), LHS_->numCols(), EigenOuterStride(LHS_->stride()) );
631 
632  if (Matrix_->upper())
633  lhsMap = lltu_.solve(rhsMap);
634  else
635  lhsMap = lltl_.solve(rhsMap);
636 
637 #else
638  this->POTRS(Matrix_->UPLO(), numRowCols_, RHS_->numCols(), AF_, LDAF_, LHS_->values(), LHS_->stride(), &INFO_);
639 #endif
640 
641  if (INFO_!=0) return(INFO_);
642  solved_ = true;
643 
644  }
645  int ierr1=0;
646  if (refineSolution_ && !inverted()) ierr1 = applyRefinement();
647  if (ierr1!=0)
648  return(ierr1);
649 
650  if (equilibrate_) ierr1 = unequilibrateLHS();
651  return(ierr1);
652 }
653 //=============================================================================
654 
655 template<typename OrdinalType, typename ScalarType>
657 {
658  TEUCHOS_TEST_FOR_EXCEPTION(!solved(), std::logic_error,
659  "SerialSpdDenseSolver<T>::applyRefinement: Must have an existing solution!");
660  TEUCHOS_TEST_FOR_EXCEPTION(A_==AF_, std::logic_error,
661  "SerialSpdDenseSolver<T>::applyRefinement: Cannot apply refinement if no original copy of A!");
662 
663 
664 #ifdef HAVE_TEUCHOSNUMERICS_EIGEN
665  // Implement templated PORFS or use Eigen.
666  return (-1);
667 #else
668  OrdinalType NRHS = RHS_->numCols();
669  FERR_.resize( NRHS );
670  BERR_.resize( NRHS );
671  allocateWORK();
672  allocateIWORK();
673 
674  INFO_ = 0;
675  std::vector<typename details::lapack_traits<ScalarType>::iwork_type> PORFS_WORK( numRowCols_ );
676  this->PORFS(Matrix_->UPLO(), numRowCols_, NRHS, A_, LDA_, AF_, LDAF_,
677  RHS_->values(), RHS_->stride(), LHS_->values(), LHS_->stride(),
678  &FERR_[0], &BERR_[0], &WORK_[0], &PORFS_WORK[0], &INFO_);
679 
680  solutionErrorsEstimated_ = true;
681  reciprocalConditionEstimated_ = true;
682  solutionRefined_ = true;
683 
684  return(INFO_);
685 #endif
686 }
687 
688 //=============================================================================
689 
690 template<typename OrdinalType, typename ScalarType>
692 {
693  if (R_.size()!=0) return(0); // Already computed
694 
695  R_.resize( numRowCols_ );
696 
697  INFO_ = 0;
698  this->POEQU (numRowCols_, AF_, LDAF_, &R_[0], &SCOND_, &AMAX_, &INFO_);
699  if (SCOND_<0.1*ScalarTraits<MagnitudeType>::one() ||
700  AMAX_ < ScalarTraits<ScalarType>::rmin() ||
702  shouldEquilibrate_ = true;
703 
704  return(INFO_);
705 }
706 
707 //=============================================================================
708 
709 template<typename OrdinalType, typename ScalarType>
711 {
712  OrdinalType i, j;
713  int ierr = 0;
714 
715  if (equilibratedA_) return(0); // Already done.
716  if (R_.size()==0) ierr = computeEquilibrateScaling(); // Compute R if needed.
717  if (ierr!=0) return(ierr); // If return value is not zero, then can't equilibrate.
718  if (Matrix_->upper()) {
719  if (A_==AF_) {
720  ScalarType * ptr;
721  for (j=0; j<numRowCols_; j++) {
722  ptr = A_ + j*LDA_;
723  ScalarType s1 = R_[j];
724  for (i=0; i<=j; i++) {
725  *ptr = *ptr*s1*R_[i];
726  ptr++;
727  }
728  }
729  }
730  else {
731  ScalarType * ptr;
732  ScalarType * ptr1;
733  for (j=0; j<numRowCols_; j++) {
734  ptr = A_ + j*LDA_;
735  ptr1 = AF_ + j*LDAF_;
736  ScalarType s1 = R_[j];
737  for (i=0; i<=j; i++) {
738  *ptr = *ptr*s1*R_[i];
739  ptr++;
740  *ptr1 = *ptr1*s1*R_[i];
741  ptr1++;
742  }
743  }
744  }
745  }
746  else {
747  if (A_==AF_) {
748  ScalarType * ptr;
749  for (j=0; j<numRowCols_; j++) {
750  ptr = A_ + j + j*LDA_;
751  ScalarType s1 = R_[j];
752  for (i=j; i<numRowCols_; i++) {
753  *ptr = *ptr*s1*R_[i];
754  ptr++;
755  }
756  }
757  }
758  else {
759  ScalarType * ptr;
760  ScalarType * ptr1;
761  for (j=0; j<numRowCols_; j++) {
762  ptr = A_ + j + j*LDA_;
763  ptr1 = AF_ + j + j*LDAF_;
764  ScalarType s1 = R_[j];
765  for (i=j; i<numRowCols_; i++) {
766  *ptr = *ptr*s1*R_[i];
767  ptr++;
768  *ptr1 = *ptr1*s1*R_[i];
769  ptr1++;
770  }
771  }
772  }
773  }
774 
775  equilibratedA_ = true;
776 
777  return(0);
778 }
779 
780 //=============================================================================
781 
782 template<typename OrdinalType, typename ScalarType>
784 {
785  OrdinalType i, j;
786  int ierr = 0;
787 
788  if (equilibratedB_) return(0); // Already done.
789  if (R_.size()==0) ierr = computeEquilibrateScaling(); // Compute R if needed.
790  if (ierr!=0) return(ierr); // Can't count on R being computed.
791 
792  OrdinalType LDB = RHS_->stride(), NRHS = RHS_->numCols();
793  ScalarType * B = RHS_->values();
794  ScalarType * ptr;
795  for (j=0; j<NRHS; j++) {
796  ptr = B + j*LDB;
797  for (i=0; i<numRowCols_; i++) {
798  *ptr = *ptr*R_[i];
799  ptr++;
800  }
801  }
802 
803  equilibratedB_ = true;
804 
805  return(0);
806 }
807 
808 //=============================================================================
809 
810 template<typename OrdinalType, typename ScalarType>
812 {
813  OrdinalType i, j;
814 
815  if (!equilibratedB_) return(0); // Nothing to do
816 
817  OrdinalType LDX = LHS_->stride(), NLHS = LHS_->numCols();
818  ScalarType * X = LHS_->values();
819  ScalarType * ptr;
820  for (j=0; j<NLHS; j++) {
821  ptr = X + j*LDX;
822  for (i=0; i<numRowCols_; i++) {
823  *ptr = *ptr*R_[i];
824  ptr++;
825  }
826  }
827 
828  return(0);
829 }
830 
831 //=============================================================================
832 
833 template<typename OrdinalType, typename ScalarType>
835 {
836 #ifdef HAVE_TEUCHOSNUMERICS_EIGEN
837  // Implement templated inverse or use Eigen.
838  return (-1);
839 #else
840  if (!factored()) factor(); // Need matrix factored.
841 
842  INFO_ = 0;
843  this->POTRI( Matrix_->UPLO(), numRowCols_, AF_, LDAF_, &INFO_);
844 
845  // Copy lower/upper triangle to upper/lower triangle to make full inverse.
846  if (Matrix_->upper())
847  Matrix_->setLower();
848  else
849  Matrix_->setUpper();
850 
851  inverted_ = true;
852  factored_ = false;
853 
854  return(INFO_);
855 #endif
856 }
857 
858 //=============================================================================
859 
860 template<typename OrdinalType, typename ScalarType>
862 {
863 #ifdef HAVE_TEUCHOSNUMERICS_EIGEN
864  // Implement templated GECON or use Eigen. Eigen currently doesn't have a condition estimation function.
865  return (-1);
866 #else
868  Value = RCOND_;
869  return(0); // Already computed, just return it.
870  }
871 
872  if ( ANORM_<ScalarTraits<MagnitudeType>::zero() ) ANORM_ = Matrix_->normOne();
873 
874  int ierr = 0;
875  if (!factored()) ierr = factor(); // Need matrix factored.
876  if (ierr!=0) return(ierr);
877 
878  allocateWORK();
879  allocateIWORK();
880 
881  // We will assume a one-norm condition number
882  INFO_ = 0;
883  std::vector<typename details::lapack_traits<ScalarType>::iwork_type> POCON_WORK( numRowCols_ );
884  this->POCON( Matrix_->UPLO(), numRowCols_, AF_, LDAF_, ANORM_, &RCOND_, &WORK_[0], &POCON_WORK[0], &INFO_);
885  reciprocalConditionEstimated_ = true;
886  Value = RCOND_;
887 
888  return(INFO_);
889 #endif
890 }
891 //=============================================================================
892 
893 template<typename OrdinalType, typename ScalarType>
895 
896  if (Matrix_!=Teuchos::null) os << "Solver Matrix" << std::endl << *Matrix_ << std::endl;
897  if (Factor_!=Teuchos::null) os << "Solver Factored Matrix" << std::endl << *Factor_ << std::endl;
898  if (LHS_ !=Teuchos::null) os << "Solver LHS" << std::endl << *LHS_ << std::endl;
899  if (RHS_ !=Teuchos::null) os << "Solver RHS" << std::endl << *RHS_ << std::endl;
900 
901 }
902 
903 } // namespace Teuchos
904 
905 #endif /* TEUCHOS_SERIALSPDDENSESOLVER_H */
MagnitudeType RCOND() const
Returns the reciprocal of the condition number of the this matrix (returns -1 if not yet computed)...
void POCON(const char UPLO, const OrdinalType n, const ScalarType *A, const OrdinalType lda, const ScalarType anorm, ScalarType *rcond, ScalarType *WORK, OrdinalType *IWORK, OrdinalType *info) const
Estimates the reciprocal of the condition number (1-norm) of a real symmetric positive definite matri...
void solveToRefinedSolution(bool flag)
Causes all solves to compute solution to best ability using iterative refinement. ...
A class for constructing and using Hermitian positive definite dense matrices.
std::vector< MagnitudeType > R() const
Returns a pointer to the row scaling vector used for equilibration.
T magnitudeType
Mandatory typedef for result of magnitude.
MagnitudeType SCOND()
Ratio of smallest to largest equilibration scale factors for the this matrix (returns -1 if not yet c...
OrdinalType numCols() const
Returns column dimension of system.
Templated serial dense matrix class.
bool transpose()
Returns true if transpose of this matrix has and will be used.
bool solutionErrorsEstimated()
Returns true if forward and backward error estimated have been computed (available via FERR() and BER...
Templated interface class to BLAS routines.
#define TEUCHOS_TEST_FOR_EXCEPTION(throw_exception_test, Exception, msg)
Macro for throwing an exception with breakpointing to ease debugging.
int factor()
Computes the in-place Cholesky factorization of the matrix using the LAPACK routine DPOTRF...
Teuchos header file which uses auto-configuration information to include necessary C++ headers...
void estimateSolutionErrors(bool flag)
Causes all solves to estimate the forward and backward solution error.
Templated BLAS wrapper.
bool factored()
Returns true if matrix is factored (factor available via AF() and LDAF()).
This class creates and provides basic support for symmetric, positive-definite dense matrices of temp...
int invert()
Inverts the this matrix.
Object for storing data and providing functionality that is common to all computational classes...
MagnitudeType AMAX() const
Returns the absolute value of the largest entry of the this matrix (returns -1 if not yet computed)...
This structure defines some basic traits for a scalar field type.
int equilibrateMatrix()
Equilibrates the this matrix.
RCP< SerialDenseMatrix< OrdinalType, ScalarType > > getLHS() const
Returns pointer to current LHS.
Templated class for solving dense linear problems.
Functionality and data that is common to all computational classes.
Templated interface class to LAPACK routines.
void POTRF(const char UPLO, const OrdinalType n, ScalarType *A, const OrdinalType lda, OrdinalType *info) const
Computes Cholesky factorization of a real symmetric positive definite matrix A.
bool shouldEquilibrate()
Returns true if the LAPACK general rules for equilibration suggest you should equilibrate the system...
void POTRS(const char UPLO, const OrdinalType n, const OrdinalType nrhs, const ScalarType *A, const OrdinalType lda, ScalarType *B, const OrdinalType ldb, OrdinalType *info) const
Solves a system of linear equations A*X=B, where A is a symmetric positive definite matrix factored b...
TEUCHOS_DEPRECATED RCP< T > rcp(T *p, Dealloc_T dealloc, bool owns_mem)
Deprecated.
int setMatrix(const RCP< SerialSymDenseMatrix< OrdinalType, ScalarType > > &A_in)
Sets the pointers for coefficient matrix.
virtual ~SerialSpdDenseSolver()
SerialSpdDenseSolver destructor.
Teuchos implementation details.
int solve()
Computes the solution X to AX = B for the this matrix and the B provided to SetVectors()..
int setVectors(const RCP< SerialDenseMatrix< OrdinalType, ScalarType > > &X, const RCP< SerialDenseMatrix< OrdinalType, ScalarType > > &B)
Sets the pointers for left and right hand side vector(s).
The base Teuchos class.
void GEMM(ETransp transa, ETransp transb, const OrdinalType m, const OrdinalType n, const OrdinalType k, const alpha_type alpha, const A_type *A, const OrdinalType lda, const B_type *B, const OrdinalType ldb, const beta_type beta, ScalarType *C, const OrdinalType ldc) const
General matrix-matrix multiply.
void PORFS(const char UPLO, const OrdinalType n, const OrdinalType nrhs, const ScalarType *A, const OrdinalType lda, const ScalarType *AF, const OrdinalType ldaf, const ScalarType *B, const OrdinalType ldb, ScalarType *X, const OrdinalType ldx, ScalarType *FERR, ScalarType *BERR, ScalarType *WORK, OrdinalType *IWORK, OrdinalType *info) const
Improves the computed solution to a system of linear equations when the coefficient matrix is symmetr...
int applyRefinement()
Apply Iterative Refinement.
RCP< SerialDenseMatrix< OrdinalType, ScalarType > > getRHS() const
Returns pointer to current RHS.
The Templated LAPACK Wrapper Class.
MagnitudeType ANORM() const
Returns the 1-Norm of the this matrix (returns -1 if not yet computed).
bool solved()
Returns true if the current set of vectors has been solved.
Templated serial, dense, symmetric matrix class.
bool equilibratedA()
Returns true if factor is equilibrated (factor available via AF() and LDAF()).
bool reciprocalConditionEstimated()
Returns true if the condition number of the this matrix has been computed (value available via Recipr...
The Teuchos namespace contains all of the classes, structs and enums used by Teuchos, as well as a number of utility routines.
int reciprocalConditionEstimate(MagnitudeType &Value)
Returns the reciprocal of the 1-norm condition number of the this matrix.
void factorWithEquilibration(bool flag)
Causes equilibration to be called just before the matrix factorization as part of the call to factor...
bool inverted()
Returns true if matrix inverse has been computed (inverse available via AF() and LDAF()).
RCP< SerialSymDenseMatrix< OrdinalType, ScalarType > > getFactoredMatrix() const
Returns pointer to factored matrix (assuming factorization has been performed).
OrdinalType numRows() const
Returns row dimension of system.
Smart reference counting pointer class for automatic garbage collection.
bool equilibratedB()
Returns true if RHS is equilibrated (RHS available via B() and LDB()).
int unequilibrateLHS()
Unscales the solution vectors if equilibration was used to solve the system.
int computeEquilibrateScaling()
Computes the scaling vector S(i) = 1/sqrt(A(i,i) of the this matrix.
std::vector< MagnitudeType > FERR() const
Returns a pointer to the forward error estimates computed by LAPACK.
std::vector< MagnitudeType > BERR() const
Returns a pointer to the backward error estimates computed by LAPACK.
bool solutionRefined()
Returns true if the current set of vectors has been refined.
RCP< SerialSymDenseMatrix< OrdinalType, ScalarType > > getMatrix() const
Returns pointer to current matrix.
void Print(std::ostream &os) const
Print service methods; defines behavior of ostream << operator.
Reference-counted pointer class and non-member templated function implementations.
void POTRI(const char UPLO, const OrdinalType n, ScalarType *A, const OrdinalType lda, OrdinalType *info) const
Computes the inverse of a real symmetric positive definite matrix A using the Cholesky factorization ...
void POEQU(const OrdinalType n, const ScalarType *A, const OrdinalType lda, MagnitudeType *S, MagnitudeType *scond, MagnitudeType *amax, OrdinalType *info) const
Computes row and column scalings intended to equilibrate a symmetric positive definite matrix A and r...
static T one()
Returns representation of one for this scalar type.
SerialSpdDenseSolver()
Default constructor; matrix should be set using setMatrix(), LHS and RHS set with setVectors()...
int equilibrateRHS()
Equilibrates the current RHS.
This class creates and provides basic support for dense rectangular matrix of templated type...