Teuchos - Trilinos Tools Package Version of the Day
Loading...
Searching...
No Matches
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_
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
97namespace 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
126 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
298 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
367template<typename OrdinalType, typename ScalarType>
368SerialQRDenseSolver<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
403template<typename OrdinalType, typename ScalarType>
406
407//=============================================================================
408
409template<typename OrdinalType, typename ScalarType>
410void 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
420template<typename OrdinalType, typename ScalarType>
421void 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
449template<typename OrdinalType, typename ScalarType>
450int 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
476template<typename OrdinalType, typename ScalarType>
477int 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
500template<typename OrdinalType, typename ScalarType>
501int 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
542template<typename OrdinalType, typename ScalarType>
543int 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
628template<typename OrdinalType, typename ScalarType>
629int 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
667template<typename OrdinalType, typename ScalarType>
668int 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
713template<typename OrdinalType, typename ScalarType>
714int 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
759template<typename OrdinalType, typename ScalarType>
760int 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
796template<typename OrdinalType, typename ScalarType>
797int 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
833template<typename OrdinalType, typename ScalarType>
834int 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
860template<typename OrdinalType, typename ScalarType>
861int 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
921template<typename OrdinalType, typename ScalarType>
922int 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
983template<typename OrdinalType, typename ScalarType>
984void 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_BLAS.hpp
Templated interface class to BLAS routines.
Teuchos_CompObject.hpp
Object for storing data and providing functionality that is common to all computational classes.
Teuchos_ConfigDefs.hpp
Teuchos header file which uses auto-configuration information to include necessary C++ headers.
Teuchos_LAPACK.hpp
Templated interface class to LAPACK routines.
Teuchos_RCP.hpp
Reference-counted pointer class and non-member templated function implementations.
Teuchos_ScalarTraits.hpp
Defines basic traits for the scalar field type.
Teuchos_SerialDenseMatrix.hpp
Templated serial dense matrix class.
Teuchos_SerialDenseSolver.hpp
Templated class for solving dense linear problems.
Teuchos::BLAS
Templated BLAS wrapper.
Definition Teuchos_BLAS.hpp:213
Teuchos::CompObject
Functionality and data that is common to all computational classes.
Definition Teuchos_CompObject.hpp:34
Teuchos::LAPACK
The Templated LAPACK Wrapper Class.
Definition Teuchos_LAPACK.hpp:65
Teuchos::Object
The base Teuchos class.
Definition Teuchos_Object.hpp:36
Teuchos::RCP
Smart reference counting pointer class for automatic garbage collection.
Definition Teuchos_RCPDecl.hpp:397
Teuchos::RCP::rcp
RCP< T > rcp(const boost::shared_ptr< T > &sptr)
Conversion function that takes in a boost::shared_ptr object and spits out a Teuchos::RCP object.
Teuchos::SerialQRDenseSolver
A class for solving dense linear problems.
Definition Teuchos_SerialQRDenseSolver.hpp:102
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::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::solved
bool solved()
Returns true if the current set of vectors has been solved.
Definition Teuchos_SerialQRDenseSolver.hpp:250
Teuchos::SerialQRDenseSolver::equilibratedB
bool equilibratedB()
Returns true if RHS is equilibrated (RHS available via getRHS()).
Definition Teuchos_SerialQRDenseSolver.hpp:244
Teuchos::SerialQRDenseSolver::getLHS
RCP< SerialDenseMatrix< OrdinalType, ScalarType > > getLHS() const
Returns pointer to current LHS.
Definition Teuchos_SerialQRDenseSolver.hpp:276
Teuchos::SerialQRDenseSolver::formedR
bool formedR()
Returns true if R has been formed explicitly.
Definition Teuchos_SerialQRDenseSolver.hpp:256
Teuchos::SerialQRDenseSolver::formR
int formR()
Explicitly forms the upper triangular matrix R.
Definition Teuchos_SerialQRDenseSolver.hpp:834
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::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::SerialQRDenseSolver::getMatrix
RCP< SerialDenseMatrix< OrdinalType, ScalarType > > getMatrix() const
Returns pointer to current matrix.
Definition Teuchos_SerialQRDenseSolver.hpp:264
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::transpose
bool transpose()
Returns true if adjoint of this matrix has and will be used.
Definition Teuchos_SerialQRDenseSolver.hpp:235
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::getQ
RCP< SerialDenseMatrix< OrdinalType, ScalarType > > getQ() const
Returns pointer to Q (assuming factorization has been performed).
Definition Teuchos_SerialQRDenseSolver.hpp:270
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::numCols
OrdinalType numCols() const
Returns column dimension of system.
Definition Teuchos_SerialQRDenseSolver.hpp:285
Teuchos::SerialQRDenseSolver::formedQ
bool formedQ()
Returns true if Q has been formed explicitly.
Definition Teuchos_SerialQRDenseSolver.hpp:253
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::numRows
OrdinalType numRows() const
Returns row dimension of system.
Definition Teuchos_SerialQRDenseSolver.hpp:282
Teuchos::SerialQRDenseSolver::~SerialQRDenseSolver
virtual ~SerialQRDenseSolver()
SerialQRDenseSolver destructor.
Definition Teuchos_SerialQRDenseSolver.hpp:404
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::getFactoredMatrix
RCP< SerialDenseMatrix< OrdinalType, ScalarType > > getFactoredMatrix() const
Returns pointer to factored matrix (assuming factorization has been performed).
Definition Teuchos_SerialQRDenseSolver.hpp:267
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::SerialQRDenseSolver::SerialQRDenseSolver
SerialQRDenseSolver()
Default constructor; matrix should be set using setMatrix(), LHS and RHS set with setVectors().
Definition Teuchos_SerialQRDenseSolver.hpp:368
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::SerialQRDenseSolver::factored
bool factored()
Returns true if matrix is factored (factor available via getFactoredMatrix()).
Definition Teuchos_SerialQRDenseSolver.hpp:238
Teuchos::SerialQRDenseSolver::Print
void Print(std::ostream &os) const
Print service methods; defines behavior of ostream << operator.
Definition Teuchos_SerialQRDenseSolver.hpp:984
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::setMatrix
int setMatrix(const RCP< SerialDenseMatrix< OrdinalType, ScalarType > > &A)
Sets the pointers for coefficient matrix.
Definition Teuchos_SerialQRDenseSolver.hpp:450
Teuchos::SerialQRDenseSolver::formQ
int formQ()
Explicitly forms the unitary matrix Q.
Definition Teuchos_SerialQRDenseSolver.hpp:797
Teuchos::SerialQRDenseSolver::equilibrateMatrix
int equilibrateMatrix()
Equilibrates the this matrix.
Definition Teuchos_SerialQRDenseSolver.hpp:668
Teuchos::SerialQRDenseSolver::getRHS
RCP< SerialDenseMatrix< OrdinalType, ScalarType > > getRHS() const
Returns pointer to current RHS.
Definition Teuchos_SerialQRDenseSolver.hpp:279
Teuchos::SerialQRDenseSolver::equilibrateRHS
int equilibrateRHS()
Equilibrates the current RHS.
Definition Teuchos_SerialQRDenseSolver.hpp:714
Teuchos::SerialQRDenseSolver::equilibratedA
bool equilibratedA()
Returns true if factor is equilibrated (factor available via getFactoredMatrix()).
Definition Teuchos_SerialQRDenseSolver.hpp:241
Teuchos::SerialQRDenseSolver::computeEquilibrateScaling
int computeEquilibrateScaling()
Determines if this matrix should be scaled.
Definition Teuchos_SerialQRDenseSolver.hpp:629
Teuchos::SerialQRDenseSolver::unequilibrateLHS
int unequilibrateLHS()
Unscales the solution vectors if equilibration was used to solve the system.
Definition Teuchos_SerialQRDenseSolver.hpp:760
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
The Teuchos namespace contains all of the classes, structs and enums used by Teuchos,...
details
Teuchos implementation details.
Teuchos::ScalarTraits
This structure defines some basic traits for a scalar field type.
Definition Teuchos_ScalarTraitsDecl.hpp:59
Teuchos::ScalarTraits::magnitude
static magnitudeType magnitude(T a)
Returns the magnitudeType of the scalar type a.
Definition Teuchos_ScalarTraitsDecl.hpp:100
Teuchos::ScalarTraits::one
static T one()
Returns representation of one for this scalar type.
Definition Teuchos_ScalarTraitsDecl.hpp:104
Teuchos::ScalarTraits::zero
static T zero()
Returns representation of zero for this scalar type.
Definition Teuchos_ScalarTraitsDecl.hpp:102
Teuchos::ScalarTraits::sfmin
static magnitudeType sfmin()
Returns safe minimum (sfmin), such that 1/sfmin does not overflow.
Definition Teuchos_ScalarTraitsDecl.hpp:82
Teuchos::ScalarTraits::prec
static magnitudeType prec()
Returns eps*base.
Definition Teuchos_ScalarTraitsDecl.hpp:86
Generated by doxygen 1.9.8