Teuchos - Trilinos Tools Package  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"
21 #include "Teuchos_SerialDenseMatrix.hpp"
22 #include "Teuchos_SerialDenseSolver.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 
106  typedef typename ScalarTraits<ScalarType>::magnitudeType MagnitudeType;
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 
125  SerialQRDenseSolver();
127 
129  virtual ~SerialQRDenseSolver();
131 
133 
134 
136 
138  int setMatrix(const RCP<SerialDenseMatrix<OrdinalType, ScalarType> >& A);
139 
141 
144  int setVectors(const RCP<SerialDenseMatrix<OrdinalType, ScalarType> >& X,
145  const RCP<SerialDenseMatrix<OrdinalType, ScalarType> >& B);
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 
183  int computeEquilibrateScaling();
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 
222  int multiplyQ (ETransp transq, SerialDenseMatrix<OrdinalType, ScalarType> &C);
223 
225 
228  int solveR (ETransp transr, SerialDenseMatrix<OrdinalType, ScalarType> &C);
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 
247  bool shouldEquilibrate() {computeEquilibrateScaling(); return(shouldEquilibrate_);}
248 
250  bool solved() {return(solved_);}
251 
253  bool formedQ() {return(formedQ_);}
254 
256  bool formedR() {return(formedR_);}
257 
259 
261 
262 
264  RCP<SerialDenseMatrix<OrdinalType, ScalarType> > getMatrix() const {return(Matrix_);}
265 
267  RCP<SerialDenseMatrix<OrdinalType, ScalarType> > getFactoredMatrix() const {return(Factor_);}
268 
270  RCP<SerialDenseMatrix<OrdinalType, ScalarType> > getQ() const {return(FactorQ_);}
271 
273  RCP<SerialDenseMatrix<OrdinalType, ScalarType> > getR() const {return(FactorR_);}
274 
276  RCP<SerialDenseMatrix<OrdinalType, ScalarType> > getLHS() const {return(LHS_);}
277 
279  RCP<SerialDenseMatrix<OrdinalType, ScalarType> > getRHS() const {return(RHS_);}
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 
307  bool equilibrate_;
308  bool shouldEquilibrate_;
309  bool equilibratedA_;
310  bool equilibratedB_;
311  OrdinalType equilibrationModeA_;
312  OrdinalType equilibrationModeB_;
313  bool transpose_;
314  bool factored_;
315  bool solved_;
316  bool matrixIsZero_;
317  bool formedQ_;
318  bool formedR_;
319 
320  Teuchos::ETransp TRANS_;
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 
334  MagnitudeType ANORM_;
335  MagnitudeType BNORM_;
336 
337  RCP<SerialDenseMatrix<OrdinalType, ScalarType> > Matrix_;
338  RCP<SerialDenseMatrix<OrdinalType, ScalarType> > LHS_;
339  RCP<SerialDenseMatrix<OrdinalType, ScalarType> > TMP_;
340  RCP<SerialDenseMatrix<OrdinalType, ScalarType> > RHS_;
341  RCP<SerialDenseMatrix<OrdinalType, ScalarType> > Factor_;
342  RCP<SerialDenseMatrix<OrdinalType, ScalarType> > FactorQ_;
343  RCP<SerialDenseMatrix<OrdinalType, ScalarType> > FactorR_;
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 
357  SerialQRDenseSolver(const SerialQRDenseSolver<OrdinalType, ScalarType>& Source);
358  SerialQRDenseSolver & operator=(const SerialQRDenseSolver<OrdinalType, ScalarType>& Source);
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>
368 SerialQRDenseSolver<OrdinalType,ScalarType>::SerialQRDenseSolver()
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>
404 SerialQRDenseSolver<OrdinalType,ScalarType>::~SerialQRDenseSolver()
405 {}
406 
407 //=============================================================================
408 
409 template<typename OrdinalType, typename ScalarType>
410 void SerialQRDenseSolver<OrdinalType,ScalarType>::resetVectors()
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>
421 void SerialQRDenseSolver<OrdinalType,ScalarType>::resetMatrix()
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;
438  ANORM_ = -ScalarTraits<MagnitudeType>::one();
439  BNORM_ = -ScalarTraits<MagnitudeType>::one();
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>
450 int SerialQRDenseSolver<OrdinalType,ScalarType>::setMatrix(const RCP<SerialDenseMatrix<OrdinalType,ScalarType> >& A)
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>
477 int SerialQRDenseSolver<OrdinalType,ScalarType>::setVectors(const RCP<SerialDenseMatrix<OrdinalType,ScalarType> >& X,
478  const RCP<SerialDenseMatrix<OrdinalType,ScalarType> >& B)
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>
501 int SerialQRDenseSolver<OrdinalType,ScalarType>::factor() {
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>
543 int SerialQRDenseSolver<OrdinalType,ScalarType>::solve() {
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>
629 int SerialQRDenseSolver<OrdinalType,ScalarType>::computeEquilibrateScaling()
630 {
631  MagnitudeType safeMin = ScalarTraits<ScalarType>::sfmin();
632  MagnitudeType prec = ScalarTraits<ScalarType>::prec();
633  ScalarType sZero = ScalarTraits<ScalarType>::zero();
634  ScalarType sOne = ScalarTraits<ScalarType>::one();
635  MagnitudeType mZero = ScalarTraits<ScalarType>::magnitude(sZero);
636 
637  MagnitudeType smlnum = ScalarTraits<ScalarType>::magnitude(safeMin/prec);
638  MagnitudeType bignum = ScalarTraits<ScalarType>::magnitude(sOne/smlnum);
639 
640  // Compute maximum absolute value of matrix entries
641  OrdinalType i, j;
642  MagnitudeType anorm = ScalarTraits<ScalarType>::magnitude(ScalarTraits<ScalarType>::zero());
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>
668 int SerialQRDenseSolver<OrdinalType,ScalarType>::equilibrateMatrix()
669 {
670  if (equilibratedA_) return(0);
671 
672  MagnitudeType safeMin = ScalarTraits<ScalarType>::sfmin();
673  MagnitudeType prec = ScalarTraits<ScalarType>::prec();
674  ScalarType sZero = ScalarTraits<ScalarType>::zero();
675  ScalarType sOne = ScalarTraits<ScalarType>::one();
676  MagnitudeType mZero = ScalarTraits<ScalarType>::magnitude(sZero);
677 
678  MagnitudeType smlnum = ScalarTraits<ScalarType>::magnitude(safeMin/prec);
679  MagnitudeType bignum = ScalarTraits<ScalarType>::magnitude(sOne/smlnum);
680  OrdinalType BW = 0;
681 
682  // Compute maximum absolute value of matrix entries
683  OrdinalType i, j;
684  MagnitudeType anorm = ScalarTraits<ScalarType>::magnitude(ScalarTraits<ScalarType>::zero());
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>
714 int SerialQRDenseSolver<OrdinalType,ScalarType>::equilibrateRHS()
715 {
716  if (equilibratedB_) return(0);
717 
718  MagnitudeType safeMin = ScalarTraits<ScalarType>::sfmin();
719  MagnitudeType prec = ScalarTraits<ScalarType>::prec();
720  ScalarType sZero = ScalarTraits<ScalarType>::zero();
721  ScalarType sOne = ScalarTraits<ScalarType>::one();
722  MagnitudeType mZero = ScalarTraits<ScalarType>::magnitude(sZero);
723 
724  MagnitudeType smlnum = ScalarTraits<ScalarType>::magnitude(safeMin/prec);
725  MagnitudeType bignum = ScalarTraits<ScalarType>::magnitude(sOne/smlnum);
726  OrdinalType BW = 0;
727 
728  // Compute maximum absolute value of rhs entries
729  OrdinalType i, j;
730  MagnitudeType bnorm = ScalarTraits<ScalarType>::magnitude(ScalarTraits<ScalarType>::zero());
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>
760 int SerialQRDenseSolver<OrdinalType,ScalarType>::unequilibrateLHS()
761 {
762  if (!equilibratedB_) return(0);
763 
764  MagnitudeType safeMin = ScalarTraits<ScalarType>::sfmin();
765  MagnitudeType prec = ScalarTraits<ScalarType>::prec();
766  ScalarType sZero = ScalarTraits<ScalarType>::zero();
767  ScalarType sOne = ScalarTraits<ScalarType>::one();
768  MagnitudeType mZero = ScalarTraits<ScalarType>::magnitude(sZero);
769  (void) mZero; // Silence "unused variable" compiler warning.
770 
771  MagnitudeType smlnum = ScalarTraits<ScalarType>::magnitude(safeMin/prec);
772  MagnitudeType bignum = ScalarTraits<ScalarType>::magnitude(sOne/smlnum);
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>
797 int SerialQRDenseSolver<OrdinalType,ScalarType>::formQ() {
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>
834 int SerialQRDenseSolver<OrdinalType,ScalarType>::formR() {
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>
861 int SerialQRDenseSolver<OrdinalType, ScalarType>::multiplyQ(ETransp transq, SerialDenseMatrix<OrdinalType, ScalarType> &C)
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
891  Teuchos::ETransp NO_TRANS = Teuchos::NO_TRANS;
892  Teuchos::ETransp TRANS = (Teuchos::ScalarTraits<ScalarType>::isComplex) ? Teuchos::CONJ_TRANS : Teuchos::TRANS;
893  Teuchos::ESide SIDE = Teuchos::LEFT_SIDE;
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>
922 int SerialQRDenseSolver<OrdinalType, ScalarType>::solveR(ETransp transr, SerialDenseMatrix<OrdinalType, ScalarType> &C)
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
954  Teuchos::ETransp NO_TRANS = Teuchos::NO_TRANS;
955  Teuchos::ETransp TRANS = (Teuchos::ScalarTraits<ScalarType>::isComplex) ? Teuchos::CONJ_TRANS : Teuchos::TRANS;
956  Teuchos::EUplo UPLO = Teuchos::UPPER_TRI;
957  Teuchos::EDiag DIAG = Teuchos::NON_UNIT_DIAG;
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>
984 void SerialQRDenseSolver<OrdinalType,ScalarType>::Print(std::ostream& os) const {
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_ */
Teuchos::SerialDenseMatrix::values
ScalarType * values() const
Data array access method.
Definition: Teuchos_SerialDenseMatrix.hpp:228
Teuchos::SerialQRDenseSolver::unequilibrateLHS
int unequilibrateLHS()
Unscales the solution vectors if equilibration was used to solve the system.
Definition: Teuchos_SerialQRDenseSolver.hpp:760
Teuchos::SerialQRDenseSolver::getMatrix
RCP< SerialDenseMatrix< OrdinalType, ScalarType > > getMatrix() const
Returns pointer to current matrix.
Definition: Teuchos_SerialQRDenseSolver.hpp:264
Teuchos::SerialQRDenseSolver::numRows
OrdinalType numRows() const
Returns row dimension of system.
Definition: Teuchos_SerialQRDenseSolver.hpp:282
Teuchos::SerialQRDenseSolver::setMatrix
int setMatrix(const RCP< SerialDenseMatrix< OrdinalType, ScalarType > > &A)
Sets the pointers for coefficient matrix.
Definition: Teuchos_SerialQRDenseSolver.hpp:450
Teuchos_SerialDenseMatrix.hpp
Templated serial dense matrix class.
Teuchos::SerialQRDenseSolver::multiplyQ
int multiplyQ(ETransp transq, SerialDenseMatrix< OrdinalType, ScalarType > &C)
Left multiply the input matrix by the unitary matrix Q or its adjoint.
Definition: Teuchos_SerialQRDenseSolver.hpp:861
Teuchos::SerialQRDenseSolver::solveR
int solveR(ETransp transr, SerialDenseMatrix< OrdinalType, ScalarType > &C)
Solve input matrix on the left with the upper triangular matrix R or its adjoint. ...
Definition: Teuchos_SerialQRDenseSolver.hpp:922
Teuchos::SerialQRDenseSolver::formQ
int formQ()
Explicitly forms the unitary matrix Q.
Definition: Teuchos_SerialQRDenseSolver.hpp:797
Teuchos_BLAS.hpp
Templated interface class to BLAS routines.
Teuchos::SerialQRDenseSolver::formedQ
bool formedQ()
Returns true if Q has been formed explicitly.
Definition: Teuchos_SerialQRDenseSolver.hpp:253
TEUCHOS_TEST_FOR_EXCEPTION
#define TEUCHOS_TEST_FOR_EXCEPTION(throw_exception_test, Exception, msg)
Macro for throwing an exception with breakpointing to ease debugging.
Definition: Teuchos_TestForException.hpp:146
Teuchos::SerialQRDenseSolver::equilibratedA
bool equilibratedA()
Returns true if factor is equilibrated (factor available via getFactoredMatrix()).
Definition: Teuchos_SerialQRDenseSolver.hpp:241
Teuchos::SerialQRDenseSolver::getRHS
RCP< SerialDenseMatrix< OrdinalType, ScalarType > > getRHS() const
Returns pointer to current RHS.
Definition: Teuchos_SerialQRDenseSolver.hpp:279
Teuchos_ConfigDefs.hpp
Teuchos header file which uses auto-configuration information to include necessary C++ headers...
Teuchos::BLAS
Templated BLAS wrapper.
Definition: Teuchos_BLAS.hpp:212
Teuchos_CompObject.hpp
Object for storing data and providing functionality that is common to all computational classes...
Teuchos::SerialQRDenseSolver::factored
bool factored()
Returns true if matrix is factored (factor available via getFactoredMatrix()).
Definition: Teuchos_SerialQRDenseSolver.hpp:238
Teuchos::SerialQRDenseSolver::solve
int solve()
Computes the solution X to AX = B for the this matrix and the B provided to SetVectors()..
Definition: Teuchos_SerialQRDenseSolver.hpp:543
Teuchos::ScalarTraits::sfmin
static magnitudeType sfmin()
Returns safe minimum (sfmin), such that 1/sfmin does not overflow.
Definition: Teuchos_ScalarTraitsDecl.hpp:82
Teuchos::ScalarTraits
This structure defines some basic traits for a scalar field type.
Definition: Teuchos_ScalarTraitsDecl.hpp:58
Teuchos::SerialQRDenseSolver
A class for solving dense linear problems.
Definition: Teuchos_SerialQRDenseSolver.hpp:100
Teuchos_SerialDenseSolver.hpp
Templated class for solving dense linear problems.
Teuchos::CompObject
Functionality and data that is common to all computational classes.
Definition: Teuchos_CompObject.hpp:33
Teuchos_LAPACK.hpp
Templated interface class to LAPACK routines.
Teuchos::SerialQRDenseSolver::solveWithTransposeFlag
void solveWithTransposeFlag(Teuchos::ETransp trans)
All subsequent function calls will work with the transpose-type set by this method (Teuchos::NO_TRANS...
Definition: Teuchos_SerialQRDenseSolver.hpp:160
Teuchos::rcp
TEUCHOS_DEPRECATED RCP< T > rcp(T *p, Dealloc_T dealloc, bool owns_mem)
Deprecated.
Definition: Teuchos_RCPDecl.hpp:1234
Teuchos::ScalarTraits::prec
static magnitudeType prec()
Returns eps*base.
Definition: Teuchos_ScalarTraitsDecl.hpp:86
Teuchos::Object
The base Teuchos class.
Definition: Teuchos_Object.hpp:36
Teuchos::SerialQRDenseSolver::Print
void Print(std::ostream &os) const
Print service methods; defines behavior of ostream &lt;&lt; operator.
Definition: Teuchos_SerialQRDenseSolver.hpp:984
Teuchos::SerialQRDenseSolver::equilibrateRHS
int equilibrateRHS()
Equilibrates the current RHS.
Definition: Teuchos_SerialQRDenseSolver.hpp:714
Teuchos::SerialQRDenseSolver::getFactoredMatrix
RCP< SerialDenseMatrix< OrdinalType, ScalarType > > getFactoredMatrix() const
Returns pointer to factored matrix (assuming factorization has been performed).
Definition: Teuchos_SerialQRDenseSolver.hpp:267
Teuchos::LAPACK
The Templated LAPACK Wrapper Class.
Definition: Teuchos_LAPACK.hpp:64
Teuchos::SerialQRDenseSolver::solveWithTranspose
void solveWithTranspose(bool flag)
If flag is true, causes all subsequent function calls to work with the adjoint of this matrix...
Definition: Teuchos_SerialQRDenseSolver.hpp:157
Teuchos::SerialQRDenseSolver::solved
bool solved()
Returns true if the current set of vectors has been solved.
Definition: Teuchos_SerialQRDenseSolver.hpp:250
Teuchos::SerialQRDenseSolver::SerialQRDenseSolver
SerialQRDenseSolver()
Default constructor; matrix should be set using setMatrix(), LHS and RHS set with setVectors()...
Definition: Teuchos_SerialQRDenseSolver.hpp:368
Teuchos::SerialQRDenseSolver::tau
std::vector< ScalarType > tau() const
Returns pointer to pivot vector (if factorization has been computed), zero otherwise.
Definition: Teuchos_SerialQRDenseSolver.hpp:288
Teuchos::SerialQRDenseSolver::getLHS
RCP< SerialDenseMatrix< OrdinalType, ScalarType > > getLHS() const
Returns pointer to current LHS.
Definition: Teuchos_SerialQRDenseSolver.hpp:276
Teuchos::SerialQRDenseSolver::equilibratedB
bool equilibratedB()
Returns true if RHS is equilibrated (RHS available via getRHS()).
Definition: Teuchos_SerialQRDenseSolver.hpp:244
Teuchos::SerialQRDenseSolver::numCols
OrdinalType numCols() const
Returns column dimension of system.
Definition: Teuchos_SerialQRDenseSolver.hpp:285
Teuchos::ScalarTraits::magnitude
static magnitudeType magnitude(T a)
Returns the magnitudeType of the scalar type a.
Definition: Teuchos_ScalarTraitsDecl.hpp:100
Teuchos::SerialDenseMatrix::numCols
OrdinalType numCols() const
Returns the column dimension of this matrix.
Definition: Teuchos_SerialDenseMatrix.hpp:326
Teuchos::SerialQRDenseSolver::shouldEquilibrate
bool shouldEquilibrate()
Returns true if the LAPACK general rules for equilibration suggest you should equilibrate the system...
Definition: Teuchos_SerialQRDenseSolver.hpp:247
Teuchos::SerialQRDenseSolver::ANORM
MagnitudeType ANORM() const
Returns the absolute value of the largest element of this matrix (returns -1 if not yet computed)...
Definition: Teuchos_SerialQRDenseSolver.hpp:291
Teuchos::SerialQRDenseSolver::factorWithEquilibration
void factorWithEquilibration(bool flag)
Causes equilibration to be called just before the matrix factorization as part of the call to factor...
Definition: Teuchos_SerialQRDenseSolver.hpp:154
Teuchos::SerialQRDenseSolver::equilibrateMatrix
int equilibrateMatrix()
Equilibrates the this matrix.
Definition: Teuchos_SerialQRDenseSolver.hpp:668
Teuchos_ScalarTraits.hpp
Defines basic traits for the scalar field type.
Teuchos::ScalarTraits::zero
static T zero()
Returns representation of zero for this scalar type.
Definition: Teuchos_ScalarTraitsDecl.hpp:102
Teuchos::SerialQRDenseSolver::~SerialQRDenseSolver
virtual ~SerialQRDenseSolver()
SerialQRDenseSolver destructor.
Definition: Teuchos_SerialQRDenseSolver.hpp:404
Teuchos::SerialQRDenseSolver::computeEquilibrateScaling
int computeEquilibrateScaling()
Determines if this matrix should be scaled.
Definition: Teuchos_SerialQRDenseSolver.hpp:629
Teuchos::RCP
Smart reference counting pointer class for automatic garbage collection.
Definition: Teuchos_RCPDecl.hpp:397
Teuchos::SerialQRDenseSolver::formR
int formR()
Explicitly forms the upper triangular matrix R.
Definition: Teuchos_SerialQRDenseSolver.hpp:834
Teuchos::SerialQRDenseSolver::getQ
RCP< SerialDenseMatrix< OrdinalType, ScalarType > > getQ() const
Returns pointer to Q (assuming factorization has been performed).
Definition: Teuchos_SerialQRDenseSolver.hpp:270
Teuchos::SerialQRDenseSolver::factor
int factor()
Computes the in-place QR factorization of the matrix using the LAPACK routine _GETRF or the Eigen cla...
Definition: Teuchos_SerialQRDenseSolver.hpp:501
Teuchos::SerialQRDenseSolver::transpose
bool transpose()
Returns true if adjoint of this matrix has and will be used.
Definition: Teuchos_SerialQRDenseSolver.hpp:235
Teuchos::SerialQRDenseSolver::formedR
bool formedR()
Returns true if R has been formed explicitly.
Definition: Teuchos_SerialQRDenseSolver.hpp:256
Teuchos_RCP.hpp
Reference-counted pointer class and non-member templated function implementations.
Teuchos::ScalarTraits::one
static T one()
Returns representation of one for this scalar type.
Definition: Teuchos_ScalarTraitsDecl.hpp:104
Teuchos::SerialQRDenseSolver::getR
RCP< SerialDenseMatrix< OrdinalType, ScalarType > > getR() const
Returns pointer to R (assuming factorization has been performed).
Definition: Teuchos_SerialQRDenseSolver.hpp:273
Teuchos::SerialQRDenseSolver::setVectors
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).
Definition: Teuchos_SerialQRDenseSolver.hpp:477
Teuchos::SerialDenseMatrix::stride
OrdinalType stride() const
Returns the stride between the columns of this matrix in memory.
Definition: Teuchos_SerialDenseMatrix.hpp:329
Teuchos::SerialDenseMatrix::numRows
OrdinalType numRows() const
Returns the row dimension of this matrix.
Definition: Teuchos_SerialDenseMatrix.hpp:323
Teuchos::SerialDenseMatrix
This class creates and provides basic support for dense rectangular matrix of templated type...
Definition: Teuchos_SerialDenseMatrix.hpp:37
Generated on Thu Nov 21 2024 09:21:49 for Teuchos - Trilinos Tools Package by   doxygen 1.8.5