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