Teuchos Package Browser (Single Doxygen Collection)  Version of the Day
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
Teuchos_SerialQRDenseSolver.hpp
Go to the documentation of this file.
1 // @HEADER
2 // *****************************************************************************
3 // Teuchos: Common Tools Package
4 //
5 // Copyright 2004 NTESS and the Teuchos contributors.
6 // SPDX-License-Identifier: BSD-3-Clause
7 // *****************************************************************************
8 // @HEADER
9 
10 #ifndef _TEUCHOS_SERIALQRDENSESOLVER_HPP_
11 #define _TEUCHOS_SERIALQRDENSESOLVER_HPP_
12 
16 #include "Teuchos_CompObject.hpp"
17 #include "Teuchos_BLAS.hpp"
18 #include "Teuchos_LAPACK.hpp"
19 #include "Teuchos_RCP.hpp"
20 #include "Teuchos_ConfigDefs.hpp"
23 #include "Teuchos_ScalarTraits.hpp"
24 
25 #ifdef HAVE_TEUCHOSNUMERICS_EIGEN
26 #include "Eigen/Dense"
27 #endif
28 
97 namespace Teuchos {
98 
99  template<typename OrdinalType, typename ScalarType>
100  class SerialQRDenseSolver : public CompObject, public Object, public BLAS<OrdinalType, ScalarType>,
101  public LAPACK<OrdinalType, ScalarType>
102  {
103 
104  public:
105 
107 #ifdef HAVE_TEUCHOSNUMERICS_EIGEN
108  typedef typename Eigen::Matrix<ScalarType,Eigen::Dynamic,Eigen::Dynamic,Eigen::ColMajor> EigenMatrix;
109  typedef typename Eigen::Matrix<ScalarType,Eigen::Dynamic,1> EigenVector;
110  typedef typename Eigen::InnerStride<Eigen::Dynamic> EigenInnerStride;
111  typedef typename Eigen::OuterStride<Eigen::Dynamic> EigenOuterStride;
112  typedef typename Eigen::Map<EigenVector,0,EigenInnerStride > EigenVectorMap;
113  typedef typename Eigen::Map<const EigenVector,0,EigenInnerStride > EigenConstVectorMap;
114  typedef typename Eigen::Map<EigenMatrix,0,EigenOuterStride > EigenMatrixMap;
115  typedef typename Eigen::Map<const EigenMatrix,0,EigenOuterStride > EigenConstMatrixMap;
116  typedef typename Eigen::PermutationMatrix<Eigen::Dynamic,Eigen::Dynamic,OrdinalType> EigenPermutationMatrix;
117  typedef typename Eigen::Array<OrdinalType,Eigen::Dynamic,1> EigenOrdinalArray;
118  typedef typename Eigen::Map<EigenOrdinalArray> EigenOrdinalArrayMap;
119  typedef typename Eigen::Array<ScalarType,Eigen::Dynamic,1> EigenScalarArray;
120  typedef typename Eigen::Map<EigenScalarArray> EigenScalarArrayMap;
121 #endif
122 
124 
127 
129  virtual ~SerialQRDenseSolver();
131 
133 
134 
136 
139 
141 
147 
149 
150 
152 
154  void factorWithEquilibration(bool flag) {equilibrate_ = flag; return;}
155 
157  void solveWithTranspose(bool flag) {transpose_ = flag; if (flag) TRANS_ = Teuchos::CONJ_TRANS; else TRANS_ = Teuchos::NO_TRANS; return;}
158 
160  void solveWithTransposeFlag(Teuchos::ETransp trans) {TRANS_ = trans; transpose_ = (trans != Teuchos::NO_TRANS) ? true : false; }
161 
163 
165 
166 
168 
171  int factor();
172 
174 
177  int solve();
178 
180 
184 
186 
190  int equilibrateMatrix();
191 
193 
197  int equilibrateRHS();
198 
200 
204  int unequilibrateLHS();
205 
207 
210  int formQ();
211 
213 
216  int formR();
217 
219 
223 
225 
230 
232 
233 
235  bool transpose() {return(transpose_);}
236 
238  bool factored() {return(factored_);}
239 
241  bool equilibratedA() {return(equilibratedA_);}
242 
244  bool equilibratedB() {return(equilibratedB_);}
245 
248 
250  bool solved() {return(solved_);}
251 
253  bool formedQ() {return(formedQ_);}
254 
256  bool formedR() {return(formedR_);}
257 
259 
261 
262 
265 
268 
271 
274 
277 
280 
282  OrdinalType numRows() const {return(M_);}
283 
285  OrdinalType numCols() const {return(N_);}
286 
288  std::vector<ScalarType> tau() const {return(TAU_);}
289 
291  MagnitudeType ANORM() const {return(ANORM_);}
292 
294 
296 
297  void Print(std::ostream& os) const;
300  protected:
301 
302  void allocateWORK() { LWORK_ = 4*N_; WORK_.resize( LWORK_ ); return;}
303  void resetMatrix();
304  void resetVectors();
305 
306 
311  OrdinalType equilibrationModeA_;
312  OrdinalType equilibrationModeB_;
314  bool factored_;
315  bool solved_;
317  bool formedQ_;
318  bool formedR_;
319 
321 
322  OrdinalType M_;
323  OrdinalType N_;
324  OrdinalType Min_MN_;
325  OrdinalType LDA_;
326  OrdinalType LDAF_;
327  OrdinalType LDQ_;
328  OrdinalType LDR_;
329  OrdinalType INFO_;
330  OrdinalType LWORK_;
331 
332  std::vector<ScalarType> TAU_;
333 
336 
344 
345  ScalarType * A_;
346  ScalarType * AF_;
347  ScalarType * Q_;
348  ScalarType * R_;
349  std::vector<ScalarType> WORK_;
350 #ifdef HAVE_TEUCHOSNUMERICS_EIGEN
351  Eigen::HouseholderQR<EigenMatrix> qr_;
352 #endif
353 
354  private:
355  // SerialQRDenseSolver copy constructor (put here because we don't want user access)
356 
359 
360  };
361 
362  // Helper traits to distinguish work arrays for real and complex-valued datatypes.
363  using namespace details;
364 
365 //=============================================================================
366 
367 template<typename OrdinalType, typename ScalarType>
369  : CompObject(),
370  equilibrate_(false),
371  shouldEquilibrate_(false),
372  equilibratedA_(false),
373  equilibratedB_(false),
374  equilibrationModeA_(0),
375  equilibrationModeB_(0),
376  transpose_(false),
377  factored_(false),
378  solved_(false),
379  matrixIsZero_(false),
380  formedQ_(false),
381  formedR_(false),
382  TRANS_(Teuchos::NO_TRANS),
383  M_(0),
384  N_(0),
385  Min_MN_(0),
386  LDA_(0),
387  LDAF_(0),
388  LDQ_(0),
389  LDR_(0),
390  INFO_(0),
391  LWORK_(0),
392  ANORM_(ScalarTraits<MagnitudeType>::zero()),
393  BNORM_(ScalarTraits<MagnitudeType>::zero()),
394  A_(0),
395  AF_(0),
396  Q_(0),
397  R_(0)
398 {
399  resetMatrix();
400 }
401 //=============================================================================
402 
403 template<typename OrdinalType, typename ScalarType>
405 {}
406 
407 //=============================================================================
408 
409 template<typename OrdinalType, typename ScalarType>
411 {
412  LHS_ = Teuchos::null;
413  TMP_ = Teuchos::null;
414  RHS_ = Teuchos::null;
415  solved_ = false;
416  equilibratedB_ = false;
417 }
418 //=============================================================================
419 
420 template<typename OrdinalType, typename ScalarType>
422 {
423  resetVectors();
424  equilibratedA_ = false;
425  equilibrationModeA_ = 0;
426  equilibrationModeB_ = 0;
427  factored_ = false;
428  matrixIsZero_ = false;
429  formedQ_ = false;
430  formedR_ = false;
431  M_ = 0;
432  N_ = 0;
433  Min_MN_ = 0;
434  LDA_ = 0;
435  LDAF_ = 0;
436  LDQ_ = 0;
437  LDR_ = 0;
440  A_ = 0;
441  AF_ = 0;
442  Q_ = 0;
443  R_ = 0;
444  INFO_ = 0;
445  LWORK_ = 0;
446 }
447 //=============================================================================
448 
449 template<typename OrdinalType, typename ScalarType>
451 {
452  TEUCHOS_TEST_FOR_EXCEPTION(A->numRows() < A->numCols(), std::invalid_argument,
453  "SerialQRDenseSolver<T>::setMatrix: the matrix A must have A.numRows() >= A.numCols()!");
454 
455  resetMatrix();
456  Matrix_ = A;
457  Factor_ = A;
458  FactorQ_ = A;
459  FactorR_ = A;
460  M_ = A->numRows();
461  N_ = A->numCols();
462  Min_MN_ = TEUCHOS_MIN(M_,N_);
463  LDA_ = A->stride();
464  LDAF_ = LDA_;
465  LDQ_ = LDA_;
466  LDR_ = N_;
467  A_ = A->values();
468  AF_ = A->values();
469  Q_ = A->values();
470  R_ = A->values();
471 
472  return(0);
473 }
474 //=============================================================================
475 
476 template<typename OrdinalType, typename ScalarType>
479 {
480  // Check that these new vectors are consistent
481  TEUCHOS_TEST_FOR_EXCEPTION(B->numCols() != X->numCols(), std::invalid_argument,
482  "SerialQRDenseSolver<T>::setVectors: X and B have different numbers of columns!");
483  TEUCHOS_TEST_FOR_EXCEPTION(B->values()==0, std::invalid_argument,
484  "SerialQRDenseSolver<T>::setVectors: B is an empty SerialDenseMatrix<T>!");
485  TEUCHOS_TEST_FOR_EXCEPTION(X->values()==0, std::invalid_argument,
486  "SerialQRDenseSolver<T>::setVectors: X is an empty SerialDenseMatrix<T>!");
487  TEUCHOS_TEST_FOR_EXCEPTION(B->stride()<1, std::invalid_argument,
488  "SerialQRDenseSolver<T>::setVectors: B has an invalid stride!");
489  TEUCHOS_TEST_FOR_EXCEPTION(X->stride()<1, std::invalid_argument,
490  "SerialQRDenseSolver<T>::setVectors: X has an invalid stride!");
491 
492  resetVectors();
493  LHS_ = X;
494  RHS_ = B;
495 
496  return(0);
497 }
498 //=============================================================================
499 
500 template<typename OrdinalType, typename ScalarType>
502 
503  if (factored()) return(0);
504 
505  // Equilibrate matrix if necessary
506  int ierr = 0;
507  if (equilibrate_) ierr = equilibrateMatrix();
508  if (ierr!=0) return(ierr);
509 
510  allocateWORK();
511  if ((int)TAU_.size()!=Min_MN_) TAU_.resize( Min_MN_ );
512 
513  INFO_ = 0;
514 
515  // Factor
516 #ifdef HAVE_TEUCHOSNUMERICS_EIGEN
517  EigenMatrixMap aMap( AF_, M_, N_, EigenOuterStride(LDAF_) );
518  EigenScalarArray tau;
519  EigenScalarArrayMap tauMap(&TAU_[0],TEUCHOS_MIN(M_,N_));
520  qr_.compute(aMap);
521  tau = qr_.hCoeffs();
522  for (OrdinalType i=0; i<tauMap.innerSize(); i++) {
523  tauMap(i) = tau(i);
524  }
525  EigenMatrix qrMat = qr_.matrixQR();
526  for (OrdinalType j=0; j<aMap.outerSize(); j++) {
527  for (OrdinalType i=0; i<aMap.innerSize(); i++) {
528  aMap(i,j) = qrMat(i,j);
529  }
530  }
531 #else
532  this->GEQRF(M_, N_, AF_, LDAF_, &TAU_[0], &WORK_[0], LWORK_, &INFO_);
533 #endif
534 
535  factored_ = true;
536 
537  return(INFO_);
538 }
539 
540 //=============================================================================
541 
542 template<typename OrdinalType, typename ScalarType>
544 
545  // Check if the matrix is zero
546  if (matrixIsZero_) {
547  LHS_->putScalar(ScalarTraits<ScalarType>::zero());
548  return(0);
549  }
550 
551  // Equilibrate RHS if necessary
552  int ierr = 0;
553  if (equilibrate_) {
554  ierr = equilibrateRHS();
555  }
556  if (ierr != 0) return(ierr);
557 
558  TEUCHOS_TEST_FOR_EXCEPTION( RHS_==Teuchos::null, std::invalid_argument,
559  "SerialQRDenseSolver<T>::solve: No right-hand side vector (RHS) has been set for the linear system!");
560  TEUCHOS_TEST_FOR_EXCEPTION( LHS_==Teuchos::null, std::invalid_argument,
561  "SerialQRDenseSolver<T>::solve: No solution vector (LHS) has been set for the linear system!");
562  if ( TRANS_ == Teuchos::NO_TRANS ) {
563  TEUCHOS_TEST_FOR_EXCEPTION( LHS_->numRows() != N_, std::invalid_argument,
564  "SerialQRDenseSolver<T>::solve: No transpose: must have LHS_->numRows() = N_");
565  } else {
566  TEUCHOS_TEST_FOR_EXCEPTION( LHS_->numRows() != M_, std::invalid_argument,
567  "SerialQRDenseSolver<T>::solve: Transpose: must have LHS_->numRows() = M_");
568  }
569 
570  if (shouldEquilibrate() && !equilibratedA_)
571  std::cout << "WARNING! SerialQRDenseSolver<T>::solve: System should be equilibrated!" << std::endl;
572 
573  // Matrix must be factored
574  if (!factored()) factor();
575 
576  TEUCHOS_TEST_FOR_EXCEPTION( (equilibratedA_ && !equilibratedB_) || (!equilibratedA_ && equilibratedB_) ,
577  std::logic_error, "SerialQRDenseSolver<T>::solve: Matrix and vectors must be similarly scaled!");
578 
579  TMP_ = rcp( new SerialDenseMatrix<OrdinalType,ScalarType>(M_, RHS_->numCols()) );
580  for (OrdinalType j=0; j<RHS_->numCols(); j++) {
581  for (OrdinalType i=0; i<RHS_->numRows(); i++) {
582  (*TMP_)(i,j) = (*RHS_)(i,j);
583  }
584  }
585 
586  INFO_ = 0;
587 
588  // Solve
589  if ( TRANS_ == Teuchos::NO_TRANS ) {
590 
591  // Solve A lhs = rhs
592  this->multiplyQ( Teuchos::CONJ_TRANS, *TMP_ );
593  this->solveR( Teuchos::NO_TRANS, *TMP_ );
594 
595  } else {
596 
597  // Solve A**H lhs = rhs
598  this->solveR( Teuchos::CONJ_TRANS, *TMP_ );
599  for (OrdinalType j = 0; j < RHS_->numCols(); j++ ) {
600  for (OrdinalType i = N_; i < M_; i++ ) {
601  (*TMP_)(i, j) = ScalarTraits<ScalarType>::zero();
602  }
603  }
604  this->multiplyQ( Teuchos::NO_TRANS, *TMP_ );
605 
606  }
607  for (OrdinalType j = 0; j < LHS_->numCols(); j++ ) {
608  for (OrdinalType i = 0; i < LHS_->numRows(); i++ ) {
609  (*LHS_)(i, j) = (*TMP_)(i,j);
610  }
611  }
612 
613  solved_ = true;
614 
615  // Unequilibrate LHS if necessary
616  if (equilibrate_) {
617  ierr = unequilibrateLHS();
618  }
619  if (ierr != 0) {
620  return ierr;
621  }
622 
623  return INFO_;
624 }
625 
626 //=============================================================================
627 
628 template<typename OrdinalType, typename ScalarType>
630 {
633  ScalarType sZero = ScalarTraits<ScalarType>::zero();
634  ScalarType sOne = ScalarTraits<ScalarType>::one();
636 
637  MagnitudeType smlnum = ScalarTraits<ScalarType>::magnitude(safeMin/prec);
639 
640  // Compute maximum absolute value of matrix entries
641  OrdinalType i, j;
643  for (j = 0; j < N_; j++) {
644  for (i = 0; i < M_; i++) {
645  anorm = TEUCHOS_MAX( anorm, ScalarTraits<ScalarType>::magnitude((*Matrix_)(i,j)) );
646  }
647  }
648 
649  ANORM_ = anorm;
650 
651  if (ANORM_ > mZero && ANORM_ < smlnum) {
652  // Scale matrix norm up to smlnum
653  shouldEquilibrate_ = true;
654  } else if (ANORM_ > bignum) {
655  // Scale matrix norm down to bignum
656  shouldEquilibrate_ = true;
657  } else if (ANORM_ == mZero) {
658  // Matrix all zero. Return zero solution
659  matrixIsZero_ = true;
660  }
661 
662  return(0);
663 }
664 
665 //=============================================================================
666 
667 template<typename OrdinalType, typename ScalarType>
669 {
670  if (equilibratedA_) return(0);
671 
674  ScalarType sZero = ScalarTraits<ScalarType>::zero();
675  ScalarType sOne = ScalarTraits<ScalarType>::one();
677 
678  MagnitudeType smlnum = ScalarTraits<ScalarType>::magnitude(safeMin/prec);
680  OrdinalType BW = 0;
681 
682  // Compute maximum absolute value of matrix entries
683  OrdinalType i, j;
685  for (j = 0; j < N_; j++) {
686  for (i = 0; i < M_; i++) {
687  anorm = TEUCHOS_MAX( anorm, ScalarTraits<ScalarType>::magnitude((*Matrix_)(i,j)) );
688  }
689  }
690 
691  ANORM_ = anorm;
692  int ierr1 = 0;
693  if (ANORM_ > mZero && ANORM_ < smlnum) {
694  // Scale matrix norm up to smlnum
695  this->LASCL(Teuchos::ETypeChar[Teuchos::FULL], BW, BW, ANORM_, smlnum, M_, N_, A_, LDA_, &INFO_);
696  equilibrationModeA_ = 1;
697  } else if (ANORM_ > bignum) {
698  // Scale matrix norm down to bignum
699  this->LASCL(Teuchos::ETypeChar[Teuchos::FULL], BW, BW, ANORM_, bignum, M_, N_, A_, LDA_, &INFO_);
700  equilibrationModeA_ = 2;
701  } else if (ANORM_ == mZero) {
702  // Matrix all zero. Return zero solution
703  matrixIsZero_ = true;
704  }
705 
706  equilibratedA_ = true;
707 
708  return(ierr1);
709 }
710 
711 //=============================================================================
712 
713 template<typename OrdinalType, typename ScalarType>
715 {
716  if (equilibratedB_) return(0);
717 
720  ScalarType sZero = ScalarTraits<ScalarType>::zero();
721  ScalarType sOne = ScalarTraits<ScalarType>::one();
723 
724  MagnitudeType smlnum = ScalarTraits<ScalarType>::magnitude(safeMin/prec);
726  OrdinalType BW = 0;
727 
728  // Compute maximum absolute value of rhs entries
729  OrdinalType i, j;
731  for (j = 0; j <RHS_->numCols(); j++) {
732  for (i = 0; i < RHS_->numRows(); i++) {
733  bnorm = TEUCHOS_MAX( bnorm, ScalarTraits<ScalarType>::magnitude((*RHS_)(i,j)) );
734  }
735  }
736 
737  BNORM_ = bnorm;
738 
739  int ierr1 = 0;
740  if (BNORM_ > mZero && BNORM_ < smlnum) {
741  // Scale RHS norm up to smlnum
742  this->LASCL(Teuchos::ETypeChar[Teuchos::FULL], BW, BW, BNORM_, smlnum, RHS_->numRows(), RHS_->numCols(),
743  RHS_->values(), RHS_->stride(), &INFO_);
744  equilibrationModeB_ = 1;
745  } else if (BNORM_ > bignum) {
746  // Scale RHS norm down to bignum
747  this->LASCL(Teuchos::ETypeChar[Teuchos::FULL], BW, BW, BNORM_, bignum, RHS_->numRows(), RHS_->numCols(),
748  RHS_->values(), RHS_->stride(), &INFO_);
749  equilibrationModeB_ = 2;
750  }
751 
752  equilibratedB_ = true;
753 
754  return(ierr1);
755 }
756 
757 //=============================================================================
758 
759 template<typename OrdinalType, typename ScalarType>
761 {
762  if (!equilibratedB_) return(0);
763 
766  ScalarType sZero = ScalarTraits<ScalarType>::zero();
767  ScalarType sOne = ScalarTraits<ScalarType>::one();
769  (void) mZero; // Silence "unused variable" compiler warning.
770 
771  MagnitudeType smlnum = ScalarTraits<ScalarType>::magnitude(safeMin/prec);
773  OrdinalType BW = 0;
774 
775  int ierr1 = 0;
776  if (equilibrationModeA_ == 1) {
777  this->LASCL(Teuchos::ETypeChar[Teuchos::FULL], BW, BW, ANORM_, smlnum, LHS_->numRows(), LHS_->numCols(),
778  LHS_->values(), LHS_->stride(), &INFO_);
779  } else if (equilibrationModeA_ == 2) {
780  this->LASCL(Teuchos::ETypeChar[Teuchos::FULL], BW, BW, ANORM_, bignum, LHS_->numRows(), LHS_->numCols(),
781  LHS_->values(), LHS_->stride(), &INFO_);
782  }
783  if (equilibrationModeB_ == 1) {
784  this->LASCL(Teuchos::ETypeChar[Teuchos::FULL], BW, BW, smlnum, BNORM_, LHS_->numRows(), LHS_->numCols(),
785  LHS_->values(), LHS_->stride(), &INFO_);
786  } else if (equilibrationModeB_ == 2) {
787  this->LASCL(Teuchos::ETypeChar[Teuchos::FULL], BW, BW, bignum, BNORM_, LHS_->numRows(), LHS_->numCols(),
788  LHS_->values(), LHS_->stride(), &INFO_);
789  }
790 
791  return(ierr1);
792 }
793 
794 //=============================================================================
795 
796 template<typename OrdinalType, typename ScalarType>
798 
799  // Matrix must be factored first
800  if (!factored()) factor();
801 
802  // Store Q separately from factored matrix
803  if (AF_ == Q_) {
804  FactorQ_ = rcp( new SerialDenseMatrix<OrdinalType,ScalarType>(*Factor_) );
805  Q_ = FactorQ_->values();
806  LDQ_ = FactorQ_->stride();
807  }
808 
809  INFO_ = 0;
810 
811  // Form Q
812 #ifdef HAVE_TEUCHOSNUMERICS_EIGEN
813  EigenMatrixMap qMap( Q_, M_, N_, EigenOuterStride(LDQ_) );
814  EigenMatrix qMat = qr_.householderQ();
815  for (OrdinalType j=0; j<qMap.outerSize(); j++) {
816  for (OrdinalType i=0; i<qMap.innerSize(); i++) {
817  qMap(i,j) = qMat(i,j);
818  }
819  }
820 #else
821  this->UNGQR(M_, N_, N_, Q_, LDQ_, &TAU_[0], &WORK_[0], LWORK_, &INFO_);
822 #endif
823 
824  if (INFO_!=0) return(INFO_);
825 
826  formedQ_ = true;
827 
828  return(INFO_);
829 }
830 
831 //=============================================================================
832 
833 template<typename OrdinalType, typename ScalarType>
835 
836  // Matrix must be factored first
837  if (!factored()) factor();
838 
839  // Store R separately from factored matrix
840  if (AF_ == R_) {
841  FactorR_ = rcp( new SerialDenseMatrix<OrdinalType,ScalarType>(N_, N_) );
842  R_ = FactorR_->values();
843  LDR_ = FactorR_->stride();
844  }
845 
846  // Form R
847  for (OrdinalType j=0; j<N_; j++) {
848  for (OrdinalType i=0; i<=j; i++) {
849  (*FactorR_)(i,j) = (*Factor_)(i,j);
850  }
851  }
852 
853  formedR_ = true;
854 
855  return(0);
856 }
857 
858 //=============================================================================
859 
860 template<typename OrdinalType, typename ScalarType>
862 {
863 
864  // Check that the matrices are consistent.
865  TEUCHOS_TEST_FOR_EXCEPTION(C.numRows()!=M_, std::invalid_argument,
866  "SerialQRDenseSolver<T>::multiplyQ: C.numRows() != M_!");
867  TEUCHOS_TEST_FOR_EXCEPTION(C.values()==0, std::invalid_argument,
868  "SerialQRDenseSolver<T>::multiplyQ: C is an empty SerialDenseMatrix<T>!");
869 
870  // Matrix must be factored
871  if (!factored()) factor();
872 
873  INFO_ = 0;
874 
875  // Multiply
876 #ifdef HAVE_TEUCHOSNUMERICS_EIGEN
877  EigenMatrixMap cMap( C.values(), C.numRows(), C.numCols(), EigenOuterStride(C.stride()) );
878  if ( transq == Teuchos::NO_TRANS ) {
879  // C = Q * C
880  cMap = qr_.householderQ() * cMap;
881  } else {
882  // C = Q**H * C
883  cMap = qr_.householderQ().adjoint() * cMap;
884  for (OrdinalType j = 0; j < C.numCols(); j++) {
885  for (OrdinalType i = N_; i < C.numRows(); i++ ) {
886  cMap(i, j) = ScalarTraits<ScalarType>::zero();
887  }
888  }
889  }
890 #else
894 
895  if ( transq == Teuchos::NO_TRANS ) {
896 
897  // C = Q * C
898  this->UNMQR(ESideChar[SIDE], ETranspChar[NO_TRANS], M_, C.numCols(), N_, AF_, LDAF_,
899  &TAU_[0], C.values(), C.stride(), &WORK_[0], LWORK_, &INFO_);
900 
901  } else {
902 
903  // C = Q**H * C
904  this->UNMQR(ESideChar[SIDE], ETranspChar[TRANS], M_, C.numCols(), N_, AF_, LDAF_,
905  &TAU_[0], C.values(), C.stride(), &WORK_[0], LWORK_, &INFO_);
906 
907  for (OrdinalType j = 0; j < C.numCols(); j++) {
908  for (OrdinalType i = N_; i < C.numRows(); i++ ) {
909  C(i, j) = ScalarTraits<ScalarType>::zero();
910  }
911  }
912  }
913 #endif
914 
915  return(INFO_);
916 
917 }
918 
919 //=============================================================================
920 
921 template<typename OrdinalType, typename ScalarType>
923 {
924 
925  // Check that the matrices are consistent.
926  TEUCHOS_TEST_FOR_EXCEPTION(C.numRows()<N_ || C.numRows()>M_, std::invalid_argument,
927  "SerialQRDenseSolver<T>::solveR: must have N_ < C.numRows() < M_!");
928  TEUCHOS_TEST_FOR_EXCEPTION(C.values()==0, std::invalid_argument,
929  "SerialQRDenseSolver<T>::solveR: C is an empty SerialDenseMatrix<T>!");
930 
931  // Matrix must be factored
932  if (!factored()) factor();
933 
934  INFO_ = 0;
935 
936  // Solve
937 #ifdef HAVE_TEUCHOSNUMERICS_EIGEN
938  EigenMatrixMap cMap( C.values(), N_, C.numCols(), EigenOuterStride(C.stride()) );
939  // Check for singularity first like TRTRS
940  for (OrdinalType j=0; j<N_; j++) {
941  if ((qr_.matrixQR())(j,j) == ScalarTraits<ScalarType>::zero()) {
942  INFO_ = j+1;
943  return(INFO_);
944  }
945  }
946  if ( transr == Teuchos::NO_TRANS ) {
947  // C = R \ C
948  qr_.matrixQR().topRows(N_).template triangularView<Eigen::Upper>().solveInPlace(cMap);
949  } else {
950  // C = R**H \ C
951  qr_.matrixQR().topRows(N_).template triangularView<Eigen::Upper>().adjoint().solveInPlace(cMap);
952  }
953 #else
958 
959  ScalarType* RMAT = (formedR_) ? R_ : AF_;
960  OrdinalType LDRMAT = (formedR_) ? LDR_ : LDAF_;
961 
962  if ( transr == Teuchos::NO_TRANS ) {
963 
964  // C = R \ C
965  this->TRTRS(EUploChar[UPLO], ETranspChar[NO_TRANS], EDiagChar[DIAG], N_, C.numCols(),
966  RMAT, LDRMAT, C.values(), C.stride(), &INFO_);
967 
968  } else {
969 
970  // C = R**H \ C
971  this->TRTRS(EUploChar[UPLO], ETranspChar[TRANS], EDiagChar[DIAG], N_, C.numCols(),
972  RMAT, LDRMAT, C.values(), C.stride(), &INFO_);
973 
974  }
975 #endif
976 
977  return(INFO_);
978 
979 }
980 
981 //=============================================================================
982 
983 template<typename OrdinalType, typename ScalarType>
985 
986  if (Matrix_!=Teuchos::null) os << "Solver Matrix" << std::endl << *Matrix_ << std::endl;
987  if (Factor_!=Teuchos::null) os << "Solver Factored Matrix" << std::endl << *Factor_ << std::endl;
988  if (Q_!=Teuchos::null) os << "Solver Factor Q" << std::endl << *Q_ << std::endl;
989  if (LHS_ !=Teuchos::null) os << "Solver LHS" << std::endl << *LHS_ << std::endl;
990  if (RHS_ !=Teuchos::null) os << "Solver RHS" << std::endl << *RHS_ << std::endl;
991 
992 }
993 
994 } // namespace Teuchos
995 
996 #endif /* _TEUCHOS_SERIALQRDENSESOLVER_HPP_ */
ScalarType * values() const
Data array access method.
RCP< SerialDenseMatrix< OrdinalType, ScalarType > > FactorR_
int unequilibrateLHS()
Unscales the solution vectors if equilibration was used to solve the system.
RCP< SerialDenseMatrix< OrdinalType, ScalarType > > getMatrix() const
Returns pointer to current matrix.
OrdinalType numRows() const
Returns row dimension of system.
int setMatrix(const RCP< SerialDenseMatrix< OrdinalType, ScalarType > > &A)
Sets the pointers for coefficient matrix.
Templated serial dense matrix class.
int multiplyQ(ETransp transq, SerialDenseMatrix< OrdinalType, ScalarType > &C)
Left multiply the input matrix by the unitary matrix Q or its adjoint.
int solveR(ETransp transr, SerialDenseMatrix< OrdinalType, ScalarType > &C)
Solve input matrix on the left with the upper triangular matrix R or its adjoint. ...
int formQ()
Explicitly forms the unitary matrix Q.
Templated interface class to BLAS routines.
bool formedQ()
Returns true if Q has been formed explicitly.
#define TEUCHOS_TEST_FOR_EXCEPTION(throw_exception_test, Exception, msg)
Macro for throwing an exception with breakpointing to ease debugging.
bool equilibratedA()
Returns true if factor is equilibrated (factor available via getFactoredMatrix()).
RCP< SerialDenseMatrix< OrdinalType, ScalarType > > getRHS() const
Returns pointer to current RHS.
RCP< SerialDenseMatrix< OrdinalType, ScalarType > > LHS_
Teuchos header file which uses auto-configuration information to include necessary C++ headers...
Templated BLAS wrapper.
RCP< SerialDenseMatrix< OrdinalType, ScalarType > > TMP_
SerialQRDenseSolver & operator=(const SerialQRDenseSolver< OrdinalType, ScalarType > &Source)
Object for storing data and providing functionality that is common to all computational classes...
bool factored()
Returns true if matrix is factored (factor available via getFactoredMatrix()).
int solve()
Computes the solution X to AX = B for the this matrix and the B provided to SetVectors()..
TEUCHOSNUMERICS_LIB_DLL_EXPORT const char EDiagChar[]
static magnitudeType sfmin()
Returns safe minimum (sfmin), such that 1/sfmin does not overflow.
This structure defines some basic traits for a scalar field type.
A class for solving dense linear problems.
Templated class for solving dense linear problems.
Functionality and data that is common to all computational classes.
RCP< SerialDenseMatrix< OrdinalType, ScalarType > > Matrix_
Templated interface class to LAPACK routines.
void solveWithTransposeFlag(Teuchos::ETransp trans)
All subsequent function calls will work with the transpose-type set by this method (Teuchos::NO_TRANS...
TEUCHOS_DEPRECATED RCP< T > rcp(T *p, Dealloc_T dealloc, bool owns_mem)
Deprecated.
static magnitudeType prec()
Returns eps*base.
RCP< SerialDenseMatrix< OrdinalType, ScalarType > > FactorQ_
TEUCHOSNUMERICS_LIB_DLL_EXPORT const char EUploChar[]
The base Teuchos class.
void Print(std::ostream &os) const
Print service methods; defines behavior of ostream &lt;&lt; operator.
TEUCHOSNUMERICS_LIB_DLL_EXPORT const char ETypeChar[]
int equilibrateRHS()
Equilibrates the current RHS.
RCP< SerialDenseMatrix< OrdinalType, ScalarType > > getFactoredMatrix() const
Returns pointer to factored matrix (assuming factorization has been performed).
ScalarTraits< ScalarType >::magnitudeType MagnitudeType
The Templated LAPACK Wrapper Class.
void solveWithTranspose(bool flag)
If flag is true, causes all subsequent function calls to work with the adjoint of this matrix...
#define TEUCHOS_MAX(x, y)
bool solved()
Returns true if the current set of vectors has been solved.
RCP< SerialDenseMatrix< OrdinalType, ScalarType > > Factor_
TEUCHOSNUMERICS_LIB_DLL_EXPORT const char ESideChar[]
SerialQRDenseSolver()
Default constructor; matrix should be set using setMatrix(), LHS and RHS set with setVectors()...
RCP< SerialDenseMatrix< OrdinalType, ScalarType > > RHS_
std::vector< ScalarType > tau() const
Returns pointer to pivot vector (if factorization has been computed), zero otherwise.
RCP< SerialDenseMatrix< OrdinalType, ScalarType > > getLHS() const
Returns pointer to current LHS.
bool equilibratedB()
Returns true if RHS is equilibrated (RHS available via getRHS()).
OrdinalType numCols() const
Returns column dimension of system.
static magnitudeType magnitude(T a)
Returns the magnitudeType of the scalar type a.
OrdinalType numCols() const
Returns the column dimension of this matrix.
bool shouldEquilibrate()
Returns true if the LAPACK general rules for equilibration suggest you should equilibrate the system...
MagnitudeType ANORM() const
Returns the absolute value of the largest element of this matrix (returns -1 if not yet computed)...
void factorWithEquilibration(bool flag)
Causes equilibration to be called just before the matrix factorization as part of the call to factor...
int equilibrateMatrix()
Equilibrates the this matrix.
Defines basic traits for the scalar field type.
static T zero()
Returns representation of zero for this scalar type.
TEUCHOSNUMERICS_LIB_DLL_EXPORT const char ETranspChar[]
virtual ~SerialQRDenseSolver()
SerialQRDenseSolver destructor.
int computeEquilibrateScaling()
Determines if this matrix should be scaled.
Smart reference counting pointer class for automatic garbage collection.
int formR()
Explicitly forms the upper triangular matrix R.
RCP< SerialDenseMatrix< OrdinalType, ScalarType > > getQ() const
Returns pointer to Q (assuming factorization has been performed).
int factor()
Computes the in-place QR factorization of the matrix using the LAPACK routine _GETRF or the Eigen cla...
#define TEUCHOS_MIN(x, y)
bool transpose()
Returns true if adjoint of this matrix has and will be used.
bool formedR()
Returns true if R has been formed explicitly.
Reference-counted pointer class and non-member templated function implementations.
static T one()
Returns representation of one for this scalar type.
RCP< SerialDenseMatrix< OrdinalType, ScalarType > > getR() const
Returns pointer to R (assuming factorization has been performed).
int setVectors(const RCP< SerialDenseMatrix< OrdinalType, ScalarType > > &X, const RCP< SerialDenseMatrix< OrdinalType, ScalarType > > &B)
Sets the pointers for left and right hand side vector(s).
OrdinalType stride() const
Returns the stride between the columns of this matrix in memory.
OrdinalType numRows() const
Returns the row dimension of this matrix.
This class creates and provides basic support for dense rectangular matrix of templated type...