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//
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_
49#include "Teuchos_BLAS.hpp"
50#include "Teuchos_LAPACK.hpp"
51#include "Teuchos_RCP.hpp"
56
57#ifdef HAVE_TEUCHOSNUMERICS_EIGEN
58#include "Eigen/Dense"
59#endif
60
129namespace 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
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
159
161 virtual ~SerialQRDenseSolver();
163
165
166
168
171
173
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
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
255
257
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
297
300
303
306
309
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
330 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
399template<typename OrdinalType, typename ScalarType>
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
435template<typename OrdinalType, typename ScalarType>
438
439//=============================================================================
440
441template<typename OrdinalType, typename ScalarType>
443{
444 LHS_ = Teuchos::null;
445 TMP_ = Teuchos::null;
446 RHS_ = Teuchos::null;
447 solved_ = false;
448 equilibratedB_ = false;
449}
450//=============================================================================
451
452template<typename OrdinalType, typename ScalarType>
453void 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;
472 A_ = 0;
473 AF_ = 0;
474 Q_ = 0;
475 R_ = 0;
476 INFO_ = 0;
477 LWORK_ = 0;
478}
479//=============================================================================
480
481template<typename OrdinalType, typename ScalarType>
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
508template<typename OrdinalType, typename ScalarType>
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
532template<typename OrdinalType, typename ScalarType>
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_) );
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
574template<typename OrdinalType, typename ScalarType>
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
660template<typename OrdinalType, typename ScalarType>
662{
668
671
672 // Compute maximum absolute value of matrix entries
673 OrdinalType i, j;
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
699template<typename OrdinalType, typename ScalarType>
701{
702 if (equilibratedA_) return(0);
703
709
712 OrdinalType BW = 0;
713
714 // Compute maximum absolute value of matrix entries
715 OrdinalType i, j;
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
745template<typename OrdinalType, typename ScalarType>
747{
748 if (equilibratedB_) return(0);
749
755
758 OrdinalType BW = 0;
759
760 // Compute maximum absolute value of rhs entries
761 OrdinalType i, j;
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
791template<typename OrdinalType, typename ScalarType>
793{
794 if (!equilibratedB_) return(0);
795
801 (void) mZero; // Silence "unused variable" compiler warning.
802
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
828template<typename OrdinalType, typename ScalarType>
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
865template<typename OrdinalType, typename ScalarType>
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
892template<typename OrdinalType, typename ScalarType>
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++ ) {
919 }
920 }
921 }
922#else
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++ ) {
942 }
943 }
944 }
945#endif
946
947 return(INFO_);
948
949}
950
951//=============================================================================
952
953template<typename OrdinalType, typename ScalarType>
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
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
1015template<typename OrdinalType, typename ScalarType>
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_ */
Templated interface class to BLAS routines.
Object for storing data and providing functionality that is common to all computational classes.
Teuchos header file which uses auto-configuration information to include necessary C++ headers.
Templated interface class to LAPACK routines.
Reference-counted pointer class and non-member templated function implementations.
Defines basic traits for the scalar field type.
Templated serial dense matrix class.
Templated class for solving dense linear problems.
Templated BLAS wrapper.
Functionality and data that is common to all computational classes.
The Templated LAPACK Wrapper Class.
The base Teuchos class.
Smart reference counting pointer class for automatic garbage collection.
RCP(ENull null_arg=null)
Initialize RCP<T> to NULL.
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.
A class for solving dense linear problems.
void factorWithEquilibration(bool flag)
Causes equilibration to be called just before the matrix factorization as part of the call to factor.
int solveR(ETransp transr, SerialDenseMatrix< OrdinalType, ScalarType > &C)
Solve input matrix on the left with the upper triangular matrix R or its adjoint.
bool solved()
Returns true if the current set of vectors has been solved.
bool equilibratedB()
Returns true if RHS is equilibrated (RHS available via getRHS()).
RCP< SerialDenseMatrix< OrdinalType, ScalarType > > getLHS() const
Returns pointer to current LHS.
bool formedR()
Returns true if R has been formed explicitly.
int formR()
Explicitly forms the upper triangular matrix R.
bool shouldEquilibrate()
Returns true if the LAPACK general rules for equilibration suggest you should equilibrate the system.
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).
RCP< SerialDenseMatrix< OrdinalType, ScalarType > > getMatrix() const
Returns pointer to current matrix.
void solveWithTranspose(bool flag)
If flag is true, causes all subsequent function calls to work with the adjoint of this matrix,...
bool transpose()
Returns true if adjoint of this matrix has and will be used.
int multiplyQ(ETransp transq, SerialDenseMatrix< OrdinalType, ScalarType > &C)
Left multiply the input matrix by the unitary matrix Q or its adjoint.
RCP< SerialDenseMatrix< OrdinalType, ScalarType > > getQ() const
Returns pointer to Q (assuming factorization has been performed).
MagnitudeType ANORM() const
Returns the absolute value of the largest element of this matrix (returns -1 if not yet computed).
OrdinalType numCols() const
Returns column dimension of system.
bool formedQ()
Returns true if Q has been formed explicitly.
std::vector< ScalarType > tau() const
Returns pointer to pivot vector (if factorization has been computed), zero otherwise.
OrdinalType numRows() const
Returns row dimension of system.
virtual ~SerialQRDenseSolver()
SerialQRDenseSolver destructor.
int factor()
Computes the in-place QR factorization of the matrix using the LAPACK routine _GETRF or the Eigen cla...
RCP< SerialDenseMatrix< OrdinalType, ScalarType > > getFactoredMatrix() const
Returns pointer to factored matrix (assuming factorization has been performed).
int solve()
Computes the solution X to AX = B for the this matrix and the B provided to SetVectors()....
SerialQRDenseSolver()
Default constructor; matrix should be set using setMatrix(), LHS and RHS set with setVectors().
void solveWithTransposeFlag(Teuchos::ETransp trans)
All subsequent function calls will work with the transpose-type set by this method (Teuchos::NO_TRANS...
bool factored()
Returns true if matrix is factored (factor available via getFactoredMatrix()).
void Print(std::ostream &os) const
Print service methods; defines behavior of ostream << operator.
RCP< SerialDenseMatrix< OrdinalType, ScalarType > > getR() const
Returns pointer to R (assuming factorization has been performed).
int setMatrix(const RCP< SerialDenseMatrix< OrdinalType, ScalarType > > &A)
Sets the pointers for coefficient matrix.
int formQ()
Explicitly forms the unitary matrix Q.
int equilibrateMatrix()
Equilibrates the this matrix.
RCP< SerialDenseMatrix< OrdinalType, ScalarType > > getRHS() const
Returns pointer to current RHS.
int equilibrateRHS()
Equilibrates the current RHS.
bool equilibratedA()
Returns true if factor is equilibrated (factor available via getFactoredMatrix()).
int computeEquilibrateScaling()
Determines if this matrix should be scaled.
int unequilibrateLHS()
Unscales the solution vectors if equilibration was used to solve the system.
#define TEUCHOS_TEST_FOR_EXCEPTION(throw_exception_test, Exception, msg)
Macro for throwing an exception with breakpointing to ease debugging.
The Teuchos namespace contains all of the classes, structs and enums used by Teuchos,...
Teuchos implementation details.
This structure defines some basic traits for a scalar field type.
static magnitudeType magnitude(T a)
Returns the magnitudeType of the scalar type a.
static T one()
Returns representation of one for this scalar type.
static T zero()
Returns representation of zero for this scalar type.
static magnitudeType sfmin()
Returns safe minimum (sfmin), such that 1/sfmin does not overflow.
static magnitudeType prec()
Returns eps*base.