Teuchos - Trilinos Tools Package 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_SerialDenseMatrix.hpp"
00057 
00125 namespace Teuchos {
00126 
00127   template<typename OrdinalType, typename ScalarType>
00128   class SerialSpdDenseSolver : public CompObject, public Object, public BLAS<OrdinalType, ScalarType>,
00129           public LAPACK<OrdinalType, ScalarType>
00130   {
00131   public:
00132     
00133     typedef typename Teuchos::ScalarTraits<ScalarType>::magnitudeType MagnitudeType;
00134     
00136 
00137 
00138     SerialSpdDenseSolver();
00139     
00140     
00142     virtual ~SerialSpdDenseSolver();
00144     
00146 
00147     
00149     int setMatrix(const RCP<SerialSymDenseMatrix<OrdinalType,ScalarType> >& A_in);
00150     
00152 
00155     int setVectors(const RCP<SerialDenseMatrix<OrdinalType, ScalarType> >& X, 
00156        const RCP<SerialDenseMatrix<OrdinalType, ScalarType> >& B);
00158     
00160 
00161     
00163 
00165     void factorWithEquilibration(bool flag) {equilibrate_ = flag; return;}
00166     
00168     void solveToRefinedSolution(bool flag) {refineSolution_ = flag; return;}
00169     
00171 
00174     void estimateSolutionErrors(bool flag);
00176     
00178 
00179     
00181 
00184     int factor();
00185     
00187 
00190     int solve();
00191     
00193 
00198     int invert();
00199     
00201 
00204     int computeEquilibrateScaling();
00205     
00207 
00210     int equilibrateMatrix();
00211     
00213 
00216     int equilibrateRHS();
00217     
00218     
00220 
00223     int applyRefinement();
00224     
00226 
00229     int unequilibrateLHS();
00230     
00232 
00238     int reciprocalConditionEstimate(MagnitudeType & Value);
00240     
00242 
00243     
00245     bool transpose() {return(transpose_);}
00246     
00248     bool factored() {return(factored_);}
00249     
00251     bool equilibratedA() {return(equilibratedA_);}
00252     
00254     bool equilibratedB() {return(equilibratedB_);}
00255     
00257     bool shouldEquilibrate() {computeEquilibrateScaling(); return(shouldEquilibrate_);}
00258     
00260     bool solutionErrorsEstimated() {return(solutionErrorsEstimated_);}
00261     
00263     bool inverted() {return(inverted_);}
00264     
00266     bool reciprocalConditionEstimated() {return(reciprocalConditionEstimated_);}
00267     
00269     bool solved() {return(solved_);}
00270     
00272     bool solutionRefined() {return(solutionRefined_);}
00274     
00276 
00277     
00279     RCP<SerialSymDenseMatrix<OrdinalType, ScalarType> > getMatrix()  const {return(Matrix_);}
00280     
00282     RCP<SerialSymDenseMatrix<OrdinalType, ScalarType> > getFactoredMatrix()  const {return(Factor_);}
00283     
00285     RCP<SerialDenseMatrix<OrdinalType, ScalarType> > getLHS() const {return(LHS_);}
00286     
00288     RCP<SerialDenseMatrix<OrdinalType, ScalarType> > getRHS() const {return(RHS_);}
00289     
00291     OrdinalType numRows()  const {return(numRowCols_);}
00292     
00294     OrdinalType numCols()  const {return(numRowCols_);}
00295     
00297     MagnitudeType ANORM()  const {return(ANORM_);}
00298     
00300     MagnitudeType RCOND()  const {return(RCOND_);}
00301     
00303 
00305     MagnitudeType SCOND() {return(SCOND_);};
00306     
00308     MagnitudeType AMAX()  const {return(AMAX_);}
00309     
00311     std::vector<ScalarType> FERR()  const {return(FERR_);}
00312     
00314     std::vector<ScalarType> BERR()  const {return(BERR_);}
00315     
00317     std::vector<ScalarType> R()  const {return(R_);}
00318     
00320 
00322 
00323 
00324     void Print(std::ostream& os) const;
00326     
00327   protected:
00328     
00329     void allocateWORK() { LWORK_ = 4*numRowCols_; WORK_.resize( LWORK_ ); return;}
00330     void allocateIWORK() { IWORK_.resize( numRowCols_ ); return;}
00331     void resetMatrix();
00332     void resetVectors();
00333     
00334     bool equilibrate_;
00335     bool shouldEquilibrate_;
00336     bool equilibratedA_;
00337     bool equilibratedB_;
00338     bool transpose_;
00339     bool factored_;
00340     bool estimateSolutionErrors_;
00341     bool solutionErrorsEstimated_;
00342     bool solved_;
00343     bool inverted_;
00344     bool reciprocalConditionEstimated_;
00345     bool refineSolution_;
00346     bool solutionRefined_;
00347     
00348     OrdinalType numRowCols_;
00349     OrdinalType LDA_;
00350     OrdinalType LDAF_;
00351     OrdinalType INFO_;
00352     OrdinalType LWORK_;
00353     
00354     std::vector<int> IWORK_;
00355     
00356     MagnitudeType ANORM_;
00357     MagnitudeType RCOND_;
00358     MagnitudeType SCOND_;
00359     MagnitudeType AMAX_;
00360     
00361     RCP<SerialSymDenseMatrix<OrdinalType, ScalarType> > Matrix_;
00362     RCP<SerialDenseMatrix<OrdinalType, ScalarType> > LHS_;
00363     RCP<SerialDenseMatrix<OrdinalType, ScalarType> > RHS_;
00364     RCP<SerialSymDenseMatrix<OrdinalType, ScalarType> > Factor_;
00365     
00366     ScalarType * A_;
00367     ScalarType * AF_;
00368     std::vector<ScalarType> FERR_;
00369     std::vector<ScalarType> BERR_;
00370     std::vector<ScalarType> WORK_;
00371     std::vector<ScalarType> R_;
00372     
00373   private:
00374     // SerialSpdDenseSolver copy constructor (put here because we don't want user access)
00375     
00376     SerialSpdDenseSolver(const SerialSpdDenseSolver<OrdinalType, ScalarType>& Source);
00377     SerialSpdDenseSolver & operator=(const SerialSpdDenseSolver<OrdinalType, ScalarType>& Source);
00378     
00379   };  
00380 
00381 //=============================================================================
00382 
00383 template<typename OrdinalType, typename ScalarType>
00384 SerialSpdDenseSolver<OrdinalType,ScalarType>::SerialSpdDenseSolver()
00385   : CompObject(),
00386     Object("Teuchos::SerialSpdDenseSolver"),
00387     equilibrate_(false),
00388     shouldEquilibrate_(false),
00389     equilibratedA_(false),
00390     equilibratedB_(false),
00391     transpose_(false),
00392     factored_(false),
00393     estimateSolutionErrors_(false),
00394     solutionErrorsEstimated_(false),
00395     solved_(false),
00396     inverted_(false),
00397     reciprocalConditionEstimated_(false),
00398     refineSolution_(false),
00399     solutionRefined_(false),
00400     numRowCols_(0),
00401     LDA_(0),
00402     LDAF_(0),
00403     INFO_(0),
00404     LWORK_(0),
00405     ANORM_(ScalarTraits<ScalarType>::zero()),
00406     RCOND_(ScalarTraits<ScalarType>::zero()),
00407     SCOND_(ScalarTraits<ScalarType>::zero()),
00408     AMAX_(ScalarTraits<ScalarType>::zero()),
00409     A_(0),
00410     AF_(0)
00411 {
00412   resetMatrix();
00413 }
00414 //=============================================================================
00415 
00416 template<typename OrdinalType, typename ScalarType>
00417 SerialSpdDenseSolver<OrdinalType,ScalarType>::~SerialSpdDenseSolver()
00418 {}
00419 
00420 //=============================================================================
00421 
00422 template<typename OrdinalType, typename ScalarType>
00423 void SerialSpdDenseSolver<OrdinalType,ScalarType>::resetVectors()
00424 {
00425   LHS_ = Teuchos::null;
00426   RHS_ = Teuchos::null;
00427   reciprocalConditionEstimated_ = false;
00428   solutionRefined_ = false;
00429   solved_ = false;
00430   solutionErrorsEstimated_ = false;
00431   equilibratedB_ = false;
00432 }
00433 //=============================================================================
00434 
00435 template<typename OrdinalType, typename ScalarType>
00436 void SerialSpdDenseSolver<OrdinalType,ScalarType>::resetMatrix()
00437 {
00438   resetVectors();
00439   equilibratedA_ = false;
00440   factored_ = false;
00441   inverted_ = false;
00442   numRowCols_ = 0;
00443   LDA_ = 0;
00444   LDAF_ = 0;
00445   ANORM_ = -ScalarTraits<MagnitudeType>::one();
00446   RCOND_ = -ScalarTraits<MagnitudeType>::one();
00447   SCOND_ = -ScalarTraits<MagnitudeType>::one();
00448   AMAX_ = -ScalarTraits<MagnitudeType>::one();
00449   A_ = 0;
00450   AF_ = 0;
00451   INFO_ = 0;
00452   LWORK_ = 0;
00453   R_.resize(0);
00454 }
00455 //=============================================================================
00456 
00457 template<typename OrdinalType, typename ScalarType>
00458 int SerialSpdDenseSolver<OrdinalType,ScalarType>::setMatrix(const RCP<SerialSymDenseMatrix<OrdinalType,ScalarType> >& A) 
00459 {
00460   resetMatrix();
00461   Matrix_ = A;
00462   Factor_ = A;
00463   numRowCols_ = A->numRows();
00464   LDA_ = A->stride();
00465   LDAF_ = LDA_;
00466   A_ = A->values();
00467   AF_ = A->values();
00468   return(0);
00469 }
00470 //=============================================================================
00471 
00472 template<typename OrdinalType, typename ScalarType>
00473 int SerialSpdDenseSolver<OrdinalType,ScalarType>::setVectors(const RCP<SerialDenseMatrix<OrdinalType,ScalarType> >& X, 
00474                    const RCP<SerialDenseMatrix<OrdinalType,ScalarType> >& B)
00475 {
00476   // Check that these new vectors are consistent.
00477   TEUCHOS_TEST_FOR_EXCEPTION(B->numRows()!=X->numRows() || B->numCols() != X->numCols(), std::invalid_argument,
00478          "SerialSpdDenseSolver<T>::setVectors: X and B are not the same size!");
00479   TEUCHOS_TEST_FOR_EXCEPTION(B->values()==0, std::invalid_argument,
00480          "SerialSpdDenseSolver<T>::setVectors: B is an empty SerialDenseMatrix<T>!");
00481   TEUCHOS_TEST_FOR_EXCEPTION(X->values()==0, std::invalid_argument,
00482          "SerialSpdDenseSolver<T>::setVectors: X is an empty SerialDenseMatrix<T>!");
00483   TEUCHOS_TEST_FOR_EXCEPTION(B->stride()<1, std::invalid_argument,
00484          "SerialSpdDenseSolver<T>::setVectors: B has an invalid stride!");
00485   TEUCHOS_TEST_FOR_EXCEPTION(X->stride()<1, std::invalid_argument,
00486          "SerialSpdDenseSolver<T>::setVectors: X has an invalid stride!");
00487 
00488   resetVectors(); 
00489   LHS_ = X;
00490   RHS_ = B;
00491   return(0);
00492 }
00493 //=============================================================================
00494 
00495 template<typename OrdinalType, typename ScalarType>
00496 void SerialSpdDenseSolver<OrdinalType,ScalarType>::estimateSolutionErrors(bool flag) 
00497 {
00498   estimateSolutionErrors_ = flag;
00499 
00500   // If the errors are estimated, this implies that the solution must be refined
00501   refineSolution_ = refineSolution_ || flag;
00502 }
00503 //=============================================================================
00504 
00505 template<typename OrdinalType, typename ScalarType>
00506 int SerialSpdDenseSolver<OrdinalType,ScalarType>::factor() {
00507 
00508   if (factored()) return(0); // Already factored
00509 
00510   TEUCHOS_TEST_FOR_EXCEPTION(inverted(), std::logic_error,
00511          "SerialSpdDenseSolver<T>::factor: Cannot factor an inverted matrix!");
00512 
00513   ANORM_ = Matrix_->normOne(); // Compute 1-Norm of A
00514 
00515 
00516   // If we want to refine the solution, then the factor must
00517   // be stored separatedly from the original matrix
00518 
00519   if (A_ == AF_)
00520     if (refineSolution_ ) {
00521       Factor_ = rcp( new SerialSymDenseMatrix<OrdinalType,ScalarType>(*Matrix_) );
00522       AF_ = Factor_->values();
00523       LDAF_ = Factor_->stride();
00524     }
00525   
00526   int ierr = 0;
00527   if (equilibrate_) ierr = equilibrateMatrix();
00528 
00529   if (ierr!=0) return(ierr);
00530   
00531   INFO_ = 0;
00532   this->POTRF(Matrix_->UPLO(), numRowCols_, AF_, LDAF_, &INFO_);
00533   factored_ = true;
00534 
00535   return(INFO_);
00536 }
00537 
00538 //=============================================================================
00539 
00540 template<typename OrdinalType, typename ScalarType>
00541 int SerialSpdDenseSolver<OrdinalType,ScalarType>::solve() {
00542 
00543   // We will call one of four routines depending on what services the user wants and 
00544   // whether or not the matrix has been inverted or factored already.
00545   //
00546   // If the matrix has been inverted, use DGEMM to compute solution.
00547   // Otherwise, if the user want the matrix to be equilibrated or wants a refined solution, we will
00548   // call the X interface.
00549   // Otherwise, if the matrix is already factored we will call the TRS interface.
00550   // Otherwise, if the matrix is unfactored we will call the SV interface.
00551 
00552   int ierr = 0;
00553   if (equilibrate_) {
00554     ierr = equilibrateRHS();
00555     equilibratedB_ = true;
00556   }
00557   if (ierr != 0) return(ierr);  // Can't equilibrate B, so return.
00558 
00559   TEUCHOS_TEST_FOR_EXCEPTION( (equilibratedA_ && !equilibratedB_) || (!equilibratedA_ && equilibratedB_) , 
00560                      std::logic_error, "SerialSpdDenseSolver<T>::solve: Matrix and vectors must be similarly scaled!");
00561   TEUCHOS_TEST_FOR_EXCEPTION( RHS_==Teuchos::null, std::invalid_argument, 
00562                      "SerialSpdDenseSolver<T>::solve: No right-hand side vector (RHS) has been set for the linear system!");
00563   TEUCHOS_TEST_FOR_EXCEPTION( LHS_==Teuchos::null, std::invalid_argument, 
00564                      "SerialSpdDenseSolver<T>::solve: No solution vector (LHS) has been set for the linear system!");
00565 
00566   if (shouldEquilibrate() && !equilibratedA_)
00567     std::cout << "WARNING!  SerialSpdDenseSolver<T>::solve: System should be equilibrated!" << std::endl;
00568 
00569   if (inverted()) {
00570 
00571     TEUCHOS_TEST_FOR_EXCEPTION( RHS_->values() == LHS_->values(), std::invalid_argument, 
00572                         "SerialSpdDenseSolver<T>::solve: X and B must be different vectors if matrix is inverted.");
00573 
00574     INFO_ = 0;
00575     this->GEMM(Teuchos::NO_TRANS, Teuchos::NO_TRANS, numRowCols_, RHS_->numCols(), 
00576          numRowCols_, 1.0, AF_, LDAF_, RHS_->values(), RHS_->stride(), 0.0, 
00577          LHS_->values(), LHS_->stride());
00578     if (INFO_!=0) return(INFO_);
00579     solved_ = true;
00580   }
00581   else {
00582 
00583     if (!factored()) factor(); // Matrix must be factored
00584     
00585     if (RHS_->values()!=LHS_->values()) {
00586        (*LHS_) = (*RHS_); // Copy B to X if needed
00587     }
00588     INFO_ = 0;
00589     this->POTRS(Matrix_->UPLO(), numRowCols_, RHS_->numCols(), AF_, LDAF_, LHS_->values(), LHS_->stride(), &INFO_);
00590     if (INFO_!=0) return(INFO_);
00591     solved_ = true;
00592 
00593   }
00594   int ierr1=0;
00595   if (refineSolution_ && !inverted()) ierr1 = applyRefinement();
00596   if (ierr1!=0) 
00597     return(ierr1);
00598   
00599   if (equilibrate_) ierr1 = unequilibrateLHS();
00600   return(ierr1);
00601 }
00602 //=============================================================================
00603 
00604 template<typename OrdinalType, typename ScalarType>
00605 int SerialSpdDenseSolver<OrdinalType,ScalarType>::applyRefinement()
00606 {
00607   TEUCHOS_TEST_FOR_EXCEPTION(!solved(), std::logic_error,
00608          "SerialSpdDenseSolver<T>::applyRefinement: Must have an existing solution!");
00609   TEUCHOS_TEST_FOR_EXCEPTION(A_==AF_, std::logic_error,
00610          "SerialSpdDenseSolver<T>::applyRefinement: Cannot apply refinement if no original copy of A!");
00611 
00612   OrdinalType NRHS = RHS_->numCols();
00613   FERR_.resize( NRHS );
00614   BERR_.resize( NRHS );
00615   allocateWORK();
00616   allocateIWORK();
00617   
00618   INFO_ = 0;
00619   this->PORFS(Matrix_->UPLO(), numRowCols_, NRHS, A_, LDA_, AF_, LDAF_,
00620         RHS_->values(), RHS_->stride(), LHS_->values(), LHS_->stride(), 
00621               &FERR_[0], &BERR_[0], &WORK_[0], &IWORK_[0], &INFO_);
00622   
00623   solutionErrorsEstimated_ = true;
00624   reciprocalConditionEstimated_ = true;
00625   solutionRefined_ = true;
00626   
00627   return(INFO_);
00628 
00629 }
00630 
00631 //=============================================================================
00632 
00633 template<typename OrdinalType, typename ScalarType>
00634 int SerialSpdDenseSolver<OrdinalType,ScalarType>::computeEquilibrateScaling() 
00635 {
00636   if (R_.size()!=0) return(0); // Already computed
00637  
00638   R_.resize( numRowCols_ );
00639   
00640   INFO_ = 0;
00641   this->POEQU (numRowCols_, AF_, LDAF_, &R_[0], &SCOND_, &AMAX_, &INFO_);
00642   if (SCOND_<0.1*ScalarTraits<MagnitudeType>::one() || 
00643       AMAX_ < ScalarTraits<ScalarType>::rmin() || 
00644       AMAX_ > ScalarTraits<ScalarType>::rmax()) 
00645     shouldEquilibrate_ = true;
00646 
00647   return(INFO_);
00648 }
00649 
00650 //=============================================================================
00651 
00652 template<typename OrdinalType, typename ScalarType>
00653 int SerialSpdDenseSolver<OrdinalType,ScalarType>::equilibrateMatrix()
00654 {
00655   OrdinalType i, j;
00656   int ierr = 0;
00657 
00658   if (equilibratedA_) return(0); // Already done.
00659   if (R_.size()==0) ierr = computeEquilibrateScaling(); // Compute R if needed.
00660   if (ierr!=0) return(ierr);     // If return value is not zero, then can't equilibrate.
00661   if (Matrix_->upper()) {
00662     if (A_==AF_) {
00663       ScalarType * ptr;
00664       for (j=0; j<numRowCols_; j++) {
00665   ptr = A_ + j*LDA_;
00666   ScalarType s1 = R_[j];
00667   for (i=0; i<=j; i++) {
00668     *ptr = *ptr*s1*R_[i];
00669     ptr++;
00670   }
00671       }
00672     }
00673     else {
00674       ScalarType * ptr;
00675       ScalarType * ptr1;
00676       for (j=0; j<numRowCols_; j++) {
00677   ptr = A_ + j*LDA_;
00678   ptr1 = AF_ + j*LDAF_;
00679   ScalarType s1 = R_[j];
00680   for (i=0; i<=j; i++) {
00681     *ptr = *ptr*s1*R_[i];
00682     ptr++;
00683     *ptr1 = *ptr1*s1*R_[i];
00684     ptr1++;
00685   }
00686       }
00687     }
00688   }
00689   else {
00690     if (A_==AF_) {
00691       ScalarType * ptr;
00692       for (j=0; j<numRowCols_; j++) {
00693   ptr = A_ + j + j*LDA_;
00694   ScalarType s1 = R_[j];
00695   for (i=j; i<numRowCols_; i++) {
00696     *ptr = *ptr*s1*R_[i];
00697     ptr++;
00698   }
00699       }
00700     }
00701     else {
00702       ScalarType * ptr;
00703       ScalarType * ptr1;
00704       for (j=0; j<numRowCols_; j++) {
00705   ptr = A_ + j + j*LDA_;
00706   ptr1 = AF_ + j + j*LDAF_;
00707   ScalarType s1 = R_[j];
00708   for (i=j; i<numRowCols_; i++) {
00709     *ptr = *ptr*s1*R_[i];
00710     ptr++;
00711     *ptr1 = *ptr1*s1*R_[i];
00712     ptr1++;
00713   }
00714       }
00715     }
00716   }
00717   
00718   equilibratedA_ = true;
00719   
00720   return(0);
00721 }
00722 
00723 //=============================================================================
00724 
00725 template<typename OrdinalType, typename ScalarType>
00726 int SerialSpdDenseSolver<OrdinalType,ScalarType>::equilibrateRHS()
00727 {
00728   OrdinalType i, j;
00729   int ierr = 0;
00730 
00731   if (equilibratedB_) return(0); // Already done.
00732   if (R_.size()==0) ierr = computeEquilibrateScaling(); // Compute R if needed.
00733   if (ierr!=0) return(ierr);     // Can't count on R being computed.
00734 
00735   OrdinalType LDB = RHS_->stride(), NRHS = RHS_->numCols();
00736   ScalarType * B = RHS_->values();
00737   ScalarType * ptr;
00738   for (j=0; j<NRHS; j++) {
00739     ptr = B + j*LDB;
00740     for (i=0; i<numRowCols_; i++) {
00741       *ptr = *ptr*R_[i];
00742       ptr++;
00743     }
00744   }
00745   
00746   equilibratedB_ = true;
00747 
00748   return(0);
00749 }
00750 
00751 //=============================================================================
00752 
00753 template<typename OrdinalType, typename ScalarType>
00754 int SerialSpdDenseSolver<OrdinalType,ScalarType>::unequilibrateLHS()
00755 {
00756   OrdinalType i, j;
00757 
00758   if (!equilibratedB_) return(0); // Nothing to do
00759 
00760   OrdinalType LDX = LHS_->stride(), NLHS = LHS_->numCols();
00761   ScalarType * X = LHS_->values();
00762   ScalarType * ptr;
00763   for (j=0; j<NLHS; j++) {
00764     ptr = X + j*LDX;
00765     for (i=0; i<numRowCols_; i++) {
00766       *ptr = *ptr*R_[i];
00767       ptr++;
00768     }
00769   }
00770 
00771   return(0);
00772 }
00773 
00774 //=============================================================================
00775 
00776 template<typename OrdinalType, typename ScalarType>
00777 int SerialSpdDenseSolver<OrdinalType,ScalarType>::invert()
00778 {
00779   if (!factored()) factor(); // Need matrix factored.
00780 
00781   INFO_ = 0;
00782   this->POTRI( Matrix_->UPLO(), numRowCols_, AF_, LDAF_, &INFO_);
00783   
00784   // Copy lower/upper triangle to upper/lower triangle to make full inverse.
00785   if (Matrix_->upper())
00786     Matrix_->setLower();
00787   else
00788     Matrix_->setUpper();
00789 
00790   inverted_ = true;
00791   factored_ = false;
00792   
00793   return(INFO_);
00794 }
00795 
00796 //=============================================================================
00797 
00798 template<typename OrdinalType, typename ScalarType>
00799 int SerialSpdDenseSolver<OrdinalType,ScalarType>::reciprocalConditionEstimate(MagnitudeType & Value)
00800 {
00801   if (reciprocalConditionEstimated()) {
00802     Value = RCOND_;
00803     return(0); // Already computed, just return it.
00804   }
00805 
00806   if ( ANORM_<ScalarTraits<MagnitudeType>::zero() ) ANORM_ = Matrix_->normOne();
00807 
00808   int ierr = 0;
00809   if (!factored()) ierr = factor(); // Need matrix factored.
00810   if (ierr!=0) return(ierr);
00811 
00812   allocateWORK();
00813   allocateIWORK();
00814 
00815   // We will assume a one-norm condition number
00816   INFO_ = 0;
00817   this->POCON( Matrix_->UPLO(), numRowCols_, AF_, LDAF_, ANORM_, &RCOND_, &WORK_[0], &IWORK_[0], &INFO_);
00818   reciprocalConditionEstimated_ = true;
00819   Value = RCOND_;
00820 
00821   return(INFO_);
00822 }
00823 //=============================================================================
00824 
00825 template<typename OrdinalType, typename ScalarType>
00826 void SerialSpdDenseSolver<OrdinalType,ScalarType>::Print(std::ostream& os) const {
00827 
00828   if (Matrix_!=Teuchos::null) os << "Solver Matrix"          << std::endl << *Matrix_ << std::endl;
00829   if (Factor_!=Teuchos::null) os << "Solver Factored Matrix" << std::endl << *Factor_ << std::endl;
00830   if (LHS_   !=Teuchos::null) os << "Solver LHS"             << std::endl << *LHS_    << std::endl;
00831   if (RHS_   !=Teuchos::null) os << "Solver RHS"             << std::endl << *RHS_    << std::endl;
00832 
00833 }
00834 
00835 } // namespace Teuchos
00836 
00837 #endif /* TEUCHOS_SERIALSPDDENSESOLVER_H */
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Defines