Teuchos Package Browser (Single Doxygen Collection) Version of the Day
Teuchos_SerialBandDenseSolver.hpp
Go to the documentation of this file.
00001 // @HEADER
00002 // ***********************************************************************
00003 //
00004 //                    Teuchos: Common Tools Package
00005 //                 Copyright (2004) Sandia Corporation
00006 //
00007 // Under terms of Contract DE-AC04-94AL85000, there is a non-exclusive
00008 // license for use of this work by or on behalf of the U.S. Government.
00009 //
00010 // Redistribution and use in source and binary forms, with or without
00011 // modification, are permitted provided that the following conditions are
00012 // met:
00013 //
00014 // 1. Redistributions of source code must retain the above copyright
00015 // notice, this list of conditions and the following disclaimer.
00016 //
00017 // 2. Redistributions in binary form must reproduce the above copyright
00018 // notice, this list of conditions and the following disclaimer in the
00019 // documentation and/or other materials provided with the distribution.
00020 //
00021 // 3. Neither the name of the Corporation nor the names of the
00022 // contributors may be used to endorse or promote products derived from
00023 // this software without specific prior written permission.
00024 //
00025 // THIS SOFTWARE IS PROVIDED BY SANDIA CORPORATION "AS IS" AND ANY
00026 // EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
00027 // IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR
00028 // PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL SANDIA CORPORATION OR THE
00029 // CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL,
00030 // EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO,
00031 // PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR
00032 // PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF
00033 // LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING
00034 // NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
00035 // SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
00036 //
00037 // Questions? Contact Michael A. Heroux (maherou@sandia.gov)
00038 //
00039 // ***********************************************************************
00040 // @HEADER
00041 
00042 #ifndef _TEUCHOS_SERIALBANDDENSESOLVER_HPP_
00043 #define _TEUCHOS_SERIALBANDDENSESOLVER_HPP_
00044 
00045 
00046 
00047 
00048 
00049 #include "Teuchos_CompObject.hpp"
00050 #include "Teuchos_BLAS.hpp"
00051 #include "Teuchos_LAPACK.hpp"
00052 #include "Teuchos_RCP.hpp"
00053 #include "Teuchos_ConfigDefs.hpp"
00054 #include "Teuchos_SerialDenseMatrix.hpp"
00055 #include "Teuchos_SerialBandDenseMatrix.hpp"
00056 #include "Teuchos_SerialDenseSolver.hpp"
00057 #include "Teuchos_ScalarTraits.hpp"
00058 
00059 #ifdef HAVE_TEUCHOSNUMERICS_EIGEN
00060 #include "Eigen/Dense"
00061 #endif
00062 
00063 namespace Teuchos {
00064 
00165   template<typename OrdinalType, typename ScalarType>
00166   class SerialBandDenseSolver : public CompObject, public Object, public BLAS<OrdinalType, ScalarType>,
00167         public LAPACK<OrdinalType, ScalarType>
00168   {
00169 
00170   public:
00171 
00172     typedef typename Teuchos::ScalarTraits<ScalarType>::magnitudeType MagnitudeType;
00173 #ifdef HAVE_TEUCHOSNUMERICS_EIGEN
00174     typedef typename Eigen::Matrix<ScalarType,Eigen::Dynamic,Eigen::Dynamic,Eigen::ColMajor> EigenMatrix;
00175     typedef typename Eigen::internal::BandMatrix<ScalarType,Eigen::Dynamic,Eigen::Dynamic,Eigen::ColMajor> EigenBandMatrix;
00176     typedef typename Eigen::Matrix<ScalarType,Eigen::Dynamic,1> EigenVector;
00177     typedef typename Eigen::InnerStride<Eigen::Dynamic> EigenInnerStride;
00178     typedef typename Eigen::OuterStride<Eigen::Dynamic> EigenOuterStride;
00179     typedef typename Eigen::Map<EigenVector,0,EigenInnerStride > EigenVectorMap;
00180     typedef typename Eigen::Map<const EigenVector,0,EigenInnerStride > EigenConstVectorMap;
00181     typedef typename Eigen::Map<EigenMatrix,0,EigenOuterStride > EigenMatrixMap;
00182     typedef typename Eigen::Map<const EigenMatrix,0,EigenOuterStride > EigenConstMatrixMap;
00183     typedef typename Eigen::PermutationMatrix<Eigen::Dynamic,Eigen::Dynamic,OrdinalType> EigenPermutationMatrix;
00184     typedef typename Eigen::Array<OrdinalType,Eigen::Dynamic,1> EigenOrdinalArray;
00185     typedef typename Eigen::Map<EigenOrdinalArray> EigenOrdinalArrayMap;
00186     typedef typename Eigen::Array<ScalarType,Eigen::Dynamic,1> EigenScalarArray;
00187     typedef typename Eigen::Map<EigenScalarArray> EigenScalarArrayMap;
00188 #endif
00189 
00191 
00192 
00194     SerialBandDenseSolver();
00195 
00197     virtual ~SerialBandDenseSolver();
00198 
00200 
00201 
00202 
00204     int setMatrix(const RCP<SerialBandDenseMatrix<OrdinalType, ScalarType> >& AB);
00205 
00207 
00210     int setVectors(const RCP<SerialDenseMatrix<OrdinalType, ScalarType> >& X,
00211        const RCP<SerialDenseMatrix<OrdinalType, ScalarType> >& B);
00212 
00214 
00215 
00216 
00219     void factorWithEquilibration(bool flag) {equilibrate_ = flag; return;}
00220 
00224     void solveWithTranspose(bool flag) {transpose_ = flag; if (flag) TRANS_ = Teuchos::TRANS; else TRANS_ = Teuchos::NO_TRANS; return;}
00225 
00228     void solveWithTransposeFlag(Teuchos::ETransp trans) {TRANS_ = trans; if (trans != Teuchos::NO_TRANS) {  transpose_ = true; } }
00229 
00231     void solveToRefinedSolution(bool flag) {refineSolution_ = flag; return;}
00232 
00234 
00237     void estimateSolutionErrors(bool flag);
00239 
00241 
00242 
00244 
00247     int factor();
00248 
00250 
00253     int solve();
00254 
00256 
00259     int computeEquilibrateScaling();
00260 
00262 
00265     int equilibrateMatrix();
00266 
00268 
00271     int equilibrateRHS();
00272 
00274 
00277     int applyRefinement();
00278 
00280 
00283     int unequilibrateLHS();
00284 
00286 
00292      int reciprocalConditionEstimate(MagnitudeType & Value);
00294 
00296 
00297 
00299     bool transpose() {return(transpose_);}
00300 
00302     bool factored() {return(factored_);}
00303 
00305     bool equilibratedA() {return(equilibratedA_);}
00306 
00308     bool equilibratedB() {return(equilibratedB_);}
00309 
00311      bool shouldEquilibrate() {computeEquilibrateScaling(); return(shouldEquilibrate_);}
00312 
00314     bool solutionErrorsEstimated() {return(solutionErrorsEstimated_);}
00315 
00317     bool reciprocalConditionEstimated() {return(reciprocalConditionEstimated_);}
00318 
00320     bool solved() {return(solved_);}
00321 
00323     bool solutionRefined() {return(solutionRefined_);}
00325 
00327 
00328 
00330     RCP<SerialBandDenseMatrix<OrdinalType, ScalarType> > getMatrix()  const {return(Matrix_);}
00331 
00333     RCP<SerialBandDenseMatrix<OrdinalType, ScalarType> > getFactoredMatrix()  const {return(Factor_);}
00334 
00336     RCP<SerialDenseMatrix<OrdinalType, ScalarType> > getLHS() const {return(LHS_);}
00337 
00339     RCP<SerialDenseMatrix<OrdinalType, ScalarType> > getRHS() const {return(RHS_);}
00340 
00342     OrdinalType numRows()  const {return(M_);}
00343 
00345     OrdinalType numCols()  const {return(N_);}
00346 
00348     OrdinalType numLower()  const {return(KL_);}
00349 
00351     OrdinalType numUpper()  const {return(KU_);}
00352 
00354     std::vector<OrdinalType> IPIV()  const {return(IPIV_);}
00355 
00357     MagnitudeType ANORM()  const {return(ANORM_);}
00358 
00360     MagnitudeType RCOND()  const {return(RCOND_);}
00361 
00363 
00365     MagnitudeType ROWCND()  const {return(ROWCND_);}
00366 
00368 
00370     MagnitudeType COLCND()  const {return(COLCND_);}
00371 
00373     MagnitudeType AMAX()  const {return(AMAX_);}
00374 
00376     std::vector<MagnitudeType> FERR()  const {return(FERR_);}
00377 
00379     std::vector<MagnitudeType> BERR()  const {return(BERR_);}
00380 
00382     std::vector<MagnitudeType> R()  const {return(R_);}
00383 
00385     std::vector<MagnitudeType> C()  const {return(C_);}
00387 
00389 
00390 
00391     void Print(std::ostream& os) const;
00393   protected:
00394 
00395     void allocateWORK() { LWORK_ = 3*N_; WORK_.resize( LWORK_ ); return;}
00396     void resetMatrix();
00397     void resetVectors();
00398 
00399 
00400     bool equilibrate_;
00401     bool shouldEquilibrate_;
00402     bool equilibratedA_;
00403     bool equilibratedB_;
00404     bool transpose_;
00405     bool factored_;
00406     bool estimateSolutionErrors_;
00407     bool solutionErrorsEstimated_;
00408     bool solved_;
00409     bool reciprocalConditionEstimated_;
00410     bool refineSolution_;
00411     bool solutionRefined_;
00412 
00413     Teuchos::ETransp TRANS_;
00414 
00415     OrdinalType M_;
00416     OrdinalType N_;
00417     OrdinalType KL_;
00418     OrdinalType KU_;
00419     OrdinalType Min_MN_;
00420     OrdinalType LDA_;
00421     OrdinalType LDAF_;
00422     OrdinalType INFO_;
00423     OrdinalType LWORK_;
00424 
00425     std::vector<OrdinalType> IPIV_;
00426     std::vector<int> IWORK_;
00427 
00428     MagnitudeType ANORM_;
00429     MagnitudeType RCOND_;
00430     MagnitudeType ROWCND_;
00431     MagnitudeType COLCND_;
00432     MagnitudeType AMAX_;
00433 
00434     RCP<SerialBandDenseMatrix<OrdinalType, ScalarType> > Matrix_;
00435     RCP<SerialDenseMatrix<OrdinalType, ScalarType> > LHS_;
00436     RCP<SerialDenseMatrix<OrdinalType, ScalarType> > RHS_;
00437     RCP<SerialBandDenseMatrix<OrdinalType, ScalarType> > Factor_;
00438 
00439     ScalarType * A_;
00440     ScalarType * AF_;
00441     std::vector<MagnitudeType> FERR_;
00442     std::vector<MagnitudeType> BERR_;
00443     std::vector<ScalarType> WORK_;
00444     std::vector<MagnitudeType> R_;
00445     std::vector<MagnitudeType> C_;
00446 #ifdef HAVE_TEUCHOSNUMERICS_EIGEN
00447     Eigen::PartialPivLU<EigenMatrix> lu_;
00448 #endif
00449 
00450   private:
00451 
00452     // SerialBandDenseSolver copy constructor (put here because we don't want user access)
00453     SerialBandDenseSolver(const SerialBandDenseSolver<OrdinalType, ScalarType>& Source);
00454     SerialBandDenseSolver & operator=(const SerialBandDenseSolver<OrdinalType, ScalarType>& Source);
00455 
00456   };
00457 
00458   // Helper traits to distinguish work arrays for real and complex-valued datatypes.
00459   using namespace details;
00460 
00461 //=============================================================================
00462 
00463 template<typename OrdinalType, typename ScalarType>
00464 SerialBandDenseSolver<OrdinalType,ScalarType>::SerialBandDenseSolver()
00465   : CompObject(),
00466     Object("Teuchos::SerialBandDenseSolver"),
00467     equilibrate_(false),
00468     shouldEquilibrate_(false),
00469     equilibratedA_(false),
00470     equilibratedB_(false),
00471     transpose_(false),
00472     factored_(false),
00473     estimateSolutionErrors_(false),
00474     solutionErrorsEstimated_(false),
00475     solved_(false),
00476     reciprocalConditionEstimated_(false),
00477     refineSolution_(false),
00478     solutionRefined_(false),
00479     TRANS_(Teuchos::NO_TRANS),
00480     M_(0),
00481     N_(0),
00482     KL_(0),
00483     KU_(0),
00484     Min_MN_(0),
00485     LDA_(0),
00486     LDAF_(0),
00487     INFO_(0),
00488     LWORK_(0),
00489     RCOND_(ScalarTraits<MagnitudeType>::zero()),
00490     ROWCND_(ScalarTraits<MagnitudeType>::zero()),
00491     COLCND_(ScalarTraits<MagnitudeType>::zero()),
00492     AMAX_(ScalarTraits<MagnitudeType>::zero()),
00493     A_(0),
00494     AF_(0)
00495 {
00496   resetMatrix();
00497 }
00498 //=============================================================================
00499 
00500 template<typename OrdinalType, typename ScalarType>
00501 SerialBandDenseSolver<OrdinalType,ScalarType>::~SerialBandDenseSolver()
00502 {}
00503 
00504 //=============================================================================
00505 
00506 template<typename OrdinalType, typename ScalarType>
00507 void SerialBandDenseSolver<OrdinalType,ScalarType>::resetVectors()
00508 {
00509   LHS_ = Teuchos::null;
00510   RHS_ = Teuchos::null;
00511   reciprocalConditionEstimated_ = false;
00512   solutionRefined_ = false;
00513   solved_ = false;
00514   solutionErrorsEstimated_ = false;
00515   equilibratedB_ = false;
00516 }
00517 //=============================================================================
00518 
00519 template<typename OrdinalType, typename ScalarType>
00520 void SerialBandDenseSolver<OrdinalType,ScalarType>::resetMatrix()
00521 {
00522   resetVectors();
00523   equilibratedA_ = false;
00524   factored_ = false;
00525   M_ = 0;
00526   N_ = 0;
00527   KL_ = 0;
00528   KU_ = 0;
00529   Min_MN_ = 0;
00530   LDA_ = 0;
00531   LDAF_ = 0;
00532   RCOND_ = -ScalarTraits<MagnitudeType>::one();
00533   ROWCND_ = -ScalarTraits<MagnitudeType>::one();
00534   COLCND_ = -ScalarTraits<MagnitudeType>::one();
00535   AMAX_ = -ScalarTraits<MagnitudeType>::one();
00536   A_ = 0;
00537   AF_ = 0;
00538   INFO_ = 0;
00539   LWORK_ = 0;
00540   R_.resize(0);
00541   C_.resize(0);
00542 }
00543 //=============================================================================
00544 
00545 template<typename OrdinalType, typename ScalarType>
00546 int SerialBandDenseSolver<OrdinalType,ScalarType>::setMatrix(const RCP<SerialBandDenseMatrix<OrdinalType,ScalarType> >& AB)
00547 {
00548 
00549   OrdinalType m = AB->numRows();
00550   OrdinalType n = AB->numCols();
00551   OrdinalType kl = AB->lowerBandwidth();
00552   OrdinalType ku = AB->upperBandwidth();
00553 
00554   // Check that the new matrix is consistent.
00555   TEUCHOS_TEST_FOR_EXCEPTION(AB->values()==0, std::invalid_argument,
00556          "SerialBandDenseSolver<T>::setMatrix: A is an empty SerialBandDenseMatrix<T>!");
00557 
00558   resetMatrix();
00559   Matrix_ = AB;
00560   Factor_ = AB;
00561   M_ = m;
00562   N_ = n;
00563   KL_ = kl;
00564   KU_ = ku-kl;
00565   Min_MN_ = TEUCHOS_MIN(M_,N_);
00566   LDA_ = AB->stride();
00567   LDAF_ = LDA_;
00568   A_ = AB->values();
00569   AF_ = AB->values();
00570 
00571   return(0);
00572 }
00573 //=============================================================================
00574 
00575 template<typename OrdinalType, typename ScalarType>
00576 int SerialBandDenseSolver<OrdinalType,ScalarType>::setVectors(const RCP<SerialDenseMatrix<OrdinalType,ScalarType> >& X,
00577                  const RCP<SerialDenseMatrix<OrdinalType,ScalarType> >& B)
00578 {
00579   // Check that these new vectors are consistent.
00580   TEUCHOS_TEST_FOR_EXCEPTION(B->numRows()!=X->numRows() || B->numCols() != X->numCols(), std::invalid_argument,
00581          "SerialBandDenseSolver<T>::setVectors: X and B are not the same size!");
00582   TEUCHOS_TEST_FOR_EXCEPTION(B->values()==0, std::invalid_argument,
00583          "SerialBandDenseSolver<T>::setVectors: B is an empty SerialDenseMatrix<T>!");
00584   TEUCHOS_TEST_FOR_EXCEPTION(X->values()==0, std::invalid_argument,
00585          "SerialBandDenseSolver<T>::setVectors: X is an empty SerialDenseMatrix<T>!");
00586   TEUCHOS_TEST_FOR_EXCEPTION(B->stride()<1, std::invalid_argument,
00587          "SerialBandDenseSolver<T>::setVectors: B has an invalid stride!");
00588   TEUCHOS_TEST_FOR_EXCEPTION(X->stride()<1, std::invalid_argument,
00589          "SerialBandDenseSolver<T>::setVectors: X has an invalid stride!");
00590 
00591   resetVectors();
00592   LHS_ = X;
00593   RHS_ = B;
00594   return(0);
00595 }
00596 //=============================================================================
00597 
00598 template<typename OrdinalType, typename ScalarType>
00599 void SerialBandDenseSolver<OrdinalType,ScalarType>::estimateSolutionErrors(bool flag)
00600 {
00601   estimateSolutionErrors_ = flag;
00602 
00603   // If the errors are estimated, this implies that the solution must be refined
00604   refineSolution_ = refineSolution_ || flag;
00605 }
00606 //=============================================================================
00607 
00608 template<typename OrdinalType, typename ScalarType>
00609 int SerialBandDenseSolver<OrdinalType,ScalarType>::factor() {
00610 
00611   if (factored()) return(0); // Already factored
00612 
00613   ANORM_ = Matrix_->normOne(); // Compute 1-Norm of A
00614 
00615   // If we want to refine the solution, then the factor must
00616   // be stored separatedly from the original matrix
00617 
00618   if (A_ == AF_)
00619     if (refineSolution_ ) {
00620       Factor_ = rcp( new SerialBandDenseMatrix<OrdinalType,ScalarType>(*Matrix_) );
00621       AF_ = Factor_->values();
00622       LDAF_ = Factor_->stride();
00623     }
00624 
00625   int ierr = 0;
00626   if (equilibrate_) ierr = equilibrateMatrix();
00627 
00628   if (ierr!=0) return(ierr);
00629 
00630   if (IPIV_.size()==0) IPIV_.resize( N_ ); // Allocated Pivot vector if not already done.
00631 
00632   INFO_ = 0;
00633 
00634 #ifdef HAVE_TEUCHOSNUMERICS_EIGEN
00635   EigenMatrixMap aMap( AF_+KL_, KL_+KU_+1, N_, EigenOuterStride(LDAF_) );
00636   EigenMatrix fullMatrix(M_, N_);
00637   for (OrdinalType j=0; j<N_; j++) {
00638     for (OrdinalType i=TEUCHOS_MAX(0,j-KU_); i<=TEUCHOS_MIN(M_-1, j+KL_); i++) {
00639       fullMatrix(i,j) = aMap(KU_+i-j, j);
00640     }
00641   }
00642 
00643   EigenPermutationMatrix p;
00644   EigenOrdinalArray ind;
00645   EigenOrdinalArrayMap ipivMap( &IPIV_[0], Min_MN_ );
00646 
00647   lu_.compute(fullMatrix);
00648   p = lu_.permutationP();
00649   ind = p.indices();
00650 
00651   for (OrdinalType i=0; i<ipivMap.innerSize(); i++) {
00652     ipivMap(i) = ind(i);
00653   }
00654 
00655 #else
00656   this->GBTRF(M_, N_, KL_, KU_, AF_, LDAF_, &IPIV_[0], &INFO_);
00657 #endif
00658 
00659   factored_ = true;
00660 
00661   return(INFO_);
00662 }
00663 
00664 //=============================================================================
00665 
00666 template<typename OrdinalType, typename ScalarType>
00667 int SerialBandDenseSolver<OrdinalType,ScalarType>::solve() {
00668 
00669   // If the user want the matrix to be equilibrated or wants a refined solution, we will
00670   // call the X interface.
00671   // Otherwise, if the matrix is already factored we will call the TRS interface.
00672   // Otherwise, if the matrix is unfactored we will call the SV interface.
00673 
00674   int ierr = 0;
00675   if (equilibrate_) {
00676     ierr = equilibrateRHS();
00677     equilibratedB_ = true;
00678   }
00679   if (ierr != 0) return(ierr);  // Can't equilibrate B, so return.
00680 
00681   TEUCHOS_TEST_FOR_EXCEPTION( (equilibratedA_ && !equilibratedB_) || (!equilibratedA_ && equilibratedB_) ,
00682          std::logic_error, "SerialBandDenseSolver<T>::solve: Matrix and vectors must be similarly scaled!");
00683   TEUCHOS_TEST_FOR_EXCEPTION( RHS_==Teuchos::null, std::invalid_argument,
00684          "SerialBandDenseSolver<T>::solve: No right-hand side vector (RHS) has been set for the linear system!");
00685   TEUCHOS_TEST_FOR_EXCEPTION( LHS_==Teuchos::null, std::invalid_argument,
00686          "SerialBandDenseSolver<T>::solve: No solution vector (LHS) has been set for the linear system!");
00687 
00688   if (shouldEquilibrate() && !equilibratedA_)
00689     std::cout << "WARNING!  SerialBandDenseSolver<T>::solve: System should be equilibrated!" << std::endl;
00690 
00691   if (!factored()) factor(); // Matrix must be factored
00692 
00693   if (RHS_->values()!=LHS_->values()) {
00694     (*LHS_) = (*RHS_); // Copy B to X if needed
00695   }
00696   INFO_ = 0;
00697 
00698 #ifdef HAVE_TEUCHOSNUMERICS_EIGEN
00699     EigenMatrixMap rhsMap( RHS_->values(), RHS_->numRows(), RHS_->numCols(), EigenOuterStride(RHS_->stride()) );
00700     EigenMatrixMap lhsMap( LHS_->values(), LHS_->numRows(), LHS_->numCols(), EigenOuterStride(LHS_->stride()) );
00701     if ( TRANS_ == Teuchos::NO_TRANS ) {
00702       lhsMap = lu_.solve(rhsMap);
00703     } else if ( TRANS_ == Teuchos::TRANS ) {
00704       lu_.matrixLU().template triangularView<Eigen::Upper>().transpose().solveInPlace(lhsMap);
00705       lu_.matrixLU().template triangularView<Eigen::UnitLower>().transpose().solveInPlace(lhsMap);
00706       EigenMatrix x = lu_.permutationP().transpose() * lhsMap;
00707       lhsMap = x;
00708     } else {
00709       lu_.matrixLU().template triangularView<Eigen::Upper>().adjoint().solveInPlace(lhsMap);
00710       lu_.matrixLU().template triangularView<Eigen::UnitLower>().adjoint().solveInPlace(lhsMap);
00711       EigenMatrix x = lu_.permutationP().transpose() * lhsMap;
00712       lhsMap = x;
00713     }
00714 #else
00715   this->GBTRS(ETranspChar[TRANS_], N_, KL_, KU_, RHS_->numCols(), AF_, LDAF_, &IPIV_[0], LHS_->values(), LHS_->stride(), &INFO_);
00716 #endif
00717 
00718   if (INFO_!=0) return(INFO_);
00719   solved_ = true;
00720 
00721   int ierr1=0;
00722   if (refineSolution_) ierr1 = applyRefinement();
00723   if (ierr1!=0)
00724     return(ierr1);
00725 
00726   if (equilibrate_) ierr1 = unequilibrateLHS();
00727 
00728   return(ierr1);
00729 }
00730 //=============================================================================
00731 
00732 template<typename OrdinalType, typename ScalarType>
00733 int SerialBandDenseSolver<OrdinalType,ScalarType>::applyRefinement()
00734 {
00735   TEUCHOS_TEST_FOR_EXCEPTION(!solved(), std::logic_error,
00736          "SerialBandDenseSolver<T>::applyRefinement: Must have an existing solution!");
00737   TEUCHOS_TEST_FOR_EXCEPTION(A_==AF_, std::logic_error,
00738          "SerialBandDenseSolver<T>::applyRefinement: Cannot apply refinement if no original copy of A!");
00739 
00740 #ifdef HAVE_TEUCHOSNUMERICS_EIGEN
00741   // Implement templated GERFS or use Eigen.
00742   return (-1);
00743 #else
00744 
00745   OrdinalType NRHS = RHS_->numCols();
00746   FERR_.resize( NRHS );
00747   BERR_.resize( NRHS );
00748   allocateWORK();
00749 
00750   INFO_ = 0;
00751   std::vector<typename details::lapack_traits<ScalarType>::iwork_type> GBRFS_WORK( N_ );
00752   this->GBRFS(ETranspChar[TRANS_], N_, KL_, KU_, NRHS, A_+KL_, LDA_, AF_, LDAF_, &IPIV_[0],
00753         RHS_->values(), RHS_->stride(), LHS_->values(), LHS_->stride(),
00754         &FERR_[0], &BERR_[0], &WORK_[0], &GBRFS_WORK[0], &INFO_);
00755 
00756   solutionErrorsEstimated_ = true;
00757   reciprocalConditionEstimated_ = true;
00758   solutionRefined_ = true;
00759 
00760   return(INFO_);
00761 #endif
00762 }
00763 
00764 //=============================================================================
00765 
00766 template<typename OrdinalType, typename ScalarType>
00767 int SerialBandDenseSolver<OrdinalType,ScalarType>::computeEquilibrateScaling()
00768 {
00769   if (R_.size()!=0) return(0); // Already computed
00770 
00771   R_.resize( M_ );
00772   C_.resize( N_ );
00773 
00774   INFO_ = 0;
00775   this->GBEQU (M_, N_, KL_, KU_, AF_+KL_, LDAF_, &R_[0], &C_[0], &ROWCND_, &COLCND_, &AMAX_, &INFO_);
00776 
00777   if (COLCND_<0.1*ScalarTraits<MagnitudeType>::one() ||
00778       ROWCND_<0.1*ScalarTraits<MagnitudeType>::one() ||
00779       AMAX_ < ScalarTraits<ScalarType>::rmin() || AMAX_ > ScalarTraits<ScalarType>::rmax() )
00780     shouldEquilibrate_ = true;
00781 
00782   return(INFO_);
00783 }
00784 
00785 //=============================================================================
00786 
00787 template<typename OrdinalType, typename ScalarType>
00788 int SerialBandDenseSolver<OrdinalType,ScalarType>::equilibrateMatrix()
00789 {
00790   OrdinalType i, j;
00791   int ierr = 0;
00792 
00793   if (equilibratedA_) return(0); // Already done.
00794   if (R_.size()==0) ierr = computeEquilibrateScaling(); // Compute R and C if needed.
00795   if (ierr!=0) return(ierr);     // If return value is not zero, then can't equilibrate.
00796 
00797   if (A_==AF_) {
00798 
00799     ScalarType * ptr;
00800     for (j=0; j<N_; j++) {
00801       ptr = A_ + KL_ + j*LDA_ + TEUCHOS_MAX(KU_-j, 0);
00802       ScalarType s1 = C_[j];
00803       for (i=TEUCHOS_MAX(0,j-KU_); i<=TEUCHOS_MIN(M_-1,j+KL_); i++) {
00804   *ptr = *ptr*s1*R_[i];
00805   ptr++;
00806       }
00807     }
00808   } else {
00809 
00810     ScalarType * ptr;
00811     ScalarType * ptrL;
00812     ScalarType * ptrU;
00813     for (j=0; j<N_; j++) {
00814       ptr = A_ + KL_ + j*LDA_ + TEUCHOS_MAX(KU_-j, 0);
00815       ScalarType s1 = C_[j];
00816       for (i=TEUCHOS_MAX(0,j-KU_); i<=TEUCHOS_MIN(M_-1,j+KL_); i++) {
00817   *ptr = *ptr*s1*R_[i];
00818   ptr++;
00819       }
00820     }
00821     for (j=0; j<N_; j++) {
00822       ptrU = AF_ + j*LDAF_ + TEUCHOS_MAX(KL_+KU_-j, 0);
00823       ptrL = AF_ + KL_ + KU_ + 1 + j*LDAF_;
00824       ScalarType s1 = C_[j];
00825       for (i=TEUCHOS_MAX(0,j-(KL_+KU_)); i<=TEUCHOS_MIN(M_-1,j); i++) {
00826   *ptrU = *ptrU*s1*R_[i];
00827   ptrU++;
00828       }
00829       for (i=TEUCHOS_MAX(0,j); i<=TEUCHOS_MIN(M_-1,j+KL_); i++) {
00830   *ptrL = *ptrL*s1*R_[i];
00831   ptrL++;
00832       }
00833     }
00834   }
00835 
00836   equilibratedA_ = true;
00837 
00838   return(0);
00839 }
00840 
00841 //=============================================================================
00842 
00843 template<typename OrdinalType, typename ScalarType>
00844 int SerialBandDenseSolver<OrdinalType,ScalarType>::equilibrateRHS()
00845 {
00846   OrdinalType i, j;
00847   int ierr = 0;
00848 
00849   if (equilibratedB_) return(0); // Already done.
00850   if (R_.size()==0) ierr = computeEquilibrateScaling(); // Compute R and C if needed.
00851   if (ierr!=0) return(ierr);     // Can't count on R and C being computed.
00852 
00853   MagnitudeType * R_tmp = &R_[0];
00854   if (transpose_) R_tmp = &C_[0];
00855 
00856   OrdinalType LDB = RHS_->stride(), NRHS = RHS_->numCols();
00857   ScalarType * B = RHS_->values();
00858   ScalarType * ptr;
00859   for (j=0; j<NRHS; j++) {
00860     ptr = B + j*LDB;
00861     for (i=0; i<M_; i++) {
00862       *ptr = *ptr*R_tmp[i];
00863       ptr++;
00864     }
00865   }
00866 
00867   equilibratedB_ = true;
00868 
00869   return(0);
00870 }
00871 
00872 
00873 //=============================================================================
00874 
00875 template<typename OrdinalType, typename ScalarType>
00876 int SerialBandDenseSolver<OrdinalType,ScalarType>::unequilibrateLHS()
00877 {
00878   OrdinalType i, j;
00879 
00880   if (!equilibratedB_) return(0); // Nothing to do
00881 
00882   MagnitudeType * C_tmp = &C_[0];
00883   if (transpose_) C_tmp = &R_[0];
00884 
00885   OrdinalType LDX = LHS_->stride(), NLHS = LHS_->numCols();
00886   ScalarType * X = LHS_->values();
00887   ScalarType * ptr;
00888   for (j=0; j<NLHS; j++) {
00889     ptr = X + j*LDX;
00890     for (i=0; i<N_; i++) {
00891       *ptr = *ptr*C_tmp[i];
00892       ptr++;
00893     }
00894   }
00895 
00896   return(0);
00897 }
00898 
00899 //=============================================================================
00900 
00901 template<typename OrdinalType, typename ScalarType>
00902 int SerialBandDenseSolver<OrdinalType,ScalarType>::reciprocalConditionEstimate(MagnitudeType & Value)
00903 {
00904 #ifdef HAVE_TEUCHOSNUMERICS_EIGEN
00905   // Implement templated GECON or use Eigen. Eigen currently doesn't have a condition estimation function.
00906   return (-1);
00907 #else
00908 
00909   if (reciprocalConditionEstimated()) {
00910     Value = RCOND_;
00911     return(0); // Already computed, just return it.
00912   }
00913 
00914   if ( ANORM_<ScalarTraits<MagnitudeType>::zero() ) ANORM_ = Matrix_->normOne();
00915 
00916   int ierr = 0;
00917   if (!factored()) ierr = factor(); // Need matrix factored.
00918   if (ierr!=0) return(ierr);
00919 
00920   allocateWORK();
00921 
00922   // We will assume a one-norm condition number
00923   INFO_ = 0;
00924   std::vector<typename details::lapack_traits<ScalarType>::iwork_type> GBCON_WORK( N_ );
00925   this->GBCON( '1', N_, KL_, KU_, AF_, LDAF_, &IPIV_[0], ANORM_, &RCOND_, &WORK_[0], &GBCON_WORK[0], &INFO_);
00926   reciprocalConditionEstimated_ = true;
00927   Value = RCOND_;
00928 
00929   return(INFO_);
00930 #endif
00931 }
00932 //=============================================================================
00933 
00934 template<typename OrdinalType, typename ScalarType>
00935 void SerialBandDenseSolver<OrdinalType,ScalarType>::Print(std::ostream& os) const {
00936 
00937   if (Matrix_!=Teuchos::null) os << "Solver Matrix"          << std::endl << *Matrix_ << std::endl;
00938   if (Factor_!=Teuchos::null) os << "Solver Factored Matrix" << std::endl << *Factor_ << std::endl;
00939   if (LHS_   !=Teuchos::null) os << "Solver LHS"             << std::endl << *LHS_    << std::endl;
00940   if (RHS_   !=Teuchos::null) os << "Solver RHS"             << std::endl << *RHS_    << std::endl;
00941 
00942 }
00943 
00944 //=============================================================================
00945 
00946 
00947 //=============================================================================
00948 
00949 
00950 } // namespace Teuchos
00951 
00952 #endif /* _TEUCHOS_SERIALBANDDENSESOLVER_HPP_ */
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Defines