Teuchos - Trilinos Tools Package Version of the Day
Teuchos_SerialDenseSolver.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_SERIALDENSESOLVER_HPP_
00043 #define _TEUCHOS_SERIALDENSESOLVER_HPP_
00044 
00048 #include "Teuchos_CompObject.hpp"
00049 #include "Teuchos_BLAS.hpp"
00050 #include "Teuchos_LAPACK.hpp"
00051 #include "Teuchos_RCP.hpp"
00052 #include "Teuchos_ConfigDefs.hpp"
00053 #include "Teuchos_SerialDenseMatrix.hpp"
00054 #include "Teuchos_ScalarTraits.hpp"
00055 
00056 #ifdef HAVE_TEUCHOSNUMERICS_EIGEN
00057 #include "Eigen/Dense"
00058 #endif
00059 
00134 namespace Teuchos {
00135 
00136   template<typename OrdinalType, typename ScalarType>
00137   class SerialDenseSolver : public CompObject, public Object, public BLAS<OrdinalType, ScalarType>,
00138           public LAPACK<OrdinalType, ScalarType>
00139   {
00140 
00141   public:
00142 
00143     typedef typename Teuchos::ScalarTraits<ScalarType>::magnitudeType MagnitudeType;
00144 #ifdef HAVE_TEUCHOSNUMERICS_EIGEN
00145     typedef typename Eigen::Matrix<ScalarType,Eigen::Dynamic,Eigen::Dynamic,Eigen::ColMajor> EigenMatrix;
00146     typedef typename Eigen::Matrix<ScalarType,Eigen::Dynamic,1> EigenVector;
00147     typedef typename Eigen::InnerStride<Eigen::Dynamic> EigenInnerStride;
00148     typedef typename Eigen::OuterStride<Eigen::Dynamic> EigenOuterStride;
00149     typedef typename Eigen::Map<EigenVector,0,EigenInnerStride > EigenVectorMap;
00150     typedef typename Eigen::Map<const EigenVector,0,EigenInnerStride > EigenConstVectorMap;
00151     typedef typename Eigen::Map<EigenMatrix,0,EigenOuterStride > EigenMatrixMap;
00152     typedef typename Eigen::Map<const EigenMatrix,0,EigenOuterStride > EigenConstMatrixMap;
00153     typedef typename Eigen::PermutationMatrix<Eigen::Dynamic,Eigen::Dynamic,OrdinalType> EigenPermutationMatrix;
00154     typedef typename Eigen::Array<OrdinalType,Eigen::Dynamic,1> EigenOrdinalArray;
00155     typedef typename Eigen::Map<EigenOrdinalArray> EigenOrdinalArrayMap;
00156     typedef typename Eigen::Array<ScalarType,Eigen::Dynamic,1> EigenScalarArray;
00157     typedef typename Eigen::Map<EigenScalarArray> EigenScalarArrayMap;
00158 #endif
00159 
00161 
00162 
00163     SerialDenseSolver();
00164 
00166     virtual ~SerialDenseSolver();
00168 
00170 
00171 
00173     int setMatrix(const RCP<SerialDenseMatrix<OrdinalType, ScalarType> >& A);
00174 
00176 
00179     int setVectors(const RCP<SerialDenseMatrix<OrdinalType, ScalarType> >& X,
00180        const RCP<SerialDenseMatrix<OrdinalType, ScalarType> >& B);
00182 
00184 
00185 
00187 
00189     void factorWithEquilibration(bool flag) {equilibrate_ = flag; return;}
00190 
00192 
00194     void solveWithTranspose(bool flag) {transpose_ = flag; if (flag) TRANS_ = Teuchos::TRANS; else TRANS_ = Teuchos::NO_TRANS; return;}
00195 
00197 
00199     void solveWithTransposeFlag(Teuchos::ETransp trans) {TRANS_ = trans; if (trans != Teuchos::NO_TRANS) {  transpose_ = true; } }
00200 
00202 
00204     void solveToRefinedSolution(bool flag) {refineSolution_ = flag; return;}
00205 
00207 
00211     void estimateSolutionErrors(bool flag);
00213 
00215 
00216 
00218 
00221     int factor();
00222 
00224 
00227     int solve();
00228 
00230 
00233     int invert();
00234 
00236 
00239     int computeEquilibrateScaling();
00240 
00242 
00246     int equilibrateMatrix();
00247 
00249 
00253     int equilibrateRHS();
00254 
00256 
00260     int applyRefinement();
00261 
00263 
00267     int unequilibrateLHS();
00268 
00270 
00276      int reciprocalConditionEstimate(MagnitudeType & Value);
00278 
00280 
00281 
00283     bool transpose() {return(transpose_);}
00284 
00286     bool factored() {return(factored_);}
00287 
00289     bool equilibratedA() {return(equilibratedA_);}
00290 
00292     bool equilibratedB() {return(equilibratedB_);}
00293 
00295      bool shouldEquilibrate() {computeEquilibrateScaling(); return(shouldEquilibrate_);}
00296 
00298     bool solutionErrorsEstimated() {return(solutionErrorsEstimated_);}
00299 
00301     bool inverted() {return(inverted_);}
00302 
00304     bool reciprocalConditionEstimated() {return(reciprocalConditionEstimated_);}
00305 
00307     bool solved() {return(solved_);}
00308 
00310     bool solutionRefined() {return(solutionRefined_);}
00312 
00314 
00315 
00317     RCP<SerialDenseMatrix<OrdinalType, ScalarType> > getMatrix()  const {return(Matrix_);}
00318 
00320     RCP<SerialDenseMatrix<OrdinalType, ScalarType> > getFactoredMatrix()  const {return(Factor_);}
00321 
00323     RCP<SerialDenseMatrix<OrdinalType, ScalarType> > getLHS() const {return(LHS_);}
00324 
00326     RCP<SerialDenseMatrix<OrdinalType, ScalarType> > getRHS() const {return(RHS_);}
00327 
00329     OrdinalType numRows()  const {return(M_);}
00330 
00332     OrdinalType numCols()  const {return(N_);}
00333 
00335     std::vector<OrdinalType> IPIV()  const {return(IPIV_);}
00336 
00338     MagnitudeType ANORM()  const {return(ANORM_);}
00339 
00341     MagnitudeType RCOND()  const {return(RCOND_);}
00342 
00344 
00346     MagnitudeType ROWCND()  const {return(ROWCND_);}
00347 
00349 
00351     MagnitudeType COLCND()  const {return(COLCND_);}
00352 
00354     MagnitudeType AMAX()  const {return(AMAX_);}
00355 
00357     std::vector<MagnitudeType> FERR()  const {return(FERR_);}
00358 
00360     std::vector<MagnitudeType> BERR()  const {return(BERR_);}
00361 
00363     std::vector<MagnitudeType> R()  const {return(R_);}
00364 
00366     std::vector<MagnitudeType> C()  const {return(C_);}
00368 
00370 
00371 
00372     void Print(std::ostream& os) const;
00374   protected:
00375 
00376     void allocateWORK() { LWORK_ = 4*N_; WORK_.resize( LWORK_ ); return;}
00377     void resetMatrix();
00378     void resetVectors();
00379 
00380 
00381     bool equilibrate_;
00382     bool shouldEquilibrate_;
00383     bool equilibratedA_;
00384     bool equilibratedB_;
00385     bool transpose_;
00386     bool factored_;
00387     bool estimateSolutionErrors_;
00388     bool solutionErrorsEstimated_;
00389     bool solved_;
00390     bool inverted_;
00391     bool reciprocalConditionEstimated_;
00392     bool refineSolution_;
00393     bool solutionRefined_;
00394 
00395     Teuchos::ETransp TRANS_;
00396 
00397     OrdinalType M_;
00398     OrdinalType N_;
00399     OrdinalType Min_MN_;
00400     OrdinalType LDA_;
00401     OrdinalType LDAF_;
00402     OrdinalType INFO_;
00403     OrdinalType LWORK_;
00404 
00405     std::vector<OrdinalType> IPIV_;
00406 
00407     MagnitudeType ANORM_;
00408     MagnitudeType RCOND_;
00409     MagnitudeType ROWCND_;
00410     MagnitudeType COLCND_;
00411     MagnitudeType AMAX_;
00412 
00413     RCP<SerialDenseMatrix<OrdinalType, ScalarType> > Matrix_;
00414     RCP<SerialDenseMatrix<OrdinalType, ScalarType> > LHS_;
00415     RCP<SerialDenseMatrix<OrdinalType, ScalarType> > RHS_;
00416     RCP<SerialDenseMatrix<OrdinalType, ScalarType> > Factor_;
00417 
00418     ScalarType * A_;
00419     ScalarType * AF_;
00420     std::vector<MagnitudeType> FERR_;
00421     std::vector<MagnitudeType> BERR_;
00422     std::vector<ScalarType> WORK_;
00423     std::vector<MagnitudeType> R_;
00424     std::vector<MagnitudeType> C_;
00425 #ifdef HAVE_TEUCHOSNUMERICS_EIGEN
00426     Eigen::PartialPivLU<EigenMatrix> lu_;
00427 #endif
00428 
00429   private:
00430     // SerialDenseSolver copy constructor (put here because we don't want user access)
00431 
00432     SerialDenseSolver(const SerialDenseSolver<OrdinalType, ScalarType>& Source);
00433     SerialDenseSolver & operator=(const SerialDenseSolver<OrdinalType, ScalarType>& Source);
00434 
00435   };
00436 
00437   namespace details {
00438 
00439     // Helper traits to distinguish work arrays for real and complex-valued datatypes.
00440     template<typename T>
00441     struct lapack_traits {
00442       typedef int iwork_type;
00443     };
00444 
00445     // Complex-valued specialization
00446     template<typename T>
00447     struct lapack_traits<std::complex<T> > {
00448       typedef typename ScalarTraits<T>::magnitudeType iwork_type;
00449     };
00450 
00451   } // end namespace details
00452 
00453 //=============================================================================
00454 
00455 template<typename OrdinalType, typename ScalarType>
00456 SerialDenseSolver<OrdinalType,ScalarType>::SerialDenseSolver()
00457   : CompObject(),
00458     Object("Teuchos::SerialDenseSolver"),
00459     equilibrate_(false),
00460     shouldEquilibrate_(false),
00461     equilibratedA_(false),
00462     equilibratedB_(false),
00463     transpose_(false),
00464     factored_(false),
00465     estimateSolutionErrors_(false),
00466     solutionErrorsEstimated_(false),
00467     solved_(false),
00468     inverted_(false),
00469     reciprocalConditionEstimated_(false),
00470     refineSolution_(false),
00471     solutionRefined_(false),
00472     TRANS_(Teuchos::NO_TRANS),
00473     M_(0),
00474     N_(0),
00475     Min_MN_(0),
00476     LDA_(0),
00477     LDAF_(0),
00478     INFO_(0),
00479     LWORK_(0),
00480     ANORM_(ScalarTraits<MagnitudeType>::zero()),
00481     RCOND_(ScalarTraits<MagnitudeType>::zero()),
00482     ROWCND_(ScalarTraits<MagnitudeType>::zero()),
00483     COLCND_(ScalarTraits<MagnitudeType>::zero()),
00484     AMAX_(ScalarTraits<MagnitudeType>::zero()),
00485     A_(0),
00486     AF_(0)
00487 {
00488   resetMatrix();
00489 }
00490 //=============================================================================
00491 
00492 template<typename OrdinalType, typename ScalarType>
00493 SerialDenseSolver<OrdinalType,ScalarType>::~SerialDenseSolver()
00494 {}
00495 
00496 //=============================================================================
00497 
00498 template<typename OrdinalType, typename ScalarType>
00499 void SerialDenseSolver<OrdinalType,ScalarType>::resetVectors()
00500 {
00501   LHS_ = Teuchos::null;
00502   RHS_ = Teuchos::null;
00503   reciprocalConditionEstimated_ = false;
00504   solutionRefined_ = false;
00505   solved_ = false;
00506   solutionErrorsEstimated_ = false;
00507   equilibratedB_ = false;
00508 }
00509 //=============================================================================
00510 
00511 template<typename OrdinalType, typename ScalarType>
00512 void SerialDenseSolver<OrdinalType,ScalarType>::resetMatrix()
00513 {
00514   resetVectors();
00515   equilibratedA_ = false;
00516   factored_ = false;
00517   inverted_ = false;
00518   M_ = 0;
00519   N_ = 0;
00520   Min_MN_ = 0;
00521   LDA_ = 0;
00522   LDAF_ = 0;
00523   ANORM_ = -ScalarTraits<MagnitudeType>::one();
00524   RCOND_ = -ScalarTraits<MagnitudeType>::one();
00525   ROWCND_ = -ScalarTraits<MagnitudeType>::one();
00526   COLCND_ = -ScalarTraits<MagnitudeType>::one();
00527   AMAX_ = -ScalarTraits<MagnitudeType>::one();
00528   A_ = 0;
00529   AF_ = 0;
00530   INFO_ = 0;
00531   LWORK_ = 0;
00532   R_.resize(0);
00533   C_.resize(0);
00534 }
00535 //=============================================================================
00536 
00537 template<typename OrdinalType, typename ScalarType>
00538 int SerialDenseSolver<OrdinalType,ScalarType>::setMatrix(const RCP<SerialDenseMatrix<OrdinalType,ScalarType> >& A)
00539 {
00540   resetMatrix();
00541   Matrix_ = A;
00542   Factor_ = A;
00543   M_ = A->numRows();
00544   N_ = A->numCols();
00545   Min_MN_ = TEUCHOS_MIN(M_,N_);
00546   LDA_ = A->stride();
00547   LDAF_ = LDA_;
00548   A_ = A->values();
00549   AF_ = A->values();
00550   return(0);
00551 }
00552 //=============================================================================
00553 
00554 template<typename OrdinalType, typename ScalarType>
00555 int SerialDenseSolver<OrdinalType,ScalarType>::setVectors(const RCP<SerialDenseMatrix<OrdinalType,ScalarType> >& X,
00556                  const RCP<SerialDenseMatrix<OrdinalType,ScalarType> >& B)
00557 {
00558   // Check that these new vectors are consistent.
00559   TEUCHOS_TEST_FOR_EXCEPTION(B->numRows()!=X->numRows() || B->numCols() != X->numCols(), std::invalid_argument,
00560          "SerialDenseSolver<T>::setVectors: X and B are not the same size!");
00561   TEUCHOS_TEST_FOR_EXCEPTION(B->values()==0, std::invalid_argument,
00562          "SerialDenseSolver<T>::setVectors: B is an empty SerialDenseMatrix<T>!");
00563   TEUCHOS_TEST_FOR_EXCEPTION(X->values()==0, std::invalid_argument,
00564          "SerialDenseSolver<T>::setVectors: X is an empty SerialDenseMatrix<T>!");
00565   TEUCHOS_TEST_FOR_EXCEPTION(B->stride()<1, std::invalid_argument,
00566          "SerialDenseSolver<T>::setVectors: B has an invalid stride!");
00567   TEUCHOS_TEST_FOR_EXCEPTION(X->stride()<1, std::invalid_argument,
00568          "SerialDenseSolver<T>::setVectors: X has an invalid stride!");
00569 
00570   resetVectors();
00571   LHS_ = X;
00572   RHS_ = B;
00573   return(0);
00574 }
00575 //=============================================================================
00576 
00577 template<typename OrdinalType, typename ScalarType>
00578 void SerialDenseSolver<OrdinalType,ScalarType>::estimateSolutionErrors(bool flag)
00579 {
00580   estimateSolutionErrors_ = flag;
00581 
00582   // If the errors are estimated, this implies that the solution must be refined
00583   refineSolution_ = refineSolution_ || flag;
00584 }
00585 //=============================================================================
00586 
00587 template<typename OrdinalType, typename ScalarType>
00588 int SerialDenseSolver<OrdinalType,ScalarType>::factor() {
00589 
00590   if (factored()) return(0); // Already factored
00591 
00592   TEUCHOS_TEST_FOR_EXCEPTION(inverted(), std::logic_error,
00593          "SerialDenseSolver<T>::factor: Cannot factor an inverted matrix!");
00594 
00595   ANORM_ = Matrix_->normOne(); // Compute 1-Norm of A
00596 
00597 
00598   // If we want to refine the solution, then the factor must
00599   // be stored separatedly from the original matrix
00600 
00601   if (A_ == AF_)
00602     if (refineSolution_ ) {
00603       Factor_ = rcp( new SerialDenseMatrix<OrdinalType,ScalarType>(*Matrix_) );
00604       AF_ = Factor_->values();
00605       LDAF_ = Factor_->stride();
00606     }
00607 
00608   int ierr = 0;
00609   if (equilibrate_) ierr = equilibrateMatrix();
00610 
00611   if (ierr!=0) return(ierr);
00612 
00613   if ((int)IPIV_.size()!=Min_MN_) IPIV_.resize( Min_MN_ ); // Allocated Pivot vector if not already done.
00614 
00615   INFO_ = 0;
00616 
00617 #ifdef HAVE_TEUCHOSNUMERICS_EIGEN
00618   EigenMatrixMap aMap( AF_, M_, N_, EigenOuterStride(LDAF_) );
00619   EigenPermutationMatrix p;
00620   EigenOrdinalArray ind;
00621   EigenOrdinalArrayMap ipivMap( &IPIV_[0], Min_MN_ );
00622 
00623   lu_.compute(aMap);
00624   p = lu_.permutationP();
00625   ind = p.indices();
00626 
00627   EigenMatrix luMat = lu_.matrixLU();
00628   for (OrdinalType j=0; j<aMap.outerSize(); j++) {
00629     for (OrdinalType i=0; i<aMap.innerSize(); i++) {
00630       aMap(i,j) = luMat(i,j);
00631     }
00632   }
00633   for (OrdinalType i=0; i<ipivMap.innerSize(); i++) {
00634     ipivMap(i) = ind(i);
00635   }
00636 #else
00637   this->GETRF(M_, N_, AF_, LDAF_, &IPIV_[0], &INFO_);
00638 #endif
00639 
00640   factored_ = true;
00641 
00642   return(INFO_);
00643 }
00644 
00645 //=============================================================================
00646 
00647 template<typename OrdinalType, typename ScalarType>
00648 int SerialDenseSolver<OrdinalType,ScalarType>::solve() {
00649 
00650   // We will call one of four routines depending on what services the user wants and
00651   // whether or not the matrix has been inverted or factored already.
00652   //
00653   // If the matrix has been inverted, use DGEMM to compute solution.
00654   // Otherwise, if the user want the matrix to be equilibrated or wants a refined solution, we will
00655   // call the X interface.
00656   // Otherwise, if the matrix is already factored we will call the TRS interface.
00657   // Otherwise, if the matrix is unfactored we will call the SV interface.
00658 
00659   int ierr = 0;
00660   if (equilibrate_) {
00661     ierr = equilibrateRHS();
00662     equilibratedB_ = true;
00663   }
00664   if (ierr != 0) return(ierr);  // Can't equilibrate B, so return.
00665 
00666   TEUCHOS_TEST_FOR_EXCEPTION( (equilibratedA_ && !equilibratedB_) || (!equilibratedA_ && equilibratedB_) ,
00667          std::logic_error, "SerialDenseSolver<T>::solve: Matrix and vectors must be similarly scaled!");
00668   TEUCHOS_TEST_FOR_EXCEPTION( RHS_==Teuchos::null, std::invalid_argument,
00669          "SerialDenseSolver<T>::solve: No right-hand side vector (RHS) has been set for the linear system!");
00670   TEUCHOS_TEST_FOR_EXCEPTION( LHS_==Teuchos::null, std::invalid_argument,
00671          "SerialDenseSolver<T>::solve: No solution vector (LHS) has been set for the linear system!");
00672 
00673   if (shouldEquilibrate() && !equilibratedA_)
00674     std::cout << "WARNING!  SerialDenseSolver<T>::solve: System should be equilibrated!" << std::endl;
00675 
00676   if (inverted()) {
00677 
00678     TEUCHOS_TEST_FOR_EXCEPTION( RHS_->values() == LHS_->values(), std::invalid_argument,
00679       "SerialDenseSolver<T>::solve: X and B must be different vectors if matrix is inverted.");
00680 
00681     INFO_ = 0;
00682     this->GEMM(TRANS_, Teuchos::NO_TRANS, N_, RHS_->numCols(), N_, 1.0, AF_, LDAF_,
00683          RHS_->values(), RHS_->stride(), 0.0, LHS_->values(), LHS_->stride());
00684     if (INFO_!=0) return(INFO_);
00685     solved_ = true;
00686   }
00687   else {
00688 
00689     if (!factored()) factor(); // Matrix must be factored
00690 
00691     if (RHS_->values()!=LHS_->values()) {
00692        (*LHS_) = (*RHS_); // Copy B to X if needed
00693     }
00694     INFO_ = 0;
00695 
00696 #ifdef HAVE_TEUCHOSNUMERICS_EIGEN
00697     EigenMatrixMap rhsMap( RHS_->values(), RHS_->numRows(), RHS_->numCols(), EigenOuterStride(RHS_->stride()) );
00698     EigenMatrixMap lhsMap( LHS_->values(), LHS_->numRows(), LHS_->numCols(), EigenOuterStride(LHS_->stride()) );
00699     if ( TRANS_ == Teuchos::NO_TRANS ) {
00700       lhsMap = lu_.solve(rhsMap);
00701     } else if ( TRANS_ == Teuchos::TRANS ) {
00702       lu_.matrixLU().template triangularView<Eigen::Upper>().transpose().solveInPlace(lhsMap);
00703       lu_.matrixLU().template triangularView<Eigen::UnitLower>().transpose().solveInPlace(lhsMap);
00704       EigenMatrix x = lu_.permutationP().transpose() * lhsMap;
00705       lhsMap = x;
00706     } else {
00707       lu_.matrixLU().template triangularView<Eigen::Upper>().adjoint().solveInPlace(lhsMap);
00708       lu_.matrixLU().template triangularView<Eigen::UnitLower>().adjoint().solveInPlace(lhsMap);
00709       EigenMatrix x = lu_.permutationP().transpose() * lhsMap;
00710       lhsMap = x;
00711     }
00712 #else
00713     this->GETRS(ETranspChar[TRANS_], N_, RHS_->numCols(), AF_, LDAF_, &IPIV_[0], LHS_->values(), LHS_->stride(), &INFO_);
00714 #endif
00715 
00716     if (INFO_!=0) return(INFO_);
00717     solved_ = true;
00718 
00719   }
00720   int ierr1=0;
00721   if (refineSolution_ && !inverted()) ierr1 = applyRefinement();
00722   if (ierr1!=0)
00723     return(ierr1);
00724 
00725   if (equilibrate_) ierr1 = unequilibrateLHS();
00726   return(ierr1);
00727 }
00728 //=============================================================================
00729 
00730 template<typename OrdinalType, typename ScalarType>
00731 int SerialDenseSolver<OrdinalType,ScalarType>::applyRefinement()
00732 {
00733   TEUCHOS_TEST_FOR_EXCEPTION(!solved(), std::logic_error,
00734          "SerialDenseSolver<T>::applyRefinement: Must have an existing solution!");
00735   TEUCHOS_TEST_FOR_EXCEPTION(A_==AF_, std::logic_error,
00736          "SerialDenseSolver<T>::applyRefinement: Cannot apply refinement if no original copy of A!");
00737 
00738 #ifdef HAVE_TEUCHOSNUMERICS_EIGEN
00739   // Implement templated GERFS or use Eigen.
00740   return (-1);
00741 #else
00742 
00743   OrdinalType NRHS = RHS_->numCols();
00744   FERR_.resize( NRHS );
00745   BERR_.resize( NRHS );
00746   allocateWORK();
00747 
00748   INFO_ = 0;
00749   std::vector<typename details::lapack_traits<ScalarType>::iwork_type> GERFS_WORK( N_ );
00750   this->GERFS(ETranspChar[TRANS_], N_, NRHS, A_, LDA_, AF_, LDAF_, &IPIV_[0],
00751         RHS_->values(), RHS_->stride(), LHS_->values(), LHS_->stride(),
00752         &FERR_[0], &BERR_[0], &WORK_[0], &GERFS_WORK[0], &INFO_);
00753 
00754   solutionErrorsEstimated_ = true;
00755   reciprocalConditionEstimated_ = true;
00756   solutionRefined_ = true;
00757 
00758   return(INFO_);
00759 #endif
00760 }
00761 
00762 //=============================================================================
00763 
00764 template<typename OrdinalType, typename ScalarType>
00765 int SerialDenseSolver<OrdinalType,ScalarType>::computeEquilibrateScaling()
00766 {
00767   if (R_.size()!=0) return(0); // Already computed
00768 
00769   R_.resize( M_ );
00770   C_.resize( N_ );
00771 
00772   INFO_ = 0;
00773   this->GEEQU (M_, N_, AF_, LDAF_, &R_[0], &C_[0], &ROWCND_, &COLCND_, &AMAX_, &INFO_);
00774 
00775   if (COLCND_<0.1*ScalarTraits<MagnitudeType>::one() ||
00776       ROWCND_<0.1*ScalarTraits<MagnitudeType>::one() ||
00777       AMAX_ < ScalarTraits<ScalarType>::rmin() || AMAX_ > ScalarTraits<ScalarType>::rmax() )
00778     shouldEquilibrate_ = true;
00779 
00780   return(INFO_);
00781 }
00782 
00783 //=============================================================================
00784 
00785 template<typename OrdinalType, typename ScalarType>
00786 int SerialDenseSolver<OrdinalType,ScalarType>::equilibrateMatrix()
00787 {
00788   OrdinalType i, j;
00789   int ierr = 0;
00790 
00791   if (equilibratedA_) return(0); // Already done.
00792   if (R_.size()==0) ierr = computeEquilibrateScaling(); // Compute R and C if needed.
00793   if (ierr!=0) return(ierr);     // If return value is not zero, then can't equilibrate.
00794   if (A_==AF_) {
00795     ScalarType * ptr;
00796     for (j=0; j<N_; j++) {
00797       ptr = A_ + j*LDA_;
00798       ScalarType s1 = C_[j];
00799       for (i=0; i<M_; i++) {
00800   *ptr = *ptr*s1*R_[i];
00801   ptr++;
00802       }
00803     }
00804   }
00805   else {
00806     ScalarType * ptr;
00807     ScalarType * ptr1;
00808     for (j=0; j<N_; j++) {
00809       ptr = A_ + j*LDA_;
00810       ptr1 = AF_ + j*LDAF_;
00811       ScalarType s1 = C_[j];
00812       for (i=0; i<M_; i++) {
00813   *ptr = *ptr*s1*R_[i];
00814   ptr++;
00815   *ptr1 = *ptr1*s1*R_[i];
00816   ptr1++;
00817       }
00818     }
00819   }
00820 
00821   equilibratedA_ = true;
00822 
00823   return(0);
00824 }
00825 
00826 //=============================================================================
00827 
00828 template<typename OrdinalType, typename ScalarType>
00829 int SerialDenseSolver<OrdinalType,ScalarType>::equilibrateRHS()
00830 {
00831   OrdinalType i, j;
00832   int ierr = 0;
00833 
00834   if (equilibratedB_) return(0); // Already done.
00835   if (R_.size()==0) ierr = computeEquilibrateScaling(); // Compute R and C if needed.
00836   if (ierr!=0) return(ierr);     // Can't count on R and C being computed.
00837 
00838   MagnitudeType * R_tmp = &R_[0];
00839   if (transpose_) R_tmp = &C_[0];
00840 
00841   OrdinalType LDB = RHS_->stride(), NRHS = RHS_->numCols();
00842   ScalarType * B = RHS_->values();
00843   ScalarType * ptr;
00844   for (j=0; j<NRHS; j++) {
00845     ptr = B + j*LDB;
00846     for (i=0; i<M_; i++) {
00847       *ptr = *ptr*R_tmp[i];
00848       ptr++;
00849     }
00850   }
00851 
00852   equilibratedB_ = true;
00853 
00854   return(0);
00855 }
00856 
00857 //=============================================================================
00858 
00859 template<typename OrdinalType, typename ScalarType>
00860 int SerialDenseSolver<OrdinalType,ScalarType>::unequilibrateLHS()
00861 {
00862   OrdinalType i, j;
00863 
00864   if (!equilibratedB_) return(0); // Nothing to do
00865 
00866   MagnitudeType * C_tmp = &C_[0];
00867   if (transpose_) C_tmp = &R_[0];
00868 
00869   OrdinalType LDX = LHS_->stride(), NLHS = LHS_->numCols();
00870   ScalarType * X = LHS_->values();
00871   ScalarType * ptr;
00872   for (j=0; j<NLHS; j++) {
00873     ptr = X + j*LDX;
00874     for (i=0; i<N_; i++) {
00875       *ptr = *ptr*C_tmp[i];
00876       ptr++;
00877     }
00878   }
00879 
00880   return(0);
00881 }
00882 
00883 //=============================================================================
00884 
00885 template<typename OrdinalType, typename ScalarType>
00886 int SerialDenseSolver<OrdinalType,ScalarType>::invert()
00887 {
00888 
00889   if (!factored()) factor(); // Need matrix factored.
00890 
00891 #ifdef HAVE_TEUCHOSNUMERICS_EIGEN
00892   EigenMatrixMap invMap( AF_, M_, N_, EigenOuterStride(LDAF_) );
00893   EigenMatrix invMat = lu_.inverse();
00894   for (OrdinalType j=0; j<invMap.outerSize(); j++) {
00895     for (OrdinalType i=0; i<invMap.innerSize(); i++) {
00896       invMap(i,j) = invMat(i,j);
00897     }
00898   }
00899 #else
00900 
00901   /* This section work with LAPACK Version 3.0 only
00902   // Setting LWORK = -1 and calling GETRI will return optimal work space size in WORK_TMP
00903   OrdinalType LWORK_TMP = -1;
00904   ScalarType WORK_TMP;
00905   GETRI( N_, AF_, LDAF_, &IPIV_[0], &WORK_TMP, &LWORK_TMP, &INFO_);
00906   LWORK_TMP = WORK_TMP; // Convert to integer
00907   if (LWORK_TMP>LWORK_) {
00908   LWORK_ = LWORK_TMP;
00909   WORK_.resize( LWORK_ );
00910   }
00911   */
00912   // This section will work with any version of LAPACK
00913   allocateWORK();
00914 
00915   INFO_ = 0;
00916   this->GETRI( N_, AF_, LDAF_, &IPIV_[0], &WORK_[0], LWORK_, &INFO_);
00917 #endif
00918 
00919   inverted_ = true;
00920   factored_ = false;
00921 
00922   return(INFO_);
00923 
00924 }
00925 
00926 //=============================================================================
00927 
00928 template<typename OrdinalType, typename ScalarType>
00929 int SerialDenseSolver<OrdinalType,ScalarType>::reciprocalConditionEstimate(MagnitudeType & Value)
00930 {
00931 #ifdef HAVE_TEUCHOSNUMERICS_EIGEN
00932   // Implement templated GECON or use Eigen. Eigen currently doesn't have a condition estimation function.
00933   return (-1);
00934 #else
00935   if (reciprocalConditionEstimated()) {
00936     Value = RCOND_;
00937     return(0); // Already computed, just return it.
00938   }
00939 
00940   if ( ANORM_<ScalarTraits<MagnitudeType>::zero() ) ANORM_ = Matrix_->normOne();
00941 
00942   int ierr = 0;
00943   if (!factored()) ierr = factor(); // Need matrix factored.
00944   if (ierr!=0) return(ierr);
00945 
00946   allocateWORK();
00947 
00948   // We will assume a one-norm condition number
00949   INFO_ = 0;
00950   std::vector<typename details::lapack_traits<ScalarType>::iwork_type> GECON_WORK( 2*N_ );
00951   this->GECON( '1', N_, AF_, LDAF_, ANORM_, &RCOND_, &WORK_[0], &GECON_WORK[0], &INFO_);
00952 
00953   reciprocalConditionEstimated_ = true;
00954   Value = RCOND_;
00955 
00956   return(INFO_);
00957 #endif
00958 }
00959 //=============================================================================
00960 
00961 template<typename OrdinalType, typename ScalarType>
00962 void SerialDenseSolver<OrdinalType,ScalarType>::Print(std::ostream& os) const {
00963 
00964   if (Matrix_!=Teuchos::null) os << "Solver Matrix"          << std::endl << *Matrix_ << std::endl;
00965   if (Factor_!=Teuchos::null) os << "Solver Factored Matrix" << std::endl << *Factor_ << std::endl;
00966   if (LHS_   !=Teuchos::null) os << "Solver LHS"             << std::endl << *LHS_    << std::endl;
00967   if (RHS_   !=Teuchos::null) os << "Solver RHS"             << std::endl << *RHS_    << std::endl;
00968 
00969 }
00970 
00971 } // namespace Teuchos
00972 
00973 #endif /* _TEUCHOS_SERIALDENSESOLVER_HPP_ */
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Defines