Teuchos - Trilinos Tools Package 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_ScalarTraits.hpp"
00055 
00056 #ifdef HAVE_TEUCHOSNUMERICS_EIGEN
00057 #include "Eigen/Dense"
00058 #endif
00059 
00128 namespace Teuchos {
00129 
00130   template<typename OrdinalType, typename ScalarType>
00131   class SerialQRDenseSolver : public CompObject, public Object, public BLAS<OrdinalType, ScalarType>,
00132             public LAPACK<OrdinalType, ScalarType>
00133   {
00134 
00135   public:
00136 
00137     typedef typename ScalarTraits<ScalarType>::magnitudeType MagnitudeType;
00138 #ifdef HAVE_TEUCHOSNUMERICS_EIGEN
00139     typedef typename Eigen::Matrix<ScalarType,Eigen::Dynamic,Eigen::Dynamic,Eigen::ColMajor> EigenMatrix;
00140     typedef typename Eigen::Matrix<ScalarType,Eigen::Dynamic,1> EigenVector;
00141     typedef typename Eigen::InnerStride<Eigen::Dynamic> EigenInnerStride;
00142     typedef typename Eigen::OuterStride<Eigen::Dynamic> EigenOuterStride;
00143     typedef typename Eigen::Map<EigenVector,0,EigenInnerStride > EigenVectorMap;
00144     typedef typename Eigen::Map<const EigenVector,0,EigenInnerStride > EigenConstVectorMap;
00145     typedef typename Eigen::Map<EigenMatrix,0,EigenOuterStride > EigenMatrixMap;
00146     typedef typename Eigen::Map<const EigenMatrix,0,EigenOuterStride > EigenConstMatrixMap;
00147     typedef typename Eigen::PermutationMatrix<Eigen::Dynamic,Eigen::Dynamic,OrdinalType> EigenPermutationMatrix;
00148     typedef typename Eigen::Array<OrdinalType,Eigen::Dynamic,1> EigenOrdinalArray;
00149     typedef typename Eigen::Map<EigenOrdinalArray> EigenOrdinalArrayMap;
00150     typedef typename Eigen::Array<ScalarType,Eigen::Dynamic,1> EigenScalarArray;
00151     typedef typename Eigen::Map<EigenScalarArray> EigenScalarArrayMap;
00152 #endif
00153 
00155 
00156 
00157     SerialQRDenseSolver();
00158 
00160     virtual ~SerialQRDenseSolver();
00162 
00164 
00165 
00167 
00169     int setMatrix(const RCP<SerialDenseMatrix<OrdinalType, ScalarType> >& A);
00170 
00172 
00175     int setVectors(const RCP<SerialDenseMatrix<OrdinalType, ScalarType> >& X,
00176        const RCP<SerialDenseMatrix<OrdinalType, ScalarType> >& B);
00178 
00180 
00181 
00183 
00185     void factorWithEquilibration(bool flag) {equilibrate_ = flag; return;}
00186 
00188     void solveWithTranspose(bool flag) {transpose_ = flag; if (flag) TRANS_ = Teuchos::CONJ_TRANS; else TRANS_ = Teuchos::NO_TRANS; return;}
00189 
00191     void solveWithTransposeFlag(Teuchos::ETransp trans) {TRANS_ = trans; if (trans != Teuchos::NO_TRANS) {  transpose_ = true; } }
00192 
00194 
00196 
00197 
00199 
00202     int factor();
00203 
00205 
00208     int solve();
00209 
00211 
00214     int computeEquilibrateScaling();
00215 
00217 
00221     int equilibrateMatrix();
00222 
00224 
00228     int equilibrateRHS();
00229 
00231 
00235     int unequilibrateLHS();
00236 
00238 
00241     int formQ();
00242 
00244 
00247     int formR();
00248 
00250 
00253     int multiplyQ (ETransp transq, SerialDenseMatrix<OrdinalType, ScalarType> &C);
00254 
00256 
00259     int solveR (ETransp transr, SerialDenseMatrix<OrdinalType, ScalarType> &C);
00261 
00263 
00264 
00266     bool transpose() {return(transpose_);}
00267 
00269     bool factored() {return(factored_);}
00270 
00272     bool equilibratedA() {return(equilibratedA_);}
00273 
00275     bool equilibratedB() {return(equilibratedB_);}
00276 
00278     bool shouldEquilibrate() {computeEquilibrateScaling(); return(shouldEquilibrate_);}
00279 
00281     bool solved() {return(solved_);}
00282 
00284     bool formedQ() {return(formedQ_);}
00285 
00287     bool formedR() {return(formedR_);}
00288 
00290 
00292 
00293 
00295     RCP<SerialDenseMatrix<OrdinalType, ScalarType> > getMatrix()  const {return(Matrix_);}
00296 
00298     RCP<SerialDenseMatrix<OrdinalType, ScalarType> > getFactoredMatrix()  const {return(Factor_);}
00299 
00301     RCP<SerialDenseMatrix<OrdinalType, ScalarType> > getQ()  const {return(FactorQ_);}
00302 
00304     RCP<SerialDenseMatrix<OrdinalType, ScalarType> > getR()  const {return(FactorR_);}
00305 
00307     RCP<SerialDenseMatrix<OrdinalType, ScalarType> > getLHS() const {return(LHS_);}
00308 
00310     RCP<SerialDenseMatrix<OrdinalType, ScalarType> > getRHS() const {return(RHS_);}
00311 
00313     OrdinalType numRows()  const {return(M_);}
00314 
00316     OrdinalType numCols()  const {return(N_);}
00317 
00319     std::vector<ScalarType> tau()  const {return(TAU_);}
00320 
00322     MagnitudeType ANORM()  const {return(ANORM_);}
00323 
00325 
00327 
00328 
00329     void Print(std::ostream& os) const;
00331   protected:
00332 
00333     void allocateWORK() { LWORK_ = 4*N_; WORK_.resize( LWORK_ ); return;}
00334     void resetMatrix();
00335     void resetVectors();
00336 
00337 
00338     bool equilibrate_;
00339     bool shouldEquilibrate_;
00340     bool equilibratedA_;
00341     bool equilibratedB_;
00342     OrdinalType equilibrationModeA_;
00343     OrdinalType equilibrationModeB_;
00344     bool transpose_;
00345     bool factored_;
00346     bool solved_;
00347     bool matrixIsZero_;
00348     bool formedQ_;
00349     bool formedR_;
00350 
00351     Teuchos::ETransp TRANS_;
00352 
00353     OrdinalType M_;
00354     OrdinalType N_;
00355     OrdinalType Min_MN_;
00356     OrdinalType LDA_;
00357     OrdinalType LDAF_;
00358     OrdinalType LDQ_;
00359     OrdinalType LDR_;
00360     OrdinalType INFO_;
00361     OrdinalType LWORK_;
00362 
00363     std::vector<ScalarType> TAU_;
00364 
00365     MagnitudeType ANORM_;
00366     MagnitudeType BNORM_;
00367 
00368     RCP<SerialDenseMatrix<OrdinalType, ScalarType> > Matrix_;
00369     RCP<SerialDenseMatrix<OrdinalType, ScalarType> > LHS_;
00370     RCP<SerialDenseMatrix<OrdinalType, ScalarType> > TMP_;
00371     RCP<SerialDenseMatrix<OrdinalType, ScalarType> > RHS_;
00372     RCP<SerialDenseMatrix<OrdinalType, ScalarType> > Factor_;
00373     RCP<SerialDenseMatrix<OrdinalType, ScalarType> > FactorQ_;
00374     RCP<SerialDenseMatrix<OrdinalType, ScalarType> > FactorR_;
00375 
00376     ScalarType * A_;
00377     ScalarType * AF_;
00378     ScalarType * Q_;
00379     ScalarType * R_;
00380     std::vector<ScalarType> WORK_;
00381 #ifdef HAVE_TEUCHOSNUMERICS_EIGEN
00382     Eigen::HouseholderQR<EigenMatrix> qr_;
00383 #endif
00384 
00385   private:
00386     // SerialQRDenseSolver copy constructor (put here because we don't want user access)
00387 
00388     SerialQRDenseSolver(const SerialQRDenseSolver<OrdinalType, ScalarType>& Source);
00389     SerialQRDenseSolver & operator=(const SerialQRDenseSolver<OrdinalType, ScalarType>& Source);
00390 
00391   };
00392 
00393   namespace details {
00394 
00395     // Helper traits to distinguish work arrays for real and complex-valued datatypes.
00396     template<typename T>
00397     struct lapack_traits {
00398       typedef int iwork_type;
00399     };
00400 
00401     // Complex-valued specialization
00402     template<typename T>
00403     struct lapack_traits<std::complex<T> > {
00404       typedef typename ScalarTraits<T>::magnitudeType iwork_type;
00405     };
00406 
00407   } // end namespace details
00408 
00409 //=============================================================================
00410 
00411 template<typename OrdinalType, typename ScalarType>
00412 SerialQRDenseSolver<OrdinalType,ScalarType>::SerialQRDenseSolver()
00413   : CompObject(),
00414     Object("Teuchos::SerialQRDenseSolver"),
00415     equilibrate_(false),
00416     shouldEquilibrate_(false),
00417     equilibratedA_(false),
00418     equilibratedB_(false),
00419     equilibrationModeA_(0),
00420     equilibrationModeB_(0),
00421     transpose_(false),
00422     factored_(false),
00423     solved_(false),
00424     matrixIsZero_(false),
00425     formedQ_(false),
00426     formedR_(false),
00427     TRANS_(Teuchos::NO_TRANS),
00428     M_(0),
00429     N_(0),
00430     Min_MN_(0),
00431     LDA_(0),
00432     LDAF_(0),
00433     LDQ_(0),
00434     LDR_(0),
00435     INFO_(0),
00436     LWORK_(0),
00437     ANORM_(ScalarTraits<MagnitudeType>::zero()),
00438     BNORM_(ScalarTraits<MagnitudeType>::zero()),
00439     A_(0),
00440     AF_(0),
00441     Q_(0),
00442     R_(0)
00443 {
00444   resetMatrix();
00445 }
00446 //=============================================================================
00447 
00448 template<typename OrdinalType, typename ScalarType>
00449 SerialQRDenseSolver<OrdinalType,ScalarType>::~SerialQRDenseSolver()
00450 {}
00451 
00452 //=============================================================================
00453 
00454 template<typename OrdinalType, typename ScalarType>
00455 void SerialQRDenseSolver<OrdinalType,ScalarType>::resetVectors()
00456 {
00457   LHS_ = Teuchos::null;
00458   TMP_ = Teuchos::null;
00459   RHS_ = Teuchos::null;
00460   solved_ = false;
00461   equilibratedB_ = false;
00462 }
00463 //=============================================================================
00464 
00465 template<typename OrdinalType, typename ScalarType>
00466 void SerialQRDenseSolver<OrdinalType,ScalarType>::resetMatrix()
00467 {
00468   resetVectors();
00469   equilibratedA_ = false;
00470   equilibrationModeA_ = 0;
00471   equilibrationModeB_ = 0;
00472   factored_ = false;
00473   matrixIsZero_ = false;
00474   formedQ_ = false;
00475   formedR_ = false;
00476   M_ = 0;
00477   N_ = 0;
00478   Min_MN_ = 0;
00479   LDA_ = 0;
00480   LDAF_ = 0;
00481   LDQ_ = 0;
00482   LDR_ = 0;
00483   ANORM_ = -ScalarTraits<MagnitudeType>::one();
00484   BNORM_ = -ScalarTraits<MagnitudeType>::one();
00485   A_ = 0;
00486   AF_ = 0;
00487   Q_ = 0;
00488   R_ = 0;
00489   INFO_ = 0;
00490   LWORK_ = 0;
00491 }
00492 //=============================================================================
00493 
00494 template<typename OrdinalType, typename ScalarType>
00495 int SerialQRDenseSolver<OrdinalType,ScalarType>::setMatrix(const RCP<SerialDenseMatrix<OrdinalType,ScalarType> >& A)
00496 {
00497   TEUCHOS_TEST_FOR_EXCEPTION(A->numRows() < A->numCols(), std::invalid_argument,
00498          "SerialQRDenseSolver<T>::setMatrix: the matrix A must have A.numRows() >= A.numCols()!");
00499 
00500   resetMatrix();
00501   Matrix_ = A;
00502   Factor_ = A;
00503   FactorQ_ = A;
00504   FactorR_ = A;
00505   M_ = A->numRows();
00506   N_ = A->numCols();
00507   Min_MN_ = TEUCHOS_MIN(M_,N_);
00508   LDA_ = A->stride();
00509   LDAF_ = LDA_;
00510   LDQ_ = LDA_;
00511   LDR_ = N_;
00512   A_ = A->values();
00513   AF_ = A->values();
00514   Q_ = A->values();
00515   R_ = A->values();
00516 
00517   return(0);
00518 }
00519 //=============================================================================
00520 
00521 template<typename OrdinalType, typename ScalarType>
00522 int SerialQRDenseSolver<OrdinalType,ScalarType>::setVectors(const RCP<SerialDenseMatrix<OrdinalType,ScalarType> >& X,
00523                  const RCP<SerialDenseMatrix<OrdinalType,ScalarType> >& B)
00524 {
00525   // Check that these new vectors are consistent
00526   TEUCHOS_TEST_FOR_EXCEPTION(B->numCols() != X->numCols(), std::invalid_argument,
00527          "SerialQRDenseSolver<T>::setVectors: X and B have different numbers of columns!");
00528   TEUCHOS_TEST_FOR_EXCEPTION(B->values()==0, std::invalid_argument,
00529          "SerialQRDenseSolver<T>::setVectors: B is an empty SerialDenseMatrix<T>!");
00530   TEUCHOS_TEST_FOR_EXCEPTION(X->values()==0, std::invalid_argument,
00531          "SerialQRDenseSolver<T>::setVectors: X is an empty SerialDenseMatrix<T>!");
00532   TEUCHOS_TEST_FOR_EXCEPTION(B->stride()<1, std::invalid_argument,
00533          "SerialQRDenseSolver<T>::setVectors: B has an invalid stride!");
00534   TEUCHOS_TEST_FOR_EXCEPTION(X->stride()<1, std::invalid_argument,
00535          "SerialQRDenseSolver<T>::setVectors: X has an invalid stride!");
00536 
00537   resetVectors();
00538   LHS_ = X;
00539   RHS_ = B;
00540 
00541   return(0);
00542 }
00543 //=============================================================================
00544 
00545 template<typename OrdinalType, typename ScalarType>
00546 int SerialQRDenseSolver<OrdinalType,ScalarType>::factor() {
00547 
00548   if (factored()) return(0);
00549 
00550   // Equilibrate matrix if necessary
00551   int ierr = 0;
00552   if (equilibrate_) ierr = equilibrateMatrix();
00553   if (ierr!=0) return(ierr);
00554 
00555   allocateWORK();
00556   if ((int)TAU_.size()!=Min_MN_) TAU_.resize( Min_MN_ );
00557 
00558   INFO_ = 0;
00559 
00560   // Factor
00561 #ifdef HAVE_TEUCHOSNUMERICS_EIGEN
00562   EigenMatrixMap aMap( AF_, M_, N_, EigenOuterStride(LDAF_) );
00563   EigenScalarArray tau;
00564   EigenScalarArrayMap tauMap(&TAU_[0],TEUCHOS_MIN(M_,N_));
00565   qr_.compute(aMap);
00566   tau = qr_.hCoeffs();
00567   for (OrdinalType i=0; i<tauMap.innerSize(); i++) {
00568     tauMap(i) = tau(i);
00569   }
00570   EigenMatrix qrMat = qr_.matrixQR();
00571   for (OrdinalType j=0; j<aMap.outerSize(); j++) {
00572     for (OrdinalType i=0; i<aMap.innerSize(); i++) {
00573       aMap(i,j) = qrMat(i,j);
00574     }
00575   }
00576 #else
00577   this->GEQRF(M_, N_, AF_, LDAF_, &TAU_[0], &WORK_[0], LWORK_, &INFO_);
00578 #endif
00579 
00580   factored_ = true;
00581 
00582   return(INFO_);
00583 }
00584 
00585 //=============================================================================
00586 
00587 template<typename OrdinalType, typename ScalarType>
00588 int SerialQRDenseSolver<OrdinalType,ScalarType>::solve() {
00589 
00590   // Check if the matrix is zero
00591   if (matrixIsZero_) {
00592     LHS_->putScalar(ScalarTraits<ScalarType>::zero());
00593     return(0);
00594   }
00595 
00596   // Equilibrate RHS if necessary
00597   int ierr = 0;
00598   if (equilibrate_) {
00599     ierr = equilibrateRHS();
00600     equilibratedB_ = true;
00601   }
00602   if (ierr != 0) return(ierr);
00603 
00604   TEUCHOS_TEST_FOR_EXCEPTION( (equilibratedA_ && !equilibratedB_) || (!equilibratedA_ && equilibratedB_) ,
00605            std::logic_error, "SerialQRDenseSolver<T>::solve: Matrix and vectors must be similarly scaled!");
00606   TEUCHOS_TEST_FOR_EXCEPTION( RHS_==Teuchos::null, std::invalid_argument,
00607          "SerialQRDenseSolver<T>::solve: No right-hand side vector (RHS) has been set for the linear system!");
00608   TEUCHOS_TEST_FOR_EXCEPTION( LHS_==Teuchos::null, std::invalid_argument,
00609          "SerialQRDenseSolver<T>::solve: No solution vector (LHS) has been set for the linear system!");
00610   if ( TRANS_ == Teuchos::NO_TRANS ) {
00611     TEUCHOS_TEST_FOR_EXCEPTION( LHS_->numRows() != N_, std::invalid_argument,
00612            "SerialQRDenseSolver<T>::solve: No transpose: must have LHS_->numRows() = N_");
00613   } else {
00614     TEUCHOS_TEST_FOR_EXCEPTION( LHS_->numRows() != M_, std::invalid_argument,
00615            "SerialQRDenseSolver<T>::solve: Transpose: must have LHS_->numRows() = M_");
00616   }
00617 
00618   if (shouldEquilibrate() && !equilibratedA_)
00619     std::cout << "WARNING!  SerialQRDenseSolver<T>::solve: System should be equilibrated!" << std::endl;
00620 
00621   // Matrix must be factored
00622   if (!factored()) factor();
00623 
00624   TMP_ = rcp( new SerialDenseMatrix<OrdinalType,ScalarType>(M_, RHS_->numCols()) );
00625   for (OrdinalType j=0; j<RHS_->numCols(); j++) {
00626     for (OrdinalType i=0; i<RHS_->numRows(); i++) {
00627       (*TMP_)(i,j) = (*RHS_)(i,j);
00628     }
00629   }
00630 
00631   INFO_ = 0;
00632 
00633   // Solve
00634   if ( TRANS_ == Teuchos::NO_TRANS ) {
00635 
00636     // Solve A lhs = rhs
00637     this->multiplyQ( Teuchos::CONJ_TRANS, *TMP_ );
00638     this->solveR( Teuchos::NO_TRANS, *TMP_ );
00639 
00640   } else {
00641 
00642     // Solve A**H lhs = rhs
00643     this->solveR( Teuchos::CONJ_TRANS, *TMP_ );
00644     for (OrdinalType j = 0; j < RHS_->numCols(); j++ ) {
00645       for (OrdinalType i = N_; i < M_; i++ ) {
00646   (*TMP_)(i, j) = ScalarTraits<ScalarType>::zero();
00647       }
00648     }
00649     this->multiplyQ( Teuchos::NO_TRANS, *TMP_ );
00650 
00651   }
00652   for (OrdinalType j = 0; j < LHS_->numCols(); j++ ) {
00653     for (OrdinalType i = 0; i < LHS_->numRows(); i++ ) {
00654       (*LHS_)(i, j) = (*TMP_)(i,j);
00655     }
00656   }
00657 
00658   solved_ = true;
00659 
00660   // Unequilibrate LHS if necessary
00661   int ierr1=0;
00662   if (equilibrate_) ierr1 = unequilibrateLHS();
00663   if (ierr != 0) return(ierr);
00664 
00665   return(INFO_);
00666 
00667 }
00668 
00669 //=============================================================================
00670 
00671 template<typename OrdinalType, typename ScalarType>
00672 int SerialQRDenseSolver<OrdinalType,ScalarType>::computeEquilibrateScaling()
00673 {
00674   MagnitudeType safeMin = ScalarTraits<ScalarType>::sfmin();
00675   MagnitudeType prec = ScalarTraits<ScalarType>::prec();
00676   ScalarType sZero = ScalarTraits<ScalarType>::zero();
00677   ScalarType sOne  = ScalarTraits<ScalarType>::one();
00678   MagnitudeType mZero = ScalarTraits<ScalarType>::magnitude(sZero);
00679 
00680   MagnitudeType smlnum = ScalarTraits<ScalarType>::magnitude(safeMin/prec);
00681   MagnitudeType bignum = ScalarTraits<ScalarType>::magnitude(sOne/smlnum);
00682 
00683   // Compute maximum absolute value of matrix entries
00684   OrdinalType i, j;
00685   MagnitudeType anorm = ScalarTraits<ScalarType>::magnitude(ScalarTraits<ScalarType>::zero());
00686   for (j = 0; j < N_; j++) {
00687     for (i = 0; i < M_; i++) {
00688       anorm = TEUCHOS_MAX( anorm, ScalarTraits<ScalarType>::magnitude((*Matrix_)(i,j)) );
00689     }
00690   }
00691 
00692   ANORM_ = anorm;
00693 
00694   if (ANORM_ > mZero && ANORM_ < smlnum) {
00695     // Scale matrix norm up to smlnum
00696     shouldEquilibrate_ = true;
00697   } else if (ANORM_ > bignum) {
00698     // Scale matrix norm down to bignum
00699     shouldEquilibrate_ = true;
00700   } else if (ANORM_ == mZero) {
00701     // Matrix all zero. Return zero solution
00702     matrixIsZero_ = true;
00703   }
00704 
00705   return(0);
00706 }
00707 
00708 //=============================================================================
00709 
00710 template<typename OrdinalType, typename ScalarType>
00711 int SerialQRDenseSolver<OrdinalType,ScalarType>::equilibrateMatrix()
00712 {
00713   if (equilibratedA_) return(0);
00714 
00715   MagnitudeType safeMin = ScalarTraits<ScalarType>::sfmin();
00716   MagnitudeType prec = ScalarTraits<ScalarType>::prec();
00717   ScalarType sZero = ScalarTraits<ScalarType>::zero();
00718   ScalarType sOne  = ScalarTraits<ScalarType>::one();
00719   MagnitudeType mZero = ScalarTraits<ScalarType>::magnitude(sZero);
00720 
00721   MagnitudeType smlnum = ScalarTraits<ScalarType>::magnitude(safeMin/prec);
00722   MagnitudeType bignum = ScalarTraits<ScalarType>::magnitude(sOne/smlnum);
00723   OrdinalType BW = 0;
00724 
00725   // Compute maximum absolute value of matrix entries
00726   OrdinalType i, j;
00727   MagnitudeType anorm = ScalarTraits<ScalarType>::magnitude(ScalarTraits<ScalarType>::zero());
00728   for (j = 0; j < N_; j++) {
00729     for (i = 0; i < M_; i++) {
00730       anorm = TEUCHOS_MAX( anorm, ScalarTraits<ScalarType>::magnitude((*Matrix_)(i,j)) );
00731     }
00732   }
00733 
00734   ANORM_ = anorm;
00735   int ierr1 = 0;
00736   if (ANORM_ > mZero && ANORM_ < smlnum) {
00737     // Scale matrix norm up to smlnum
00738     this->LASCL(Teuchos::ETypeChar[Teuchos::FULL], BW, BW, ANORM_, smlnum, M_, N_, A_, LDA_, &INFO_);
00739     equilibrationModeA_ = 1;
00740   } else if (ANORM_ > bignum) {
00741     // Scale matrix norm down to bignum
00742     this->LASCL(Teuchos::ETypeChar[Teuchos::FULL], BW, BW, ANORM_, bignum, M_, N_, A_, LDA_, &INFO_);
00743     equilibrationModeA_ = 2;
00744   } else if (ANORM_ == mZero) {
00745     // Matrix all zero. Return zero solution
00746     matrixIsZero_ = true;
00747   }
00748 
00749   equilibratedA_ = true;
00750 
00751   return(ierr1);
00752 }
00753 
00754 //=============================================================================
00755 
00756 template<typename OrdinalType, typename ScalarType>
00757 int SerialQRDenseSolver<OrdinalType,ScalarType>::equilibrateRHS()
00758 {
00759   if (equilibratedB_) return(0);
00760 
00761   MagnitudeType safeMin = ScalarTraits<ScalarType>::sfmin();
00762   MagnitudeType prec = ScalarTraits<ScalarType>::prec();
00763   ScalarType sZero = ScalarTraits<ScalarType>::zero();
00764   ScalarType sOne  = ScalarTraits<ScalarType>::one();
00765   MagnitudeType mZero = ScalarTraits<ScalarType>::magnitude(sZero);
00766 
00767   MagnitudeType smlnum = ScalarTraits<ScalarType>::magnitude(safeMin/prec);
00768   MagnitudeType bignum = ScalarTraits<ScalarType>::magnitude(sOne/smlnum);
00769   OrdinalType BW = 0;
00770 
00771   // Compute maximum absolute value of rhs entries
00772   OrdinalType i, j;
00773   MagnitudeType bnorm = ScalarTraits<ScalarType>::magnitude(ScalarTraits<ScalarType>::zero());
00774   for (j = 0; j <RHS_->numCols(); j++) {
00775     for (i = 0; i < RHS_->numRows(); i++) {
00776       bnorm = TEUCHOS_MAX( bnorm, ScalarTraits<ScalarType>::magnitude((*RHS_)(i,j)) );
00777     }
00778   }
00779 
00780   BNORM_ = bnorm;
00781 
00782   int ierr1 = 0;
00783   if (BNORM_ > mZero && BNORM_ < smlnum) {
00784     // Scale RHS norm up to smlnum
00785     this->LASCL(Teuchos::ETypeChar[Teuchos::FULL], BW, BW, BNORM_, smlnum, RHS_->numRows(), RHS_->numCols(),
00786         RHS_->values(), RHS_->stride(), &INFO_);
00787     equilibrationModeB_ = 1;
00788   } else if (BNORM_ > bignum) {
00789     // Scale RHS norm down to bignum
00790     this->LASCL(Teuchos::ETypeChar[Teuchos::FULL], BW, BW, BNORM_, bignum, RHS_->numRows(), RHS_->numCols(),
00791         RHS_->values(), RHS_->stride(), &INFO_);
00792     equilibrationModeB_ = 2;
00793   }
00794 
00795   equilibratedB_ = true;
00796 
00797   return(ierr1);
00798 }
00799 
00800 //=============================================================================
00801 
00802 template<typename OrdinalType, typename ScalarType>
00803 int SerialQRDenseSolver<OrdinalType,ScalarType>::unequilibrateLHS()
00804 {
00805   if (!equilibratedB_) return(0);
00806 
00807   MagnitudeType safeMin = ScalarTraits<ScalarType>::sfmin();
00808   MagnitudeType prec = ScalarTraits<ScalarType>::prec();
00809   ScalarType sZero = ScalarTraits<ScalarType>::zero();
00810   ScalarType sOne  = ScalarTraits<ScalarType>::one();
00811   MagnitudeType mZero = ScalarTraits<ScalarType>::magnitude(sZero);
00812 
00813   MagnitudeType smlnum = ScalarTraits<ScalarType>::magnitude(safeMin/prec);
00814   MagnitudeType bignum = ScalarTraits<ScalarType>::magnitude(sOne/smlnum);
00815   OrdinalType BW = 0;
00816 
00817   int ierr1 = 0;
00818   if (equilibrationModeA_ == 1) {
00819     this->LASCL(Teuchos::ETypeChar[Teuchos::FULL], BW, BW, ANORM_, smlnum, LHS_->numRows(), LHS_->numCols(),
00820     LHS_->values(), LHS_->stride(), &INFO_);
00821   } else if (equilibrationModeA_ == 2) {
00822     this->LASCL(Teuchos::ETypeChar[Teuchos::FULL], BW, BW, ANORM_, bignum, LHS_->numRows(), LHS_->numCols(),
00823     LHS_->values(), LHS_->stride(), &INFO_);
00824   }
00825   if (equilibrationModeB_ == 1) {
00826     this->LASCL(Teuchos::ETypeChar[Teuchos::FULL], BW, BW, smlnum, BNORM_, LHS_->numRows(), LHS_->numCols(),
00827     LHS_->values(), LHS_->stride(), &INFO_);
00828   } else if (equilibrationModeB_ == 2) {
00829     this->LASCL(Teuchos::ETypeChar[Teuchos::FULL], BW, BW, bignum, BNORM_, LHS_->numRows(), LHS_->numCols(),
00830     LHS_->values(), LHS_->stride(), &INFO_);
00831   }
00832 
00833   return(ierr1);
00834 }
00835 
00836 //=============================================================================
00837 
00838 template<typename OrdinalType, typename ScalarType>
00839 int SerialQRDenseSolver<OrdinalType,ScalarType>::formQ() {
00840 
00841   // Matrix must be factored first
00842   if (!factored()) factor();
00843 
00844   // Store Q separately from factored matrix
00845   if (AF_ == Q_) {
00846     FactorQ_ = rcp( new SerialDenseMatrix<OrdinalType,ScalarType>(*Factor_) );
00847     Q_ = FactorQ_->values();
00848     LDQ_ = FactorQ_->stride();
00849   }
00850 
00851   INFO_ = 0;
00852 
00853   // Form Q
00854 #ifdef HAVE_TEUCHOSNUMERICS_EIGEN
00855   EigenMatrixMap qMap( Q_, M_, N_, EigenOuterStride(LDQ_) );
00856   EigenMatrix qMat = qr_.householderQ();
00857   for (OrdinalType j=0; j<qMap.outerSize(); j++) {
00858     for (OrdinalType i=0; i<qMap.innerSize(); i++) {
00859       qMap(i,j) = qMat(i,j);
00860     }
00861   }
00862 #else
00863   this->UNGQR(M_, N_, N_, Q_, LDQ_, &TAU_[0], &WORK_[0], LWORK_, &INFO_);
00864 #endif
00865 
00866   if (INFO_!=0) return(INFO_);
00867 
00868   formedQ_ = true;
00869 
00870   return(INFO_);
00871 }
00872 
00873 //=============================================================================
00874 
00875 template<typename OrdinalType, typename ScalarType>
00876 int SerialQRDenseSolver<OrdinalType,ScalarType>::formR() {
00877 
00878   // Matrix must be factored first
00879   if (!factored()) factor();
00880 
00881   // Store R separately from factored matrix
00882   if (AF_ == R_) {
00883     FactorR_ = rcp( new SerialDenseMatrix<OrdinalType,ScalarType>(N_, N_) );
00884     R_ = FactorR_->values();
00885     LDR_ = FactorR_->stride();
00886   }
00887 
00888   // Form R
00889   for (OrdinalType j=0; j<N_; j++) {
00890     for (OrdinalType i=0; i<=j; i++) {
00891       (*FactorR_)(i,j) = (*Factor_)(i,j);
00892     }
00893   }
00894 
00895   formedR_ = true;
00896 
00897   return(0);
00898 }
00899 
00900 //=============================================================================
00901 
00902 template<typename OrdinalType, typename ScalarType>
00903 int  SerialQRDenseSolver<OrdinalType, ScalarType>::multiplyQ(ETransp transq, SerialDenseMatrix<OrdinalType, ScalarType> &C)
00904 {
00905 
00906   // Check that the matrices are consistent.
00907   TEUCHOS_TEST_FOR_EXCEPTION(C.numRows()!=M_, std::invalid_argument,
00908          "SerialQRDenseSolver<T>::multiplyQ: C.numRows() != M_!");
00909   TEUCHOS_TEST_FOR_EXCEPTION(C.values()==0, std::invalid_argument,
00910          "SerialQRDenseSolver<T>::multiplyQ: C is an empty SerialDenseMatrix<T>!");
00911 
00912   // Matrix must be factored
00913   if (!factored()) factor();
00914 
00915   INFO_ = 0;
00916 
00917   // Multiply
00918 #ifdef HAVE_TEUCHOSNUMERICS_EIGEN
00919   EigenMatrixMap cMap( C.values(), C.numRows(), C.numCols(), EigenOuterStride(C.stride()) );
00920   if ( transq == Teuchos::NO_TRANS ) {
00921     // C = Q * C
00922     cMap = qr_.householderQ() * cMap;
00923   } else {
00924     // C = Q**H * C
00925     cMap = qr_.householderQ().adjoint() * cMap;
00926     for (OrdinalType j = 0; j < C.numCols(); j++) {
00927       for (OrdinalType i = N_; i < C.numRows(); i++ ) {
00928   cMap(i, j) = ScalarTraits<ScalarType>::zero();
00929       }
00930     }
00931   }
00932 #else
00933   Teuchos::ETransp NO_TRANS = Teuchos::NO_TRANS;
00934   Teuchos::ETransp TRANS = (Teuchos::ScalarTraits<ScalarType>::isComplex) ? Teuchos::CONJ_TRANS : Teuchos::TRANS;
00935   Teuchos::ESide SIDE = Teuchos::LEFT_SIDE;
00936 
00937   if ( transq == Teuchos::NO_TRANS ) {
00938 
00939     // C = Q * C
00940     this->UNMQR(ESideChar[SIDE], ETranspChar[NO_TRANS], M_, C.numCols(), N_, AF_, LDAF_,
00941     &TAU_[0], C.values(), C.stride(), &WORK_[0], LWORK_, &INFO_);
00942 
00943   } else {
00944 
00945     // C = Q**H * C
00946     this->UNMQR(ESideChar[SIDE], ETranspChar[TRANS], M_, C.numCols(), N_, AF_, LDAF_,
00947     &TAU_[0], C.values(), C.stride(), &WORK_[0], LWORK_, &INFO_);
00948 
00949     for (OrdinalType j = 0; j < C.numCols(); j++) {
00950       for (OrdinalType i = N_; i < C.numRows(); i++ ) {
00951   C(i, j) = ScalarTraits<ScalarType>::zero();
00952       }
00953     }
00954   }
00955 #endif
00956 
00957   return(INFO_);
00958 
00959 }
00960 
00961 //=============================================================================
00962 
00963 template<typename OrdinalType, typename ScalarType>
00964 int  SerialQRDenseSolver<OrdinalType, ScalarType>::solveR(ETransp transr, SerialDenseMatrix<OrdinalType, ScalarType> &C)
00965 {
00966 
00967   // Check that the matrices are consistent.
00968   TEUCHOS_TEST_FOR_EXCEPTION(C.numRows()<N_ || C.numRows()>M_, std::invalid_argument,
00969          "SerialQRDenseSolver<T>::solveR: must have N_ < C.numRows() < M_!");
00970   TEUCHOS_TEST_FOR_EXCEPTION(C.values()==0, std::invalid_argument,
00971          "SerialQRDenseSolver<T>::solveR: C is an empty SerialDenseMatrix<T>!");
00972 
00973   // Matrix must be factored
00974   if (!factored()) factor();
00975 
00976   INFO_ = 0;
00977 
00978   // Solve
00979 #ifdef HAVE_TEUCHOSNUMERICS_EIGEN
00980   EigenMatrixMap cMap( C.values(), N_, C.numCols(), EigenOuterStride(C.stride()) );
00981   // Check for singularity first like TRTRS
00982   for (OrdinalType j=0; j<N_; j++) {
00983     if ((qr_.matrixQR())(j,j) == ScalarTraits<ScalarType>::zero()) {
00984       INFO_ = j+1;
00985       return(INFO_);
00986     }
00987   }
00988   if ( transr == Teuchos::NO_TRANS ) {
00989     // C = R \ C
00990     qr_.matrixQR().topRows(N_).template triangularView<Eigen::Upper>().solveInPlace(cMap);
00991   } else {
00992     // C = R**H \ C
00993     qr_.matrixQR().topRows(N_).template triangularView<Eigen::Upper>().adjoint().solveInPlace(cMap);
00994   }
00995 #else
00996   Teuchos::ETransp NO_TRANS = Teuchos::NO_TRANS;
00997   Teuchos::ETransp TRANS = (Teuchos::ScalarTraits<ScalarType>::isComplex) ? Teuchos::CONJ_TRANS : Teuchos::TRANS;
00998   Teuchos::EUplo UPLO = Teuchos::UPPER_TRI;
00999   Teuchos::EDiag DIAG = Teuchos::NON_UNIT_DIAG;
01000 
01001   ScalarType* RMAT = (formedR_) ? R_ : AF_;
01002   OrdinalType LDRMAT = (formedR_) ? LDR_ : LDAF_;
01003 
01004   if ( transr == Teuchos::NO_TRANS ) {
01005 
01006     // C = R \ C
01007     this->TRTRS(EUploChar[UPLO], ETranspChar[NO_TRANS], EDiagChar[DIAG], N_, C.numCols(),
01008     RMAT, LDRMAT, C.values(), C.stride(), &INFO_);
01009 
01010   } else {
01011 
01012     // C = R**H \ C
01013     this->TRTRS(EUploChar[UPLO], ETranspChar[TRANS], EDiagChar[DIAG], N_, C.numCols(),
01014     RMAT, LDRMAT, C.values(), C.stride(), &INFO_);
01015 
01016   }
01017 #endif
01018 
01019   return(INFO_);
01020 
01021 }
01022 
01023 //=============================================================================
01024 
01025 template<typename OrdinalType, typename ScalarType>
01026 void SerialQRDenseSolver<OrdinalType,ScalarType>::Print(std::ostream& os) const {
01027 
01028   if (Matrix_!=Teuchos::null) os << "Solver Matrix"          << std::endl << *Matrix_ << std::endl;
01029   if (Factor_!=Teuchos::null) os << "Solver Factored Matrix" << std::endl << *Factor_ << std::endl;
01030   if (Q_!=Teuchos::null) os << "Solver Factor Q" << std::endl << *Q_ << std::endl;
01031   if (LHS_   !=Teuchos::null) os << "Solver LHS"             << std::endl << *LHS_    << std::endl;
01032   if (RHS_   !=Teuchos::null) os << "Solver RHS"             << std::endl << *RHS_    << std::endl;
01033 
01034 }
01035 
01036 } // namespace Teuchos
01037 
01038 #endif /* _TEUCHOS_SERIALQRDENSESOLVER_HPP_ */
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Defines