Teuchos Package Browser (Single Doxygen Collection) Version of the Day
Teuchos_SerialQRDenseSolver.hpp
Go to the documentation of this file.
00001 // @HEADER
00002 // ***********************************************************************
00003 //
00004 //                    Teuchos: Common Tools Package
00005 //                 Copyright (2004) Sandia Corporation
00006 //
00007 // Under terms of Contract DE-AC04-94AL85000, there is a non-exclusive
00008 // license for use of this work by or on behalf of the U.S. Government.
00009 //
00010 // Redistribution and use in source and binary forms, with or without
00011 // modification, are permitted provided that the following conditions are
00012 // met:
00013 //
00014 // 1. Redistributions of source code must retain the above copyright
00015 // notice, this list of conditions and the following disclaimer.
00016 //
00017 // 2. Redistributions in binary form must reproduce the above copyright
00018 // notice, this list of conditions and the following disclaimer in the
00019 // documentation and/or other materials provided with the distribution.
00020 //
00021 // 3. Neither the name of the Corporation nor the names of the
00022 // contributors may be used to endorse or promote products derived from
00023 // this software without specific prior written permission.
00024 //
00025 // THIS SOFTWARE IS PROVIDED BY SANDIA CORPORATION "AS IS" AND ANY
00026 // EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
00027 // IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR
00028 // PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL SANDIA CORPORATION OR THE
00029 // CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL,
00030 // EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO,
00031 // PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR
00032 // PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF
00033 // LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING
00034 // NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
00035 // SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
00036 //
00037 // Questions? Contact Michael A. Heroux (maherou@sandia.gov)
00038 //
00039 // ***********************************************************************
00040 // @HEADER
00041 
00042 #ifndef _TEUCHOS_SERIALQRDENSESOLVER_HPP_
00043 #define _TEUCHOS_SERIALQRDENSESOLVER_HPP_
00044 
00048 #include "Teuchos_CompObject.hpp"
00049 #include "Teuchos_BLAS.hpp"
00050 #include "Teuchos_LAPACK.hpp"
00051 #include "Teuchos_RCP.hpp"
00052 #include "Teuchos_ConfigDefs.hpp"
00053 #include "Teuchos_SerialDenseMatrix.hpp"
00054 #include "Teuchos_SerialDenseSolver.hpp"
00055 #include "Teuchos_ScalarTraits.hpp"
00056 
00057 #ifdef HAVE_TEUCHOSNUMERICS_EIGEN
00058 #include "Eigen/Dense"
00059 #endif
00060 
00129 namespace Teuchos {
00130 
00131   template<typename OrdinalType, typename ScalarType>
00132   class SerialQRDenseSolver : public CompObject, public Object, public BLAS<OrdinalType, ScalarType>,
00133                               public LAPACK<OrdinalType, ScalarType>
00134   {
00135 
00136   public:
00137 
00138     typedef typename 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     SerialQRDenseSolver();
00159 
00161     virtual ~SerialQRDenseSolver();
00163 
00165 
00166 
00168 
00170     int setMatrix(const RCP<SerialDenseMatrix<OrdinalType, ScalarType> >& A);
00171 
00173 
00176     int setVectors(const RCP<SerialDenseMatrix<OrdinalType, ScalarType> >& X,
00177                    const RCP<SerialDenseMatrix<OrdinalType, ScalarType> >& B);
00179 
00181 
00182 
00184 
00186     void factorWithEquilibration(bool flag) {equilibrate_ = flag; return;}
00187 
00189     void solveWithTranspose(bool flag) {transpose_ = flag; if (flag) TRANS_ = Teuchos::CONJ_TRANS; else TRANS_ = Teuchos::NO_TRANS; return;}
00190 
00192     void solveWithTransposeFlag(Teuchos::ETransp trans) {TRANS_ = trans; if (trans != Teuchos::NO_TRANS) {  transpose_ = true; } }
00193 
00195 
00197 
00198 
00200 
00203     int factor();
00204 
00206 
00209     int solve();
00210 
00212 
00215     int computeEquilibrateScaling();
00216 
00218 
00222     int equilibrateMatrix();
00223 
00225 
00229     int equilibrateRHS();
00230 
00232 
00236     int unequilibrateLHS();
00237 
00239 
00242     int formQ();
00243 
00245 
00248     int formR();
00249 
00251 
00254     int multiplyQ (ETransp transq, SerialDenseMatrix<OrdinalType, ScalarType> &C);
00255 
00257 
00260     int solveR (ETransp transr, SerialDenseMatrix<OrdinalType, ScalarType> &C);
00262 
00264 
00265 
00267     bool transpose() {return(transpose_);}
00268 
00270     bool factored() {return(factored_);}
00271 
00273     bool equilibratedA() {return(equilibratedA_);}
00274 
00276     bool equilibratedB() {return(equilibratedB_);}
00277 
00279     bool shouldEquilibrate() {computeEquilibrateScaling(); return(shouldEquilibrate_);}
00280 
00282     bool solved() {return(solved_);}
00283 
00285     bool formedQ() {return(formedQ_);}
00286 
00288     bool formedR() {return(formedR_);}
00289 
00291 
00293 
00294 
00296     RCP<SerialDenseMatrix<OrdinalType, ScalarType> > getMatrix()  const {return(Matrix_);}
00297 
00299     RCP<SerialDenseMatrix<OrdinalType, ScalarType> > getFactoredMatrix()  const {return(Factor_);}
00300 
00302     RCP<SerialDenseMatrix<OrdinalType, ScalarType> > getQ()  const {return(FactorQ_);}
00303 
00305     RCP<SerialDenseMatrix<OrdinalType, ScalarType> > getR()  const {return(FactorR_);}
00306 
00308     RCP<SerialDenseMatrix<OrdinalType, ScalarType> > getLHS() const {return(LHS_);}
00309 
00311     RCP<SerialDenseMatrix<OrdinalType, ScalarType> > getRHS() const {return(RHS_);}
00312 
00314     OrdinalType numRows()  const {return(M_);}
00315 
00317     OrdinalType numCols()  const {return(N_);}
00318 
00320     std::vector<ScalarType> tau()  const {return(TAU_);}
00321 
00323     MagnitudeType ANORM()  const {return(ANORM_);}
00324 
00326 
00328 
00329 
00330     void Print(std::ostream& os) const;
00332   protected:
00333 
00334     void allocateWORK() { LWORK_ = 4*N_; WORK_.resize( LWORK_ ); return;}
00335     void resetMatrix();
00336     void resetVectors();
00337 
00338 
00339     bool equilibrate_;
00340     bool shouldEquilibrate_;
00341     bool equilibratedA_;
00342     bool equilibratedB_;
00343     OrdinalType equilibrationModeA_;
00344     OrdinalType equilibrationModeB_;
00345     bool transpose_;
00346     bool factored_;
00347     bool solved_;
00348     bool matrixIsZero_;
00349     bool formedQ_;
00350     bool formedR_;
00351 
00352     Teuchos::ETransp TRANS_;
00353 
00354     OrdinalType M_;
00355     OrdinalType N_;
00356     OrdinalType Min_MN_;
00357     OrdinalType LDA_;
00358     OrdinalType LDAF_;
00359     OrdinalType LDQ_;
00360     OrdinalType LDR_;
00361     OrdinalType INFO_;
00362     OrdinalType LWORK_;
00363 
00364     std::vector<ScalarType> TAU_;
00365 
00366     MagnitudeType ANORM_;
00367     MagnitudeType BNORM_;
00368 
00369     RCP<SerialDenseMatrix<OrdinalType, ScalarType> > Matrix_;
00370     RCP<SerialDenseMatrix<OrdinalType, ScalarType> > LHS_;
00371     RCP<SerialDenseMatrix<OrdinalType, ScalarType> > TMP_;
00372     RCP<SerialDenseMatrix<OrdinalType, ScalarType> > RHS_;
00373     RCP<SerialDenseMatrix<OrdinalType, ScalarType> > Factor_;
00374     RCP<SerialDenseMatrix<OrdinalType, ScalarType> > FactorQ_;
00375     RCP<SerialDenseMatrix<OrdinalType, ScalarType> > FactorR_;
00376 
00377     ScalarType * A_;
00378     ScalarType * AF_;
00379     ScalarType * Q_;
00380     ScalarType * R_;
00381     std::vector<ScalarType> WORK_;
00382 #ifdef HAVE_TEUCHOSNUMERICS_EIGEN
00383     Eigen::HouseholderQR<EigenMatrix> qr_;
00384 #endif
00385 
00386   private:
00387     // SerialQRDenseSolver copy constructor (put here because we don't want user access)
00388 
00389     SerialQRDenseSolver(const SerialQRDenseSolver<OrdinalType, ScalarType>& Source);
00390     SerialQRDenseSolver & operator=(const SerialQRDenseSolver<OrdinalType, ScalarType>& Source);
00391 
00392   };
00393 
00394   // Helper traits to distinguish work arrays for real and complex-valued datatypes.
00395   using namespace details;
00396 
00397 //=============================================================================
00398 
00399 template<typename OrdinalType, typename ScalarType>
00400 SerialQRDenseSolver<OrdinalType,ScalarType>::SerialQRDenseSolver()
00401   : CompObject(),
00402     Object("Teuchos::SerialQRDenseSolver"),
00403     equilibrate_(false),
00404     shouldEquilibrate_(false),
00405     equilibratedA_(false),
00406     equilibratedB_(false),
00407     equilibrationModeA_(0),
00408     equilibrationModeB_(0),
00409     transpose_(false),
00410     factored_(false),
00411     solved_(false),
00412     matrixIsZero_(false),
00413     formedQ_(false),
00414     formedR_(false),
00415     TRANS_(Teuchos::NO_TRANS),
00416     M_(0),
00417     N_(0),
00418     Min_MN_(0),
00419     LDA_(0),
00420     LDAF_(0),
00421     LDQ_(0),
00422     LDR_(0),
00423     INFO_(0),
00424     LWORK_(0),
00425     ANORM_(ScalarTraits<MagnitudeType>::zero()),
00426     BNORM_(ScalarTraits<MagnitudeType>::zero()),
00427     A_(0),
00428     AF_(0),
00429     Q_(0),
00430     R_(0)
00431 {
00432   resetMatrix();
00433 }
00434 //=============================================================================
00435 
00436 template<typename OrdinalType, typename ScalarType>
00437 SerialQRDenseSolver<OrdinalType,ScalarType>::~SerialQRDenseSolver()
00438 {}
00439 
00440 //=============================================================================
00441 
00442 template<typename OrdinalType, typename ScalarType>
00443 void SerialQRDenseSolver<OrdinalType,ScalarType>::resetVectors()
00444 {
00445   LHS_ = Teuchos::null;
00446   TMP_ = Teuchos::null;
00447   RHS_ = Teuchos::null;
00448   solved_ = false;
00449   equilibratedB_ = false;
00450 }
00451 //=============================================================================
00452 
00453 template<typename OrdinalType, typename ScalarType>
00454 void SerialQRDenseSolver<OrdinalType,ScalarType>::resetMatrix()
00455 {
00456   resetVectors();
00457   equilibratedA_ = false;
00458   equilibrationModeA_ = 0;
00459   equilibrationModeB_ = 0;
00460   factored_ = false;
00461   matrixIsZero_ = false;
00462   formedQ_ = false;
00463   formedR_ = false;
00464   M_ = 0;
00465   N_ = 0;
00466   Min_MN_ = 0;
00467   LDA_ = 0;
00468   LDAF_ = 0;
00469   LDQ_ = 0;
00470   LDR_ = 0;
00471   ANORM_ = -ScalarTraits<MagnitudeType>::one();
00472   BNORM_ = -ScalarTraits<MagnitudeType>::one();
00473   A_ = 0;
00474   AF_ = 0;
00475   Q_ = 0;
00476   R_ = 0;
00477   INFO_ = 0;
00478   LWORK_ = 0;
00479 }
00480 //=============================================================================
00481 
00482 template<typename OrdinalType, typename ScalarType>
00483 int SerialQRDenseSolver<OrdinalType,ScalarType>::setMatrix(const RCP<SerialDenseMatrix<OrdinalType,ScalarType> >& A)
00484 {
00485   TEUCHOS_TEST_FOR_EXCEPTION(A->numRows() < A->numCols(), std::invalid_argument,
00486                      "SerialQRDenseSolver<T>::setMatrix: the matrix A must have A.numRows() >= A.numCols()!");
00487 
00488   resetMatrix();
00489   Matrix_ = A;
00490   Factor_ = A;
00491   FactorQ_ = A;
00492   FactorR_ = A;
00493   M_ = A->numRows();
00494   N_ = A->numCols();
00495   Min_MN_ = TEUCHOS_MIN(M_,N_);
00496   LDA_ = A->stride();
00497   LDAF_ = LDA_;
00498   LDQ_ = LDA_;
00499   LDR_ = N_;
00500   A_ = A->values();
00501   AF_ = A->values();
00502   Q_ = A->values();
00503   R_ = A->values();
00504 
00505   return(0);
00506 }
00507 //=============================================================================
00508 
00509 template<typename OrdinalType, typename ScalarType>
00510 int SerialQRDenseSolver<OrdinalType,ScalarType>::setVectors(const RCP<SerialDenseMatrix<OrdinalType,ScalarType> >& X,
00511                                                            const RCP<SerialDenseMatrix<OrdinalType,ScalarType> >& B)
00512 {
00513   // Check that these new vectors are consistent
00514   TEUCHOS_TEST_FOR_EXCEPTION(B->numCols() != X->numCols(), std::invalid_argument,
00515                      "SerialQRDenseSolver<T>::setVectors: X and B have different numbers of columns!");
00516   TEUCHOS_TEST_FOR_EXCEPTION(B->values()==0, std::invalid_argument,
00517                      "SerialQRDenseSolver<T>::setVectors: B is an empty SerialDenseMatrix<T>!");
00518   TEUCHOS_TEST_FOR_EXCEPTION(X->values()==0, std::invalid_argument,
00519                      "SerialQRDenseSolver<T>::setVectors: X is an empty SerialDenseMatrix<T>!");
00520   TEUCHOS_TEST_FOR_EXCEPTION(B->stride()<1, std::invalid_argument,
00521                      "SerialQRDenseSolver<T>::setVectors: B has an invalid stride!");
00522   TEUCHOS_TEST_FOR_EXCEPTION(X->stride()<1, std::invalid_argument,
00523                      "SerialQRDenseSolver<T>::setVectors: X has an invalid stride!");
00524 
00525   resetVectors();
00526   LHS_ = X;
00527   RHS_ = B;
00528 
00529   return(0);
00530 }
00531 //=============================================================================
00532 
00533 template<typename OrdinalType, typename ScalarType>
00534 int SerialQRDenseSolver<OrdinalType,ScalarType>::factor() {
00535 
00536   if (factored()) return(0);
00537 
00538   // Equilibrate matrix if necessary
00539   int ierr = 0;
00540   if (equilibrate_) ierr = equilibrateMatrix();
00541   if (ierr!=0) return(ierr);
00542 
00543   allocateWORK();
00544   if ((int)TAU_.size()!=Min_MN_) TAU_.resize( Min_MN_ );
00545 
00546   INFO_ = 0;
00547 
00548   // Factor
00549 #ifdef HAVE_TEUCHOSNUMERICS_EIGEN
00550   EigenMatrixMap aMap( AF_, M_, N_, EigenOuterStride(LDAF_) );
00551   EigenScalarArray tau;
00552   EigenScalarArrayMap tauMap(&TAU_[0],TEUCHOS_MIN(M_,N_));
00553   qr_.compute(aMap);
00554   tau = qr_.hCoeffs();
00555   for (OrdinalType i=0; i<tauMap.innerSize(); i++) {
00556     tauMap(i) = tau(i);
00557   }
00558   EigenMatrix qrMat = qr_.matrixQR();
00559   for (OrdinalType j=0; j<aMap.outerSize(); j++) {
00560     for (OrdinalType i=0; i<aMap.innerSize(); i++) {
00561       aMap(i,j) = qrMat(i,j);
00562     }
00563   }
00564 #else
00565   this->GEQRF(M_, N_, AF_, LDAF_, &TAU_[0], &WORK_[0], LWORK_, &INFO_);
00566 #endif
00567 
00568   factored_ = true;
00569 
00570   return(INFO_);
00571 }
00572 
00573 //=============================================================================
00574 
00575 template<typename OrdinalType, typename ScalarType>
00576 int SerialQRDenseSolver<OrdinalType,ScalarType>::solve() {
00577 
00578   // Check if the matrix is zero
00579   if (matrixIsZero_) {
00580     LHS_->putScalar(ScalarTraits<ScalarType>::zero());
00581     return(0);
00582   }
00583 
00584   // Equilibrate RHS if necessary
00585   int ierr = 0;
00586   if (equilibrate_) {
00587     ierr = equilibrateRHS();
00588     equilibratedB_ = true;
00589   }
00590   if (ierr != 0) return(ierr);
00591 
00592   TEUCHOS_TEST_FOR_EXCEPTION( (equilibratedA_ && !equilibratedB_) || (!equilibratedA_ && equilibratedB_) ,
00593                      std::logic_error, "SerialQRDenseSolver<T>::solve: Matrix and vectors must be similarly scaled!");
00594   TEUCHOS_TEST_FOR_EXCEPTION( RHS_==Teuchos::null, std::invalid_argument,
00595                      "SerialQRDenseSolver<T>::solve: No right-hand side vector (RHS) has been set for the linear system!");
00596   TEUCHOS_TEST_FOR_EXCEPTION( LHS_==Teuchos::null, std::invalid_argument,
00597                      "SerialQRDenseSolver<T>::solve: No solution vector (LHS) has been set for the linear system!");
00598   if ( TRANS_ == Teuchos::NO_TRANS ) {
00599     TEUCHOS_TEST_FOR_EXCEPTION( LHS_->numRows() != N_, std::invalid_argument,
00600                      "SerialQRDenseSolver<T>::solve: No transpose: must have LHS_->numRows() = N_");
00601   } else {
00602     TEUCHOS_TEST_FOR_EXCEPTION( LHS_->numRows() != M_, std::invalid_argument,
00603                      "SerialQRDenseSolver<T>::solve: Transpose: must have LHS_->numRows() = M_");
00604   }
00605 
00606   if (shouldEquilibrate() && !equilibratedA_)
00607     std::cout << "WARNING!  SerialQRDenseSolver<T>::solve: System should be equilibrated!" << std::endl;
00608 
00609   // Matrix must be factored
00610   if (!factored()) factor();
00611 
00612   TMP_ = rcp( new SerialDenseMatrix<OrdinalType,ScalarType>(M_, RHS_->numCols()) );
00613   for (OrdinalType j=0; j<RHS_->numCols(); j++) {
00614     for (OrdinalType i=0; i<RHS_->numRows(); i++) {
00615       (*TMP_)(i,j) = (*RHS_)(i,j);
00616     }
00617   }
00618 
00619   INFO_ = 0;
00620 
00621   // Solve
00622   if ( TRANS_ == Teuchos::NO_TRANS ) {
00623 
00624     // Solve A lhs = rhs
00625     this->multiplyQ( Teuchos::CONJ_TRANS, *TMP_ );
00626     this->solveR( Teuchos::NO_TRANS, *TMP_ );
00627 
00628   } else {
00629 
00630     // Solve A**H lhs = rhs
00631     this->solveR( Teuchos::CONJ_TRANS, *TMP_ );
00632     for (OrdinalType j = 0; j < RHS_->numCols(); j++ ) {
00633       for (OrdinalType i = N_; i < M_; i++ ) {
00634         (*TMP_)(i, j) = ScalarTraits<ScalarType>::zero();
00635       }
00636     }
00637     this->multiplyQ( Teuchos::NO_TRANS, *TMP_ );
00638 
00639   }
00640   for (OrdinalType j = 0; j < LHS_->numCols(); j++ ) {
00641     for (OrdinalType i = 0; i < LHS_->numRows(); i++ ) {
00642       (*LHS_)(i, j) = (*TMP_)(i,j);
00643     }
00644   }
00645 
00646   solved_ = true;
00647 
00648   // Unequilibrate LHS if necessary
00649   if (equilibrate_) {
00650     ierr = unequilibrateLHS();
00651   }
00652   if (ierr != 0) {
00653     return ierr;
00654   }
00655 
00656   return INFO_;
00657 }
00658 
00659 //=============================================================================
00660 
00661 template<typename OrdinalType, typename ScalarType>
00662 int SerialQRDenseSolver<OrdinalType,ScalarType>::computeEquilibrateScaling()
00663 {
00664   MagnitudeType safeMin = ScalarTraits<ScalarType>::sfmin();
00665   MagnitudeType prec = ScalarTraits<ScalarType>::prec();
00666   ScalarType sZero = ScalarTraits<ScalarType>::zero();
00667   ScalarType sOne  = ScalarTraits<ScalarType>::one();
00668   MagnitudeType mZero = ScalarTraits<ScalarType>::magnitude(sZero);
00669 
00670   MagnitudeType smlnum = ScalarTraits<ScalarType>::magnitude(safeMin/prec);
00671   MagnitudeType bignum = ScalarTraits<ScalarType>::magnitude(sOne/smlnum);
00672 
00673   // Compute maximum absolute value of matrix entries
00674   OrdinalType i, j;
00675   MagnitudeType anorm = ScalarTraits<ScalarType>::magnitude(ScalarTraits<ScalarType>::zero());
00676   for (j = 0; j < N_; j++) {
00677     for (i = 0; i < M_; i++) {
00678       anorm = TEUCHOS_MAX( anorm, ScalarTraits<ScalarType>::magnitude((*Matrix_)(i,j)) );
00679     }
00680   }
00681 
00682   ANORM_ = anorm;
00683 
00684   if (ANORM_ > mZero && ANORM_ < smlnum) {
00685     // Scale matrix norm up to smlnum
00686     shouldEquilibrate_ = true;
00687   } else if (ANORM_ > bignum) {
00688     // Scale matrix norm down to bignum
00689     shouldEquilibrate_ = true;
00690   } else if (ANORM_ == mZero) {
00691     // Matrix all zero. Return zero solution
00692     matrixIsZero_ = true;
00693   }
00694 
00695   return(0);
00696 }
00697 
00698 //=============================================================================
00699 
00700 template<typename OrdinalType, typename ScalarType>
00701 int SerialQRDenseSolver<OrdinalType,ScalarType>::equilibrateMatrix()
00702 {
00703   if (equilibratedA_) return(0);
00704 
00705   MagnitudeType safeMin = ScalarTraits<ScalarType>::sfmin();
00706   MagnitudeType prec = ScalarTraits<ScalarType>::prec();
00707   ScalarType sZero = ScalarTraits<ScalarType>::zero();
00708   ScalarType sOne  = ScalarTraits<ScalarType>::one();
00709   MagnitudeType mZero = ScalarTraits<ScalarType>::magnitude(sZero);
00710 
00711   MagnitudeType smlnum = ScalarTraits<ScalarType>::magnitude(safeMin/prec);
00712   MagnitudeType bignum = ScalarTraits<ScalarType>::magnitude(sOne/smlnum);
00713   OrdinalType BW = 0;
00714 
00715   // Compute maximum absolute value of matrix entries
00716   OrdinalType i, j;
00717   MagnitudeType anorm = ScalarTraits<ScalarType>::magnitude(ScalarTraits<ScalarType>::zero());
00718   for (j = 0; j < N_; j++) {
00719     for (i = 0; i < M_; i++) {
00720       anorm = TEUCHOS_MAX( anorm, ScalarTraits<ScalarType>::magnitude((*Matrix_)(i,j)) );
00721     }
00722   }
00723 
00724   ANORM_ = anorm;
00725   int ierr1 = 0;
00726   if (ANORM_ > mZero && ANORM_ < smlnum) {
00727     // Scale matrix norm up to smlnum
00728     this->LASCL(Teuchos::ETypeChar[Teuchos::FULL], BW, BW, ANORM_, smlnum, M_, N_, A_, LDA_, &INFO_);
00729     equilibrationModeA_ = 1;
00730   } else if (ANORM_ > bignum) {
00731     // Scale matrix norm down to bignum
00732     this->LASCL(Teuchos::ETypeChar[Teuchos::FULL], BW, BW, ANORM_, bignum, M_, N_, A_, LDA_, &INFO_);
00733     equilibrationModeA_ = 2;
00734   } else if (ANORM_ == mZero) {
00735     // Matrix all zero. Return zero solution
00736     matrixIsZero_ = true;
00737   }
00738 
00739   equilibratedA_ = true;
00740 
00741   return(ierr1);
00742 }
00743 
00744 //=============================================================================
00745 
00746 template<typename OrdinalType, typename ScalarType>
00747 int SerialQRDenseSolver<OrdinalType,ScalarType>::equilibrateRHS()
00748 {
00749   if (equilibratedB_) return(0);
00750 
00751   MagnitudeType safeMin = ScalarTraits<ScalarType>::sfmin();
00752   MagnitudeType prec = ScalarTraits<ScalarType>::prec();
00753   ScalarType sZero = ScalarTraits<ScalarType>::zero();
00754   ScalarType sOne  = ScalarTraits<ScalarType>::one();
00755   MagnitudeType mZero = ScalarTraits<ScalarType>::magnitude(sZero);
00756 
00757   MagnitudeType smlnum = ScalarTraits<ScalarType>::magnitude(safeMin/prec);
00758   MagnitudeType bignum = ScalarTraits<ScalarType>::magnitude(sOne/smlnum);
00759   OrdinalType BW = 0;
00760 
00761   // Compute maximum absolute value of rhs entries
00762   OrdinalType i, j;
00763   MagnitudeType bnorm = ScalarTraits<ScalarType>::magnitude(ScalarTraits<ScalarType>::zero());
00764   for (j = 0; j <RHS_->numCols(); j++) {
00765     for (i = 0; i < RHS_->numRows(); i++) {
00766       bnorm = TEUCHOS_MAX( bnorm, ScalarTraits<ScalarType>::magnitude((*RHS_)(i,j)) );
00767     }
00768   }
00769 
00770   BNORM_ = bnorm;
00771 
00772   int ierr1 = 0;
00773   if (BNORM_ > mZero && BNORM_ < smlnum) {
00774     // Scale RHS norm up to smlnum
00775     this->LASCL(Teuchos::ETypeChar[Teuchos::FULL], BW, BW, BNORM_, smlnum, RHS_->numRows(), RHS_->numCols(),
00776                 RHS_->values(), RHS_->stride(), &INFO_);
00777     equilibrationModeB_ = 1;
00778   } else if (BNORM_ > bignum) {
00779     // Scale RHS norm down to bignum
00780     this->LASCL(Teuchos::ETypeChar[Teuchos::FULL], BW, BW, BNORM_, bignum, RHS_->numRows(), RHS_->numCols(),
00781                 RHS_->values(), RHS_->stride(), &INFO_);
00782     equilibrationModeB_ = 2;
00783   }
00784 
00785   equilibratedB_ = true;
00786 
00787   return(ierr1);
00788 }
00789 
00790 //=============================================================================
00791 
00792 template<typename OrdinalType, typename ScalarType>
00793 int SerialQRDenseSolver<OrdinalType,ScalarType>::unequilibrateLHS()
00794 {
00795   if (!equilibratedB_) return(0);
00796 
00797   MagnitudeType safeMin = ScalarTraits<ScalarType>::sfmin();
00798   MagnitudeType prec = ScalarTraits<ScalarType>::prec();
00799   ScalarType sZero = ScalarTraits<ScalarType>::zero();
00800   ScalarType sOne  = ScalarTraits<ScalarType>::one();
00801   MagnitudeType mZero = ScalarTraits<ScalarType>::magnitude(sZero);
00802   (void) mZero; // Silence "unused variable" compiler warning.
00803 
00804   MagnitudeType smlnum = ScalarTraits<ScalarType>::magnitude(safeMin/prec);
00805   MagnitudeType bignum = ScalarTraits<ScalarType>::magnitude(sOne/smlnum);
00806   OrdinalType BW = 0;
00807 
00808   int ierr1 = 0;
00809   if (equilibrationModeA_ == 1) {
00810     this->LASCL(Teuchos::ETypeChar[Teuchos::FULL], BW, BW, ANORM_, smlnum, LHS_->numRows(), LHS_->numCols(),
00811                 LHS_->values(), LHS_->stride(), &INFO_);
00812   } else if (equilibrationModeA_ == 2) {
00813     this->LASCL(Teuchos::ETypeChar[Teuchos::FULL], BW, BW, ANORM_, bignum, LHS_->numRows(), LHS_->numCols(),
00814                 LHS_->values(), LHS_->stride(), &INFO_);
00815   }
00816   if (equilibrationModeB_ == 1) {
00817     this->LASCL(Teuchos::ETypeChar[Teuchos::FULL], BW, BW, smlnum, BNORM_, LHS_->numRows(), LHS_->numCols(),
00818                 LHS_->values(), LHS_->stride(), &INFO_);
00819   } else if (equilibrationModeB_ == 2) {
00820     this->LASCL(Teuchos::ETypeChar[Teuchos::FULL], BW, BW, bignum, BNORM_, LHS_->numRows(), LHS_->numCols(),
00821                 LHS_->values(), LHS_->stride(), &INFO_);
00822   }
00823 
00824   return(ierr1);
00825 }
00826 
00827 //=============================================================================
00828 
00829 template<typename OrdinalType, typename ScalarType>
00830 int SerialQRDenseSolver<OrdinalType,ScalarType>::formQ() {
00831 
00832   // Matrix must be factored first
00833   if (!factored()) factor();
00834 
00835   // Store Q separately from factored matrix
00836   if (AF_ == Q_) {
00837     FactorQ_ = rcp( new SerialDenseMatrix<OrdinalType,ScalarType>(*Factor_) );
00838     Q_ = FactorQ_->values();
00839     LDQ_ = FactorQ_->stride();
00840   }
00841 
00842   INFO_ = 0;
00843 
00844   // Form Q
00845 #ifdef HAVE_TEUCHOSNUMERICS_EIGEN
00846   EigenMatrixMap qMap( Q_, M_, N_, EigenOuterStride(LDQ_) );
00847   EigenMatrix qMat = qr_.householderQ();
00848   for (OrdinalType j=0; j<qMap.outerSize(); j++) {
00849     for (OrdinalType i=0; i<qMap.innerSize(); i++) {
00850       qMap(i,j) = qMat(i,j);
00851     }
00852   }
00853 #else
00854   this->UNGQR(M_, N_, N_, Q_, LDQ_, &TAU_[0], &WORK_[0], LWORK_, &INFO_);
00855 #endif
00856 
00857   if (INFO_!=0) return(INFO_);
00858 
00859   formedQ_ = true;
00860 
00861   return(INFO_);
00862 }
00863 
00864 //=============================================================================
00865 
00866 template<typename OrdinalType, typename ScalarType>
00867 int SerialQRDenseSolver<OrdinalType,ScalarType>::formR() {
00868 
00869   // Matrix must be factored first
00870   if (!factored()) factor();
00871 
00872   // Store R separately from factored matrix
00873   if (AF_ == R_) {
00874     FactorR_ = rcp( new SerialDenseMatrix<OrdinalType,ScalarType>(N_, N_) );
00875     R_ = FactorR_->values();
00876     LDR_ = FactorR_->stride();
00877   }
00878 
00879   // Form R
00880   for (OrdinalType j=0; j<N_; j++) {
00881     for (OrdinalType i=0; i<=j; i++) {
00882       (*FactorR_)(i,j) = (*Factor_)(i,j);
00883     }
00884   }
00885 
00886   formedR_ = true;
00887 
00888   return(0);
00889 }
00890 
00891 //=============================================================================
00892 
00893 template<typename OrdinalType, typename ScalarType>
00894 int  SerialQRDenseSolver<OrdinalType, ScalarType>::multiplyQ(ETransp transq, SerialDenseMatrix<OrdinalType, ScalarType> &C)
00895 {
00896 
00897   // Check that the matrices are consistent.
00898   TEUCHOS_TEST_FOR_EXCEPTION(C.numRows()!=M_, std::invalid_argument,
00899                      "SerialQRDenseSolver<T>::multiplyQ: C.numRows() != M_!");
00900   TEUCHOS_TEST_FOR_EXCEPTION(C.values()==0, std::invalid_argument,
00901                      "SerialQRDenseSolver<T>::multiplyQ: C is an empty SerialDenseMatrix<T>!");
00902 
00903   // Matrix must be factored
00904   if (!factored()) factor();
00905 
00906   INFO_ = 0;
00907 
00908   // Multiply
00909 #ifdef HAVE_TEUCHOSNUMERICS_EIGEN
00910   EigenMatrixMap cMap( C.values(), C.numRows(), C.numCols(), EigenOuterStride(C.stride()) );
00911   if ( transq == Teuchos::NO_TRANS ) {
00912     // C = Q * C
00913     cMap = qr_.householderQ() * cMap;
00914   } else {
00915     // C = Q**H * C
00916     cMap = qr_.householderQ().adjoint() * cMap;
00917     for (OrdinalType j = 0; j < C.numCols(); j++) {
00918       for (OrdinalType i = N_; i < C.numRows(); i++ ) {
00919         cMap(i, j) = ScalarTraits<ScalarType>::zero();
00920       }
00921     }
00922   }
00923 #else
00924   Teuchos::ETransp NO_TRANS = Teuchos::NO_TRANS;
00925   Teuchos::ETransp TRANS = (Teuchos::ScalarTraits<ScalarType>::isComplex) ? Teuchos::CONJ_TRANS : Teuchos::TRANS;
00926   Teuchos::ESide SIDE = Teuchos::LEFT_SIDE;
00927 
00928   if ( transq == Teuchos::NO_TRANS ) {
00929 
00930     // C = Q * C
00931     this->UNMQR(ESideChar[SIDE], ETranspChar[NO_TRANS], M_, C.numCols(), N_, AF_, LDAF_,
00932                 &TAU_[0], C.values(), C.stride(), &WORK_[0], LWORK_, &INFO_);
00933 
00934   } else {
00935 
00936     // C = Q**H * C
00937     this->UNMQR(ESideChar[SIDE], ETranspChar[TRANS], M_, C.numCols(), N_, AF_, LDAF_,
00938                 &TAU_[0], C.values(), C.stride(), &WORK_[0], LWORK_, &INFO_);
00939 
00940     for (OrdinalType j = 0; j < C.numCols(); j++) {
00941       for (OrdinalType i = N_; i < C.numRows(); i++ ) {
00942         C(i, j) = ScalarTraits<ScalarType>::zero();
00943       }
00944     }
00945   }
00946 #endif
00947 
00948   return(INFO_);
00949 
00950 }
00951 
00952 //=============================================================================
00953 
00954 template<typename OrdinalType, typename ScalarType>
00955 int  SerialQRDenseSolver<OrdinalType, ScalarType>::solveR(ETransp transr, SerialDenseMatrix<OrdinalType, ScalarType> &C)
00956 {
00957 
00958   // Check that the matrices are consistent.
00959   TEUCHOS_TEST_FOR_EXCEPTION(C.numRows()<N_ || C.numRows()>M_, std::invalid_argument,
00960                      "SerialQRDenseSolver<T>::solveR: must have N_ < C.numRows() < M_!");
00961   TEUCHOS_TEST_FOR_EXCEPTION(C.values()==0, std::invalid_argument,
00962                      "SerialQRDenseSolver<T>::solveR: C is an empty SerialDenseMatrix<T>!");
00963 
00964   // Matrix must be factored
00965   if (!factored()) factor();
00966 
00967   INFO_ = 0;
00968 
00969   // Solve
00970 #ifdef HAVE_TEUCHOSNUMERICS_EIGEN
00971   EigenMatrixMap cMap( C.values(), N_, C.numCols(), EigenOuterStride(C.stride()) );
00972   // Check for singularity first like TRTRS
00973   for (OrdinalType j=0; j<N_; j++) {
00974     if ((qr_.matrixQR())(j,j) == ScalarTraits<ScalarType>::zero()) {
00975       INFO_ = j+1;
00976       return(INFO_);
00977     }
00978   }
00979   if ( transr == Teuchos::NO_TRANS ) {
00980     // C = R \ C
00981     qr_.matrixQR().topRows(N_).template triangularView<Eigen::Upper>().solveInPlace(cMap);
00982   } else {
00983     // C = R**H \ C
00984     qr_.matrixQR().topRows(N_).template triangularView<Eigen::Upper>().adjoint().solveInPlace(cMap);
00985   }
00986 #else
00987   Teuchos::ETransp NO_TRANS = Teuchos::NO_TRANS;
00988   Teuchos::ETransp TRANS = (Teuchos::ScalarTraits<ScalarType>::isComplex) ? Teuchos::CONJ_TRANS : Teuchos::TRANS;
00989   Teuchos::EUplo UPLO = Teuchos::UPPER_TRI;
00990   Teuchos::EDiag DIAG = Teuchos::NON_UNIT_DIAG;
00991 
00992   ScalarType* RMAT = (formedR_) ? R_ : AF_;
00993   OrdinalType LDRMAT = (formedR_) ? LDR_ : LDAF_;
00994 
00995   if ( transr == Teuchos::NO_TRANS ) {
00996 
00997     // C = R \ C
00998     this->TRTRS(EUploChar[UPLO], ETranspChar[NO_TRANS], EDiagChar[DIAG], N_, C.numCols(),
00999                 RMAT, LDRMAT, C.values(), C.stride(), &INFO_);
01000 
01001   } else {
01002 
01003     // C = R**H \ C
01004     this->TRTRS(EUploChar[UPLO], ETranspChar[TRANS], EDiagChar[DIAG], N_, C.numCols(),
01005                 RMAT, LDRMAT, C.values(), C.stride(), &INFO_);
01006 
01007   }
01008 #endif
01009 
01010   return(INFO_);
01011 
01012 }
01013 
01014 //=============================================================================
01015 
01016 template<typename OrdinalType, typename ScalarType>
01017 void SerialQRDenseSolver<OrdinalType,ScalarType>::Print(std::ostream& os) const {
01018 
01019   if (Matrix_!=Teuchos::null) os << "Solver Matrix"          << std::endl << *Matrix_ << std::endl;
01020   if (Factor_!=Teuchos::null) os << "Solver Factored Matrix" << std::endl << *Factor_ << std::endl;
01021   if (Q_!=Teuchos::null) os << "Solver Factor Q" << std::endl << *Q_ << std::endl;
01022   if (LHS_   !=Teuchos::null) os << "Solver LHS"             << std::endl << *LHS_    << std::endl;
01023   if (RHS_   !=Teuchos::null) os << "Solver RHS"             << std::endl << *RHS_    << std::endl;
01024 
01025 }
01026 
01027 } // namespace Teuchos
01028 
01029 #endif /* _TEUCHOS_SERIALQRDENSESOLVER_HPP_ */
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Defines