Teuchos Package Browser (Single Doxygen Collection) Version of the Day
Loading...
Searching...
No Matches
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
51#include "Teuchos_BLAS.hpp"
52#include "Teuchos_LAPACK.hpp"
53#include "Teuchos_RCP.hpp"
58
59#ifdef HAVE_TEUCHOSNUMERICS_EIGEN
60#include "Eigen/Dense"
61#endif
62
130namespace 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
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
157
159
160
162 virtual ~SerialSpdDenseSolver();
164
166
167
170
172
178
180
181
183
186
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
239
242 int applyRefinement();
243
245
248 int unequilibrateLHS();
249
251
259
261
262
264 bool transpose() {return(transpose_);}
265
267 bool factored() {return(factored_);}
268
271
274
277
280
282 bool inverted() {return(inverted_);}
283
286
288 bool solved() {return(solved_);}
289
293
295
296
299
302
305
308
310 OrdinalType numRows() const {return(numRowCols_);}
311
313 OrdinalType numCols() const {return(numRowCols_);}
314
316 MagnitudeType ANORM() const {return(ANORM_);}
317
319 MagnitudeType RCOND() const {return(RCOND_);}
320
322
325
327 MagnitudeType AMAX() const {return(AMAX_);}
328
330 std::vector<MagnitudeType> FERR() const {return(FERR_);}
331
333 std::vector<MagnitudeType> BERR() const {return(BERR_);}
334
336 std::vector<MagnitudeType> R() const {return(R_);}
337
339
341
342
343 void Print(std::ostream& os) const;
345
346 protected:
347
348 void allocateWORK() { LWORK_ = 4*numRowCols_; WORK_.resize( LWORK_ ); return;}
349 void allocateIWORK() { IWORK_.resize( numRowCols_ ); return;}
350 void resetMatrix();
351 void resetVectors();
352
366
372
373 std::vector<int> IWORK_;
374
379
384
387 std::vector<MagnitudeType> FERR_;
388 std::vector<MagnitudeType> BERR_;
389 std::vector<ScalarType> WORK_;
390 std::vector<MagnitudeType> R_;
391#ifdef HAVE_TEUCHOSNUMERICS_EIGEN
392 Eigen::LLT<EigenMatrix,Eigen::Upper> lltu_;
393 Eigen::LLT<EigenMatrix,Eigen::Lower> lltl_;
394#endif
395
396 private:
397 // SerialSpdDenseSolver copy constructor (put here because we don't want user access)
398
401
402 };
403
404 // Helper traits to distinguish work arrays for real and complex-valued datatypes.
405 using namespace details;
406
407//=============================================================================
408
409template<typename OrdinalType, typename ScalarType>
411 : CompObject(),
412 equilibrate_(false),
413 shouldEquilibrate_(false),
414 equilibratedA_(false),
415 equilibratedB_(false),
416 transpose_(false),
417 factored_(false),
418 estimateSolutionErrors_(false),
419 solutionErrorsEstimated_(false),
420 solved_(false),
421 inverted_(false),
422 reciprocalConditionEstimated_(false),
423 refineSolution_(false),
424 solutionRefined_(false),
425 numRowCols_(0),
426 LDA_(0),
427 LDAF_(0),
428 INFO_(0),
429 LWORK_(0),
430 ANORM_(ScalarTraits<MagnitudeType>::zero()),
431 RCOND_(ScalarTraits<MagnitudeType>::zero()),
432 SCOND_(ScalarTraits<MagnitudeType>::zero()),
433 AMAX_(ScalarTraits<MagnitudeType>::zero()),
434 A_(0),
435 AF_(0)
436{
437 resetMatrix();
438}
439//=============================================================================
440
441template<typename OrdinalType, typename ScalarType>
444
445//=============================================================================
446
447template<typename OrdinalType, typename ScalarType>
449{
450 LHS_ = Teuchos::null;
451 RHS_ = Teuchos::null;
452 reciprocalConditionEstimated_ = false;
453 solutionRefined_ = false;
454 solved_ = false;
455 solutionErrorsEstimated_ = false;
456 equilibratedB_ = false;
457}
458//=============================================================================
459
460template<typename OrdinalType, typename ScalarType>
462{
463 resetVectors();
464 equilibratedA_ = false;
465 factored_ = false;
466 inverted_ = false;
467 numRowCols_ = 0;
468 LDA_ = 0;
469 LDAF_ = 0;
474 A_ = 0;
475 AF_ = 0;
476 INFO_ = 0;
477 LWORK_ = 0;
478 R_.resize(0);
479}
480//=============================================================================
481
482template<typename OrdinalType, typename ScalarType>
484{
485 resetMatrix();
486 Matrix_ = A;
487 Factor_ = A;
488 numRowCols_ = A->numRows();
489 LDA_ = A->stride();
490 LDAF_ = LDA_;
491 A_ = A->values();
492 AF_ = A->values();
493 return(0);
494}
495//=============================================================================
496
497template<typename OrdinalType, typename ScalarType>
500{
501 // Check that these new vectors are consistent.
502 TEUCHOS_TEST_FOR_EXCEPTION(B->numRows()!=X->numRows() || B->numCols() != X->numCols(), std::invalid_argument,
503 "SerialSpdDenseSolver<T>::setVectors: X and B are not the same size!");
504 TEUCHOS_TEST_FOR_EXCEPTION(B->values()==0, std::invalid_argument,
505 "SerialSpdDenseSolver<T>::setVectors: B is an empty SerialDenseMatrix<T>!");
506 TEUCHOS_TEST_FOR_EXCEPTION(X->values()==0, std::invalid_argument,
507 "SerialSpdDenseSolver<T>::setVectors: X is an empty SerialDenseMatrix<T>!");
508 TEUCHOS_TEST_FOR_EXCEPTION(B->stride()<1, std::invalid_argument,
509 "SerialSpdDenseSolver<T>::setVectors: B has an invalid stride!");
510 TEUCHOS_TEST_FOR_EXCEPTION(X->stride()<1, std::invalid_argument,
511 "SerialSpdDenseSolver<T>::setVectors: X has an invalid stride!");
512
513 resetVectors();
514 LHS_ = X;
515 RHS_ = B;
516 return(0);
517}
518//=============================================================================
519
520template<typename OrdinalType, typename ScalarType>
522{
523 estimateSolutionErrors_ = flag;
524
525 // If the errors are estimated, this implies that the solution must be refined
526 refineSolution_ = refineSolution_ || flag;
527}
528//=============================================================================
529
530template<typename OrdinalType, typename ScalarType>
532
533 if (factored()) return(0); // Already factored
534
535 TEUCHOS_TEST_FOR_EXCEPTION(inverted(), std::logic_error,
536 "SerialSpdDenseSolver<T>::factor: Cannot factor an inverted matrix!");
537
538 ANORM_ = Matrix_->normOne(); // Compute 1-Norm of A
539
540
541 // If we want to refine the solution, then the factor must
542 // be stored separatedly from the original matrix
543
544 if (A_ == AF_)
545 if (refineSolution_ ) {
546 Factor_ = rcp( new SerialSymDenseMatrix<OrdinalType,ScalarType>(*Matrix_) );
547 AF_ = Factor_->values();
548 LDAF_ = Factor_->stride();
549 }
550
551 int ierr = 0;
552 if (equilibrate_) ierr = equilibrateMatrix();
553
554 if (ierr!=0) return(ierr);
555
556 INFO_ = 0;
557
558#ifdef HAVE_TEUCHOSNUMERICS_EIGEN
559 EigenMatrixMap aMap( AF_, numRowCols_, numRowCols_, EigenOuterStride(LDAF_) );
560
561 if (Matrix_->upper())
562 lltu_.compute(aMap);
563 else
564 lltl_.compute(aMap);
565
566#else
567 this->POTRF(Matrix_->UPLO(), numRowCols_, AF_, LDAF_, &INFO_);
568#endif
569
570 factored_ = true;
571
572 return(INFO_);
573}
574
575//=============================================================================
576
577template<typename OrdinalType, typename ScalarType>
579
580 // We will call one of four routines depending on what services the user wants and
581 // whether or not the matrix has been inverted or factored already.
582 //
583 // If the matrix has been inverted, use DGEMM to compute solution.
584 // Otherwise, if the user want the matrix to be equilibrated or wants a refined solution, we will
585 // call the X interface.
586 // Otherwise, if the matrix is already factored we will call the TRS interface.
587 // Otherwise, if the matrix is unfactored we will call the SV interface.
588
589 int ierr = 0;
590 if (equilibrate_) {
591 ierr = equilibrateRHS();
592 }
593 if (ierr != 0) return(ierr); // Can't equilibrate B, so return.
594
595 TEUCHOS_TEST_FOR_EXCEPTION( RHS_==Teuchos::null, std::invalid_argument,
596 "SerialSpdDenseSolver<T>::solve: No right-hand side vector (RHS) has been set for the linear system!");
597 TEUCHOS_TEST_FOR_EXCEPTION( LHS_==Teuchos::null, std::invalid_argument,
598 "SerialSpdDenseSolver<T>::solve: No solution vector (LHS) has been set for the linear system!");
599
600 if (inverted()) {
601
602 TEUCHOS_TEST_FOR_EXCEPTION( RHS_->values() == LHS_->values(), std::invalid_argument,
603 "SerialSpdDenseSolver<T>::solve: X and B must be different vectors if matrix is inverted.");
604 TEUCHOS_TEST_FOR_EXCEPTION( (equilibratedA_ && !equilibratedB_) || (!equilibratedA_ && equilibratedB_) ,
605 std::logic_error, "SerialSpdDenseSolver<T>::solve: Matrix and vectors must be similarly scaled!");
606
607 INFO_ = 0;
608 this->GEMM(Teuchos::NO_TRANS, Teuchos::NO_TRANS, numRowCols_, RHS_->numCols(),
609 numRowCols_, 1.0, AF_, LDAF_, RHS_->values(), RHS_->stride(), 0.0,
610 LHS_->values(), LHS_->stride());
611 if (INFO_!=0) return(INFO_);
612 solved_ = true;
613 }
614 else {
615
616 if (!factored()) factor(); // Matrix must be factored
617
618 TEUCHOS_TEST_FOR_EXCEPTION( (equilibratedA_ && !equilibratedB_) || (!equilibratedA_ && equilibratedB_) ,
619 std::logic_error, "SerialSpdDenseSolver<T>::solve: Matrix and vectors must be similarly scaled!");
620
621 if (RHS_->values()!=LHS_->values()) {
622 (*LHS_) = (*RHS_); // Copy B to X if needed
623 }
624 INFO_ = 0;
625
626#ifdef HAVE_TEUCHOSNUMERICS_EIGEN
627 EigenMatrixMap rhsMap( RHS_->values(), RHS_->numRows(), RHS_->numCols(), EigenOuterStride(RHS_->stride()) );
628 EigenMatrixMap lhsMap( LHS_->values(), LHS_->numRows(), LHS_->numCols(), EigenOuterStride(LHS_->stride()) );
629
630 if (Matrix_->upper())
631 lhsMap = lltu_.solve(rhsMap);
632 else
633 lhsMap = lltl_.solve(rhsMap);
634
635#else
636 this->POTRS(Matrix_->UPLO(), numRowCols_, RHS_->numCols(), AF_, LDAF_, LHS_->values(), LHS_->stride(), &INFO_);
637#endif
638
639 if (INFO_!=0) return(INFO_);
640 solved_ = true;
641
642 }
643
644 // Warn users that the system should be equilibrated, but it wasn't.
645 if (shouldEquilibrate() && !equilibratedA_)
646 std::cout << "WARNING! SerialSpdDenseSolver<T>::solve: System should be equilibrated!" << std::endl;
647
648 int ierr1=0;
649 if (refineSolution_ && !inverted()) ierr1 = applyRefinement();
650 if (ierr1!=0)
651 return(ierr1);
652
653 if (equilibrate_) ierr1 = unequilibrateLHS();
654 return(ierr1);
655}
656//=============================================================================
657
658template<typename OrdinalType, typename ScalarType>
660{
661 TEUCHOS_TEST_FOR_EXCEPTION(!solved(), std::logic_error,
662 "SerialSpdDenseSolver<T>::applyRefinement: Must have an existing solution!");
663 TEUCHOS_TEST_FOR_EXCEPTION(A_==AF_, std::logic_error,
664 "SerialSpdDenseSolver<T>::applyRefinement: Cannot apply refinement if no original copy of A!");
665
666
667#ifdef HAVE_TEUCHOSNUMERICS_EIGEN
668 // Implement templated PORFS or use Eigen.
669 return (-1);
670#else
671 OrdinalType NRHS = RHS_->numCols();
672 FERR_.resize( NRHS );
673 BERR_.resize( NRHS );
674 allocateWORK();
675 allocateIWORK();
676
677 INFO_ = 0;
678 std::vector<typename details::lapack_traits<ScalarType>::iwork_type> PORFS_WORK( numRowCols_ );
679 this->PORFS(Matrix_->UPLO(), numRowCols_, NRHS, A_, LDA_, AF_, LDAF_,
680 RHS_->values(), RHS_->stride(), LHS_->values(), LHS_->stride(),
681 &FERR_[0], &BERR_[0], &WORK_[0], &PORFS_WORK[0], &INFO_);
682
683 solutionErrorsEstimated_ = true;
684 reciprocalConditionEstimated_ = true;
685 solutionRefined_ = true;
686
687 return(INFO_);
688#endif
689}
690
691//=============================================================================
692
693template<typename OrdinalType, typename ScalarType>
695{
696 if (R_.size()!=0) return(0); // Already computed
697
698 R_.resize( numRowCols_ );
699
700 INFO_ = 0;
701 this->POEQU (numRowCols_, AF_, LDAF_, &R_[0], &SCOND_, &AMAX_, &INFO_);
702 if (SCOND_<0.1*ScalarTraits<MagnitudeType>::one() ||
705 shouldEquilibrate_ = true;
706
707 return(INFO_);
708}
709
710//=============================================================================
711
712template<typename OrdinalType, typename ScalarType>
714{
715 OrdinalType i, j;
716 int ierr = 0;
717
718 if (equilibratedA_) return(0); // Already done.
719 if (R_.size()==0) ierr = computeEquilibrateScaling(); // Compute R if needed.
720 if (ierr!=0) return(ierr); // If return value is not zero, then can't equilibrate.
721 if (Matrix_->upper()) {
722 if (A_==AF_) {
723 ScalarType * ptr;
724 for (j=0; j<numRowCols_; j++) {
725 ptr = A_ + j*LDA_;
726 ScalarType s1 = R_[j];
727 for (i=0; i<=j; i++) {
728 *ptr = *ptr*s1*R_[i];
729 ptr++;
730 }
731 }
732 }
733 else {
734 ScalarType * ptr;
736 for (j=0; j<numRowCols_; j++) {
737 ptr = A_ + j*LDA_;
738 ptr1 = AF_ + j*LDAF_;
739 ScalarType s1 = R_[j];
740 for (i=0; i<=j; i++) {
741 *ptr = *ptr*s1*R_[i];
742 ptr++;
743 *ptr1 = *ptr1*s1*R_[i];
744 ptr1++;
745 }
746 }
747 }
748 }
749 else {
750 if (A_==AF_) {
751 ScalarType * ptr;
752 for (j=0; j<numRowCols_; j++) {
753 ptr = A_ + j + j*LDA_;
754 ScalarType s1 = R_[j];
755 for (i=j; i<numRowCols_; i++) {
756 *ptr = *ptr*s1*R_[i];
757 ptr++;
758 }
759 }
760 }
761 else {
762 ScalarType * ptr;
764 for (j=0; j<numRowCols_; j++) {
765 ptr = A_ + j + j*LDA_;
766 ptr1 = AF_ + j + j*LDAF_;
767 ScalarType s1 = R_[j];
768 for (i=j; i<numRowCols_; i++) {
769 *ptr = *ptr*s1*R_[i];
770 ptr++;
771 *ptr1 = *ptr1*s1*R_[i];
772 ptr1++;
773 }
774 }
775 }
776 }
777
778 equilibratedA_ = true;
779
780 return(0);
781}
782
783//=============================================================================
784
785template<typename OrdinalType, typename ScalarType>
787{
788 OrdinalType i, j;
789 int ierr = 0;
790
791 if (equilibratedB_) return(0); // Already done.
792 if (R_.size()==0) ierr = computeEquilibrateScaling(); // Compute R if needed.
793 if (ierr!=0) return(ierr); // Can't count on R being computed.
794
795 OrdinalType LDB = RHS_->stride(), NRHS = RHS_->numCols();
796 ScalarType * B = RHS_->values();
797 ScalarType * ptr;
798 for (j=0; j<NRHS; j++) {
799 ptr = B + j*LDB;
800 for (i=0; i<numRowCols_; i++) {
801 *ptr = *ptr*R_[i];
802 ptr++;
803 }
804 }
805
806 equilibratedB_ = true;
807
808 return(0);
809}
810
811//=============================================================================
812
813template<typename OrdinalType, typename ScalarType>
815{
816 OrdinalType i, j;
817
818 if (!equilibratedB_) return(0); // Nothing to do
819
820 OrdinalType LDX = LHS_->stride(), NLHS = LHS_->numCols();
821 ScalarType * X = LHS_->values();
822 ScalarType * ptr;
823 for (j=0; j<NLHS; j++) {
824 ptr = X + j*LDX;
825 for (i=0; i<numRowCols_; i++) {
826 *ptr = *ptr*R_[i];
827 ptr++;
828 }
829 }
830
831 return(0);
832}
833
834//=============================================================================
835
836template<typename OrdinalType, typename ScalarType>
838{
839#ifdef HAVE_TEUCHOSNUMERICS_EIGEN
840 // Implement templated inverse or use Eigen.
841 return (-1);
842#else
843 if (!factored()) factor(); // Need matrix factored.
844
845 INFO_ = 0;
846 this->POTRI( Matrix_->UPLO(), numRowCols_, AF_, LDAF_, &INFO_);
847
848 // Copy lower/upper triangle to upper/lower triangle to make full inverse.
849 if (Matrix_->upper())
850 Matrix_->setLower();
851 else
852 Matrix_->setUpper();
853
854 inverted_ = true;
855 factored_ = false;
856
857 return(INFO_);
858#endif
859}
860
861//=============================================================================
862
863template<typename OrdinalType, typename ScalarType>
865{
866#ifdef HAVE_TEUCHOSNUMERICS_EIGEN
867 // Implement templated GECON or use Eigen. Eigen currently doesn't have a condition estimation function.
868 return (-1);
869#else
870 if (reciprocalConditionEstimated()) {
871 Value = RCOND_;
872 return(0); // Already computed, just return it.
873 }
874
875 if ( ANORM_<ScalarTraits<MagnitudeType>::zero() ) ANORM_ = Matrix_->normOne();
876
877 int ierr = 0;
878 if (!factored()) ierr = factor(); // Need matrix factored.
879 if (ierr!=0) return(ierr);
880
881 allocateWORK();
882 allocateIWORK();
883
884 // We will assume a one-norm condition number
885 INFO_ = 0;
886 std::vector<typename details::lapack_traits<ScalarType>::iwork_type> POCON_WORK( numRowCols_ );
887 this->POCON( Matrix_->UPLO(), numRowCols_, AF_, LDAF_, ANORM_, &RCOND_, &WORK_[0], &POCON_WORK[0], &INFO_);
888 reciprocalConditionEstimated_ = true;
889 Value = RCOND_;
890
891 return(INFO_);
892#endif
893}
894//=============================================================================
895
896template<typename OrdinalType, typename ScalarType>
898
899 if (Matrix_!=Teuchos::null) os << "Solver Matrix" << std::endl << *Matrix_ << std::endl;
900 if (Factor_!=Teuchos::null) os << "Solver Factored Matrix" << std::endl << *Factor_ << std::endl;
901 if (LHS_ !=Teuchos::null) os << "Solver LHS" << std::endl << *LHS_ << std::endl;
902 if (RHS_ !=Teuchos::null) os << "Solver RHS" << std::endl << *RHS_ << std::endl;
903
904}
905
906} // namespace Teuchos
907
908#endif /* TEUCHOS_SERIALSPDDENSESOLVER_H */
Templated interface class to BLAS routines.
Object for storing data and providing functionality that is common to all computational classes.
Teuchos header file which uses auto-configuration information to include necessary C++ headers.
Templated interface class to LAPACK routines.
Reference-counted pointer class and non-member templated function implementations.
Templated serial dense matrix class.
Templated class for solving dense linear problems.
Templated serial, dense, symmetric matrix class.
Templated BLAS wrapper.
Functionality and data that is common to all computational classes.
The Templated LAPACK Wrapper Class.
The base Teuchos class.
Smart reference counting pointer class for automatic garbage collection.
Concrete serial communicator subclass.
A class for constructing and using Hermitian positive definite dense matrices.
SerialSpdDenseSolver(const SerialSpdDenseSolver< OrdinalType, ScalarType > &Source)
MagnitudeType SCOND()
Ratio of smallest to largest equilibration scale factors for the this matrix (returns -1 if not yet c...
int setMatrix(const RCP< SerialSymDenseMatrix< OrdinalType, ScalarType > > &A_in)
Sets the pointers for coefficient matrix.
OrdinalType numCols() const
Returns column dimension of system.
Teuchos::ScalarTraits< ScalarType >::magnitudeType MagnitudeType
RCP< SerialSymDenseMatrix< OrdinalType, ScalarType > > Factor_
SerialSpdDenseSolver & operator=(const SerialSpdDenseSolver< OrdinalType, ScalarType > &Source)
std::vector< MagnitudeType > FERR() const
Returns a pointer to the forward error estimates computed by LAPACK.
int applyRefinement()
Apply Iterative Refinement.
bool transpose()
Returns true if transpose of this matrix has and will be used.
int solve()
Computes the solution X to AX = B for the this matrix and the B provided to SetVectors()....
std::vector< MagnitudeType > R() const
Returns a pointer to the row scaling vector used for equilibration.
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).
int factor()
Computes the in-place Cholesky factorization of the matrix using the LAPACK routine DPOTRF.
RCP< SerialDenseMatrix< OrdinalType, ScalarType > > getRHS() const
Returns pointer to current RHS.
void estimateSolutionErrors(bool flag)
Causes all solves to estimate the forward and backward solution error.
bool solutionRefined()
Returns true if the current set of vectors has been refined.
virtual ~SerialSpdDenseSolver()
SerialSpdDenseSolver destructor.
RCP< SerialSymDenseMatrix< OrdinalType, ScalarType > > getMatrix() const
Returns pointer to current matrix.
bool solved()
Returns true if the current set of vectors has been solved.
bool reciprocalConditionEstimated()
Returns true if the condition number of the this matrix has been computed (value available via Recipr...
RCP< SerialSymDenseMatrix< OrdinalType, ScalarType > > getFactoredMatrix() const
Returns pointer to factored matrix (assuming factorization has been performed).
void solveToRefinedSolution(bool flag)
Causes all solves to compute solution to best ability using iterative refinement.
OrdinalType numRows() const
Returns row dimension of system.
RCP< SerialDenseMatrix< OrdinalType, ScalarType > > RHS_
RCP< SerialDenseMatrix< OrdinalType, ScalarType > > LHS_
MagnitudeType ANORM() const
Returns the 1-Norm of the this matrix (returns -1 if not yet computed).
bool solutionErrorsEstimated()
Returns true if forward and backward error estimated have been computed (available via FERR() and BER...
MagnitudeType AMAX() const
Returns the absolute value of the largest entry of the this matrix (returns -1 if not yet computed).
bool factored()
Returns true if matrix is factored (factor available via AF() and LDAF()).
void Print(std::ostream &os) const
Print service methods; defines behavior of ostream << operator.
SerialSpdDenseSolver()
Default constructor; matrix should be set using setMatrix(), LHS and RHS set with setVectors().
std::vector< MagnitudeType > BERR() const
Returns a pointer to the backward error estimates computed by LAPACK.
int unequilibrateLHS()
Unscales the solution vectors if equilibration was used to solve the system.
void factorWithEquilibration(bool flag)
Causes equilibration to be called just before the matrix factorization as part of the call to factor.
int equilibrateMatrix()
Equilibrates the this matrix.
int reciprocalConditionEstimate(MagnitudeType &Value)
Returns the reciprocal of the 1-norm condition number of the this matrix.
bool equilibratedB()
Returns true if RHS is equilibrated (RHS available via B() and LDB()).
RCP< SerialSymDenseMatrix< OrdinalType, ScalarType > > Matrix_
RCP< SerialDenseMatrix< OrdinalType, ScalarType > > getLHS() const
Returns pointer to current LHS.
bool inverted()
Returns true if matrix inverse has been computed (inverse available via AF() and LDAF()).
int computeEquilibrateScaling()
Computes the scaling vector S(i) = 1/sqrt(A(i,i) of the this matrix.
MagnitudeType RCOND() const
Returns the reciprocal of the condition number of the this matrix (returns -1 if not yet computed).
int equilibrateRHS()
Equilibrates the current RHS.
bool shouldEquilibrate()
Returns true if the LAPACK general rules for equilibration suggest you should equilibrate the system.
bool equilibratedA()
Returns true if factor is equilibrated (factor available via AF() and LDAF()).
#define TEUCHOS_TEST_FOR_EXCEPTION(throw_exception_test, Exception, msg)
Macro for throwing an exception with breakpointing to ease debugging.
Definition PackageA.cpp:3
Definition PackageB.cpp:3
TEUCHOS_DEPRECATED RCP< T > rcp(T *p, Dealloc_T dealloc, bool owns_mem)
Deprecated.
Teuchos implementation details.
This structure defines some basic traits for a scalar field type.
static T one()
Returns representation of one for this scalar type.
T magnitudeType
Mandatory typedef for result of magnitude.