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