Teuchos Package Browser (Single Doxygen Collection) Version of the Day
Loading...
Searching...
No Matches
Teuchos_SerialBandDenseSolver.hpp
Go to the documentation of this file.
1// @HEADER
2// ***********************************************************************
3//
4// Teuchos: Common Tools Package
5// Copyright (2004) 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 Michael A. Heroux (maherou@sandia.gov)
38//
39// ***********************************************************************
40// @HEADER
41
42#ifndef _TEUCHOS_SERIALBANDDENSESOLVER_HPP_
43#define _TEUCHOS_SERIALBANDDENSESOLVER_HPP_
48
50#include "Teuchos_BLAS.hpp"
51#include "Teuchos_LAPACK.hpp"
52#include "Teuchos_RCP.hpp"
58
59#ifdef HAVE_TEUCHOSNUMERICS_EIGEN
60#include "Eigen/Dense"
61#endif
62
63namespace Teuchos {
64
165 template<typename OrdinalType, typename ScalarType>
166 class SerialBandDenseSolver : public CompObject, public Object, public BLAS<OrdinalType, ScalarType>,
167 public LAPACK<OrdinalType, ScalarType>
168 {
169
170 public:
171
173#ifdef HAVE_TEUCHOSNUMERICS_EIGEN
174 typedef typename Eigen::Matrix<ScalarType,Eigen::Dynamic,Eigen::Dynamic,Eigen::ColMajor> EigenMatrix;
175 typedef typename Eigen::internal::BandMatrix<ScalarType,Eigen::Dynamic,Eigen::Dynamic,Eigen::ColMajor> EigenBandMatrix;
176 typedef typename Eigen::Matrix<ScalarType,Eigen::Dynamic,1> EigenVector;
177 typedef typename Eigen::InnerStride<Eigen::Dynamic> EigenInnerStride;
178 typedef typename Eigen::OuterStride<Eigen::Dynamic> EigenOuterStride;
179 typedef typename Eigen::Map<EigenVector,0,EigenInnerStride > EigenVectorMap;
180 typedef typename Eigen::Map<const EigenVector,0,EigenInnerStride > EigenConstVectorMap;
181 typedef typename Eigen::Map<EigenMatrix,0,EigenOuterStride > EigenMatrixMap;
182 typedef typename Eigen::Map<const EigenMatrix,0,EigenOuterStride > EigenConstMatrixMap;
183 typedef typename Eigen::PermutationMatrix<Eigen::Dynamic,Eigen::Dynamic,OrdinalType> EigenPermutationMatrix;
184 typedef typename Eigen::Array<OrdinalType,Eigen::Dynamic,1> EigenOrdinalArray;
185 typedef typename Eigen::Map<EigenOrdinalArray> EigenOrdinalArrayMap;
186 typedef typename Eigen::Array<ScalarType,Eigen::Dynamic,1> EigenScalarArray;
187 typedef typename Eigen::Map<EigenScalarArray> EigenScalarArrayMap;
188#endif
189
191
192
195
197 virtual ~SerialBandDenseSolver();
198
200
202
205
207
212
214
216
220
225
229
232
234
237 void estimateSolutionErrors(bool flag);
239
241
242
244
247 int factor();
248
250
253 int solve();
254
256
260
262
265 int equilibrateMatrix();
266
268
271 int equilibrateRHS();
272
274
277 int applyRefinement();
278
280
283 int unequilibrateLHS();
284
286
294
296
297
299 bool transpose() {return(transpose_);}
300
302 bool factored() {return(factored_);}
303
306
309
312
315
318
320 bool solved() {return(solved_);}
321
325
327
328
331
334
337
340
342 OrdinalType numRows() const {return(M_);}
343
345 OrdinalType numCols() const {return(N_);}
346
348 OrdinalType numLower() const {return(KL_);}
349
351 OrdinalType numUpper() const {return(KU_);}
352
354 std::vector<OrdinalType> IPIV() const {return(IPIV_);}
355
357 MagnitudeType ANORM() const {return(ANORM_);}
358
360 MagnitudeType RCOND() const {return(RCOND_);}
361
363
365 MagnitudeType ROWCND() const {return(ROWCND_);}
366
368
370 MagnitudeType COLCND() const {return(COLCND_);}
371
373 MagnitudeType AMAX() const {return(AMAX_);}
374
376 std::vector<MagnitudeType> FERR() const {return(FERR_);}
377
379 std::vector<MagnitudeType> BERR() const {return(BERR_);}
380
382 std::vector<MagnitudeType> R() const {return(R_);}
383
385 std::vector<MagnitudeType> C() const {return(C_);}
387
389
390
391 void Print(std::ostream& os) const;
393 protected:
394
395 void allocateWORK() { LWORK_ = 3*N_; WORK_.resize( LWORK_ ); return;}
396 void resetMatrix();
397 void resetVectors();
398
399
412
414
424
425 std::vector<OrdinalType> IPIV_;
426 std::vector<int> IWORK_;
427
433
438
441 std::vector<MagnitudeType> FERR_;
442 std::vector<MagnitudeType> BERR_;
443 std::vector<ScalarType> WORK_;
444 std::vector<MagnitudeType> R_;
445 std::vector<MagnitudeType> C_;
446#ifdef HAVE_TEUCHOSNUMERICS_EIGEN
447 Eigen::PartialPivLU<EigenMatrix> lu_;
448#endif
449
450 private:
451
452 // SerialBandDenseSolver copy constructor (put here because we don't want user access)
455
456 };
457
458 // Helper traits to distinguish work arrays for real and complex-valued datatypes.
459 using namespace details;
460
461//=============================================================================
462
463template<typename OrdinalType, typename ScalarType>
465 : CompObject(),
466 equilibrate_(false),
467 shouldEquilibrate_(false),
468 equilibratedA_(false),
469 equilibratedB_(false),
470 transpose_(false),
471 factored_(false),
472 estimateSolutionErrors_(false),
473 solutionErrorsEstimated_(false),
474 solved_(false),
475 reciprocalConditionEstimated_(false),
476 refineSolution_(false),
477 solutionRefined_(false),
478 TRANS_(Teuchos::NO_TRANS),
479 M_(0),
480 N_(0),
481 KL_(0),
482 KU_(0),
483 Min_MN_(0),
484 LDA_(0),
485 LDAF_(0),
486 INFO_(0),
487 LWORK_(0),
488 ANORM_(0.0),
489 RCOND_(ScalarTraits<MagnitudeType>::zero()),
490 ROWCND_(ScalarTraits<MagnitudeType>::zero()),
491 COLCND_(ScalarTraits<MagnitudeType>::zero()),
492 AMAX_(ScalarTraits<MagnitudeType>::zero()),
493 A_(NULL),
494 AF_(NULL)
495{
496 resetMatrix();
497}
498//=============================================================================
499
500template<typename OrdinalType, typename ScalarType>
503
504//=============================================================================
505
506template<typename OrdinalType, typename ScalarType>
508{
509 LHS_ = Teuchos::null;
510 RHS_ = Teuchos::null;
511 reciprocalConditionEstimated_ = false;
512 solutionRefined_ = false;
513 solved_ = false;
514 solutionErrorsEstimated_ = false;
515 equilibratedB_ = false;
516}
517//=============================================================================
518
519template<typename OrdinalType, typename ScalarType>
521{
522 resetVectors();
523 equilibratedA_ = false;
524 factored_ = false;
525 M_ = 0;
526 N_ = 0;
527 KL_ = 0;
528 KU_ = 0;
529 Min_MN_ = 0;
530 LDA_ = 0;
531 LDAF_ = 0;
536 A_ = 0;
537 AF_ = 0;
538 INFO_ = 0;
539 LWORK_ = 0;
540 R_.resize(0);
541 C_.resize(0);
542}
543//=============================================================================
544
545template<typename OrdinalType, typename ScalarType>
547{
548
549 OrdinalType m = AB->numRows();
550 OrdinalType n = AB->numCols();
551 OrdinalType kl = AB->lowerBandwidth();
552 OrdinalType ku = AB->upperBandwidth();
553
554 // Check that the new matrix is consistent.
555 TEUCHOS_TEST_FOR_EXCEPTION(AB->values()==0, std::invalid_argument,
556 "SerialBandDenseSolver<T>::setMatrix: A is an empty SerialBandDenseMatrix<T>!");
557
558 resetMatrix();
559 Matrix_ = AB;
560 Factor_ = AB;
561 M_ = m;
562 N_ = n;
563 KL_ = kl;
564 KU_ = ku-kl;
565 Min_MN_ = TEUCHOS_MIN(M_,N_);
566 LDA_ = AB->stride();
567 LDAF_ = LDA_;
568 A_ = AB->values();
569 AF_ = AB->values();
570
571 return(0);
572}
573//=============================================================================
574
575template<typename OrdinalType, typename ScalarType>
578{
579 // Check that these new vectors are consistent.
580 TEUCHOS_TEST_FOR_EXCEPTION(B->numRows()!=X->numRows() || B->numCols() != X->numCols(), std::invalid_argument,
581 "SerialBandDenseSolver<T>::setVectors: X and B are not the same size!");
582 TEUCHOS_TEST_FOR_EXCEPTION(B->values()==0, std::invalid_argument,
583 "SerialBandDenseSolver<T>::setVectors: B is an empty SerialDenseMatrix<T>!");
584 TEUCHOS_TEST_FOR_EXCEPTION(X->values()==0, std::invalid_argument,
585 "SerialBandDenseSolver<T>::setVectors: X is an empty SerialDenseMatrix<T>!");
586 TEUCHOS_TEST_FOR_EXCEPTION(B->stride()<1, std::invalid_argument,
587 "SerialBandDenseSolver<T>::setVectors: B has an invalid stride!");
588 TEUCHOS_TEST_FOR_EXCEPTION(X->stride()<1, std::invalid_argument,
589 "SerialBandDenseSolver<T>::setVectors: X has an invalid stride!");
590
591 resetVectors();
592 LHS_ = X;
593 RHS_ = B;
594 return(0);
595}
596//=============================================================================
597
598template<typename OrdinalType, typename ScalarType>
600{
601 estimateSolutionErrors_ = flag;
602
603 // If the errors are estimated, this implies that the solution must be refined
604 refineSolution_ = refineSolution_ || flag;
605}
606//=============================================================================
607
608template<typename OrdinalType, typename ScalarType>
610
611 if (factored()) return(0); // Already factored
612
613 ANORM_ = Matrix_->normOne(); // Compute 1-Norm of A
614
615 // If we want to refine the solution, then the factor must
616 // be stored separatedly from the original matrix
617
618 if (A_ == AF_)
619 if (refineSolution_ ) {
620 Factor_ = rcp( new SerialBandDenseMatrix<OrdinalType,ScalarType>(*Matrix_) );
621 AF_ = Factor_->values();
622 LDAF_ = Factor_->stride();
623 }
624
625 int ierr = 0;
626 if (equilibrate_) ierr = equilibrateMatrix();
627
628 if (ierr!=0) return(ierr);
629
630 if (IPIV_.size()==0) IPIV_.resize( N_ ); // Allocated Pivot vector if not already done.
631
632 INFO_ = 0;
633
634#ifdef HAVE_TEUCHOSNUMERICS_EIGEN
635 EigenMatrixMap aMap( AF_+KL_, KL_+KU_+1, N_, EigenOuterStride(LDAF_) );
636 EigenMatrix fullMatrix(M_, N_);
637 for (OrdinalType j=0; j<N_; j++) {
638 for (OrdinalType i=TEUCHOS_MAX(0,j-KU_); i<=TEUCHOS_MIN(M_-1, j+KL_); i++) {
639 fullMatrix(i,j) = aMap(KU_+i-j, j);
640 }
641 }
642
645 EigenOrdinalArrayMap ipivMap( &IPIV_[0], Min_MN_ );
646
647 lu_.compute(fullMatrix);
648 p = lu_.permutationP();
649 ind = p.indices();
650
651 for (OrdinalType i=0; i<ipivMap.innerSize(); i++) {
652 ipivMap(i) = ind(i);
653 }
654
655#else
656 this->GBTRF(M_, N_, KL_, KU_, AF_, LDAF_, &IPIV_[0], &INFO_);
657#endif
658
659 factored_ = true;
660
661 return(INFO_);
662}
663
664//=============================================================================
665
666template<typename OrdinalType, typename ScalarType>
668
669 // If the user want the matrix to be equilibrated or wants a refined solution, we will
670 // call the X interface.
671 // Otherwise, if the matrix is already factored we will call the TRS interface.
672 // Otherwise, if the matrix is unfactored we will call the SV interface.
673
674 int ierr = 0;
675 if (equilibrate_) {
676 ierr = equilibrateRHS();
677 }
678 if (ierr != 0) return(ierr); // Can't equilibrate B, so return.
679
680 TEUCHOS_TEST_FOR_EXCEPTION( RHS_==Teuchos::null, std::invalid_argument,
681 "SerialBandDenseSolver<T>::solve: No right-hand side vector (RHS) has been set for the linear system!");
682 TEUCHOS_TEST_FOR_EXCEPTION( LHS_==Teuchos::null, std::invalid_argument,
683 "SerialBandDenseSolver<T>::solve: No solution vector (LHS) has been set for the linear system!");
684
685 if (!factored()) factor(); // Matrix must be factored
686
687 TEUCHOS_TEST_FOR_EXCEPTION( (equilibratedA_ && !equilibratedB_) || (!equilibratedA_ && equilibratedB_) ,
688 std::logic_error, "SerialBandDenseSolver<T>::solve: Matrix and vectors must be similarly scaled!");
689
690 if (shouldEquilibrate() && !equilibratedA_)
691 std::cout << "WARNING! SerialBandDenseSolver<T>::solve: System should be equilibrated!" << std::endl;
692
693 if (RHS_->values()!=LHS_->values()) {
694 (*LHS_) = (*RHS_); // Copy B to X if needed
695 }
696 INFO_ = 0;
697
698#ifdef HAVE_TEUCHOSNUMERICS_EIGEN
699 EigenMatrixMap rhsMap( RHS_->values(), RHS_->numRows(), RHS_->numCols(), EigenOuterStride(RHS_->stride()) );
700 EigenMatrixMap lhsMap( LHS_->values(), LHS_->numRows(), LHS_->numCols(), EigenOuterStride(LHS_->stride()) );
701 if ( TRANS_ == Teuchos::NO_TRANS ) {
702 lhsMap = lu_.solve(rhsMap);
703 } else if ( TRANS_ == Teuchos::TRANS ) {
704 lu_.matrixLU().template triangularView<Eigen::Upper>().transpose().solveInPlace(lhsMap);
705 lu_.matrixLU().template triangularView<Eigen::UnitLower>().transpose().solveInPlace(lhsMap);
706 EigenMatrix x = lu_.permutationP().transpose() * lhsMap;
707 lhsMap = x;
708 } else {
709 lu_.matrixLU().template triangularView<Eigen::Upper>().adjoint().solveInPlace(lhsMap);
710 lu_.matrixLU().template triangularView<Eigen::UnitLower>().adjoint().solveInPlace(lhsMap);
711 EigenMatrix x = lu_.permutationP().transpose() * lhsMap;
712 lhsMap = x;
713 }
714#else
715 this->GBTRS(ETranspChar[TRANS_], N_, KL_, KU_, RHS_->numCols(), AF_, LDAF_, &IPIV_[0], LHS_->values(), LHS_->stride(), &INFO_);
716#endif
717
718 if (INFO_!=0) return(INFO_);
719 solved_ = true;
720
721 int ierr1=0;
722 if (refineSolution_) ierr1 = applyRefinement();
723 if (ierr1!=0)
724 return(ierr1);
725
726 if (equilibrate_) ierr1 = unequilibrateLHS();
727
728 return(ierr1);
729}
730//=============================================================================
731
732template<typename OrdinalType, typename ScalarType>
734{
735 TEUCHOS_TEST_FOR_EXCEPTION(!solved(), std::logic_error,
736 "SerialBandDenseSolver<T>::applyRefinement: Must have an existing solution!");
737 TEUCHOS_TEST_FOR_EXCEPTION(A_==AF_, std::logic_error,
738 "SerialBandDenseSolver<T>::applyRefinement: Cannot apply refinement if no original copy of A!");
739
740#ifdef HAVE_TEUCHOSNUMERICS_EIGEN
741 // Implement templated GERFS or use Eigen.
742 return (-1);
743#else
744
745 OrdinalType NRHS = RHS_->numCols();
746 FERR_.resize( NRHS );
747 BERR_.resize( NRHS );
748 allocateWORK();
749
750 INFO_ = 0;
751 std::vector<typename details::lapack_traits<ScalarType>::iwork_type> GBRFS_WORK( N_ );
752 this->GBRFS(ETranspChar[TRANS_], N_, KL_, KU_, NRHS, A_+KL_, LDA_, AF_, LDAF_, &IPIV_[0],
753 RHS_->values(), RHS_->stride(), LHS_->values(), LHS_->stride(),
754 &FERR_[0], &BERR_[0], &WORK_[0], &GBRFS_WORK[0], &INFO_);
755
756 solutionErrorsEstimated_ = true;
757 reciprocalConditionEstimated_ = true;
758 solutionRefined_ = true;
759
760 return(INFO_);
761#endif
762}
763
764//=============================================================================
765
766template<typename OrdinalType, typename ScalarType>
768{
769 if (R_.size()!=0) return(0); // Already computed
770
771 R_.resize( M_ );
772 C_.resize( N_ );
773
774 INFO_ = 0;
775 this->GBEQU (M_, N_, KL_, KU_, AF_+KL_, LDAF_, &R_[0], &C_[0], &ROWCND_, &COLCND_, &AMAX_, &INFO_);
776
777 if (COLCND_<0.1*ScalarTraits<MagnitudeType>::one() ||
778 ROWCND_<0.1*ScalarTraits<MagnitudeType>::one() ||
780 shouldEquilibrate_ = true;
781
782 return(INFO_);
783}
784
785//=============================================================================
786
787template<typename OrdinalType, typename ScalarType>
789{
790 OrdinalType i, j;
791 int ierr = 0;
792
793 if (equilibratedA_) return(0); // Already done.
794 if (R_.size()==0) ierr = computeEquilibrateScaling(); // Compute R and C if needed.
795 if (ierr!=0) return(ierr); // If return value is not zero, then can't equilibrate.
796
797 if (A_==AF_) {
798
799 ScalarType * ptr;
800 for (j=0; j<N_; j++) {
801 ptr = A_ + KL_ + j*LDA_ + TEUCHOS_MAX(KU_-j, 0);
802 ScalarType s1 = C_[j];
803 for (i=TEUCHOS_MAX(0,j-KU_); i<=TEUCHOS_MIN(M_-1,j+KL_); i++) {
804 *ptr = *ptr*s1*R_[i];
805 ptr++;
806 }
807 }
808 } else {
809
810 ScalarType * ptr;
813 for (j=0; j<N_; j++) {
814 ptr = A_ + KL_ + j*LDA_ + TEUCHOS_MAX(KU_-j, 0);
815 ScalarType s1 = C_[j];
816 for (i=TEUCHOS_MAX(0,j-KU_); i<=TEUCHOS_MIN(M_-1,j+KL_); i++) {
817 *ptr = *ptr*s1*R_[i];
818 ptr++;
819 }
820 }
821 for (j=0; j<N_; j++) {
822 ptrU = AF_ + j*LDAF_ + TEUCHOS_MAX(KL_+KU_-j, 0);
823 ptrL = AF_ + KL_ + KU_ + 1 + j*LDAF_;
824 ScalarType s1 = C_[j];
825 for (i=TEUCHOS_MAX(0,j-(KL_+KU_)); i<=TEUCHOS_MIN(M_-1,j); i++) {
826 *ptrU = *ptrU*s1*R_[i];
827 ptrU++;
828 }
829 for (i=TEUCHOS_MAX(0,j); i<=TEUCHOS_MIN(M_-1,j+KL_); i++) {
830 *ptrL = *ptrL*s1*R_[i];
831 ptrL++;
832 }
833 }
834 }
835
836 equilibratedA_ = true;
837
838 return(0);
839}
840
841//=============================================================================
842
843template<typename OrdinalType, typename ScalarType>
845{
846 OrdinalType i, j;
847 int ierr = 0;
848
849 if (equilibratedB_) return(0); // Already done.
850 if (R_.size()==0) ierr = computeEquilibrateScaling(); // Compute R and C if needed.
851 if (ierr!=0) return(ierr); // Can't count on R and C being computed.
852
853 MagnitudeType * R_tmp = &R_[0];
854 if (transpose_) R_tmp = &C_[0];
855
856 OrdinalType LDB = RHS_->stride(), NRHS = RHS_->numCols();
857 ScalarType * B = RHS_->values();
858 ScalarType * ptr;
859 for (j=0; j<NRHS; j++) {
860 ptr = B + j*LDB;
861 for (i=0; i<M_; i++) {
862 *ptr = *ptr*R_tmp[i];
863 ptr++;
864 }
865 }
866
867 equilibratedB_ = true;
868
869 return(0);
870}
871
872
873//=============================================================================
874
875template<typename OrdinalType, typename ScalarType>
877{
878 OrdinalType i, j;
879
880 if (!equilibratedB_) return(0); // Nothing to do
881
882 MagnitudeType * C_tmp = &C_[0];
883 if (transpose_) C_tmp = &R_[0];
884
885 OrdinalType LDX = LHS_->stride(), NLHS = LHS_->numCols();
886 ScalarType * X = LHS_->values();
887 ScalarType * ptr;
888 for (j=0; j<NLHS; j++) {
889 ptr = X + j*LDX;
890 for (i=0; i<N_; i++) {
891 *ptr = *ptr*C_tmp[i];
892 ptr++;
893 }
894 }
895
896 return(0);
897}
898
899//=============================================================================
900
901template<typename OrdinalType, typename ScalarType>
903{
904#ifdef HAVE_TEUCHOSNUMERICS_EIGEN
905 // Implement templated GECON or use Eigen. Eigen currently doesn't have a condition estimation function.
906 return (-1);
907#else
908
909 if (reciprocalConditionEstimated()) {
910 Value = RCOND_;
911 return(0); // Already computed, just return it.
912 }
913
914 if ( ANORM_<ScalarTraits<MagnitudeType>::zero() ) ANORM_ = Matrix_->normOne();
915
916 int ierr = 0;
917 if (!factored()) ierr = factor(); // Need matrix factored.
918 if (ierr!=0) return(ierr);
919
920 allocateWORK();
921
922 // We will assume a one-norm condition number
923 INFO_ = 0;
924 std::vector<typename details::lapack_traits<ScalarType>::iwork_type> GBCON_WORK( N_ );
925 this->GBCON( '1', N_, KL_, KU_, AF_, LDAF_, &IPIV_[0], ANORM_, &RCOND_, &WORK_[0], &GBCON_WORK[0], &INFO_);
926 reciprocalConditionEstimated_ = true;
927 Value = RCOND_;
928
929 return(INFO_);
930#endif
931}
932//=============================================================================
933
934template<typename OrdinalType, typename ScalarType>
936
937 if (Matrix_!=Teuchos::null) os << "Solver Matrix" << std::endl << *Matrix_ << std::endl;
938 if (Factor_!=Teuchos::null) os << "Solver Factored Matrix" << std::endl << *Factor_ << std::endl;
939 if (LHS_ !=Teuchos::null) os << "Solver LHS" << std::endl << *LHS_ << std::endl;
940 if (RHS_ !=Teuchos::null) os << "Solver RHS" << std::endl << *RHS_ << std::endl;
941
942}
943
944//=============================================================================
945
946
947//=============================================================================
948
949
950} // namespace Teuchos
951
952#endif /* _TEUCHOS_SERIALBANDDENSESOLVER_HPP_ */
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.
#define TEUCHOS_MIN(x, y)
#define TEUCHOS_MAX(x, y)
Templated interface class to LAPACK routines.
Reference-counted pointer class and non-member templated function implementations.
Defines basic traits for the scalar field type.
Templated serial dense matrix class.
Templated serial dense matrix class.
Templated class for solving dense linear problems.
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.
A class for representing and solving banded dense linear systems.
bool equilibratedA()
Returns true if factor is equilibrated (factor available via AF() and LDAF()).
void estimateSolutionErrors(bool flag)
Causes all solves to estimate the forward and backward solution error.
MagnitudeType ANORM() const
Returns the 1-Norm of the this matrix (returns -1 if not yet computed).
bool shouldEquilibrate()
Returns true if the LAPACK general rules for equilibration suggest you should equilibrate the system.
SerialBandDenseSolver(const SerialBandDenseSolver< OrdinalType, ScalarType > &Source)
MagnitudeType ROWCND() const
Ratio of smallest to largest row scale factors for the this matrix (returns -1 if not yet computed).
RCP< SerialDenseMatrix< OrdinalType, ScalarType > > getRHS() const
Returns pointer to current RHS.
OrdinalType numRows() const
Returns row dimension of system.
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).
bool factored()
Returns true if matrix is factored (factor available via AF() and LDAF()).
RCP< SerialDenseMatrix< OrdinalType, ScalarType > > RHS_
std::vector< OrdinalType > IPIV() const
Returns pointer to pivot vector (if factorization has been computed), zero otherwise.
int applyRefinement()
Apply Iterative Refinement.
bool solutionErrorsEstimated()
Returns true if forward and backward error estimated have been computed (available via FERR() and BER...
OrdinalType numLower() const
Returns lower bandwidth of system.
bool solved()
Returns true if the current set of vectors has been solved.
MagnitudeType RCOND() const
Returns the reciprocal of the condition number of the this matrix (returns -1 if not yet computed).
OrdinalType numUpper() const
Returns upper bandwidth of system.
virtual ~SerialBandDenseSolver()
SerialBandDenseSolver destructor.
RCP< SerialDenseMatrix< OrdinalType, ScalarType > > LHS_
int setMatrix(const RCP< SerialBandDenseMatrix< OrdinalType, ScalarType > > &AB)
Sets the pointer for coefficient matrix.
RCP< SerialDenseMatrix< OrdinalType, ScalarType > > getLHS() const
Returns pointer to current LHS.
std::vector< MagnitudeType > FERR() const
Returns a pointer to the forward error estimates computed by LAPACK.
int factor()
Computes the in-place LU factorization of the matrix using the LAPACK routine _GBTRF.
std::vector< MagnitudeType > BERR() const
Returns a pointer to the backward error estimates computed by LAPACK.
SerialBandDenseSolver & operator=(const SerialBandDenseSolver< OrdinalType, ScalarType > &Source)
std::vector< MagnitudeType > R() const
Returns a pointer to the row scaling vector used for equilibration.
RCP< SerialBandDenseMatrix< OrdinalType, ScalarType > > Matrix_
bool transpose()
Returns true if transpose of this matrix has and will be used.
Teuchos::ScalarTraits< ScalarType >::magnitudeType MagnitudeType
MagnitudeType AMAX() const
Returns the absolute value of the largest entry of the this matrix (returns -1 if not yet computed).
std::vector< MagnitudeType > C() const
Returns a pointer to the column scale vector used for equilibration.
MagnitudeType COLCND() const
Ratio of smallest to largest column scale factors for the this matrix (returns -1 if not yet computed...
RCP< SerialBandDenseMatrix< OrdinalType, ScalarType > > Factor_
int computeEquilibrateScaling()
Computes the scaling vector S(i) = 1/sqrt(A(i,i)) of the this matrix.
RCP< SerialBandDenseMatrix< OrdinalType, ScalarType > > getFactoredMatrix() const
Returns pointer to factored matrix (assuming factorization has been performed).
bool reciprocalConditionEstimated()
Returns true if the condition number of the this matrix has been computed (value available via Recipr...
int reciprocalConditionEstimate(MagnitudeType &Value)
Returns the reciprocal of the 1-norm condition number of the this matrix.
int unequilibrateLHS()
Unscales the solution vectors if equilibration was used to solve the system.
bool equilibratedB()
Returns true if RHS is equilibrated (RHS available via B() and LDB()).
void Print(std::ostream &os) const
Print service methods; defines behavior of ostream << operator.
int equilibrateMatrix()
Equilibrates the this matrix.
bool solutionRefined()
Returns true if the current set of vectors has been refined.
OrdinalType numCols() const
Returns column dimension of system.
int equilibrateRHS()
Equilibrates the current RHS.
void solveWithTransposeFlag(Teuchos::ETransp trans)
RCP< SerialBandDenseMatrix< OrdinalType, ScalarType > > getMatrix() const
Returns pointer to current matrix.
int solve()
Computes the solution X to AX = B for the this matrix and the B provided to SetVectors()....
void solveToRefinedSolution(bool flag)
Set whether or not to use iterative refinement to improve solutions to linear systems.
SerialBandDenseSolver()
Default constructor; matrix should be set using setMatrix(), LHS and RHS set with setVectors().
Concrete serial communicator subclass.
#define TEUCHOS_TEST_FOR_EXCEPTION(throw_exception_test, Exception, msg)
Macro for throwing an exception with breakpointing to ease debugging.
Definition PackageB.cpp:3
TEUCHOSNUMERICS_LIB_DLL_EXPORT const char ETranspChar[]
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.