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