Teuchos Package Browser (Single Doxygen Collection)  Version of the Day
Teuchos_SerialQRDenseSolver.hpp
Go to the documentation of this file.
1 // @HEADER
2 // ***********************************************************************
3 //
4 // Teuchos: Common Tools Package
5 // Copyright (2004) Sandia Corporation
6 //
7 // Under terms of Contract DE-AC04-94AL85000, there is a non-exclusive
8 // license for use of this work by or on behalf of the U.S. Government.
9 //
10 // Redistribution and use in source and binary forms, with or without
11 // modification, are permitted provided that the following conditions are
12 // met:
13 //
14 // 1. Redistributions of source code must retain the above copyright
15 // notice, this list of conditions and the following disclaimer.
16 //
17 // 2. Redistributions in binary form must reproduce the above copyright
18 // notice, this list of conditions and the following disclaimer in the
19 // documentation and/or other materials provided with the distribution.
20 //
21 // 3. Neither the name of the Corporation nor the names of the
22 // contributors may be used to endorse or promote products derived from
23 // this software without specific prior written permission.
24 //
25 // THIS SOFTWARE IS PROVIDED BY SANDIA CORPORATION "AS IS" AND ANY
26 // EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
27 // IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR
28 // PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL SANDIA CORPORATION OR THE
29 // CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL,
30 // EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO,
31 // PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR
32 // PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF
33 // LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING
34 // NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
35 // SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
36 //
37 // Questions? Contact Michael A. Heroux (maherou@sandia.gov)
38 //
39 // ***********************************************************************
40 // @HEADER
41 
42 #ifndef _TEUCHOS_SERIALQRDENSESOLVER_HPP_
43 #define _TEUCHOS_SERIALQRDENSESOLVER_HPP_
44 
48 #include "Teuchos_CompObject.hpp"
49 #include "Teuchos_BLAS.hpp"
50 #include "Teuchos_LAPACK.hpp"
51 #include "Teuchos_RCP.hpp"
52 #include "Teuchos_ConfigDefs.hpp"
53 #include "Teuchos_SerialDenseMatrix.hpp"
54 #include "Teuchos_SerialDenseSolver.hpp"
55 #include "Teuchos_ScalarTraits.hpp"
56 
57 #ifdef HAVE_TEUCHOSNUMERICS_EIGEN
58 #include "Eigen/Dense"
59 #endif
60 
129 namespace Teuchos {
130 
131  template<typename OrdinalType, typename ScalarType>
132  class SerialQRDenseSolver : public CompObject, public Object, public BLAS<OrdinalType, ScalarType>,
133  public LAPACK<OrdinalType, ScalarType>
134  {
135 
136  public:
137 
138  typedef typename ScalarTraits<ScalarType>::magnitudeType MagnitudeType;
139 #ifdef HAVE_TEUCHOSNUMERICS_EIGEN
140  typedef typename Eigen::Matrix<ScalarType,Eigen::Dynamic,Eigen::Dynamic,Eigen::ColMajor> EigenMatrix;
141  typedef typename Eigen::Matrix<ScalarType,Eigen::Dynamic,1> EigenVector;
142  typedef typename Eigen::InnerStride<Eigen::Dynamic> EigenInnerStride;
143  typedef typename Eigen::OuterStride<Eigen::Dynamic> EigenOuterStride;
144  typedef typename Eigen::Map<EigenVector,0,EigenInnerStride > EigenVectorMap;
145  typedef typename Eigen::Map<const EigenVector,0,EigenInnerStride > EigenConstVectorMap;
146  typedef typename Eigen::Map<EigenMatrix,0,EigenOuterStride > EigenMatrixMap;
147  typedef typename Eigen::Map<const EigenMatrix,0,EigenOuterStride > EigenConstMatrixMap;
148  typedef typename Eigen::PermutationMatrix<Eigen::Dynamic,Eigen::Dynamic,OrdinalType> EigenPermutationMatrix;
149  typedef typename Eigen::Array<OrdinalType,Eigen::Dynamic,1> EigenOrdinalArray;
150  typedef typename Eigen::Map<EigenOrdinalArray> EigenOrdinalArrayMap;
151  typedef typename Eigen::Array<ScalarType,Eigen::Dynamic,1> EigenScalarArray;
152  typedef typename Eigen::Map<EigenScalarArray> EigenScalarArrayMap;
153 #endif
154 
156 
157  SerialQRDenseSolver();
159 
161  virtual ~SerialQRDenseSolver();
163 
165 
166 
168 
170  int setMatrix(const RCP<SerialDenseMatrix<OrdinalType, ScalarType> >& A);
171 
173 
176  int setVectors(const RCP<SerialDenseMatrix<OrdinalType, ScalarType> >& X,
177  const RCP<SerialDenseMatrix<OrdinalType, ScalarType> >& B);
179 
181 
182 
184 
186  void factorWithEquilibration(bool flag) {equilibrate_ = flag; return;}
187 
189  void solveWithTranspose(bool flag) {transpose_ = flag; if (flag) TRANS_ = Teuchos::CONJ_TRANS; else TRANS_ = Teuchos::NO_TRANS; return;}
190 
192  void solveWithTransposeFlag(Teuchos::ETransp trans) {TRANS_ = trans; transpose_ = (trans != Teuchos::NO_TRANS) ? true : false; }
193 
195 
197 
198 
200 
203  int factor();
204 
206 
209  int solve();
210 
212 
215  int computeEquilibrateScaling();
216 
218 
222  int equilibrateMatrix();
223 
225 
229  int equilibrateRHS();
230 
232 
236  int unequilibrateLHS();
237 
239 
242  int formQ();
243 
245 
248  int formR();
249 
251 
254  int multiplyQ (ETransp transq, SerialDenseMatrix<OrdinalType, ScalarType> &C);
255 
257 
260  int solveR (ETransp transr, SerialDenseMatrix<OrdinalType, ScalarType> &C);
262 
264 
265 
267  bool transpose() {return(transpose_);}
268 
270  bool factored() {return(factored_);}
271 
273  bool equilibratedA() {return(equilibratedA_);}
274 
276  bool equilibratedB() {return(equilibratedB_);}
277 
279  bool shouldEquilibrate() {computeEquilibrateScaling(); return(shouldEquilibrate_);}
280 
282  bool solved() {return(solved_);}
283 
285  bool formedQ() {return(formedQ_);}
286 
288  bool formedR() {return(formedR_);}
289 
291 
293 
294 
296  RCP<SerialDenseMatrix<OrdinalType, ScalarType> > getMatrix() const {return(Matrix_);}
297 
299  RCP<SerialDenseMatrix<OrdinalType, ScalarType> > getFactoredMatrix() const {return(Factor_);}
300 
302  RCP<SerialDenseMatrix<OrdinalType, ScalarType> > getQ() const {return(FactorQ_);}
303 
305  RCP<SerialDenseMatrix<OrdinalType, ScalarType> > getR() const {return(FactorR_);}
306 
308  RCP<SerialDenseMatrix<OrdinalType, ScalarType> > getLHS() const {return(LHS_);}
309 
311  RCP<SerialDenseMatrix<OrdinalType, ScalarType> > getRHS() const {return(RHS_);}
312 
314  OrdinalType numRows() const {return(M_);}
315 
317  OrdinalType numCols() const {return(N_);}
318 
320  std::vector<ScalarType> tau() const {return(TAU_);}
321 
323  MagnitudeType ANORM() const {return(ANORM_);}
324 
326 
328 
329  void Print(std::ostream& os) const;
332  protected:
333 
334  void allocateWORK() { LWORK_ = 4*N_; WORK_.resize( LWORK_ ); return;}
335  void resetMatrix();
336  void resetVectors();
337 
338 
339  bool equilibrate_;
340  bool shouldEquilibrate_;
341  bool equilibratedA_;
342  bool equilibratedB_;
343  OrdinalType equilibrationModeA_;
344  OrdinalType equilibrationModeB_;
345  bool transpose_;
346  bool factored_;
347  bool solved_;
348  bool matrixIsZero_;
349  bool formedQ_;
350  bool formedR_;
351 
352  Teuchos::ETransp TRANS_;
353 
354  OrdinalType M_;
355  OrdinalType N_;
356  OrdinalType Min_MN_;
357  OrdinalType LDA_;
358  OrdinalType LDAF_;
359  OrdinalType LDQ_;
360  OrdinalType LDR_;
361  OrdinalType INFO_;
362  OrdinalType LWORK_;
363 
364  std::vector<ScalarType> TAU_;
365 
366  MagnitudeType ANORM_;
367  MagnitudeType BNORM_;
368 
369  RCP<SerialDenseMatrix<OrdinalType, ScalarType> > Matrix_;
370  RCP<SerialDenseMatrix<OrdinalType, ScalarType> > LHS_;
371  RCP<SerialDenseMatrix<OrdinalType, ScalarType> > TMP_;
372  RCP<SerialDenseMatrix<OrdinalType, ScalarType> > RHS_;
373  RCP<SerialDenseMatrix<OrdinalType, ScalarType> > Factor_;
374  RCP<SerialDenseMatrix<OrdinalType, ScalarType> > FactorQ_;
375  RCP<SerialDenseMatrix<OrdinalType, ScalarType> > FactorR_;
376 
377  ScalarType * A_;
378  ScalarType * AF_;
379  ScalarType * Q_;
380  ScalarType * R_;
381  std::vector<ScalarType> WORK_;
382 #ifdef HAVE_TEUCHOSNUMERICS_EIGEN
383  Eigen::HouseholderQR<EigenMatrix> qr_;
384 #endif
385 
386  private:
387  // SerialQRDenseSolver copy constructor (put here because we don't want user access)
388 
389  SerialQRDenseSolver(const SerialQRDenseSolver<OrdinalType, ScalarType>& Source);
390  SerialQRDenseSolver & operator=(const SerialQRDenseSolver<OrdinalType, ScalarType>& Source);
391 
392  };
393 
394  // Helper traits to distinguish work arrays for real and complex-valued datatypes.
395  using namespace details;
396 
397 //=============================================================================
398 
399 template<typename OrdinalType, typename ScalarType>
400 SerialQRDenseSolver<OrdinalType,ScalarType>::SerialQRDenseSolver()
401  : CompObject(),
402  equilibrate_(false),
403  shouldEquilibrate_(false),
404  equilibratedA_(false),
405  equilibratedB_(false),
406  equilibrationModeA_(0),
407  equilibrationModeB_(0),
408  transpose_(false),
409  factored_(false),
410  solved_(false),
411  matrixIsZero_(false),
412  formedQ_(false),
413  formedR_(false),
414  TRANS_(Teuchos::NO_TRANS),
415  M_(0),
416  N_(0),
417  Min_MN_(0),
418  LDA_(0),
419  LDAF_(0),
420  LDQ_(0),
421  LDR_(0),
422  INFO_(0),
423  LWORK_(0),
424  ANORM_(ScalarTraits<MagnitudeType>::zero()),
425  BNORM_(ScalarTraits<MagnitudeType>::zero()),
426  A_(0),
427  AF_(0),
428  Q_(0),
429  R_(0)
430 {
431  resetMatrix();
432 }
433 //=============================================================================
434 
435 template<typename OrdinalType, typename ScalarType>
436 SerialQRDenseSolver<OrdinalType,ScalarType>::~SerialQRDenseSolver()
437 {}
438 
439 //=============================================================================
440 
441 template<typename OrdinalType, typename ScalarType>
442 void SerialQRDenseSolver<OrdinalType,ScalarType>::resetVectors()
443 {
444  LHS_ = Teuchos::null;
445  TMP_ = Teuchos::null;
446  RHS_ = Teuchos::null;
447  solved_ = false;
448  equilibratedB_ = false;
449 }
450 //=============================================================================
451 
452 template<typename OrdinalType, typename ScalarType>
453 void SerialQRDenseSolver<OrdinalType,ScalarType>::resetMatrix()
454 {
455  resetVectors();
456  equilibratedA_ = false;
457  equilibrationModeA_ = 0;
458  equilibrationModeB_ = 0;
459  factored_ = false;
460  matrixIsZero_ = false;
461  formedQ_ = false;
462  formedR_ = false;
463  M_ = 0;
464  N_ = 0;
465  Min_MN_ = 0;
466  LDA_ = 0;
467  LDAF_ = 0;
468  LDQ_ = 0;
469  LDR_ = 0;
470  ANORM_ = -ScalarTraits<MagnitudeType>::one();
471  BNORM_ = -ScalarTraits<MagnitudeType>::one();
472  A_ = 0;
473  AF_ = 0;
474  Q_ = 0;
475  R_ = 0;
476  INFO_ = 0;
477  LWORK_ = 0;
478 }
479 //=============================================================================
480 
481 template<typename OrdinalType, typename ScalarType>
482 int SerialQRDenseSolver<OrdinalType,ScalarType>::setMatrix(const RCP<SerialDenseMatrix<OrdinalType,ScalarType> >& A)
483 {
484  TEUCHOS_TEST_FOR_EXCEPTION(A->numRows() < A->numCols(), std::invalid_argument,
485  "SerialQRDenseSolver<T>::setMatrix: the matrix A must have A.numRows() >= A.numCols()!");
486 
487  resetMatrix();
488  Matrix_ = A;
489  Factor_ = A;
490  FactorQ_ = A;
491  FactorR_ = A;
492  M_ = A->numRows();
493  N_ = A->numCols();
494  Min_MN_ = TEUCHOS_MIN(M_,N_);
495  LDA_ = A->stride();
496  LDAF_ = LDA_;
497  LDQ_ = LDA_;
498  LDR_ = N_;
499  A_ = A->values();
500  AF_ = A->values();
501  Q_ = A->values();
502  R_ = A->values();
503 
504  return(0);
505 }
506 //=============================================================================
507 
508 template<typename OrdinalType, typename ScalarType>
509 int SerialQRDenseSolver<OrdinalType,ScalarType>::setVectors(const RCP<SerialDenseMatrix<OrdinalType,ScalarType> >& X,
510  const RCP<SerialDenseMatrix<OrdinalType,ScalarType> >& B)
511 {
512  // Check that these new vectors are consistent
513  TEUCHOS_TEST_FOR_EXCEPTION(B->numCols() != X->numCols(), std::invalid_argument,
514  "SerialQRDenseSolver<T>::setVectors: X and B have different numbers of columns!");
515  TEUCHOS_TEST_FOR_EXCEPTION(B->values()==0, std::invalid_argument,
516  "SerialQRDenseSolver<T>::setVectors: B is an empty SerialDenseMatrix<T>!");
517  TEUCHOS_TEST_FOR_EXCEPTION(X->values()==0, std::invalid_argument,
518  "SerialQRDenseSolver<T>::setVectors: X is an empty SerialDenseMatrix<T>!");
519  TEUCHOS_TEST_FOR_EXCEPTION(B->stride()<1, std::invalid_argument,
520  "SerialQRDenseSolver<T>::setVectors: B has an invalid stride!");
521  TEUCHOS_TEST_FOR_EXCEPTION(X->stride()<1, std::invalid_argument,
522  "SerialQRDenseSolver<T>::setVectors: X has an invalid stride!");
523 
524  resetVectors();
525  LHS_ = X;
526  RHS_ = B;
527 
528  return(0);
529 }
530 //=============================================================================
531 
532 template<typename OrdinalType, typename ScalarType>
533 int SerialQRDenseSolver<OrdinalType,ScalarType>::factor() {
534 
535  if (factored()) return(0);
536 
537  // Equilibrate matrix if necessary
538  int ierr = 0;
539  if (equilibrate_) ierr = equilibrateMatrix();
540  if (ierr!=0) return(ierr);
541 
542  allocateWORK();
543  if ((int)TAU_.size()!=Min_MN_) TAU_.resize( Min_MN_ );
544 
545  INFO_ = 0;
546 
547  // Factor
548 #ifdef HAVE_TEUCHOSNUMERICS_EIGEN
549  EigenMatrixMap aMap( AF_, M_, N_, EigenOuterStride(LDAF_) );
550  EigenScalarArray tau;
551  EigenScalarArrayMap tauMap(&TAU_[0],TEUCHOS_MIN(M_,N_));
552  qr_.compute(aMap);
553  tau = qr_.hCoeffs();
554  for (OrdinalType i=0; i<tauMap.innerSize(); i++) {
555  tauMap(i) = tau(i);
556  }
557  EigenMatrix qrMat = qr_.matrixQR();
558  for (OrdinalType j=0; j<aMap.outerSize(); j++) {
559  for (OrdinalType i=0; i<aMap.innerSize(); i++) {
560  aMap(i,j) = qrMat(i,j);
561  }
562  }
563 #else
564  this->GEQRF(M_, N_, AF_, LDAF_, &TAU_[0], &WORK_[0], LWORK_, &INFO_);
565 #endif
566 
567  factored_ = true;
568 
569  return(INFO_);
570 }
571 
572 //=============================================================================
573 
574 template<typename OrdinalType, typename ScalarType>
575 int SerialQRDenseSolver<OrdinalType,ScalarType>::solve() {
576 
577  // Check if the matrix is zero
578  if (matrixIsZero_) {
579  LHS_->putScalar(ScalarTraits<ScalarType>::zero());
580  return(0);
581  }
582 
583  // Equilibrate RHS if necessary
584  int ierr = 0;
585  if (equilibrate_) {
586  ierr = equilibrateRHS();
587  }
588  if (ierr != 0) return(ierr);
589 
590  TEUCHOS_TEST_FOR_EXCEPTION( RHS_==Teuchos::null, std::invalid_argument,
591  "SerialQRDenseSolver<T>::solve: No right-hand side vector (RHS) has been set for the linear system!");
592  TEUCHOS_TEST_FOR_EXCEPTION( LHS_==Teuchos::null, std::invalid_argument,
593  "SerialQRDenseSolver<T>::solve: No solution vector (LHS) has been set for the linear system!");
594  if ( TRANS_ == Teuchos::NO_TRANS ) {
595  TEUCHOS_TEST_FOR_EXCEPTION( LHS_->numRows() != N_, std::invalid_argument,
596  "SerialQRDenseSolver<T>::solve: No transpose: must have LHS_->numRows() = N_");
597  } else {
598  TEUCHOS_TEST_FOR_EXCEPTION( LHS_->numRows() != M_, std::invalid_argument,
599  "SerialQRDenseSolver<T>::solve: Transpose: must have LHS_->numRows() = M_");
600  }
601 
602  if (shouldEquilibrate() && !equilibratedA_)
603  std::cout << "WARNING! SerialQRDenseSolver<T>::solve: System should be equilibrated!" << std::endl;
604 
605  // Matrix must be factored
606  if (!factored()) factor();
607 
608  TEUCHOS_TEST_FOR_EXCEPTION( (equilibratedA_ && !equilibratedB_) || (!equilibratedA_ && equilibratedB_) ,
609  std::logic_error, "SerialQRDenseSolver<T>::solve: Matrix and vectors must be similarly scaled!");
610 
611  TMP_ = rcp( new SerialDenseMatrix<OrdinalType,ScalarType>(M_, RHS_->numCols()) );
612  for (OrdinalType j=0; j<RHS_->numCols(); j++) {
613  for (OrdinalType i=0; i<RHS_->numRows(); i++) {
614  (*TMP_)(i,j) = (*RHS_)(i,j);
615  }
616  }
617 
618  INFO_ = 0;
619 
620  // Solve
621  if ( TRANS_ == Teuchos::NO_TRANS ) {
622 
623  // Solve A lhs = rhs
624  this->multiplyQ( Teuchos::CONJ_TRANS, *TMP_ );
625  this->solveR( Teuchos::NO_TRANS, *TMP_ );
626 
627  } else {
628 
629  // Solve A**H lhs = rhs
630  this->solveR( Teuchos::CONJ_TRANS, *TMP_ );
631  for (OrdinalType j = 0; j < RHS_->numCols(); j++ ) {
632  for (OrdinalType i = N_; i < M_; i++ ) {
633  (*TMP_)(i, j) = ScalarTraits<ScalarType>::zero();
634  }
635  }
636  this->multiplyQ( Teuchos::NO_TRANS, *TMP_ );
637 
638  }
639  for (OrdinalType j = 0; j < LHS_->numCols(); j++ ) {
640  for (OrdinalType i = 0; i < LHS_->numRows(); i++ ) {
641  (*LHS_)(i, j) = (*TMP_)(i,j);
642  }
643  }
644 
645  solved_ = true;
646 
647  // Unequilibrate LHS if necessary
648  if (equilibrate_) {
649  ierr = unequilibrateLHS();
650  }
651  if (ierr != 0) {
652  return ierr;
653  }
654 
655  return INFO_;
656 }
657 
658 //=============================================================================
659 
660 template<typename OrdinalType, typename ScalarType>
661 int SerialQRDenseSolver<OrdinalType,ScalarType>::computeEquilibrateScaling()
662 {
663  MagnitudeType safeMin = ScalarTraits<ScalarType>::sfmin();
664  MagnitudeType prec = ScalarTraits<ScalarType>::prec();
665  ScalarType sZero = ScalarTraits<ScalarType>::zero();
666  ScalarType sOne = ScalarTraits<ScalarType>::one();
667  MagnitudeType mZero = ScalarTraits<ScalarType>::magnitude(sZero);
668 
669  MagnitudeType smlnum = ScalarTraits<ScalarType>::magnitude(safeMin/prec);
670  MagnitudeType bignum = ScalarTraits<ScalarType>::magnitude(sOne/smlnum);
671 
672  // Compute maximum absolute value of matrix entries
673  OrdinalType i, j;
674  MagnitudeType anorm = ScalarTraits<ScalarType>::magnitude(ScalarTraits<ScalarType>::zero());
675  for (j = 0; j < N_; j++) {
676  for (i = 0; i < M_; i++) {
677  anorm = TEUCHOS_MAX( anorm, ScalarTraits<ScalarType>::magnitude((*Matrix_)(i,j)) );
678  }
679  }
680 
681  ANORM_ = anorm;
682 
683  if (ANORM_ > mZero && ANORM_ < smlnum) {
684  // Scale matrix norm up to smlnum
685  shouldEquilibrate_ = true;
686  } else if (ANORM_ > bignum) {
687  // Scale matrix norm down to bignum
688  shouldEquilibrate_ = true;
689  } else if (ANORM_ == mZero) {
690  // Matrix all zero. Return zero solution
691  matrixIsZero_ = true;
692  }
693 
694  return(0);
695 }
696 
697 //=============================================================================
698 
699 template<typename OrdinalType, typename ScalarType>
700 int SerialQRDenseSolver<OrdinalType,ScalarType>::equilibrateMatrix()
701 {
702  if (equilibratedA_) return(0);
703 
704  MagnitudeType safeMin = ScalarTraits<ScalarType>::sfmin();
705  MagnitudeType prec = ScalarTraits<ScalarType>::prec();
706  ScalarType sZero = ScalarTraits<ScalarType>::zero();
707  ScalarType sOne = ScalarTraits<ScalarType>::one();
708  MagnitudeType mZero = ScalarTraits<ScalarType>::magnitude(sZero);
709 
710  MagnitudeType smlnum = ScalarTraits<ScalarType>::magnitude(safeMin/prec);
711  MagnitudeType bignum = ScalarTraits<ScalarType>::magnitude(sOne/smlnum);
712  OrdinalType BW = 0;
713 
714  // Compute maximum absolute value of matrix entries
715  OrdinalType i, j;
716  MagnitudeType anorm = ScalarTraits<ScalarType>::magnitude(ScalarTraits<ScalarType>::zero());
717  for (j = 0; j < N_; j++) {
718  for (i = 0; i < M_; i++) {
719  anorm = TEUCHOS_MAX( anorm, ScalarTraits<ScalarType>::magnitude((*Matrix_)(i,j)) );
720  }
721  }
722 
723  ANORM_ = anorm;
724  int ierr1 = 0;
725  if (ANORM_ > mZero && ANORM_ < smlnum) {
726  // Scale matrix norm up to smlnum
727  this->LASCL(Teuchos::ETypeChar[Teuchos::FULL], BW, BW, ANORM_, smlnum, M_, N_, A_, LDA_, &INFO_);
728  equilibrationModeA_ = 1;
729  } else if (ANORM_ > bignum) {
730  // Scale matrix norm down to bignum
731  this->LASCL(Teuchos::ETypeChar[Teuchos::FULL], BW, BW, ANORM_, bignum, M_, N_, A_, LDA_, &INFO_);
732  equilibrationModeA_ = 2;
733  } else if (ANORM_ == mZero) {
734  // Matrix all zero. Return zero solution
735  matrixIsZero_ = true;
736  }
737 
738  equilibratedA_ = true;
739 
740  return(ierr1);
741 }
742 
743 //=============================================================================
744 
745 template<typename OrdinalType, typename ScalarType>
746 int SerialQRDenseSolver<OrdinalType,ScalarType>::equilibrateRHS()
747 {
748  if (equilibratedB_) return(0);
749 
750  MagnitudeType safeMin = ScalarTraits<ScalarType>::sfmin();
751  MagnitudeType prec = ScalarTraits<ScalarType>::prec();
752  ScalarType sZero = ScalarTraits<ScalarType>::zero();
753  ScalarType sOne = ScalarTraits<ScalarType>::one();
754  MagnitudeType mZero = ScalarTraits<ScalarType>::magnitude(sZero);
755 
756  MagnitudeType smlnum = ScalarTraits<ScalarType>::magnitude(safeMin/prec);
757  MagnitudeType bignum = ScalarTraits<ScalarType>::magnitude(sOne/smlnum);
758  OrdinalType BW = 0;
759 
760  // Compute maximum absolute value of rhs entries
761  OrdinalType i, j;
762  MagnitudeType bnorm = ScalarTraits<ScalarType>::magnitude(ScalarTraits<ScalarType>::zero());
763  for (j = 0; j <RHS_->numCols(); j++) {
764  for (i = 0; i < RHS_->numRows(); i++) {
765  bnorm = TEUCHOS_MAX( bnorm, ScalarTraits<ScalarType>::magnitude((*RHS_)(i,j)) );
766  }
767  }
768 
769  BNORM_ = bnorm;
770 
771  int ierr1 = 0;
772  if (BNORM_ > mZero && BNORM_ < smlnum) {
773  // Scale RHS norm up to smlnum
774  this->LASCL(Teuchos::ETypeChar[Teuchos::FULL], BW, BW, BNORM_, smlnum, RHS_->numRows(), RHS_->numCols(),
775  RHS_->values(), RHS_->stride(), &INFO_);
776  equilibrationModeB_ = 1;
777  } else if (BNORM_ > bignum) {
778  // Scale RHS norm down to bignum
779  this->LASCL(Teuchos::ETypeChar[Teuchos::FULL], BW, BW, BNORM_, bignum, RHS_->numRows(), RHS_->numCols(),
780  RHS_->values(), RHS_->stride(), &INFO_);
781  equilibrationModeB_ = 2;
782  }
783 
784  equilibratedB_ = true;
785 
786  return(ierr1);
787 }
788 
789 //=============================================================================
790 
791 template<typename OrdinalType, typename ScalarType>
792 int SerialQRDenseSolver<OrdinalType,ScalarType>::unequilibrateLHS()
793 {
794  if (!equilibratedB_) return(0);
795 
796  MagnitudeType safeMin = ScalarTraits<ScalarType>::sfmin();
797  MagnitudeType prec = ScalarTraits<ScalarType>::prec();
798  ScalarType sZero = ScalarTraits<ScalarType>::zero();
799  ScalarType sOne = ScalarTraits<ScalarType>::one();
800  MagnitudeType mZero = ScalarTraits<ScalarType>::magnitude(sZero);
801  (void) mZero; // Silence "unused variable" compiler warning.
802 
803  MagnitudeType smlnum = ScalarTraits<ScalarType>::magnitude(safeMin/prec);
804  MagnitudeType bignum = ScalarTraits<ScalarType>::magnitude(sOne/smlnum);
805  OrdinalType BW = 0;
806 
807  int ierr1 = 0;
808  if (equilibrationModeA_ == 1) {
809  this->LASCL(Teuchos::ETypeChar[Teuchos::FULL], BW, BW, ANORM_, smlnum, LHS_->numRows(), LHS_->numCols(),
810  LHS_->values(), LHS_->stride(), &INFO_);
811  } else if (equilibrationModeA_ == 2) {
812  this->LASCL(Teuchos::ETypeChar[Teuchos::FULL], BW, BW, ANORM_, bignum, LHS_->numRows(), LHS_->numCols(),
813  LHS_->values(), LHS_->stride(), &INFO_);
814  }
815  if (equilibrationModeB_ == 1) {
816  this->LASCL(Teuchos::ETypeChar[Teuchos::FULL], BW, BW, smlnum, BNORM_, LHS_->numRows(), LHS_->numCols(),
817  LHS_->values(), LHS_->stride(), &INFO_);
818  } else if (equilibrationModeB_ == 2) {
819  this->LASCL(Teuchos::ETypeChar[Teuchos::FULL], BW, BW, bignum, BNORM_, LHS_->numRows(), LHS_->numCols(),
820  LHS_->values(), LHS_->stride(), &INFO_);
821  }
822 
823  return(ierr1);
824 }
825 
826 //=============================================================================
827 
828 template<typename OrdinalType, typename ScalarType>
829 int SerialQRDenseSolver<OrdinalType,ScalarType>::formQ() {
830 
831  // Matrix must be factored first
832  if (!factored()) factor();
833 
834  // Store Q separately from factored matrix
835  if (AF_ == Q_) {
836  FactorQ_ = rcp( new SerialDenseMatrix<OrdinalType,ScalarType>(*Factor_) );
837  Q_ = FactorQ_->values();
838  LDQ_ = FactorQ_->stride();
839  }
840 
841  INFO_ = 0;
842 
843  // Form Q
844 #ifdef HAVE_TEUCHOSNUMERICS_EIGEN
845  EigenMatrixMap qMap( Q_, M_, N_, EigenOuterStride(LDQ_) );
846  EigenMatrix qMat = qr_.householderQ();
847  for (OrdinalType j=0; j<qMap.outerSize(); j++) {
848  for (OrdinalType i=0; i<qMap.innerSize(); i++) {
849  qMap(i,j) = qMat(i,j);
850  }
851  }
852 #else
853  this->UNGQR(M_, N_, N_, Q_, LDQ_, &TAU_[0], &WORK_[0], LWORK_, &INFO_);
854 #endif
855 
856  if (INFO_!=0) return(INFO_);
857 
858  formedQ_ = true;
859 
860  return(INFO_);
861 }
862 
863 //=============================================================================
864 
865 template<typename OrdinalType, typename ScalarType>
866 int SerialQRDenseSolver<OrdinalType,ScalarType>::formR() {
867 
868  // Matrix must be factored first
869  if (!factored()) factor();
870 
871  // Store R separately from factored matrix
872  if (AF_ == R_) {
873  FactorR_ = rcp( new SerialDenseMatrix<OrdinalType,ScalarType>(N_, N_) );
874  R_ = FactorR_->values();
875  LDR_ = FactorR_->stride();
876  }
877 
878  // Form R
879  for (OrdinalType j=0; j<N_; j++) {
880  for (OrdinalType i=0; i<=j; i++) {
881  (*FactorR_)(i,j) = (*Factor_)(i,j);
882  }
883  }
884 
885  formedR_ = true;
886 
887  return(0);
888 }
889 
890 //=============================================================================
891 
892 template<typename OrdinalType, typename ScalarType>
893 int SerialQRDenseSolver<OrdinalType, ScalarType>::multiplyQ(ETransp transq, SerialDenseMatrix<OrdinalType, ScalarType> &C)
894 {
895 
896  // Check that the matrices are consistent.
897  TEUCHOS_TEST_FOR_EXCEPTION(C.numRows()!=M_, std::invalid_argument,
898  "SerialQRDenseSolver<T>::multiplyQ: C.numRows() != M_!");
899  TEUCHOS_TEST_FOR_EXCEPTION(C.values()==0, std::invalid_argument,
900  "SerialQRDenseSolver<T>::multiplyQ: C is an empty SerialDenseMatrix<T>!");
901 
902  // Matrix must be factored
903  if (!factored()) factor();
904 
905  INFO_ = 0;
906 
907  // Multiply
908 #ifdef HAVE_TEUCHOSNUMERICS_EIGEN
909  EigenMatrixMap cMap( C.values(), C.numRows(), C.numCols(), EigenOuterStride(C.stride()) );
910  if ( transq == Teuchos::NO_TRANS ) {
911  // C = Q * C
912  cMap = qr_.householderQ() * cMap;
913  } else {
914  // C = Q**H * C
915  cMap = qr_.householderQ().adjoint() * cMap;
916  for (OrdinalType j = 0; j < C.numCols(); j++) {
917  for (OrdinalType i = N_; i < C.numRows(); i++ ) {
918  cMap(i, j) = ScalarTraits<ScalarType>::zero();
919  }
920  }
921  }
922 #else
923  Teuchos::ETransp NO_TRANS = Teuchos::NO_TRANS;
924  Teuchos::ETransp TRANS = (Teuchos::ScalarTraits<ScalarType>::isComplex) ? Teuchos::CONJ_TRANS : Teuchos::TRANS;
925  Teuchos::ESide SIDE = Teuchos::LEFT_SIDE;
926 
927  if ( transq == Teuchos::NO_TRANS ) {
928 
929  // C = Q * C
930  this->UNMQR(ESideChar[SIDE], ETranspChar[NO_TRANS], M_, C.numCols(), N_, AF_, LDAF_,
931  &TAU_[0], C.values(), C.stride(), &WORK_[0], LWORK_, &INFO_);
932 
933  } else {
934 
935  // C = Q**H * C
936  this->UNMQR(ESideChar[SIDE], ETranspChar[TRANS], M_, C.numCols(), N_, AF_, LDAF_,
937  &TAU_[0], C.values(), C.stride(), &WORK_[0], LWORK_, &INFO_);
938 
939  for (OrdinalType j = 0; j < C.numCols(); j++) {
940  for (OrdinalType i = N_; i < C.numRows(); i++ ) {
941  C(i, j) = ScalarTraits<ScalarType>::zero();
942  }
943  }
944  }
945 #endif
946 
947  return(INFO_);
948 
949 }
950 
951 //=============================================================================
952 
953 template<typename OrdinalType, typename ScalarType>
954 int SerialQRDenseSolver<OrdinalType, ScalarType>::solveR(ETransp transr, SerialDenseMatrix<OrdinalType, ScalarType> &C)
955 {
956 
957  // Check that the matrices are consistent.
958  TEUCHOS_TEST_FOR_EXCEPTION(C.numRows()<N_ || C.numRows()>M_, std::invalid_argument,
959  "SerialQRDenseSolver<T>::solveR: must have N_ < C.numRows() < M_!");
960  TEUCHOS_TEST_FOR_EXCEPTION(C.values()==0, std::invalid_argument,
961  "SerialQRDenseSolver<T>::solveR: C is an empty SerialDenseMatrix<T>!");
962 
963  // Matrix must be factored
964  if (!factored()) factor();
965 
966  INFO_ = 0;
967 
968  // Solve
969 #ifdef HAVE_TEUCHOSNUMERICS_EIGEN
970  EigenMatrixMap cMap( C.values(), N_, C.numCols(), EigenOuterStride(C.stride()) );
971  // Check for singularity first like TRTRS
972  for (OrdinalType j=0; j<N_; j++) {
973  if ((qr_.matrixQR())(j,j) == ScalarTraits<ScalarType>::zero()) {
974  INFO_ = j+1;
975  return(INFO_);
976  }
977  }
978  if ( transr == Teuchos::NO_TRANS ) {
979  // C = R \ C
980  qr_.matrixQR().topRows(N_).template triangularView<Eigen::Upper>().solveInPlace(cMap);
981  } else {
982  // C = R**H \ C
983  qr_.matrixQR().topRows(N_).template triangularView<Eigen::Upper>().adjoint().solveInPlace(cMap);
984  }
985 #else
986  Teuchos::ETransp NO_TRANS = Teuchos::NO_TRANS;
987  Teuchos::ETransp TRANS = (Teuchos::ScalarTraits<ScalarType>::isComplex) ? Teuchos::CONJ_TRANS : Teuchos::TRANS;
988  Teuchos::EUplo UPLO = Teuchos::UPPER_TRI;
989  Teuchos::EDiag DIAG = Teuchos::NON_UNIT_DIAG;
990 
991  ScalarType* RMAT = (formedR_) ? R_ : AF_;
992  OrdinalType LDRMAT = (formedR_) ? LDR_ : LDAF_;
993 
994  if ( transr == Teuchos::NO_TRANS ) {
995 
996  // C = R \ C
997  this->TRTRS(EUploChar[UPLO], ETranspChar[NO_TRANS], EDiagChar[DIAG], N_, C.numCols(),
998  RMAT, LDRMAT, C.values(), C.stride(), &INFO_);
999 
1000  } else {
1001 
1002  // C = R**H \ C
1003  this->TRTRS(EUploChar[UPLO], ETranspChar[TRANS], EDiagChar[DIAG], N_, C.numCols(),
1004  RMAT, LDRMAT, C.values(), C.stride(), &INFO_);
1005 
1006  }
1007 #endif
1008 
1009  return(INFO_);
1010 
1011 }
1012 
1013 //=============================================================================
1014 
1015 template<typename OrdinalType, typename ScalarType>
1016 void SerialQRDenseSolver<OrdinalType,ScalarType>::Print(std::ostream& os) const {
1017 
1018  if (Matrix_!=Teuchos::null) os << "Solver Matrix" << std::endl << *Matrix_ << std::endl;
1019  if (Factor_!=Teuchos::null) os << "Solver Factored Matrix" << std::endl << *Factor_ << std::endl;
1020  if (Q_!=Teuchos::null) os << "Solver Factor Q" << std::endl << *Q_ << std::endl;
1021  if (LHS_ !=Teuchos::null) os << "Solver LHS" << std::endl << *LHS_ << std::endl;
1022  if (RHS_ !=Teuchos::null) os << "Solver RHS" << std::endl << *RHS_ << std::endl;
1023 
1024 }
1025 
1026 } // namespace Teuchos
1027 
1028 #endif /* _TEUCHOS_SERIALQRDENSESOLVER_HPP_ */
C
Definition: PackageC.cpp:3
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:533
Teuchos::SerialQRDenseSolver::solved
bool solved()
Returns true if the current set of vectors has been solved.
Definition: Teuchos_SerialQRDenseSolver.hpp:282
Teuchos::SerialQRDenseSolver::Factor_
RCP< SerialDenseMatrix< OrdinalType, ScalarType > > Factor_
Definition: Teuchos_SerialQRDenseSolver.hpp:373
Teuchos_RCP.hpp
Reference-counted pointer class and non-member templated function implementations.
Teuchos::SerialQRDenseSolver::BNORM_
MagnitudeType BNORM_
Definition: Teuchos_SerialQRDenseSolver.hpp:367
Teuchos::SerialQRDenseSolver::Q_
ScalarType * Q_
Definition: Teuchos_SerialQRDenseSolver.hpp:379
Teuchos::ScalarTraits::sfmin
static magnitudeType sfmin()
Returns safe minimum (sfmin), such that 1/sfmin does not overflow.
Definition: Teuchos_ScalarTraitsDecl.hpp:112
Teuchos::SerialQRDenseSolver::transpose
bool transpose()
Returns true if adjoint of this matrix has and will be used.
Definition: Teuchos_SerialQRDenseSolver.hpp:267
Teuchos::SerialQRDenseSolver::equilibratedB
bool equilibratedB()
Returns true if RHS is equilibrated (RHS available via getRHS()).
Definition: Teuchos_SerialQRDenseSolver.hpp:276
Teuchos::SerialQRDenseSolver::SerialQRDenseSolver
SerialQRDenseSolver()
Default constructor; matrix should be set using setMatrix(), LHS and RHS set with setVectors().
Definition: Teuchos_SerialQRDenseSolver.hpp:400
Teuchos::SerialQRDenseSolver::N_
OrdinalType N_
Definition: Teuchos_SerialQRDenseSolver.hpp:355
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:189
Teuchos::SerialQRDenseSolver::solved_
bool solved_
Definition: Teuchos_SerialQRDenseSolver.hpp:347
Teuchos::SerialQRDenseSolver::TRANS_
Teuchos::ETransp TRANS_
Definition: Teuchos_SerialQRDenseSolver.hpp:352
Teuchos::SerialQRDenseSolver::equilibrationModeA_
OrdinalType equilibrationModeA_
Definition: Teuchos_SerialQRDenseSolver.hpp:343
Teuchos::SerialQRDenseSolver::matrixIsZero_
bool matrixIsZero_
Definition: Teuchos_SerialQRDenseSolver.hpp:348
Teuchos::ScalarTraits::magnitude
static magnitudeType magnitude(T a)
Returns the magnitudeType of the scalar type a.
Definition: Teuchos_ScalarTraitsDecl.hpp:130
Teuchos::Scalar
Definition: Teuchos_YamlParser.cpp:149
Teuchos::EDiagChar
const TEUCHOSNUMERICS_LIB_DLL_EXPORT char EDiagChar[]
Definition: Teuchos_BLAS.cpp:69
Teuchos::ScalarTraits::zero
static T zero()
Returns representation of zero for this scalar type.
Definition: Teuchos_ScalarTraitsDecl.hpp:132
Teuchos::SerialQRDenseSolver::LDA_
OrdinalType LDA_
Definition: Teuchos_SerialQRDenseSolver.hpp:357
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:575
Teuchos::SerialQRDenseSolver::resetVectors
void resetVectors()
Definition: Teuchos_SerialQRDenseSolver.hpp:442
Teuchos::NON_UNIT_DIAG
Definition: Teuchos_BLAS_types.hpp:107
Teuchos::SerialQRDenseSolver::formedQ_
bool formedQ_
Definition: Teuchos_SerialQRDenseSolver.hpp:349
Teuchos::SerialQRDenseSolver::getLHS
RCP< SerialDenseMatrix< OrdinalType, ScalarType > > getLHS() const
Returns pointer to current LHS.
Definition: Teuchos_SerialQRDenseSolver.hpp:308
Teuchos::SerialQRDenseSolver::shouldEquilibrate_
bool shouldEquilibrate_
Definition: Teuchos_SerialQRDenseSolver.hpp:340
Teuchos::SerialQRDenseSolver::getMatrix
RCP< SerialDenseMatrix< OrdinalType, ScalarType > > getMatrix() const
Returns pointer to current matrix.
Definition: Teuchos_SerialQRDenseSolver.hpp:296
Teuchos::NO_TRANS
Definition: Teuchos_BLAS_types.hpp:94
Teuchos::SerialQRDenseSolver::equilibratedB_
bool equilibratedB_
Definition: Teuchos_SerialQRDenseSolver.hpp:342
Teuchos::SerialQRDenseSolver::FactorQ_
RCP< SerialDenseMatrix< OrdinalType, ScalarType > > FactorQ_
Definition: Teuchos_SerialQRDenseSolver.hpp:374
Teuchos::SerialQRDenseSolver::equilibrateMatrix
int equilibrateMatrix()
Equilibrates the this matrix.
Definition: Teuchos_SerialQRDenseSolver.hpp:700
Teuchos::SerialQRDenseSolver::formedR_
bool formedR_
Definition: Teuchos_SerialQRDenseSolver.hpp:350
Teuchos::rcp
TEUCHOS_DEPRECATED RCP< T > rcp(T *p, Dealloc_T dealloc, bool owns_mem)
Deprecated.
Definition: Teuchos_RCPDecl.hpp:1224
Teuchos::BLAS
Templated BLAS wrapper.
Definition: Teuchos_BLAS.hpp:244
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:192
Teuchos::ScalarTraits::prec
static magnitudeType prec()
Returns eps*base.
Definition: Teuchos_ScalarTraitsDecl.hpp:116
Teuchos::SerialQRDenseSolver::computeEquilibrateScaling
int computeEquilibrateScaling()
Determines if this matrix should be scaled.
Definition: Teuchos_SerialQRDenseSolver.hpp:661
Teuchos::SerialQRDenseSolver::equilibrationModeB_
OrdinalType equilibrationModeB_
Definition: Teuchos_SerialQRDenseSolver.hpp:344
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:893
Teuchos::SerialQRDenseSolver::formQ
int formQ()
Explicitly forms the unitary matrix Q.
Definition: Teuchos_SerialQRDenseSolver.hpp:829
Teuchos::SerialQRDenseSolver::~SerialQRDenseSolver
virtual ~SerialQRDenseSolver()
SerialQRDenseSolver destructor.
Definition: Teuchos_SerialQRDenseSolver.hpp:436
Teuchos::SerialQRDenseSolver::resetMatrix
void resetMatrix()
Definition: Teuchos_SerialQRDenseSolver.hpp:453
Teuchos::SerialQRDenseSolver::getFactoredMatrix
RCP< SerialDenseMatrix< OrdinalType, ScalarType > > getFactoredMatrix() const
Returns pointer to factored matrix (assuming factorization has been performed).
Definition: Teuchos_SerialQRDenseSolver.hpp:299
Teuchos::SerialQRDenseSolver::numRows
OrdinalType numRows() const
Returns row dimension of system.
Definition: Teuchos_SerialQRDenseSolver.hpp:314
Teuchos_LAPACK.hpp
Templated interface class to LAPACK routines.
Teuchos::SerialQRDenseSolver
A class for solving dense linear problems.
Definition: Teuchos_SerialQRDenseSolver.hpp:132
Teuchos_CompObject.hpp
Object for storing data and providing functionality that is common to all computational classes.
Teuchos::RCP
Smart reference counting pointer class for automatic garbage collection.
Definition: Teuchos_RCPDecl.hpp:429
Teuchos::SerialQRDenseSolver::A_
ScalarType * A_
Definition: Teuchos_SerialQRDenseSolver.hpp:377
Teuchos::EUploChar
const TEUCHOSNUMERICS_LIB_DLL_EXPORT char EUploChar[]
Definition: Teuchos_BLAS.cpp:68
Teuchos::SerialQRDenseSolver::formedR
bool formedR()
Returns true if R has been formed explicitly.
Definition: Teuchos_SerialQRDenseSolver.hpp:288
details
Teuchos implementation details.
Teuchos::SerialQRDenseSolver::allocateWORK
void allocateWORK()
Definition: Teuchos_SerialQRDenseSolver.hpp:334
Teuchos::TRANS
Definition: Teuchos_BLAS_types.hpp:95
Teuchos::SerialQRDenseSolver::operator=
SerialQRDenseSolver & operator=(const SerialQRDenseSolver< OrdinalType, ScalarType > &Source)
B
Definition: PackageB.cpp:3
Teuchos::SerialQRDenseSolver::setMatrix
int setMatrix(const RCP< SerialDenseMatrix< OrdinalType, ScalarType > > &A)
Sets the pointers for coefficient matrix.
Definition: Teuchos_SerialQRDenseSolver.hpp:482
Teuchos::ScalarTraits::one
static T one()
Returns representation of one for this scalar type.
Definition: Teuchos_ScalarTraitsDecl.hpp:134
Teuchos::LEFT_SIDE
Definition: Teuchos_BLAS_types.hpp:89
Teuchos::SerialQRDenseSolver::Min_MN_
OrdinalType Min_MN_
Definition: Teuchos_SerialQRDenseSolver.hpp:356
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:186
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:323
Teuchos::ScalarTraits
This structure defines some basic traits for a scalar field type.
Definition: Teuchos_ScalarTraitsDecl.hpp:90
Teuchos::SerialQRDenseSolver::getRHS
RCP< SerialDenseMatrix< OrdinalType, ScalarType > > getRHS() const
Returns pointer to current RHS.
Definition: Teuchos_SerialQRDenseSolver.hpp:311
Teuchos::SerialQRDenseSolver::WORK_
std::vector< ScalarType > WORK_
Definition: Teuchos_SerialQRDenseSolver.hpp:381
Teuchos::SerialQRDenseSolver::Matrix_
RCP< SerialDenseMatrix< OrdinalType, ScalarType > > Matrix_
Definition: Teuchos_SerialQRDenseSolver.hpp:369
Teuchos::ESideChar
const TEUCHOSNUMERICS_LIB_DLL_EXPORT char ESideChar[]
Definition: Teuchos_BLAS.cpp:66
Teuchos::CompObject
Functionality and data that is common to all computational classes.
Definition: Teuchos_CompObject.hpp:65
Teuchos_ConfigDefs.hpp
Teuchos header file which uses auto-configuration information to include necessary C++ headers.
TEUCHOS_MIN
#define TEUCHOS_MIN(x, y)
Definition: Teuchos_ConfigDefs.hpp:162
Teuchos_SerialDenseMatrix.hpp
Templated serial dense matrix class.
Teuchos::SerialQRDenseSolver::TAU_
std::vector< ScalarType > TAU_
Definition: Teuchos_SerialQRDenseSolver.hpp:364
Teuchos::SerialQRDenseSolver::LDR_
OrdinalType LDR_
Definition: Teuchos_SerialQRDenseSolver.hpp:360
Teuchos::SerialQRDenseSolver::AF_
ScalarType * AF_
Definition: Teuchos_SerialQRDenseSolver.hpp:378
Teuchos::SerialQRDenseSolver::LWORK_
OrdinalType LWORK_
Definition: Teuchos_SerialQRDenseSolver.hpp:362
Teuchos::SerialQRDenseSolver::formR
int formR()
Explicitly forms the upper triangular matrix R.
Definition: Teuchos_SerialQRDenseSolver.hpp:866
Teuchos::SerialQRDenseSolver::formedQ
bool formedQ()
Returns true if Q has been formed explicitly.
Definition: Teuchos_SerialQRDenseSolver.hpp:285
Teuchos::SerialQRDenseSolver::R_
ScalarType * R_
Definition: Teuchos_SerialQRDenseSolver.hpp:380
Teuchos::SerialQRDenseSolver::LDQ_
OrdinalType LDQ_
Definition: Teuchos_SerialQRDenseSolver.hpp:359
Teuchos::SerialQRDenseSolver::factored_
bool factored_
Definition: Teuchos_SerialQRDenseSolver.hpp:346
Teuchos::SerialQRDenseSolver::equilibrateRHS
int equilibrateRHS()
Equilibrates the current RHS.
Definition: Teuchos_SerialQRDenseSolver.hpp:746
Teuchos::ETypeChar
const TEUCHOSNUMERICS_LIB_DLL_EXPORT char ETypeChar[]
Definition: Teuchos_BLAS.cpp:70
A
Definition: PackageA.cpp:3
Teuchos::SerialQRDenseSolver::LHS_
RCP< SerialDenseMatrix< OrdinalType, ScalarType > > LHS_
Definition: Teuchos_SerialQRDenseSolver.hpp:370
Teuchos::SerialQRDenseSolver::tau
std::vector< ScalarType > tau() const
Returns pointer to pivot vector (if factorization has been computed), zero otherwise.
Definition: Teuchos_SerialQRDenseSolver.hpp:320
Teuchos::null
Definition: Teuchos_ENull.hpp:54
Teuchos::CONJ_TRANS
Definition: Teuchos_BLAS_types.hpp:96
Teuchos::SerialQRDenseSolver::LDAF_
OrdinalType LDAF_
Definition: Teuchos_SerialQRDenseSolver.hpp:358
Teuchos::SerialQRDenseSolver::ANORM_
MagnitudeType ANORM_
Definition: Teuchos_SerialQRDenseSolver.hpp:366
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:509
Teuchos::SerialQRDenseSolver::MagnitudeType
ScalarTraits< ScalarType >::magnitudeType MagnitudeType
Definition: Teuchos_SerialQRDenseSolver.hpp:138
Teuchos_BLAS.hpp
Templated interface class to BLAS routines.
Teuchos::ETranspChar
const TEUCHOSNUMERICS_LIB_DLL_EXPORT char ETranspChar[]
Definition: Teuchos_BLAS.cpp:67
TEUCHOS_MAX
#define TEUCHOS_MAX(x, y)
Definition: Teuchos_ConfigDefs.hpp:161
Teuchos::SerialQRDenseSolver::equilibratedA
bool equilibratedA()
Returns true if factor is equilibrated (factor available via getFactoredMatrix()).
Definition: Teuchos_SerialQRDenseSolver.hpp:273
Teuchos::SerialQRDenseSolver::TMP_
RCP< SerialDenseMatrix< OrdinalType, ScalarType > > TMP_
Definition: Teuchos_SerialQRDenseSolver.hpp:371
Teuchos::Object
The base Teuchos class.
Definition: Teuchos_Object.hpp:68
Teuchos::SerialQRDenseSolver::FactorR_
RCP< SerialDenseMatrix< OrdinalType, ScalarType > > FactorR_
Definition: Teuchos_SerialQRDenseSolver.hpp:375
Teuchos::ESide
ESide
Definition: Teuchos_BLAS_types.hpp:88
Teuchos::SerialQRDenseSolver::RHS_
RCP< SerialDenseMatrix< OrdinalType, ScalarType > > RHS_
Definition: Teuchos_SerialQRDenseSolver.hpp:372
Teuchos::SerialQRDenseSolver::unequilibrateLHS
int unequilibrateLHS()
Unscales the solution vectors if equilibration was used to solve the system.
Definition: Teuchos_SerialQRDenseSolver.hpp:792
Teuchos::SerialDenseMatrix
This class creates and provides basic support for dense rectangular matrix of templated type.
Definition: Teuchos_SerialDenseMatrix.hpp:67
Teuchos::SerialQRDenseSolver::equilibrate_
bool equilibrate_
Definition: Teuchos_SerialQRDenseSolver.hpp:339
Teuchos::SerialQRDenseSolver::shouldEquilibrate
bool shouldEquilibrate()
Returns true if the LAPACK general rules for equilibration suggest you should equilibrate the system.
Definition: Teuchos_SerialQRDenseSolver.hpp:279
Teuchos::SerialQRDenseSolver::factored
bool factored()
Returns true if matrix is factored (factor available via getFactoredMatrix()).
Definition: Teuchos_SerialQRDenseSolver.hpp:270
Teuchos_ScalarTraits.hpp
Defines basic traits for the scalar field type.
Teuchos::SerialQRDenseSolver::numCols
OrdinalType numCols() const
Returns column dimension of system.
Definition: Teuchos_SerialQRDenseSolver.hpp:317
Teuchos
Definition: Teuchos_AbstractFactory.hpp:47
Teuchos::FULL
Definition: Teuchos_BLAS_types.hpp:111
Teuchos::LAPACK
The Templated LAPACK Wrapper Class.
Definition: Teuchos_LAPACK.hpp:96
Teuchos::EUplo
EUplo
Definition: Teuchos_BLAS_types.hpp:99
Teuchos::SerialQRDenseSolver::getR
RCP< SerialDenseMatrix< OrdinalType, ScalarType > > getR() const
Returns pointer to R (assuming factorization has been performed).
Definition: Teuchos_SerialQRDenseSolver.hpp:305
Teuchos::EDiag
EDiag
Definition: Teuchos_BLAS_types.hpp:105
Teuchos::SerialQRDenseSolver::Print
void Print(std::ostream &os) const
Print service methods; defines behavior of ostream << operator.
Definition: Teuchos_SerialQRDenseSolver.hpp:1016
Teuchos::SerialQRDenseSolver::transpose_
bool transpose_
Definition: Teuchos_SerialQRDenseSolver.hpp:345
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:170
Teuchos_SerialDenseSolver.hpp
Templated class for solving dense linear problems.
Teuchos::SerialQRDenseSolver::M_
OrdinalType M_
Definition: Teuchos_SerialQRDenseSolver.hpp:354
Teuchos::ETransp
ETransp
Definition: Teuchos_BLAS_types.hpp:93
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:954
Teuchos::SerialQRDenseSolver::equilibratedA_
bool equilibratedA_
Definition: Teuchos_SerialQRDenseSolver.hpp:341
Teuchos::SerialQRDenseSolver::getQ
RCP< SerialDenseMatrix< OrdinalType, ScalarType > > getQ() const
Returns pointer to Q (assuming factorization has been performed).
Definition: Teuchos_SerialQRDenseSolver.hpp:302
Teuchos::UPPER_TRI
Definition: Teuchos_BLAS_types.hpp:100
Teuchos::SerialQRDenseSolver::INFO_
OrdinalType INFO_
Definition: Teuchos_SerialQRDenseSolver.hpp:361
Generated by   doxygen 1.8.16