Teuchos - Trilinos Tools Package Version of the Day
Loading...
Searching...
No Matches
Teuchos_SerialDenseSolver.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_SERIALDENSESOLVER_HPP_
43#define _TEUCHOS_SERIALDENSESOLVER_HPP_
49#include "Teuchos_BLAS.hpp"
50#include "Teuchos_LAPACK.hpp"
51#include "Teuchos_RCP.hpp"
55
56#ifdef HAVE_TEUCHOSNUMERICS_EIGEN
57#include "Eigen/Dense"
58#endif
59
134namespace Teuchos {
135
136 template<typename OrdinalType, typename ScalarType>
137 class SerialDenseSolver : public CompObject, public Object, public BLAS<OrdinalType, ScalarType>,
138 public LAPACK<OrdinalType, ScalarType>
139 {
140
141 public:
142
143 typedef typename Teuchos::ScalarTraits<ScalarType>::magnitudeType MagnitudeType;
144#ifdef HAVE_TEUCHOSNUMERICS_EIGEN
145 typedef typename Eigen::Matrix<ScalarType,Eigen::Dynamic,Eigen::Dynamic,Eigen::ColMajor> EigenMatrix;
146 typedef typename Eigen::Matrix<ScalarType,Eigen::Dynamic,1> EigenVector;
147 typedef typename Eigen::InnerStride<Eigen::Dynamic> EigenInnerStride;
148 typedef typename Eigen::OuterStride<Eigen::Dynamic> EigenOuterStride;
149 typedef typename Eigen::Map<EigenVector,0,EigenInnerStride > EigenVectorMap;
150 typedef typename Eigen::Map<const EigenVector,0,EigenInnerStride > EigenConstVectorMap;
151 typedef typename Eigen::Map<EigenMatrix,0,EigenOuterStride > EigenMatrixMap;
152 typedef typename Eigen::Map<const EigenMatrix,0,EigenOuterStride > EigenConstMatrixMap;
153 typedef typename Eigen::PermutationMatrix<Eigen::Dynamic,Eigen::Dynamic,OrdinalType> EigenPermutationMatrix;
154 typedef typename Eigen::Array<OrdinalType,Eigen::Dynamic,1> EigenOrdinalArray;
155 typedef typename Eigen::Map<EigenOrdinalArray> EigenOrdinalArrayMap;
156 typedef typename Eigen::Array<ScalarType,Eigen::Dynamic,1> EigenScalarArray;
157 typedef typename Eigen::Map<EigenScalarArray> EigenScalarArrayMap;
158#endif
159
161
162
164
166 virtual ~SerialDenseSolver();
168
170
171
174
176
182
184
185
187
189 void factorWithEquilibration(bool flag) {equilibrate_ = flag; return;}
190
192
194 void solveWithTranspose(bool flag) {transpose_ = flag; if (flag) TRANS_ = Teuchos::TRANS; else TRANS_ = Teuchos::NO_TRANS; return;}
195
197
199 void solveWithTransposeFlag(Teuchos::ETransp trans) {TRANS_ = trans; transpose_ = (trans != Teuchos::NO_TRANS) ? true : false; }
200
202
204 void solveToRefinedSolution(bool flag) {refineSolution_ = flag; return;}
205
207
211 void estimateSolutionErrors(bool flag);
213
215
216
218
221 int factor();
222
224
227 int solve();
228
230
233 int invert();
234
236
240
242
246 int equilibrateMatrix();
247
249
253 int equilibrateRHS();
254
256
260 int applyRefinement();
261
263
267 int unequilibrateLHS();
268
270
276 int reciprocalConditionEstimate(MagnitudeType & Value);
278
280
281
283 bool transpose() {return(transpose_);}
284
286 bool factored() {return(factored_);}
287
289 bool equilibratedA() {return(equilibratedA_);}
290
292 bool equilibratedB() {return(equilibratedB_);}
293
295 bool shouldEquilibrate() {computeEquilibrateScaling(); return(shouldEquilibrate_);}
296
298 bool solutionErrorsEstimated() {return(solutionErrorsEstimated_);}
299
301 bool inverted() {return(inverted_);}
302
304 bool reciprocalConditionEstimated() {return(reciprocalConditionEstimated_);}
305
307 bool solved() {return(solved_);}
308
310 bool solutionRefined() {return(solutionRefined_);}
312
314
315
318
321
324
327
329 OrdinalType numRows() const {return(M_);}
330
332 OrdinalType numCols() const {return(N_);}
333
335 std::vector<OrdinalType> IPIV() const {return(IPIV_);}
336
338 MagnitudeType ANORM() const {return(ANORM_);}
339
341 MagnitudeType RCOND() const {return(RCOND_);}
342
344
346 MagnitudeType ROWCND() const {return(ROWCND_);}
347
349
351 MagnitudeType COLCND() const {return(COLCND_);}
352
354 MagnitudeType AMAX() const {return(AMAX_);}
355
357 std::vector<MagnitudeType> FERR() const {return(FERR_);}
358
360 std::vector<MagnitudeType> BERR() const {return(BERR_);}
361
363 std::vector<MagnitudeType> R() const {return(R_);}
364
366 std::vector<MagnitudeType> C() const {return(C_);}
368
370
371
372 void Print(std::ostream& os) const;
374 protected:
375
376 void allocateWORK() { LWORK_ = 4*N_; WORK_.resize( LWORK_ ); return;}
377 void resetMatrix();
378 void resetVectors();
379
380
381 bool equilibrate_;
382 bool shouldEquilibrate_;
383 bool equilibratedA_;
384 bool equilibratedB_;
385 bool transpose_;
386 bool factored_;
387 bool estimateSolutionErrors_;
388 bool solutionErrorsEstimated_;
389 bool solved_;
390 bool inverted_;
391 bool reciprocalConditionEstimated_;
392 bool refineSolution_;
393 bool solutionRefined_;
394
395 Teuchos::ETransp TRANS_;
396
397 OrdinalType M_;
398 OrdinalType N_;
399 OrdinalType Min_MN_;
400 OrdinalType LDA_;
401 OrdinalType LDAF_;
402 OrdinalType INFO_;
403 OrdinalType LWORK_;
404
405 std::vector<OrdinalType> IPIV_;
406
407 MagnitudeType ANORM_;
408 MagnitudeType RCOND_;
409 MagnitudeType ROWCND_;
410 MagnitudeType COLCND_;
411 MagnitudeType AMAX_;
412
413 RCP<SerialDenseMatrix<OrdinalType, ScalarType> > Matrix_;
414 RCP<SerialDenseMatrix<OrdinalType, ScalarType> > LHS_;
415 RCP<SerialDenseMatrix<OrdinalType, ScalarType> > RHS_;
416 RCP<SerialDenseMatrix<OrdinalType, ScalarType> > Factor_;
417
418 ScalarType * A_;
419 ScalarType * AF_;
420 std::vector<MagnitudeType> FERR_;
421 std::vector<MagnitudeType> BERR_;
422 std::vector<ScalarType> WORK_;
423 std::vector<MagnitudeType> R_;
424 std::vector<MagnitudeType> C_;
425#ifdef HAVE_TEUCHOSNUMERICS_EIGEN
426 Eigen::PartialPivLU<EigenMatrix> lu_;
427#endif
428
429 private:
430 // SerialDenseSolver copy constructor (put here because we don't want user access)
431
432 SerialDenseSolver(const SerialDenseSolver<OrdinalType, ScalarType>& Source);
433 SerialDenseSolver & operator=(const SerialDenseSolver<OrdinalType, ScalarType>& Source);
434
435 };
436
437 namespace details {
438
439 // Helper traits to distinguish work arrays for real and complex-valued datatypes.
440 template<typename T>
441 struct lapack_traits {
442 typedef int iwork_type;
443 };
444
445 // Complex-valued specialization
446 template<typename T>
447 struct lapack_traits<std::complex<T> > {
448 typedef typename ScalarTraits<T>::magnitudeType iwork_type;
449 };
450
451 } // end namespace details
452
453//=============================================================================
454
455template<typename OrdinalType, typename ScalarType>
457 : CompObject(),
458 equilibrate_(false),
459 shouldEquilibrate_(false),
460 equilibratedA_(false),
461 equilibratedB_(false),
462 transpose_(false),
463 factored_(false),
464 estimateSolutionErrors_(false),
465 solutionErrorsEstimated_(false),
466 solved_(false),
467 inverted_(false),
468 reciprocalConditionEstimated_(false),
469 refineSolution_(false),
470 solutionRefined_(false),
471 TRANS_(Teuchos::NO_TRANS),
472 M_(0),
473 N_(0),
474 Min_MN_(0),
475 LDA_(0),
476 LDAF_(0),
477 INFO_(0),
478 LWORK_(0),
479 ANORM_(ScalarTraits<MagnitudeType>::zero()),
480 RCOND_(ScalarTraits<MagnitudeType>::zero()),
481 ROWCND_(ScalarTraits<MagnitudeType>::zero()),
482 COLCND_(ScalarTraits<MagnitudeType>::zero()),
483 AMAX_(ScalarTraits<MagnitudeType>::zero()),
484 A_(0),
485 AF_(0)
486{
487 resetMatrix();
488}
489//=============================================================================
490
491template<typename OrdinalType, typename ScalarType>
494
495//=============================================================================
496
497template<typename OrdinalType, typename ScalarType>
499{
500 LHS_ = Teuchos::null;
501 RHS_ = Teuchos::null;
502 reciprocalConditionEstimated_ = false;
503 solutionRefined_ = false;
504 solved_ = false;
505 solutionErrorsEstimated_ = false;
506 equilibratedB_ = false;
507}
508//=============================================================================
509
510template<typename OrdinalType, typename ScalarType>
511void SerialDenseSolver<OrdinalType,ScalarType>::resetMatrix()
512{
513 resetVectors();
514 equilibratedA_ = false;
515 factored_ = false;
516 inverted_ = false;
517 M_ = 0;
518 N_ = 0;
519 Min_MN_ = 0;
520 LDA_ = 0;
521 LDAF_ = 0;
527 A_ = 0;
528 AF_ = 0;
529 INFO_ = 0;
530 LWORK_ = 0;
531 R_.resize(0);
532 C_.resize(0);
533}
534//=============================================================================
535
536template<typename OrdinalType, typename ScalarType>
538{
539 resetMatrix();
540 Matrix_ = A;
541 Factor_ = A;
542 M_ = A->numRows();
543 N_ = A->numCols();
544 Min_MN_ = TEUCHOS_MIN(M_,N_);
545 LDA_ = A->stride();
546 LDAF_ = LDA_;
547 A_ = A->values();
548 AF_ = A->values();
549 return(0);
550}
551//=============================================================================
552
553template<typename OrdinalType, typename ScalarType>
556{
557 // Check that these new vectors are consistent.
558 TEUCHOS_TEST_FOR_EXCEPTION(B->numRows()!=X->numRows() || B->numCols() != X->numCols(), std::invalid_argument,
559 "SerialDenseSolver<T>::setVectors: X and B are not the same size!");
560 TEUCHOS_TEST_FOR_EXCEPTION(B->values()==0, std::invalid_argument,
561 "SerialDenseSolver<T>::setVectors: B is an empty SerialDenseMatrix<T>!");
562 TEUCHOS_TEST_FOR_EXCEPTION(X->values()==0, std::invalid_argument,
563 "SerialDenseSolver<T>::setVectors: X is an empty SerialDenseMatrix<T>!");
564 TEUCHOS_TEST_FOR_EXCEPTION(B->stride()<1, std::invalid_argument,
565 "SerialDenseSolver<T>::setVectors: B has an invalid stride!");
566 TEUCHOS_TEST_FOR_EXCEPTION(X->stride()<1, std::invalid_argument,
567 "SerialDenseSolver<T>::setVectors: X has an invalid stride!");
568
569 resetVectors();
570 LHS_ = X;
571 RHS_ = B;
572 return(0);
573}
574//=============================================================================
575
576template<typename OrdinalType, typename ScalarType>
578{
579 estimateSolutionErrors_ = flag;
580
581 // If the errors are estimated, this implies that the solution must be refined
582 refineSolution_ = refineSolution_ || flag;
583}
584//=============================================================================
585
586template<typename OrdinalType, typename ScalarType>
588
589 if (factored()) return(0); // Already factored
590
591 TEUCHOS_TEST_FOR_EXCEPTION(inverted(), std::logic_error,
592 "SerialDenseSolver<T>::factor: Cannot factor an inverted matrix!");
593
594 ANORM_ = Matrix_->normOne(); // Compute 1-Norm of A
595
596
597 // If we want to refine the solution, then the factor must
598 // be stored separatedly from the original matrix
599
600 if (A_ == AF_)
601 if (refineSolution_ ) {
603 AF_ = Factor_->values();
604 LDAF_ = Factor_->stride();
605 }
606
607 int ierr = 0;
608 if (equilibrate_) ierr = equilibrateMatrix();
609
610 if (ierr!=0) return(ierr);
611
612 if ((int)IPIV_.size()!=Min_MN_) IPIV_.resize( Min_MN_ ); // Allocated Pivot vector if not already done.
613
614 INFO_ = 0;
615
616#ifdef HAVE_TEUCHOSNUMERICS_EIGEN
617 EigenMatrixMap aMap( AF_, M_, N_, EigenOuterStride(LDAF_) );
620 EigenOrdinalArrayMap ipivMap( &IPIV_[0], Min_MN_ );
621
622 lu_.compute(aMap);
623 p = lu_.permutationP();
624 ind = p.indices();
625
626 EigenMatrix luMat = lu_.matrixLU();
627 for (OrdinalType j=0; j<aMap.outerSize(); j++) {
628 for (OrdinalType i=0; i<aMap.innerSize(); i++) {
629 aMap(i,j) = luMat(i,j);
630 }
631 }
632 for (OrdinalType i=0; i<ipivMap.innerSize(); i++) {
633 ipivMap(i) = ind(i);
634 }
635#else
636 this->GETRF(M_, N_, AF_, LDAF_, &IPIV_[0], &INFO_);
637#endif
638
639 factored_ = true;
640
641 return(INFO_);
642}
643
644//=============================================================================
645
646template<typename OrdinalType, typename ScalarType>
648
649 // We will call one of four routines depending on what services the user wants and
650 // whether or not the matrix has been inverted or factored already.
651 //
652 // If the matrix has been inverted, use DGEMM to compute solution.
653 // Otherwise, if the user want the matrix to be equilibrated or wants a refined solution, we will
654 // call the X interface.
655 // Otherwise, if the matrix is already factored we will call the TRS interface.
656 // Otherwise, if the matrix is unfactored we will call the SV interface.
657
658 int ierr = 0;
659 if (equilibrate_) {
660 ierr = equilibrateRHS();
661 }
662 if (ierr != 0) return(ierr); // Can't equilibrate B, so return.
663
664 TEUCHOS_TEST_FOR_EXCEPTION( RHS_==Teuchos::null, std::invalid_argument,
665 "SerialDenseSolver<T>::solve: No right-hand side vector (RHS) has been set for the linear system!");
666 TEUCHOS_TEST_FOR_EXCEPTION( LHS_==Teuchos::null, std::invalid_argument,
667 "SerialDenseSolver<T>::solve: No solution vector (LHS) has been set for the linear system!");
668
669 if (inverted()) {
670
671 TEUCHOS_TEST_FOR_EXCEPTION( RHS_->values() == LHS_->values(), std::invalid_argument,
672 "SerialDenseSolver<T>::solve: X and B must be different vectors if matrix is inverted.");
673 TEUCHOS_TEST_FOR_EXCEPTION( (equilibratedA_ && !equilibratedB_) || (!equilibratedA_ && equilibratedB_) ,
674 std::logic_error, "SerialDenseSolver<T>::solve: Matrix and vectors must be similarly scaled!");
675
676 INFO_ = 0;
677 this->GEMM(TRANS_, Teuchos::NO_TRANS, N_, RHS_->numCols(), N_, 1.0, AF_, LDAF_,
678 RHS_->values(), RHS_->stride(), 0.0, LHS_->values(), LHS_->stride());
679 if (INFO_!=0) return(INFO_);
680 solved_ = true;
681 }
682 else {
683
684 if (!factored()) factor(); // Matrix must be factored
685
686 TEUCHOS_TEST_FOR_EXCEPTION( (equilibratedA_ && !equilibratedB_) || (!equilibratedA_ && equilibratedB_) ,
687 std::logic_error, "SerialDenseSolver<T>::solve: Matrix and vectors must be similarly scaled!");
688
689 if (RHS_->values()!=LHS_->values()) {
690 (*LHS_) = (*RHS_); // Copy B to X if needed
691 }
692 INFO_ = 0;
693
694#ifdef HAVE_TEUCHOSNUMERICS_EIGEN
695 EigenMatrixMap rhsMap( RHS_->values(), RHS_->numRows(), RHS_->numCols(), EigenOuterStride(RHS_->stride()) );
696 EigenMatrixMap lhsMap( LHS_->values(), LHS_->numRows(), LHS_->numCols(), EigenOuterStride(LHS_->stride()) );
697 if ( TRANS_ == Teuchos::NO_TRANS ) {
698 lhsMap = lu_.solve(rhsMap);
699 } else if ( TRANS_ == Teuchos::TRANS ) {
700 lu_.matrixLU().template triangularView<Eigen::Upper>().transpose().solveInPlace(lhsMap);
701 lu_.matrixLU().template triangularView<Eigen::UnitLower>().transpose().solveInPlace(lhsMap);
702 EigenMatrix x = lu_.permutationP().transpose() * lhsMap;
703 lhsMap = x;
704 } else {
705 lu_.matrixLU().template triangularView<Eigen::Upper>().adjoint().solveInPlace(lhsMap);
706 lu_.matrixLU().template triangularView<Eigen::UnitLower>().adjoint().solveInPlace(lhsMap);
707 EigenMatrix x = lu_.permutationP().transpose() * lhsMap;
708 lhsMap = x;
709 }
710#else
711 this->GETRS(ETranspChar[TRANS_], N_, RHS_->numCols(), AF_, LDAF_, &IPIV_[0], LHS_->values(), LHS_->stride(), &INFO_);
712#endif
713
714 if (INFO_!=0) return(INFO_);
715 solved_ = true;
716
717 }
718
719 // Warn the user that the matrix should be equilibrated, if it isn't being done already.
720 if (shouldEquilibrate() && !equilibratedA_)
721 std::cout << "WARNING! SerialDenseSolver<T>::solve: System should be equilibrated!" << std::endl;
722
723 int ierr1=0;
724 if (refineSolution_ && !inverted()) ierr1 = applyRefinement();
725 if (ierr1!=0)
726 return(ierr1);
727
728 if (equilibrate_) ierr1 = unequilibrateLHS();
729 return(ierr1);
730}
731//=============================================================================
732
733template<typename OrdinalType, typename ScalarType>
735{
736 TEUCHOS_TEST_FOR_EXCEPTION(!solved(), std::logic_error,
737 "SerialDenseSolver<T>::applyRefinement: Must have an existing solution!");
738 TEUCHOS_TEST_FOR_EXCEPTION(A_==AF_, std::logic_error,
739 "SerialDenseSolver<T>::applyRefinement: Cannot apply refinement if no original copy of A!");
740
741#ifdef HAVE_TEUCHOSNUMERICS_EIGEN
742 // Implement templated GERFS or use Eigen.
743 return (-1);
744#else
745
746 OrdinalType NRHS = RHS_->numCols();
747 FERR_.resize( NRHS );
748 BERR_.resize( NRHS );
749 allocateWORK();
750
751 INFO_ = 0;
752 std::vector<typename details::lapack_traits<ScalarType>::iwork_type> GERFS_WORK( N_ );
753 this->GERFS(ETranspChar[TRANS_], N_, NRHS, A_, LDA_, AF_, LDAF_, &IPIV_[0],
754 RHS_->values(), RHS_->stride(), LHS_->values(), LHS_->stride(),
755 &FERR_[0], &BERR_[0], &WORK_[0], &GERFS_WORK[0], &INFO_);
756
757 solutionErrorsEstimated_ = true;
758 reciprocalConditionEstimated_ = true;
759 solutionRefined_ = true;
760
761 return(INFO_);
762#endif
763}
764
765//=============================================================================
766
767template<typename OrdinalType, typename ScalarType>
769{
770 if (R_.size()!=0) return(0); // Already computed
771
772 R_.resize( M_ );
773 C_.resize( N_ );
774
775 INFO_ = 0;
776 this->GEEQU (M_, N_, AF_, LDAF_, &R_[0], &C_[0], &ROWCND_, &COLCND_, &AMAX_, &INFO_);
777
778 if (COLCND_<0.1*ScalarTraits<MagnitudeType>::one() ||
779 ROWCND_<0.1*ScalarTraits<MagnitudeType>::one() ||
781 shouldEquilibrate_ = true;
782
783 return(INFO_);
784}
785
786//=============================================================================
787
788template<typename OrdinalType, typename ScalarType>
790{
791 OrdinalType i, j;
792 int ierr = 0;
793
794 if (equilibratedA_) return(0); // Already done.
795 if (R_.size()==0) ierr = computeEquilibrateScaling(); // Compute R and C if needed.
796 if (ierr!=0) return(ierr); // If return value is not zero, then can't equilibrate.
797 if (A_==AF_) {
798 ScalarType * ptr;
799 for (j=0; j<N_; j++) {
800 ptr = A_ + j*LDA_;
801 ScalarType s1 = C_[j];
802 for (i=0; i<M_; i++) {
803 *ptr = *ptr*s1*R_[i];
804 ptr++;
805 }
806 }
807 }
808 else {
809 ScalarType * ptr;
811 for (j=0; j<N_; j++) {
812 ptr = A_ + j*LDA_;
813 ptr1 = AF_ + j*LDAF_;
814 ScalarType s1 = C_[j];
815 for (i=0; i<M_; i++) {
816 *ptr = *ptr*s1*R_[i];
817 ptr++;
818 *ptr1 = *ptr1*s1*R_[i];
819 ptr1++;
820 }
821 }
822 }
823
824 equilibratedA_ = true;
825
826 return(0);
827}
828
829//=============================================================================
830
831template<typename OrdinalType, typename ScalarType>
833{
834 OrdinalType i, j;
835 int ierr = 0;
836
837 if (equilibratedB_) return(0); // Already done.
838 if (R_.size()==0) ierr = computeEquilibrateScaling(); // Compute R and C if needed.
839 if (ierr!=0) return(ierr); // Can't count on R and C being computed.
840
841 MagnitudeType * R_tmp = &R_[0];
842 if (transpose_) R_tmp = &C_[0];
843
844 OrdinalType LDB = RHS_->stride(), NRHS = RHS_->numCols();
845 ScalarType * B = RHS_->values();
846 ScalarType * ptr;
847 for (j=0; j<NRHS; j++) {
848 ptr = B + j*LDB;
849 for (i=0; i<M_; i++) {
850 *ptr = *ptr*R_tmp[i];
851 ptr++;
852 }
853 }
854
855 equilibratedB_ = true;
856
857 return(0);
858}
859
860//=============================================================================
861
862template<typename OrdinalType, typename ScalarType>
864{
865 OrdinalType i, j;
866
867 if (!equilibratedB_) return(0); // Nothing to do
868
869 MagnitudeType * C_tmp = &C_[0];
870 if (transpose_) C_tmp = &R_[0];
871
872 OrdinalType LDX = LHS_->stride(), NLHS = LHS_->numCols();
873 ScalarType * X = LHS_->values();
874 ScalarType * ptr;
875 for (j=0; j<NLHS; j++) {
876 ptr = X + j*LDX;
877 for (i=0; i<N_; i++) {
878 *ptr = *ptr*C_tmp[i];
879 ptr++;
880 }
881 }
882
883 return(0);
884}
885
886//=============================================================================
887
888template<typename OrdinalType, typename ScalarType>
890{
891
892 if (!factored()) factor(); // Need matrix factored.
893
894#ifdef HAVE_TEUCHOSNUMERICS_EIGEN
895 EigenMatrixMap invMap( AF_, M_, N_, EigenOuterStride(LDAF_) );
896 EigenMatrix invMat = lu_.inverse();
897 for (OrdinalType j=0; j<invMap.outerSize(); j++) {
898 for (OrdinalType i=0; i<invMap.innerSize(); i++) {
899 invMap(i,j) = invMat(i,j);
900 }
901 }
902#else
903
904 /* This section work with LAPACK Version 3.0 only
905 // Setting LWORK = -1 and calling GETRI will return optimal work space size in WORK_TMP
906 OrdinalType LWORK_TMP = -1;
907 ScalarType WORK_TMP;
908 GETRI( N_, AF_, LDAF_, &IPIV_[0], &WORK_TMP, &LWORK_TMP, &INFO_);
909 LWORK_TMP = WORK_TMP; // Convert to integer
910 if (LWORK_TMP>LWORK_) {
911 LWORK_ = LWORK_TMP;
912 WORK_.resize( LWORK_ );
913 }
914 */
915 // This section will work with any version of LAPACK
916 allocateWORK();
917
918 INFO_ = 0;
919 this->GETRI( N_, AF_, LDAF_, &IPIV_[0], &WORK_[0], LWORK_, &INFO_);
920#endif
921
922 inverted_ = true;
923 factored_ = false;
924
925 return(INFO_);
926
927}
928
929//=============================================================================
930
931template<typename OrdinalType, typename ScalarType>
933{
934#ifdef HAVE_TEUCHOSNUMERICS_EIGEN
935 // Implement templated GECON or use Eigen. Eigen currently doesn't have a condition estimation function.
936 return (-1);
937#else
938 if (reciprocalConditionEstimated()) {
939 Value = RCOND_;
940 return(0); // Already computed, just return it.
941 }
942
943 if ( ANORM_<ScalarTraits<MagnitudeType>::zero() ) ANORM_ = Matrix_->normOne();
944
945 int ierr = 0;
946 if (!factored()) ierr = factor(); // Need matrix factored.
947 if (ierr!=0) return(ierr);
948
949 allocateWORK();
950
951 // We will assume a one-norm condition number
952 INFO_ = 0;
953 std::vector<typename details::lapack_traits<ScalarType>::iwork_type> GECON_WORK( 2*N_ );
954 this->GECON( '1', N_, AF_, LDAF_, ANORM_, &RCOND_, &WORK_[0], &GECON_WORK[0], &INFO_);
955
956 reciprocalConditionEstimated_ = true;
957 Value = RCOND_;
958
959 return(INFO_);
960#endif
961}
962//=============================================================================
963
964template<typename OrdinalType, typename ScalarType>
966
967 if (Matrix_!=Teuchos::null) os << "Solver Matrix" << std::endl << *Matrix_ << std::endl;
968 if (Factor_!=Teuchos::null) os << "Solver Factored Matrix" << std::endl << *Factor_ << std::endl;
969 if (LHS_ !=Teuchos::null) os << "Solver LHS" << std::endl << *LHS_ << std::endl;
970 if (RHS_ !=Teuchos::null) os << "Solver RHS" << std::endl << *RHS_ << std::endl;
971
972}
973
974} // namespace Teuchos
975
976#endif /* _TEUCHOS_SERIALDENSESOLVER_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 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.
Ptr< T > ptr() const
Get a safer wrapper raw C++ pointer to the underlying object.
A class for solving dense linear problems.
std::vector< MagnitudeType > BERR() const
Returns a pointer to the backward error estimates computed by LAPACK.
MagnitudeType ROWCND() const
Ratio of smallest to largest row scale factors for the this matrix (returns -1 if not yet computed).
bool solutionRefined()
Returns true if the current set of vectors has been refined.
void factorWithEquilibration(bool flag)
Causes equilibration to be called just before the matrix factorization as part of the call to factor.
OrdinalType numRows() const
Returns row dimension of system.
bool factored()
Returns true if matrix is factored (factor available via getFactoredMatrix()).
RCP< SerialDenseMatrix< OrdinalType, ScalarType > > getRHS() const
Returns pointer to current RHS.
MagnitudeType RCOND() const
Returns the reciprocal of the condition number of the this matrix (returns -1 if not yet computed).
RCP< SerialDenseMatrix< OrdinalType, ScalarType > > getLHS() const
Returns pointer to current LHS.
std::vector< MagnitudeType > C() const
Returns a pointer to the column scale vector used for equilibration.
int setMatrix(const RCP< SerialDenseMatrix< OrdinalType, ScalarType > > &A)
Sets the pointers for coefficient matrix.
int computeEquilibrateScaling()
Computes the scaling vector S(i) = 1/sqrt(A(i,i)) of the this matrix.
int equilibrateRHS()
Equilibrates the current RHS.
bool inverted()
Returns true if matrix inverse has been computed (inverse available via getFactoredMatrix()).
bool equilibratedA()
Returns true if factor is equilibrated (factor available via getFactoredMatrix()).
MagnitudeType COLCND() const
Ratio of smallest to largest column scale factors for the this matrix (returns -1 if not yet computed...
RCP< SerialDenseMatrix< OrdinalType, ScalarType > > getMatrix() const
Returns pointer to current matrix.
bool shouldEquilibrate()
Returns true if the LAPACK general rules for equilibration suggest you should equilibrate the system.
bool solved()
Returns true if the current set of vectors has been solved.
std::vector< MagnitudeType > FERR() const
Returns a pointer to the forward error estimates computed by LAPACK.
RCP< SerialDenseMatrix< OrdinalType, ScalarType > > getFactoredMatrix() const
Returns pointer to factored matrix (assuming factorization has been performed).
void estimateSolutionErrors(bool flag)
Causes all solves to estimate the forward and backward solution error.
int invert()
Inverts the this matrix.
int factor()
Computes the in-place LU factorization of the matrix using the LAPACK routine _GETRF.
int equilibrateMatrix()
Equilibrates the this matrix.
void solveWithTransposeFlag(Teuchos::ETransp trans)
All subsequent function calls will work with the transpose-type set by this method (Teuchos::NO_TRANS...
bool reciprocalConditionEstimated()
Returns true if the condition number of the this matrix has been computed (value available via Recipr...
int reciprocalConditionEstimate(MagnitudeType &Value)
Returns the reciprocal of the 1-norm condition number of the this matrix.
void Print(std::ostream &os) const
Print service methods; defines behavior of ostream << operator.
int applyRefinement()
Apply Iterative Refinement.
int unequilibrateLHS()
Unscales the solution vectors if equilibration was used to solve the system.
MagnitudeType AMAX() const
Returns the absolute value of the largest entry of the this matrix (returns -1 if not yet computed).
void solveToRefinedSolution(bool flag)
Causes all solves to compute solution to best ability using iterative refinement.
void solveWithTranspose(bool flag)
If flag is true, causes all subsequent function calls to work with the transpose of this matrix,...
virtual ~SerialDenseSolver()
SerialDenseSolver destructor.
SerialDenseSolver()
Default constructor; matrix should be set using setMatrix(), LHS and RHS set with 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).
OrdinalType numCols() const
Returns column dimension of system.
bool transpose()
Returns true if transpose of this matrix has and will be used.
std::vector< OrdinalType > IPIV() const
Returns pointer to pivot vector (if factorization has been computed), zero otherwise.
MagnitudeType ANORM() const
Returns the 1-Norm of the this matrix (returns -1 if not yet computed).
bool equilibratedB()
Returns true if RHS is equilibrated (RHS available via getRHS()).
int solve()
Computes the solution X to AX = B for the this matrix and the B provided to SetVectors()....
std::vector< MagnitudeType > R() const
Returns a pointer to the row scaling vector used for equilibration.
bool solutionErrorsEstimated()
Returns true if forward and backward error estimated have been computed (available via FERR() and BER...
#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 T one()
Returns representation of one for this scalar type.
T magnitudeType
Mandatory typedef for result of magnitude.