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 // This library is free software; you can redistribute it and/or modify
00011 // it under the terms of the GNU Lesser General Public License as
00012 // published by the Free Software Foundation; either version 2.1 of the
00013 // License, or (at your option) any later version.
00014 //  
00015 // This library is distributed in the hope that it will be useful, but
00016 // WITHOUT ANY WARRANTY; without even the implied warranty of
00017 // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
00018 // Lesser General Public License for more details.
00019 //  
00020 // You should have received a copy of the GNU Lesser General Public
00021 // License along with this library; if not, write to the Free Software
00022 // Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307
00023 // USA
00024 // Questions? Contact Michael A. Heroux (maherou@sandia.gov) 
00025 // 
00026 // ***********************************************************************
00027 // @HEADER
00028 
00029 #ifndef _TEUCHOS_SERIALDENSESOLVER_HPP_
00030 #define _TEUCHOS_SERIALDENSESOLVER_HPP_
00031 
00035 #include "Teuchos_CompObject.hpp"
00036 #include "Teuchos_BLAS.hpp"
00037 #include "Teuchos_LAPACK.hpp"
00038 #include "Teuchos_RCP.hpp"
00039 #include "Teuchos_ConfigDefs.hpp"
00040 #include "Teuchos_SerialDenseMatrix.hpp"
00041 
00116 namespace Teuchos {
00117 
00118   template<typename OrdinalType, typename ScalarType>
00119   class SerialDenseSolver : public CompObject, public Object, public BLAS<OrdinalType, ScalarType>,
00120           public LAPACK<OrdinalType, ScalarType>
00121   {
00122 
00123   public:
00124 
00125     typedef typename Teuchos::ScalarTraits<ScalarType>::magnitudeType MagnitudeType;
00126   
00128 
00129 
00130     SerialDenseSolver();
00131     
00133     virtual ~SerialDenseSolver();
00135     
00137 
00138     
00140     int setMatrix(const RCP<SerialDenseMatrix<OrdinalType, ScalarType> >& A);
00141 
00143 
00146     int setVectors(const RCP<SerialDenseMatrix<OrdinalType, ScalarType> >& X, 
00147        const RCP<SerialDenseMatrix<OrdinalType, ScalarType> >& B);
00149     
00151 
00152     
00154 
00156     void factorWithEquilibration(bool flag) {equilibrate_ = flag; return;}
00157     
00159     void solveWithTranspose(bool flag) {transpose_ = flag; if (flag) TRANS_ = Teuchos::TRANS; else TRANS_ = Teuchos::NO_TRANS; return;}
00160     
00162     void solveToRefinedSolution(bool flag) {refineSolution_ = flag; return;}
00163     
00165 
00168     void estimateSolutionErrors(bool flag);
00170     
00172 
00173     
00175 
00178     int factor();
00179     
00181 
00184     int solve();
00185     
00187 
00190     int invert();
00191     
00193 
00196     int computeEquilibrateScaling();
00197     
00199 
00202     int equilibrateMatrix();
00203     
00205 
00208     int equilibrateRHS();
00209         
00211 
00214     int applyRefinement();
00215     
00217 
00220     int unequilibrateLHS();
00221     
00223 
00229      int reciprocalConditionEstimate(MagnitudeType & Value);
00231     
00233 
00234     
00236     bool transpose() {return(transpose_);}
00237     
00239     bool factored() {return(factored_);}
00240     
00242     bool equilibratedA() {return(equilibratedA_);}
00243     
00245     bool equilibratedB() {return(equilibratedB_);}
00246     
00248      bool shouldEquilibrate() {computeEquilibrateScaling(); return(shouldEquilibrate_);}
00249     
00251     bool solutionErrorsEstimated() {return(solutionErrorsEstimated_);}
00252     
00254     bool inverted() {return(inverted_);}
00255     
00257     bool reciprocalConditionEstimated() {return(reciprocalConditionEstimated_);}
00258     
00260     bool solved() {return(solved_);}
00261     
00263     bool solutionRefined() {return(solutionRefined_);}
00265     
00267 
00268     
00270     RCP<SerialDenseMatrix<OrdinalType, ScalarType> > getMatrix()  const {return(Matrix_);}
00271     
00273     RCP<SerialDenseMatrix<OrdinalType, ScalarType> > getFactoredMatrix()  const {return(Factor_);}
00274     
00276     RCP<SerialDenseMatrix<OrdinalType, ScalarType> > getLHS() const {return(LHS_);}
00277     
00279     RCP<SerialDenseMatrix<OrdinalType, ScalarType> > getRHS() const {return(RHS_);}
00280     
00282     OrdinalType numRows()  const {return(M_);}
00283     
00285     OrdinalType numCols()  const {return(N_);}
00286     
00288     std::vector<OrdinalType> IPIV()  const {return(IPIV_);}
00289     
00291     MagnitudeType ANORM()  const {return(ANORM_);}
00292     
00294     MagnitudeType RCOND()  const {return(RCOND_);}
00295     
00297 
00299     MagnitudeType ROWCND()  const {return(ROWCND_);}
00300     
00302 
00304     MagnitudeType COLCND()  const {return(COLCND_);}
00305     
00307     MagnitudeType AMAX()  const {return(AMAX_);}
00308     
00310     std::vector<MagnitudeType> FERR()  const {return(FERR_);}
00311     
00313     std::vector<MagnitudeType> BERR()  const {return(BERR_);}
00314     
00316     std::vector<MagnitudeType> R()  const {return(R_);}
00317     
00319     std::vector<MagnitudeType> C()  const {return(C_);}
00321     
00323 
00324 
00325     void Print(std::ostream& os) const;
00327   protected:
00328     
00329     void allocateWORK() { LWORK_ = 4*N_; WORK_.resize( LWORK_ ); return;}
00330     void allocateIWORK() { IWORK_.resize( N_ ); return;}
00331     void resetMatrix();
00332     void resetVectors();
00333     
00334     
00335     bool equilibrate_;
00336     bool shouldEquilibrate_;
00337     bool equilibratedA_;
00338     bool equilibratedB_;
00339     bool transpose_;
00340     bool factored_;
00341     bool estimateSolutionErrors_;
00342     bool solutionErrorsEstimated_;
00343     bool solved_;
00344     bool inverted_;
00345     bool reciprocalConditionEstimated_;
00346     bool refineSolution_;
00347     bool solutionRefined_;
00348     
00349     Teuchos::ETransp TRANS_;
00350     
00351     OrdinalType M_;
00352     OrdinalType N_;
00353     OrdinalType Min_MN_;
00354     OrdinalType LDA_;
00355     OrdinalType LDAF_;
00356     OrdinalType INFO_;
00357     OrdinalType LWORK_;
00358     
00359     std::vector<OrdinalType> IPIV_;
00360     std::vector<int> IWORK_;
00361     
00362     MagnitudeType ANORM_;
00363     MagnitudeType RCOND_;
00364     MagnitudeType ROWCND_;
00365     MagnitudeType COLCND_;
00366     MagnitudeType AMAX_;
00367     
00368     RCP<SerialDenseMatrix<OrdinalType, ScalarType> > Matrix_;
00369     RCP<SerialDenseMatrix<OrdinalType, ScalarType> > LHS_;
00370     RCP<SerialDenseMatrix<OrdinalType, ScalarType> > RHS_;
00371     RCP<SerialDenseMatrix<OrdinalType, ScalarType> > Factor_;
00372     
00373     ScalarType * A_;
00374     ScalarType * AF_;
00375     std::vector<MagnitudeType> FERR_;
00376     std::vector<MagnitudeType> BERR_;
00377     std::vector<ScalarType> WORK_;
00378     std::vector<MagnitudeType> R_;
00379     std::vector<MagnitudeType> C_;
00380     
00381   private:
00382     // SerialDenseSolver copy constructor (put here because we don't want user access)
00383     
00384     SerialDenseSolver(const SerialDenseSolver<OrdinalType, ScalarType>& Source);
00385     SerialDenseSolver & operator=(const SerialDenseSolver<OrdinalType, ScalarType>& Source);
00386     
00387   };  
00388 
00389 //=============================================================================
00390 
00391 template<typename OrdinalType, typename ScalarType>
00392 SerialDenseSolver<OrdinalType,ScalarType>::SerialDenseSolver()
00393   : CompObject(),
00394     Object("Teuchos::SerialDenseSolver"),
00395     equilibrate_(false),
00396     shouldEquilibrate_(false),
00397     equilibratedA_(false),
00398     equilibratedB_(false),
00399     transpose_(false),
00400     factored_(false),
00401     estimateSolutionErrors_(false),
00402     solutionErrorsEstimated_(false),
00403     solved_(false),
00404     inverted_(false),
00405     reciprocalConditionEstimated_(false),
00406     refineSolution_(false),
00407     solutionRefined_(false),
00408     TRANS_(Teuchos::NO_TRANS),
00409     M_(0),
00410     N_(0),
00411     Min_MN_(0),
00412     LDA_(0),
00413     LDAF_(0),
00414     INFO_(0),
00415     LWORK_(0),
00416     ANORM_(ScalarTraits<ScalarType>::zero()),
00417     RCOND_(ScalarTraits<ScalarType>::zero()),
00418     ROWCND_(ScalarTraits<ScalarType>::zero()),
00419     COLCND_(ScalarTraits<ScalarType>::zero()),
00420     AMAX_(ScalarTraits<ScalarType>::zero()),
00421     A_(0),
00422     AF_(0)
00423 {
00424   resetMatrix();
00425 }
00426 //=============================================================================
00427 
00428 template<typename OrdinalType, typename ScalarType>
00429 SerialDenseSolver<OrdinalType,ScalarType>::~SerialDenseSolver()
00430 {}
00431 
00432 //=============================================================================
00433 
00434 template<typename OrdinalType, typename ScalarType>
00435 void SerialDenseSolver<OrdinalType,ScalarType>::resetVectors()
00436 {
00437   LHS_ = Teuchos::null;
00438   RHS_ = Teuchos::null;
00439   reciprocalConditionEstimated_ = false;
00440   solutionRefined_ = false;
00441   solved_ = false;
00442   solutionErrorsEstimated_ = false;
00443   equilibratedB_ = false;
00444 }
00445 //=============================================================================
00446 
00447 template<typename OrdinalType, typename ScalarType>
00448 void SerialDenseSolver<OrdinalType,ScalarType>::resetMatrix()
00449 {
00450   resetVectors();
00451   equilibratedA_ = false;
00452   factored_ = false;
00453   inverted_ = false;
00454   M_ = 0;
00455   N_ = 0;
00456   Min_MN_ = 0;
00457   LDA_ = 0;
00458   LDAF_ = 0;
00459   ANORM_ = -ScalarTraits<MagnitudeType>::one();
00460   RCOND_ = -ScalarTraits<MagnitudeType>::one();
00461   ROWCND_ = -ScalarTraits<MagnitudeType>::one();
00462   COLCND_ = -ScalarTraits<MagnitudeType>::one();
00463   AMAX_ = -ScalarTraits<MagnitudeType>::one();
00464   A_ = 0;
00465   AF_ = 0;
00466   INFO_ = 0;
00467   LWORK_ = 0;
00468   R_.resize(0);
00469   C_.resize(0);
00470 }
00471 //=============================================================================
00472 
00473 template<typename OrdinalType, typename ScalarType>
00474 int SerialDenseSolver<OrdinalType,ScalarType>::setMatrix(const RCP<SerialDenseMatrix<OrdinalType,ScalarType> >& A) 
00475 {
00476   resetMatrix();
00477   Matrix_ = A;
00478   Factor_ = A;
00479   M_ = A->numRows();
00480   N_ = A->numCols();
00481   Min_MN_ = TEUCHOS_MIN(M_,N_);
00482   LDA_ = A->stride();
00483   LDAF_ = LDA_;
00484   A_ = A->values();
00485   AF_ = A->values();
00486   return(0);
00487 }
00488 //=============================================================================
00489 
00490 template<typename OrdinalType, typename ScalarType>
00491 int SerialDenseSolver<OrdinalType,ScalarType>::setVectors(const RCP<SerialDenseMatrix<OrdinalType,ScalarType> >& X, 
00492                  const RCP<SerialDenseMatrix<OrdinalType,ScalarType> >& B)
00493 {
00494   // Check that these new vectors are consistent.
00495   TEST_FOR_EXCEPTION(B->numRows()!=X->numRows() || B->numCols() != X->numCols(), std::invalid_argument,
00496          "SerialDenseSolver<T>::setVectors: X and B are not the same size!");
00497   TEST_FOR_EXCEPTION(B->values()==0, std::invalid_argument,
00498          "SerialDenseSolver<T>::setVectors: B is an empty SerialDenseMatrix<T>!");
00499   TEST_FOR_EXCEPTION(X->values()==0, std::invalid_argument,
00500          "SerialDenseSolver<T>::setVectors: X is an empty SerialDenseMatrix<T>!");
00501   TEST_FOR_EXCEPTION(B->stride()<1, std::invalid_argument,
00502          "SerialDenseSolver<T>::setVectors: B has an invalid stride!");
00503   TEST_FOR_EXCEPTION(X->stride()<1, std::invalid_argument,
00504          "SerialDenseSolver<T>::setVectors: X has an invalid stride!");
00505 
00506   resetVectors(); 
00507   LHS_ = X;
00508   RHS_ = B;
00509   return(0);
00510 }
00511 //=============================================================================
00512 
00513 template<typename OrdinalType, typename ScalarType>
00514 void SerialDenseSolver<OrdinalType,ScalarType>::estimateSolutionErrors(bool flag) 
00515 {
00516   estimateSolutionErrors_ = flag;
00517 
00518   // If the errors are estimated, this implies that the solution must be refined
00519   refineSolution_ = refineSolution_ || flag;
00520 }
00521 //=============================================================================
00522 
00523 template<typename OrdinalType, typename ScalarType>
00524 int SerialDenseSolver<OrdinalType,ScalarType>::factor() {
00525 
00526   if (factored()) return(0); // Already factored
00527 
00528   TEST_FOR_EXCEPTION(inverted(), std::logic_error,
00529          "SerialDenseSolver<T>::factor: Cannot factor an inverted matrix!");
00530 
00531   ANORM_ = Matrix_->normOne(); // Compute 1-Norm of A
00532 
00533 
00534   // If we want to refine the solution, then the factor must
00535   // be stored separatedly from the original matrix
00536 
00537   if (A_ == AF_)
00538     if (refineSolution_ ) {
00539       Factor_ = rcp( new SerialDenseMatrix<OrdinalType,ScalarType>(*Matrix_) );
00540       AF_ = Factor_->values();
00541       LDAF_ = Factor_->stride();
00542     }
00543   
00544   int ierr = 0;
00545   if (equilibrate_) ierr = equilibrateMatrix();
00546 
00547   if (ierr!=0) return(ierr);
00548   
00549   if (IPIV_.size()==0) IPIV_.resize( Min_MN_ ); // Allocated Pivot vector if not already done.
00550 
00551   INFO_ = 0;
00552   this->GETRF(M_, N_, AF_, LDAF_, &IPIV_[0], &INFO_);
00553   factored_ = true;
00554 
00555   return(INFO_);
00556 }
00557 
00558 //=============================================================================
00559 
00560 template<typename OrdinalType, typename ScalarType>
00561 int SerialDenseSolver<OrdinalType,ScalarType>::solve() {
00562 
00563   // We will call one of four routines depending on what services the user wants and 
00564   // whether or not the matrix has been inverted or factored already.
00565   //
00566   // If the matrix has been inverted, use DGEMM to compute solution.
00567   // Otherwise, if the user want the matrix to be equilibrated or wants a refined solution, we will
00568   // call the X interface.
00569   // Otherwise, if the matrix is already factored we will call the TRS interface.
00570   // Otherwise, if the matrix is unfactored we will call the SV interface.
00571 
00572   int ierr = 0;
00573   if (equilibrate_) {
00574     ierr = equilibrateRHS();
00575     equilibratedB_ = true;
00576   }
00577   if (ierr != 0) return(ierr);  // Can't equilibrate B, so return.
00578 
00579   TEST_FOR_EXCEPTION( (equilibratedA_ && !equilibratedB_) || (!equilibratedA_ && equilibratedB_) , 
00580                      std::logic_error, "SerialDenseSolver<T>::solve: Matrix and vectors must be similarly scaled!");
00581   TEST_FOR_EXCEPTION( RHS_==Teuchos::null, std::invalid_argument, 
00582                      "SerialDenseSolver<T>::solve: No right-hand side vector (RHS) has been set for the linear system!");
00583   TEST_FOR_EXCEPTION( LHS_==Teuchos::null, std::invalid_argument, 
00584                      "SerialDenseSolver<T>::solve: No solution vector (LHS) has been set for the linear system!");
00585 
00586   if (shouldEquilibrate() && !equilibratedA_)
00587     std::cout << "WARNING!  SerialDenseSolver<T>::solve: System should be equilibrated!" << std::endl;
00588 
00589   if (inverted()) {
00590 
00591     TEST_FOR_EXCEPTION( RHS_->values() == LHS_->values(), std::invalid_argument, 
00592                         "SerialDenseSolver<T>::solve: X and B must be different vectors if matrix is inverted.");
00593 
00594     INFO_ = 0;
00595     this->GEMM(TRANS_, Teuchos::NO_TRANS, N_, RHS_->numCols(), N_, 1.0, AF_, LDAF_,
00596          RHS_->values(), RHS_->stride(), 0.0, LHS_->values(), LHS_->stride());
00597     if (INFO_!=0) return(INFO_);
00598     solved_ = true;
00599   }
00600   else {
00601 
00602     if (!factored()) factor(); // Matrix must be factored
00603     
00604     if (RHS_->values()!=LHS_->values()) {
00605        (*LHS_) = (*RHS_); // Copy B to X if needed
00606     }
00607     INFO_ = 0;
00608     this->GETRS(ETranspChar[TRANS_], N_, RHS_->numCols(), AF_, LDAF_, &IPIV_[0], LHS_->values(), LHS_->stride(), &INFO_);
00609     if (INFO_!=0) return(INFO_);
00610     solved_ = true;
00611 
00612   }
00613   int ierr1=0;
00614   if (refineSolution_ && !inverted()) ierr1 = applyRefinement();
00615   if (ierr1!=0) 
00616     return(ierr1);
00617   else
00618     return(ierr);
00619   
00620   if (equilibrate_) ierr1 = unequilibrateLHS();
00621   return(ierr1);
00622 }
00623 //=============================================================================
00624 
00625 template<typename OrdinalType, typename ScalarType>
00626 int SerialDenseSolver<OrdinalType,ScalarType>::applyRefinement()
00627 {
00628   TEST_FOR_EXCEPTION(!solved(), std::logic_error,
00629          "SerialDenseSolver<T>::applyRefinement: Must have an existing solution!");
00630   TEST_FOR_EXCEPTION(A_==AF_, std::logic_error,
00631          "SerialDenseSolver<T>::applyRefinement: Cannot apply refinement if no original copy of A!");
00632 
00633   OrdinalType NRHS = RHS_->numCols();
00634   FERR_.resize( NRHS );
00635   BERR_.resize( NRHS );
00636   allocateWORK();
00637   allocateIWORK();
00638   
00639   INFO_ = 0;
00640   this->GERFS(ETranspChar[TRANS_], N_, NRHS, A_, LDA_, AF_, LDAF_, &IPIV_[0],
00641         RHS_->values(), RHS_->stride(), LHS_->values(), LHS_->stride(), 
00642               &FERR_[0], &BERR_[0], &WORK_[0], &IWORK_[0], &INFO_);
00643   
00644   solutionErrorsEstimated_ = true;
00645   reciprocalConditionEstimated_ = true;
00646   solutionRefined_ = true;
00647   
00648   return(INFO_);
00649 
00650 }
00651 
00652 //=============================================================================
00653 
00654 template<typename OrdinalType, typename ScalarType>
00655 int SerialDenseSolver<OrdinalType,ScalarType>::computeEquilibrateScaling() 
00656 {
00657   if (R_.size()!=0) return(0); // Already computed
00658  
00659   R_.resize( M_ );
00660   C_.resize( N_ );
00661   
00662   INFO_ = 0;
00663   this->GEEQU (M_, N_, AF_, LDAF_, &R_[0], &C_[0], &ROWCND_, &COLCND_, &AMAX_, &INFO_);
00664 
00665   if (COLCND_<0.1*ScalarTraits<MagnitudeType>::one() || 
00666       ROWCND_<0.1*ScalarTraits<MagnitudeType>::one() || 
00667       AMAX_ < ScalarTraits<ScalarType>::rmin() || AMAX_ > ScalarTraits<ScalarType>::rmax() ) 
00668     shouldEquilibrate_ = true;
00669 
00670   return(INFO_);
00671 }
00672 
00673 //=============================================================================
00674 
00675 template<typename OrdinalType, typename ScalarType>
00676 int SerialDenseSolver<OrdinalType,ScalarType>::equilibrateMatrix()
00677 {
00678   OrdinalType i, j;
00679   int ierr = 0;
00680 
00681   if (equilibratedA_) return(0); // Already done.
00682   if (R_.size()==0) ierr = computeEquilibrateScaling(); // Compute R and C if needed.
00683   if (ierr!=0) return(ierr);     // If return value is not zero, then can't equilibrate.
00684   if (A_==AF_) {
00685     ScalarType * ptr;
00686     for (j=0; j<N_; j++) {
00687       ptr = A_ + j*LDA_;
00688       ScalarType s1 = C_[j];
00689       for (i=0; i<M_; i++) {
00690   *ptr = *ptr*s1*R_[i];
00691   ptr++;
00692       }
00693     }
00694   }
00695   else {
00696     ScalarType * ptr;
00697     ScalarType * ptr1;
00698     for (j=0; j<N_; j++) {
00699       ptr = A_ + j*LDA_;
00700       ptr1 = AF_ + j*LDAF_;
00701       ScalarType s1 = C_[j];
00702       for (i=0; i<M_; i++) {
00703   *ptr = *ptr*s1*R_[i];
00704   ptr++;
00705   *ptr1 = *ptr1*s1*R_[i];
00706   ptr1++;
00707       }
00708     }
00709   }
00710   
00711   equilibratedA_ = true;
00712   
00713   return(0);
00714 }
00715 
00716 //=============================================================================
00717 
00718 template<typename OrdinalType, typename ScalarType>
00719 int SerialDenseSolver<OrdinalType,ScalarType>::equilibrateRHS()
00720 {
00721   OrdinalType i, j;
00722   int ierr = 0;
00723 
00724   if (equilibratedB_) return(0); // Already done.
00725   if (R_.size()==0) ierr = computeEquilibrateScaling(); // Compute R and C if needed.
00726   if (ierr!=0) return(ierr);     // Can't count on R and C being computed.
00727 
00728   ScalarType * R_tmp = &R_[0];
00729   if (transpose_) R_tmp = &C_[0];
00730 
00731   OrdinalType LDB = RHS_->stride(), NRHS = RHS_->numCols();
00732   ScalarType * B = RHS_->values();
00733   ScalarType * ptr;
00734   for (j=0; j<NRHS; j++) {
00735     ptr = B + j*LDB;
00736     for (i=0; i<M_; i++) {
00737       *ptr = *ptr*R_tmp[i];
00738       ptr++;
00739     }
00740   }
00741   
00742   equilibratedB_ = true;
00743 
00744   return(0);
00745 }
00746 
00747 //=============================================================================
00748 
00749 template<typename OrdinalType, typename ScalarType>
00750 int SerialDenseSolver<OrdinalType,ScalarType>::unequilibrateLHS()
00751 {
00752   OrdinalType i, j;
00753 
00754   if (!equilibratedB_) return(0); // Nothing to do
00755 
00756   ScalarType * C_tmp = &C_[0];
00757   if (transpose_) C_tmp = &R_[0];
00758 
00759   OrdinalType LDX = RHS_->stride(), NRHS = RHS_->numCols();
00760   ScalarType * X = RHS_->values();
00761   ScalarType * ptr;
00762   for (j=0; j<NRHS; j++) {
00763     ptr = X + j*LDX;
00764     for (i=0; i<N_; i++) {
00765       *ptr = *ptr*C_tmp[i];
00766       ptr++;
00767     }
00768   }
00769 
00770   return(0);
00771 }
00772 
00773 //=============================================================================
00774 
00775 template<typename OrdinalType, typename ScalarType>
00776 int SerialDenseSolver<OrdinalType,ScalarType>::invert()
00777 {
00778   if (!factored()) factor(); // Need matrix factored.
00779 
00780   /* This section work with LAPACK Version 3.0 only 
00781   // Setting LWORK = -1 and calling GETRI will return optimal work space size in WORK_TMP
00782   OrdinalType LWORK_TMP = -1;
00783   ScalarType WORK_TMP;
00784   GETRI( N_, AF_, LDAF_, &IPIV_[0], &WORK_TMP, &LWORK_TMP, &INFO_);
00785   LWORK_TMP = WORK_TMP; // Convert to integer
00786   if (LWORK_TMP>LWORK_) {
00787   LWORK_ = LWORK_TMP;
00788   WORK_.resize( LWORK_ );
00789   }
00790   */
00791   // This section will work with any version of LAPACK 
00792   allocateWORK();
00793 
00794   INFO_ = 0;
00795   this->GETRI( N_, AF_, LDAF_, &IPIV_[0], &WORK_[0], LWORK_, &INFO_);
00796 
00797   inverted_ = true;
00798   factored_ = false;
00799   
00800   return(INFO_);
00801 }
00802 
00803 //=============================================================================
00804 
00805 template<typename OrdinalType, typename ScalarType>
00806 int SerialDenseSolver<OrdinalType,ScalarType>::reciprocalConditionEstimate(MagnitudeType & Value)
00807 {
00808   if (reciprocalConditionEstimated()) {
00809     Value = RCOND_;
00810     return(0); // Already computed, just return it.
00811   }
00812 
00813   if ( ANORM_<ScalarTraits<MagnitudeType>::zero() ) ANORM_ = Matrix_->normOne();
00814 
00815   int ierr = 0;
00816   if (!factored()) ierr = factor(); // Need matrix factored.
00817   if (ierr!=0) return(ierr);
00818 
00819   allocateWORK();
00820   allocateIWORK();
00821 
00822   // We will assume a one-norm condition number
00823   INFO_ = 0;
00824   this->GECON( '1', N_, AF_, LDAF_, ANORM_, &RCOND_, &WORK_[0], &IWORK_[0], &INFO_);
00825   reciprocalConditionEstimated_ = true;
00826   Value = RCOND_;
00827 
00828   return(INFO_);
00829 }
00830 //=============================================================================
00831 
00832 template<typename OrdinalType, typename ScalarType>
00833 void SerialDenseSolver<OrdinalType,ScalarType>::Print(std::ostream& os) const {
00834 
00835   if (Matrix_!=Teuchos::null) os << "Solver Matrix"          << std::endl << *Matrix_ << std::endl;
00836   if (Factor_!=Teuchos::null) os << "Solver Factored Matrix" << std::endl << *Factor_ << std::endl;
00837   if (LHS_   !=Teuchos::null) os << "Solver LHS"             << std::endl << *LHS_    << std::endl;
00838   if (RHS_   !=Teuchos::null) os << "Solver RHS"             << std::endl << *RHS_    << std::endl;
00839 
00840 }
00841 
00842 } // namespace Teuchos
00843 
00844 #endif /* _TEUCHOS_SERIALDENSESOLVER_HPP_ */
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Defines

Generated on Tue Oct 20 10:14:00 2009 for Teuchos Package Browser (Single Doxygen Collection) by  doxygen 1.6.1