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 
00129 namespace Teuchos {
00130 
00131   template<typename OrdinalType, typename ScalarType>
00132   class SerialDenseSolver : public CompObject, public Object, public BLAS<OrdinalType, ScalarType>,
00133           public LAPACK<OrdinalType, ScalarType>
00134   {
00135 
00136   public:
00137 
00138     typedef typename Teuchos::ScalarTraits<ScalarType>::magnitudeType MagnitudeType;
00139   
00141 
00142 
00143     SerialDenseSolver();
00144     
00146     virtual ~SerialDenseSolver();
00148     
00150 
00151     
00153     int setMatrix(const RCP<SerialDenseMatrix<OrdinalType, ScalarType> >& A);
00154 
00156 
00159     int setVectors(const RCP<SerialDenseMatrix<OrdinalType, ScalarType> >& X, 
00160        const RCP<SerialDenseMatrix<OrdinalType, ScalarType> >& B);
00162     
00164 
00165     
00167 
00169     void factorWithEquilibration(bool flag) {equilibrate_ = flag; return;}
00170     
00172     void solveWithTranspose(bool flag) {transpose_ = flag; if (flag) TRANS_ = Teuchos::TRANS; else TRANS_ = Teuchos::NO_TRANS; return;}
00173     
00175     void solveToRefinedSolution(bool flag) {refineSolution_ = flag; return;}
00176     
00178 
00181     void estimateSolutionErrors(bool flag);
00183     
00185 
00186     
00188 
00191     int factor();
00192     
00194 
00197     int solve();
00198     
00200 
00203     int invert();
00204     
00206 
00209     int computeEquilibrateScaling();
00210     
00212 
00215     int equilibrateMatrix();
00216     
00218 
00221     int equilibrateRHS();
00222         
00224 
00227     int applyRefinement();
00228     
00230 
00233     int unequilibrateLHS();
00234     
00236 
00242      int reciprocalConditionEstimate(MagnitudeType & Value);
00244     
00246 
00247     
00249     bool transpose() {return(transpose_);}
00250     
00252     bool factored() {return(factored_);}
00253     
00255     bool equilibratedA() {return(equilibratedA_);}
00256     
00258     bool equilibratedB() {return(equilibratedB_);}
00259     
00261      bool shouldEquilibrate() {computeEquilibrateScaling(); return(shouldEquilibrate_);}
00262     
00264     bool solutionErrorsEstimated() {return(solutionErrorsEstimated_);}
00265     
00267     bool inverted() {return(inverted_);}
00268     
00270     bool reciprocalConditionEstimated() {return(reciprocalConditionEstimated_);}
00271     
00273     bool solved() {return(solved_);}
00274     
00276     bool solutionRefined() {return(solutionRefined_);}
00278     
00280 
00281     
00283     RCP<SerialDenseMatrix<OrdinalType, ScalarType> > getMatrix()  const {return(Matrix_);}
00284     
00286     RCP<SerialDenseMatrix<OrdinalType, ScalarType> > getFactoredMatrix()  const {return(Factor_);}
00287     
00289     RCP<SerialDenseMatrix<OrdinalType, ScalarType> > getLHS() const {return(LHS_);}
00290     
00292     RCP<SerialDenseMatrix<OrdinalType, ScalarType> > getRHS() const {return(RHS_);}
00293     
00295     OrdinalType numRows()  const {return(M_);}
00296     
00298     OrdinalType numCols()  const {return(N_);}
00299     
00301     std::vector<OrdinalType> IPIV()  const {return(IPIV_);}
00302     
00304     MagnitudeType ANORM()  const {return(ANORM_);}
00305     
00307     MagnitudeType RCOND()  const {return(RCOND_);}
00308     
00310 
00312     MagnitudeType ROWCND()  const {return(ROWCND_);}
00313     
00315 
00317     MagnitudeType COLCND()  const {return(COLCND_);}
00318     
00320     MagnitudeType AMAX()  const {return(AMAX_);}
00321     
00323     std::vector<MagnitudeType> FERR()  const {return(FERR_);}
00324     
00326     std::vector<MagnitudeType> BERR()  const {return(BERR_);}
00327     
00329     std::vector<MagnitudeType> R()  const {return(R_);}
00330     
00332     std::vector<MagnitudeType> C()  const {return(C_);}
00334     
00336 
00337 
00338     void Print(std::ostream& os) const;
00340   protected:
00341     
00342     void allocateWORK() { LWORK_ = 4*N_; WORK_.resize( LWORK_ ); return;}
00343     void allocateIWORK() { IWORK_.resize( N_ ); return;}
00344     void resetMatrix();
00345     void resetVectors();
00346     
00347     
00348     bool equilibrate_;
00349     bool shouldEquilibrate_;
00350     bool equilibratedA_;
00351     bool equilibratedB_;
00352     bool transpose_;
00353     bool factored_;
00354     bool estimateSolutionErrors_;
00355     bool solutionErrorsEstimated_;
00356     bool solved_;
00357     bool inverted_;
00358     bool reciprocalConditionEstimated_;
00359     bool refineSolution_;
00360     bool solutionRefined_;
00361     
00362     Teuchos::ETransp TRANS_;
00363     
00364     OrdinalType M_;
00365     OrdinalType N_;
00366     OrdinalType Min_MN_;
00367     OrdinalType LDA_;
00368     OrdinalType LDAF_;
00369     OrdinalType INFO_;
00370     OrdinalType LWORK_;
00371     
00372     std::vector<OrdinalType> IPIV_;
00373     std::vector<int> IWORK_;
00374     
00375     MagnitudeType ANORM_;
00376     MagnitudeType RCOND_;
00377     MagnitudeType ROWCND_;
00378     MagnitudeType COLCND_;
00379     MagnitudeType AMAX_;
00380     
00381     RCP<SerialDenseMatrix<OrdinalType, ScalarType> > Matrix_;
00382     RCP<SerialDenseMatrix<OrdinalType, ScalarType> > LHS_;
00383     RCP<SerialDenseMatrix<OrdinalType, ScalarType> > RHS_;
00384     RCP<SerialDenseMatrix<OrdinalType, ScalarType> > Factor_;
00385     
00386     ScalarType * A_;
00387     ScalarType * AF_;
00388     std::vector<MagnitudeType> FERR_;
00389     std::vector<MagnitudeType> BERR_;
00390     std::vector<ScalarType> WORK_;
00391     std::vector<MagnitudeType> R_;
00392     std::vector<MagnitudeType> C_;
00393     
00394   private:
00395     // SerialDenseSolver copy constructor (put here because we don't want user access)
00396     
00397     SerialDenseSolver(const SerialDenseSolver<OrdinalType, ScalarType>& Source);
00398     SerialDenseSolver & operator=(const SerialDenseSolver<OrdinalType, ScalarType>& Source);
00399     
00400   };  
00401 
00402 //=============================================================================
00403 
00404 template<typename OrdinalType, typename ScalarType>
00405 SerialDenseSolver<OrdinalType,ScalarType>::SerialDenseSolver()
00406   : CompObject(),
00407     Object("Teuchos::SerialDenseSolver"),
00408     equilibrate_(false),
00409     shouldEquilibrate_(false),
00410     equilibratedA_(false),
00411     equilibratedB_(false),
00412     transpose_(false),
00413     factored_(false),
00414     estimateSolutionErrors_(false),
00415     solutionErrorsEstimated_(false),
00416     solved_(false),
00417     inverted_(false),
00418     reciprocalConditionEstimated_(false),
00419     refineSolution_(false),
00420     solutionRefined_(false),
00421     TRANS_(Teuchos::NO_TRANS),
00422     M_(0),
00423     N_(0),
00424     Min_MN_(0),
00425     LDA_(0),
00426     LDAF_(0),
00427     INFO_(0),
00428     LWORK_(0),
00429     ANORM_(ScalarTraits<MagnitudeType>::zero()),
00430     RCOND_(ScalarTraits<MagnitudeType>::zero()),
00431     ROWCND_(ScalarTraits<MagnitudeType>::zero()),
00432     COLCND_(ScalarTraits<MagnitudeType>::zero()),
00433     AMAX_(ScalarTraits<MagnitudeType>::zero()),
00434     A_(0),
00435     AF_(0)
00436 {
00437   resetMatrix();
00438 }
00439 //=============================================================================
00440 
00441 template<typename OrdinalType, typename ScalarType>
00442 SerialDenseSolver<OrdinalType,ScalarType>::~SerialDenseSolver()
00443 {}
00444 
00445 //=============================================================================
00446 
00447 template<typename OrdinalType, typename ScalarType>
00448 void SerialDenseSolver<OrdinalType,ScalarType>::resetVectors()
00449 {
00450   LHS_ = Teuchos::null;
00451   RHS_ = Teuchos::null;
00452   reciprocalConditionEstimated_ = false;
00453   solutionRefined_ = false;
00454   solved_ = false;
00455   solutionErrorsEstimated_ = false;
00456   equilibratedB_ = false;
00457 }
00458 //=============================================================================
00459 
00460 template<typename OrdinalType, typename ScalarType>
00461 void SerialDenseSolver<OrdinalType,ScalarType>::resetMatrix()
00462 {
00463   resetVectors();
00464   equilibratedA_ = false;
00465   factored_ = false;
00466   inverted_ = false;
00467   M_ = 0;
00468   N_ = 0;
00469   Min_MN_ = 0;
00470   LDA_ = 0;
00471   LDAF_ = 0;
00472   ANORM_ = -ScalarTraits<MagnitudeType>::one();
00473   RCOND_ = -ScalarTraits<MagnitudeType>::one();
00474   ROWCND_ = -ScalarTraits<MagnitudeType>::one();
00475   COLCND_ = -ScalarTraits<MagnitudeType>::one();
00476   AMAX_ = -ScalarTraits<MagnitudeType>::one();
00477   A_ = 0;
00478   AF_ = 0;
00479   INFO_ = 0;
00480   LWORK_ = 0;
00481   R_.resize(0);
00482   C_.resize(0);
00483 }
00484 //=============================================================================
00485 
00486 template<typename OrdinalType, typename ScalarType>
00487 int SerialDenseSolver<OrdinalType,ScalarType>::setMatrix(const RCP<SerialDenseMatrix<OrdinalType,ScalarType> >& A) 
00488 {
00489   resetMatrix();
00490   Matrix_ = A;
00491   Factor_ = A;
00492   M_ = A->numRows();
00493   N_ = A->numCols();
00494   Min_MN_ = TEUCHOS_MIN(M_,N_);
00495   LDA_ = A->stride();
00496   LDAF_ = LDA_;
00497   A_ = A->values();
00498   AF_ = A->values();
00499   return(0);
00500 }
00501 //=============================================================================
00502 
00503 template<typename OrdinalType, typename ScalarType>
00504 int SerialDenseSolver<OrdinalType,ScalarType>::setVectors(const RCP<SerialDenseMatrix<OrdinalType,ScalarType> >& X, 
00505                  const RCP<SerialDenseMatrix<OrdinalType,ScalarType> >& B)
00506 {
00507   // Check that these new vectors are consistent.
00508   TEST_FOR_EXCEPTION(B->numRows()!=X->numRows() || B->numCols() != X->numCols(), std::invalid_argument,
00509          "SerialDenseSolver<T>::setVectors: X and B are not the same size!");
00510   TEST_FOR_EXCEPTION(B->values()==0, std::invalid_argument,
00511          "SerialDenseSolver<T>::setVectors: B is an empty SerialDenseMatrix<T>!");
00512   TEST_FOR_EXCEPTION(X->values()==0, std::invalid_argument,
00513          "SerialDenseSolver<T>::setVectors: X is an empty SerialDenseMatrix<T>!");
00514   TEST_FOR_EXCEPTION(B->stride()<1, std::invalid_argument,
00515          "SerialDenseSolver<T>::setVectors: B has an invalid stride!");
00516   TEST_FOR_EXCEPTION(X->stride()<1, std::invalid_argument,
00517          "SerialDenseSolver<T>::setVectors: X has an invalid stride!");
00518 
00519   resetVectors(); 
00520   LHS_ = X;
00521   RHS_ = B;
00522   return(0);
00523 }
00524 //=============================================================================
00525 
00526 template<typename OrdinalType, typename ScalarType>
00527 void SerialDenseSolver<OrdinalType,ScalarType>::estimateSolutionErrors(bool flag) 
00528 {
00529   estimateSolutionErrors_ = flag;
00530 
00531   // If the errors are estimated, this implies that the solution must be refined
00532   refineSolution_ = refineSolution_ || flag;
00533 }
00534 //=============================================================================
00535 
00536 template<typename OrdinalType, typename ScalarType>
00537 int SerialDenseSolver<OrdinalType,ScalarType>::factor() {
00538 
00539   if (factored()) return(0); // Already factored
00540 
00541   TEST_FOR_EXCEPTION(inverted(), std::logic_error,
00542          "SerialDenseSolver<T>::factor: Cannot factor an inverted matrix!");
00543 
00544   ANORM_ = Matrix_->normOne(); // Compute 1-Norm of A
00545 
00546 
00547   // If we want to refine the solution, then the factor must
00548   // be stored separatedly from the original matrix
00549 
00550   if (A_ == AF_)
00551     if (refineSolution_ ) {
00552       Factor_ = rcp( new SerialDenseMatrix<OrdinalType,ScalarType>(*Matrix_) );
00553       AF_ = Factor_->values();
00554       LDAF_ = Factor_->stride();
00555     }
00556   
00557   int ierr = 0;
00558   if (equilibrate_) ierr = equilibrateMatrix();
00559 
00560   if (ierr!=0) return(ierr);
00561   
00562   if (IPIV_.size()==0) IPIV_.resize( Min_MN_ ); // Allocated Pivot vector if not already done.
00563 
00564   INFO_ = 0;
00565   this->GETRF(M_, N_, AF_, LDAF_, &IPIV_[0], &INFO_);
00566   factored_ = true;
00567 
00568   return(INFO_);
00569 }
00570 
00571 //=============================================================================
00572 
00573 template<typename OrdinalType, typename ScalarType>
00574 int SerialDenseSolver<OrdinalType,ScalarType>::solve() {
00575 
00576   // We will call one of four routines depending on what services the user wants and 
00577   // whether or not the matrix has been inverted or factored already.
00578   //
00579   // If the matrix has been inverted, use DGEMM to compute solution.
00580   // Otherwise, if the user want the matrix to be equilibrated or wants a refined solution, we will
00581   // call the X interface.
00582   // Otherwise, if the matrix is already factored we will call the TRS interface.
00583   // Otherwise, if the matrix is unfactored we will call the SV interface.
00584 
00585   int ierr = 0;
00586   if (equilibrate_) {
00587     ierr = equilibrateRHS();
00588     equilibratedB_ = true;
00589   }
00590   if (ierr != 0) return(ierr);  // Can't equilibrate B, so return.
00591 
00592   TEST_FOR_EXCEPTION( (equilibratedA_ && !equilibratedB_) || (!equilibratedA_ && equilibratedB_) , 
00593                      std::logic_error, "SerialDenseSolver<T>::solve: Matrix and vectors must be similarly scaled!");
00594   TEST_FOR_EXCEPTION( RHS_==Teuchos::null, std::invalid_argument, 
00595                      "SerialDenseSolver<T>::solve: No right-hand side vector (RHS) has been set for the linear system!");
00596   TEST_FOR_EXCEPTION( LHS_==Teuchos::null, std::invalid_argument, 
00597                      "SerialDenseSolver<T>::solve: No solution vector (LHS) has been set for the linear system!");
00598 
00599   if (shouldEquilibrate() && !equilibratedA_)
00600     std::cout << "WARNING!  SerialDenseSolver<T>::solve: System should be equilibrated!" << std::endl;
00601 
00602   if (inverted()) {
00603 
00604     TEST_FOR_EXCEPTION( RHS_->values() == LHS_->values(), std::invalid_argument, 
00605                         "SerialDenseSolver<T>::solve: X and B must be different vectors if matrix is inverted.");
00606 
00607     INFO_ = 0;
00608     this->GEMM(TRANS_, Teuchos::NO_TRANS, N_, RHS_->numCols(), N_, 1.0, AF_, LDAF_,
00609          RHS_->values(), RHS_->stride(), 0.0, LHS_->values(), LHS_->stride());
00610     if (INFO_!=0) return(INFO_);
00611     solved_ = true;
00612   }
00613   else {
00614 
00615     if (!factored()) factor(); // Matrix must be factored
00616     
00617     if (RHS_->values()!=LHS_->values()) {
00618        (*LHS_) = (*RHS_); // Copy B to X if needed
00619     }
00620     INFO_ = 0;
00621     this->GETRS(ETranspChar[TRANS_], N_, RHS_->numCols(), AF_, LDAF_, &IPIV_[0], LHS_->values(), LHS_->stride(), &INFO_);
00622     if (INFO_!=0) return(INFO_);
00623     solved_ = true;
00624 
00625   }
00626   int ierr1=0;
00627   if (refineSolution_ && !inverted()) ierr1 = applyRefinement();
00628   if (ierr1!=0) 
00629     return(ierr1);
00630   
00631   if (equilibrate_) ierr1 = unequilibrateLHS();
00632   return(ierr1);
00633 }
00634 //=============================================================================
00635 
00636 template<typename OrdinalType, typename ScalarType>
00637 int SerialDenseSolver<OrdinalType,ScalarType>::applyRefinement()
00638 {
00639   TEST_FOR_EXCEPTION(!solved(), std::logic_error,
00640          "SerialDenseSolver<T>::applyRefinement: Must have an existing solution!");
00641   TEST_FOR_EXCEPTION(A_==AF_, std::logic_error,
00642          "SerialDenseSolver<T>::applyRefinement: Cannot apply refinement if no original copy of A!");
00643 
00644   OrdinalType NRHS = RHS_->numCols();
00645   FERR_.resize( NRHS );
00646   BERR_.resize( NRHS );
00647   allocateWORK();
00648   allocateIWORK();
00649   
00650   INFO_ = 0;
00651   this->GERFS(ETranspChar[TRANS_], N_, NRHS, A_, LDA_, AF_, LDAF_, &IPIV_[0],
00652         RHS_->values(), RHS_->stride(), LHS_->values(), LHS_->stride(), 
00653               &FERR_[0], &BERR_[0], &WORK_[0], &IWORK_[0], &INFO_);
00654   
00655   solutionErrorsEstimated_ = true;
00656   reciprocalConditionEstimated_ = true;
00657   solutionRefined_ = true;
00658   
00659   return(INFO_);
00660 
00661 }
00662 
00663 //=============================================================================
00664 
00665 template<typename OrdinalType, typename ScalarType>
00666 int SerialDenseSolver<OrdinalType,ScalarType>::computeEquilibrateScaling() 
00667 {
00668   if (R_.size()!=0) return(0); // Already computed
00669  
00670   R_.resize( M_ );
00671   C_.resize( N_ );
00672   
00673   INFO_ = 0;
00674   this->GEEQU (M_, N_, AF_, LDAF_, &R_[0], &C_[0], &ROWCND_, &COLCND_, &AMAX_, &INFO_);
00675 
00676   if (COLCND_<0.1*ScalarTraits<MagnitudeType>::one() || 
00677       ROWCND_<0.1*ScalarTraits<MagnitudeType>::one() || 
00678       AMAX_ < ScalarTraits<ScalarType>::rmin() || AMAX_ > ScalarTraits<ScalarType>::rmax() ) 
00679     shouldEquilibrate_ = true;
00680 
00681   return(INFO_);
00682 }
00683 
00684 //=============================================================================
00685 
00686 template<typename OrdinalType, typename ScalarType>
00687 int SerialDenseSolver<OrdinalType,ScalarType>::equilibrateMatrix()
00688 {
00689   OrdinalType i, j;
00690   int ierr = 0;
00691 
00692   if (equilibratedA_) return(0); // Already done.
00693   if (R_.size()==0) ierr = computeEquilibrateScaling(); // Compute R and C if needed.
00694   if (ierr!=0) return(ierr);     // If return value is not zero, then can't equilibrate.
00695   if (A_==AF_) {
00696     ScalarType * ptr;
00697     for (j=0; j<N_; j++) {
00698       ptr = A_ + j*LDA_;
00699       ScalarType s1 = C_[j];
00700       for (i=0; i<M_; i++) {
00701   *ptr = *ptr*s1*R_[i];
00702   ptr++;
00703       }
00704     }
00705   }
00706   else {
00707     ScalarType * ptr;
00708     ScalarType * ptr1;
00709     for (j=0; j<N_; j++) {
00710       ptr = A_ + j*LDA_;
00711       ptr1 = AF_ + j*LDAF_;
00712       ScalarType s1 = C_[j];
00713       for (i=0; i<M_; i++) {
00714   *ptr = *ptr*s1*R_[i];
00715   ptr++;
00716   *ptr1 = *ptr1*s1*R_[i];
00717   ptr1++;
00718       }
00719     }
00720   }
00721   
00722   equilibratedA_ = true;
00723   
00724   return(0);
00725 }
00726 
00727 //=============================================================================
00728 
00729 template<typename OrdinalType, typename ScalarType>
00730 int SerialDenseSolver<OrdinalType,ScalarType>::equilibrateRHS()
00731 {
00732   OrdinalType i, j;
00733   int ierr = 0;
00734 
00735   if (equilibratedB_) return(0); // Already done.
00736   if (R_.size()==0) ierr = computeEquilibrateScaling(); // Compute R and C if needed.
00737   if (ierr!=0) return(ierr);     // Can't count on R and C being computed.
00738 
00739   MagnitudeType * R_tmp = &R_[0];
00740   if (transpose_) R_tmp = &C_[0];
00741 
00742   OrdinalType LDB = RHS_->stride(), NRHS = RHS_->numCols();
00743   ScalarType * B = RHS_->values();
00744   ScalarType * ptr;
00745   for (j=0; j<NRHS; j++) {
00746     ptr = B + j*LDB;
00747     for (i=0; i<M_; i++) {
00748       *ptr = *ptr*R_tmp[i];
00749       ptr++;
00750     }
00751   }
00752   
00753   equilibratedB_ = true;
00754 
00755   return(0);
00756 }
00757 
00758 //=============================================================================
00759 
00760 template<typename OrdinalType, typename ScalarType>
00761 int SerialDenseSolver<OrdinalType,ScalarType>::unequilibrateLHS()
00762 {
00763   OrdinalType i, j;
00764 
00765   if (!equilibratedB_) return(0); // Nothing to do
00766 
00767   MagnitudeType * C_tmp = &C_[0];
00768   if (transpose_) C_tmp = &R_[0];
00769 
00770   OrdinalType LDX = LHS_->stride(), NLHS = LHS_->numCols();
00771   ScalarType * X = LHS_->values();
00772   ScalarType * ptr;
00773   for (j=0; j<NLHS; j++) {
00774     ptr = X + j*LDX;
00775     for (i=0; i<N_; i++) {
00776       *ptr = *ptr*C_tmp[i];
00777       ptr++;
00778     }
00779   }
00780 
00781   return(0);
00782 }
00783 
00784 //=============================================================================
00785 
00786 template<typename OrdinalType, typename ScalarType>
00787 int SerialDenseSolver<OrdinalType,ScalarType>::invert()
00788 {
00789   if (!factored()) factor(); // Need matrix factored.
00790 
00791   /* This section work with LAPACK Version 3.0 only 
00792   // Setting LWORK = -1 and calling GETRI will return optimal work space size in WORK_TMP
00793   OrdinalType LWORK_TMP = -1;
00794   ScalarType WORK_TMP;
00795   GETRI( N_, AF_, LDAF_, &IPIV_[0], &WORK_TMP, &LWORK_TMP, &INFO_);
00796   LWORK_TMP = WORK_TMP; // Convert to integer
00797   if (LWORK_TMP>LWORK_) {
00798   LWORK_ = LWORK_TMP;
00799   WORK_.resize( LWORK_ );
00800   }
00801   */
00802   // This section will work with any version of LAPACK 
00803   allocateWORK();
00804 
00805   INFO_ = 0;
00806   this->GETRI( N_, AF_, LDAF_, &IPIV_[0], &WORK_[0], LWORK_, &INFO_);
00807 
00808   inverted_ = true;
00809   factored_ = false;
00810   
00811   return(INFO_);
00812 }
00813 
00814 //=============================================================================
00815 
00816 template<typename OrdinalType, typename ScalarType>
00817 int SerialDenseSolver<OrdinalType,ScalarType>::reciprocalConditionEstimate(MagnitudeType & Value)
00818 {
00819   if (reciprocalConditionEstimated()) {
00820     Value = RCOND_;
00821     return(0); // Already computed, just return it.
00822   }
00823 
00824   if ( ANORM_<ScalarTraits<MagnitudeType>::zero() ) ANORM_ = Matrix_->normOne();
00825 
00826   int ierr = 0;
00827   if (!factored()) ierr = factor(); // Need matrix factored.
00828   if (ierr!=0) return(ierr);
00829 
00830   allocateWORK();
00831   allocateIWORK();
00832 
00833   // We will assume a one-norm condition number
00834   INFO_ = 0;
00835   this->GECON( '1', N_, AF_, LDAF_, ANORM_, &RCOND_, &WORK_[0], &IWORK_[0], &INFO_);
00836   reciprocalConditionEstimated_ = true;
00837   Value = RCOND_;
00838 
00839   return(INFO_);
00840 }
00841 //=============================================================================
00842 
00843 template<typename OrdinalType, typename ScalarType>
00844 void SerialDenseSolver<OrdinalType,ScalarType>::Print(std::ostream& os) const {
00845 
00846   if (Matrix_!=Teuchos::null) os << "Solver Matrix"          << std::endl << *Matrix_ << std::endl;
00847   if (Factor_!=Teuchos::null) os << "Solver Factored Matrix" << std::endl << *Factor_ << std::endl;
00848   if (LHS_   !=Teuchos::null) os << "Solver LHS"             << std::endl << *LHS_    << std::endl;
00849   if (RHS_   !=Teuchos::null) os << "Solver RHS"             << std::endl << *RHS_    << std::endl;
00850 
00851 }
00852 
00853 } // namespace Teuchos
00854 
00855 #endif /* _TEUCHOS_SERIALDENSESOLVER_HPP_ */
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Defines