Teuchos_SerialSpdDenseSolver.hpp

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

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