00001
00002
00003
00004
00005
00006
00007
00008
00009
00010
00011
00012
00013
00014
00015
00016
00017
00018
00019
00020
00021
00022
00023
00024
00025
00026
00027
00028
00029 #ifndef _TEUCHOS_SERIALDENSESOLVER_HPP_
00030 #define _TEUCHOS_SERIALDENSESOLVER_HPP_
00031
00035 #include "Teuchos_CompObject.hpp"
00036 #include "Teuchos_BLAS.hpp"
00037 #include "Teuchos_LAPACK.hpp"
00038 #include "Teuchos_RCP.hpp"
00039 #include "Teuchos_ConfigDefs.hpp"
00040 #include "Teuchos_SerialDenseMatrix.hpp"
00041
00116 namespace Teuchos {
00117
00118 template<typename OrdinalType, typename ScalarType>
00119 class SerialDenseSolver : public CompObject, public Object, public BLAS<OrdinalType, ScalarType>,
00120 public LAPACK<OrdinalType, ScalarType>
00121 {
00122
00123 public:
00124
00125 typedef typename Teuchos::ScalarTraits<ScalarType>::magnitudeType MagnitudeType;
00126
00128
00129
00130 SerialDenseSolver();
00131
00133 virtual ~SerialDenseSolver();
00135
00137
00138
00140 int setMatrix(const RCP<SerialDenseMatrix<OrdinalType, ScalarType> >& A);
00141
00143
00146 int setVectors(const RCP<SerialDenseMatrix<OrdinalType, ScalarType> >& X,
00147 const RCP<SerialDenseMatrix<OrdinalType, ScalarType> >& B);
00149
00151
00152
00154
00156 void factorWithEquilibration(bool flag) {equilibrate_ = flag; return;}
00157
00159 void solveWithTranspose(bool flag) {transpose_ = flag; if (flag) TRANS_ = Teuchos::TRANS; else TRANS_ = Teuchos::NO_TRANS; return;}
00160
00162 void solveToRefinedSolution(bool flag) {refineSolution_ = flag; return;}
00163
00165
00168 void estimateSolutionErrors(bool flag);
00170
00172
00173
00175
00178 int factor();
00179
00181
00184 int solve();
00185
00187
00190 int invert();
00191
00193
00196 int computeEquilibrateScaling();
00197
00199
00202 int equilibrateMatrix();
00203
00205
00208 int equilibrateRHS();
00209
00211
00214 int applyRefinement();
00215
00217
00220 int unequilibrateLHS();
00221
00223
00229 int reciprocalConditionEstimate(MagnitudeType & Value);
00231
00233
00234
00236 bool transpose() {return(transpose_);}
00237
00239 bool factored() {return(factored_);}
00240
00242 bool equilibratedA() {return(equilibratedA_);}
00243
00245 bool equilibratedB() {return(equilibratedB_);}
00246
00248 bool shouldEquilibrate() {computeEquilibrateScaling(); return(shouldEquilibrate_);}
00249
00251 bool solutionErrorsEstimated() {return(solutionErrorsEstimated_);}
00252
00254 bool inverted() {return(inverted_);}
00255
00257 bool reciprocalConditionEstimated() {return(reciprocalConditionEstimated_);}
00258
00260 bool solved() {return(solved_);}
00261
00263 bool solutionRefined() {return(solutionRefined_);}
00265
00267
00268
00270 RCP<SerialDenseMatrix<OrdinalType, ScalarType> > getMatrix() const {return(Matrix_);}
00271
00273 RCP<SerialDenseMatrix<OrdinalType, ScalarType> > getFactoredMatrix() const {return(Factor_);}
00274
00276 RCP<SerialDenseMatrix<OrdinalType, ScalarType> > getLHS() const {return(LHS_);}
00277
00279 RCP<SerialDenseMatrix<OrdinalType, ScalarType> > getRHS() const {return(RHS_);}
00280
00282 OrdinalType numRows() const {return(M_);}
00283
00285 OrdinalType numCols() const {return(N_);}
00286
00288 std::vector<OrdinalType> IPIV() const {return(IPIV_);}
00289
00291 MagnitudeType ANORM() const {return(ANORM_);}
00292
00294 MagnitudeType RCOND() const {return(RCOND_);}
00295
00297
00299 MagnitudeType ROWCND() const {return(ROWCND_);}
00300
00302
00304 MagnitudeType COLCND() const {return(COLCND_);}
00305
00307 MagnitudeType AMAX() const {return(AMAX_);}
00308
00310 std::vector<MagnitudeType> FERR() const {return(FERR_);}
00311
00313 std::vector<MagnitudeType> BERR() const {return(BERR_);}
00314
00316 std::vector<MagnitudeType> R() const {return(R_);}
00317
00319 std::vector<MagnitudeType> C() const {return(C_);}
00321
00323
00324
00325 void Print(std::ostream& os) const;
00327 protected:
00328
00329 void allocateWORK() { LWORK_ = 4*N_; WORK_.resize( LWORK_ ); return;}
00330 void allocateIWORK() { IWORK_.resize( N_ ); return;}
00331 void resetMatrix();
00332 void resetVectors();
00333
00334
00335 bool equilibrate_;
00336 bool shouldEquilibrate_;
00337 bool equilibratedA_;
00338 bool equilibratedB_;
00339 bool transpose_;
00340 bool factored_;
00341 bool estimateSolutionErrors_;
00342 bool solutionErrorsEstimated_;
00343 bool solved_;
00344 bool inverted_;
00345 bool reciprocalConditionEstimated_;
00346 bool refineSolution_;
00347 bool solutionRefined_;
00348
00349 Teuchos::ETransp TRANS_;
00350
00351 OrdinalType M_;
00352 OrdinalType N_;
00353 OrdinalType Min_MN_;
00354 OrdinalType LDA_;
00355 OrdinalType LDAF_;
00356 OrdinalType INFO_;
00357 OrdinalType LWORK_;
00358
00359 std::vector<OrdinalType> IPIV_;
00360 std::vector<int> IWORK_;
00361
00362 MagnitudeType ANORM_;
00363 MagnitudeType RCOND_;
00364 MagnitudeType ROWCND_;
00365 MagnitudeType COLCND_;
00366 MagnitudeType AMAX_;
00367
00368 RCP<SerialDenseMatrix<OrdinalType, ScalarType> > Matrix_;
00369 RCP<SerialDenseMatrix<OrdinalType, ScalarType> > LHS_;
00370 RCP<SerialDenseMatrix<OrdinalType, ScalarType> > RHS_;
00371 RCP<SerialDenseMatrix<OrdinalType, ScalarType> > Factor_;
00372
00373 ScalarType * A_;
00374 ScalarType * AF_;
00375 std::vector<MagnitudeType> FERR_;
00376 std::vector<MagnitudeType> BERR_;
00377 std::vector<ScalarType> WORK_;
00378 std::vector<MagnitudeType> R_;
00379 std::vector<MagnitudeType> C_;
00380
00381 private:
00382
00383
00384 SerialDenseSolver(const SerialDenseSolver<OrdinalType, ScalarType>& Source);
00385 SerialDenseSolver & operator=(const SerialDenseSolver<OrdinalType, ScalarType>& Source);
00386
00387 };
00388
00389
00390
00391 template<typename OrdinalType, typename ScalarType>
00392 SerialDenseSolver<OrdinalType,ScalarType>::SerialDenseSolver()
00393 : CompObject(),
00394 Object("Teuchos::SerialDenseSolver"),
00395 equilibrate_(false),
00396 shouldEquilibrate_(false),
00397 equilibratedA_(false),
00398 equilibratedB_(false),
00399 transpose_(false),
00400 factored_(false),
00401 estimateSolutionErrors_(false),
00402 solutionErrorsEstimated_(false),
00403 solved_(false),
00404 inverted_(false),
00405 reciprocalConditionEstimated_(false),
00406 refineSolution_(false),
00407 solutionRefined_(false),
00408 TRANS_(Teuchos::NO_TRANS),
00409 M_(0),
00410 N_(0),
00411 Min_MN_(0),
00412 LDA_(0),
00413 LDAF_(0),
00414 INFO_(0),
00415 LWORK_(0),
00416 ANORM_(ScalarTraits<ScalarType>::zero()),
00417 RCOND_(ScalarTraits<ScalarType>::zero()),
00418 ROWCND_(ScalarTraits<ScalarType>::zero()),
00419 COLCND_(ScalarTraits<ScalarType>::zero()),
00420 AMAX_(ScalarTraits<ScalarType>::zero()),
00421 A_(0),
00422 AF_(0)
00423 {
00424 resetMatrix();
00425 }
00426
00427
00428 template<typename OrdinalType, typename ScalarType>
00429 SerialDenseSolver<OrdinalType,ScalarType>::~SerialDenseSolver()
00430 {}
00431
00432
00433
00434 template<typename OrdinalType, typename ScalarType>
00435 void SerialDenseSolver<OrdinalType,ScalarType>::resetVectors()
00436 {
00437 LHS_ = Teuchos::null;
00438 RHS_ = Teuchos::null;
00439 reciprocalConditionEstimated_ = false;
00440 solutionRefined_ = false;
00441 solved_ = false;
00442 solutionErrorsEstimated_ = false;
00443 equilibratedB_ = false;
00444 }
00445
00446
00447 template<typename OrdinalType, typename ScalarType>
00448 void SerialDenseSolver<OrdinalType,ScalarType>::resetMatrix()
00449 {
00450 resetVectors();
00451 equilibratedA_ = false;
00452 factored_ = false;
00453 inverted_ = false;
00454 M_ = 0;
00455 N_ = 0;
00456 Min_MN_ = 0;
00457 LDA_ = 0;
00458 LDAF_ = 0;
00459 ANORM_ = -ScalarTraits<MagnitudeType>::one();
00460 RCOND_ = -ScalarTraits<MagnitudeType>::one();
00461 ROWCND_ = -ScalarTraits<MagnitudeType>::one();
00462 COLCND_ = -ScalarTraits<MagnitudeType>::one();
00463 AMAX_ = -ScalarTraits<MagnitudeType>::one();
00464 A_ = 0;
00465 AF_ = 0;
00466 INFO_ = 0;
00467 LWORK_ = 0;
00468 R_.resize(0);
00469 C_.resize(0);
00470 }
00471
00472
00473 template<typename OrdinalType, typename ScalarType>
00474 int SerialDenseSolver<OrdinalType,ScalarType>::setMatrix(const RCP<SerialDenseMatrix<OrdinalType,ScalarType> >& A)
00475 {
00476 resetMatrix();
00477 Matrix_ = A;
00478 Factor_ = A;
00479 M_ = A->numRows();
00480 N_ = A->numCols();
00481 Min_MN_ = TEUCHOS_MIN(M_,N_);
00482 LDA_ = A->stride();
00483 LDAF_ = LDA_;
00484 A_ = A->values();
00485 AF_ = A->values();
00486 return(0);
00487 }
00488
00489
00490 template<typename OrdinalType, typename ScalarType>
00491 int SerialDenseSolver<OrdinalType,ScalarType>::setVectors(const RCP<SerialDenseMatrix<OrdinalType,ScalarType> >& X,
00492 const RCP<SerialDenseMatrix<OrdinalType,ScalarType> >& B)
00493 {
00494
00495 TEST_FOR_EXCEPTION(B->numRows()!=X->numRows() || B->numCols() != X->numCols(), std::invalid_argument,
00496 "SerialDenseSolver<T>::setVectors: X and B are not the same size!");
00497 TEST_FOR_EXCEPTION(B->values()==0, std::invalid_argument,
00498 "SerialDenseSolver<T>::setVectors: B is an empty SerialDenseMatrix<T>!");
00499 TEST_FOR_EXCEPTION(X->values()==0, std::invalid_argument,
00500 "SerialDenseSolver<T>::setVectors: X is an empty SerialDenseMatrix<T>!");
00501 TEST_FOR_EXCEPTION(B->stride()<1, std::invalid_argument,
00502 "SerialDenseSolver<T>::setVectors: B has an invalid stride!");
00503 TEST_FOR_EXCEPTION(X->stride()<1, std::invalid_argument,
00504 "SerialDenseSolver<T>::setVectors: X has an invalid stride!");
00505
00506 resetVectors();
00507 LHS_ = X;
00508 RHS_ = B;
00509 return(0);
00510 }
00511
00512
00513 template<typename OrdinalType, typename ScalarType>
00514 void SerialDenseSolver<OrdinalType,ScalarType>::estimateSolutionErrors(bool flag)
00515 {
00516 estimateSolutionErrors_ = flag;
00517
00518
00519 refineSolution_ = refineSolution_ || flag;
00520 }
00521
00522
00523 template<typename OrdinalType, typename ScalarType>
00524 int SerialDenseSolver<OrdinalType,ScalarType>::factor() {
00525
00526 if (factored()) return(0);
00527
00528 TEST_FOR_EXCEPTION(inverted(), std::logic_error,
00529 "SerialDenseSolver<T>::factor: Cannot factor an inverted matrix!");
00530
00531 ANORM_ = Matrix_->normOne();
00532
00533
00534
00535
00536
00537 if (A_ == AF_)
00538 if (refineSolution_ ) {
00539 Factor_ = rcp( new SerialDenseMatrix<OrdinalType,ScalarType>(*Matrix_) );
00540 AF_ = Factor_->values();
00541 LDAF_ = Factor_->stride();
00542 }
00543
00544 int ierr = 0;
00545 if (equilibrate_) ierr = equilibrateMatrix();
00546
00547 if (ierr!=0) return(ierr);
00548
00549 if (IPIV_.size()==0) IPIV_.resize( Min_MN_ );
00550
00551 INFO_ = 0;
00552 this->GETRF(M_, N_, AF_, LDAF_, &IPIV_[0], &INFO_);
00553 factored_ = true;
00554
00555 return(INFO_);
00556 }
00557
00558
00559
00560 template<typename OrdinalType, typename ScalarType>
00561 int SerialDenseSolver<OrdinalType,ScalarType>::solve() {
00562
00563
00564
00565
00566
00567
00568
00569
00570
00571
00572 int ierr = 0;
00573 if (equilibrate_) {
00574 ierr = equilibrateRHS();
00575 equilibratedB_ = true;
00576 }
00577 if (ierr != 0) return(ierr);
00578
00579 TEST_FOR_EXCEPTION( (equilibratedA_ && !equilibratedB_) || (!equilibratedA_ && equilibratedB_) ,
00580 std::logic_error, "SerialDenseSolver<T>::solve: Matrix and vectors must be similarly scaled!");
00581 TEST_FOR_EXCEPTION( RHS_==Teuchos::null, std::invalid_argument,
00582 "SerialDenseSolver<T>::solve: No right-hand side vector (RHS) has been set for the linear system!");
00583 TEST_FOR_EXCEPTION( LHS_==Teuchos::null, std::invalid_argument,
00584 "SerialDenseSolver<T>::solve: No solution vector (LHS) has been set for the linear system!");
00585
00586 if (shouldEquilibrate() && !equilibratedA_)
00587 std::cout << "WARNING! SerialDenseSolver<T>::solve: System should be equilibrated!" << std::endl;
00588
00589 if (inverted()) {
00590
00591 TEST_FOR_EXCEPTION( RHS_->values() == LHS_->values(), std::invalid_argument,
00592 "SerialDenseSolver<T>::solve: X and B must be different vectors if matrix is inverted.");
00593
00594 INFO_ = 0;
00595 this->GEMM(TRANS_, Teuchos::NO_TRANS, N_, RHS_->numCols(), N_, 1.0, AF_, LDAF_,
00596 RHS_->values(), RHS_->stride(), 0.0, LHS_->values(), LHS_->stride());
00597 if (INFO_!=0) return(INFO_);
00598 solved_ = true;
00599 }
00600 else {
00601
00602 if (!factored()) factor();
00603
00604 if (RHS_->values()!=LHS_->values()) {
00605 (*LHS_) = (*RHS_);
00606 }
00607 INFO_ = 0;
00608 this->GETRS(ETranspChar[TRANS_], N_, RHS_->numCols(), AF_, LDAF_, &IPIV_[0], LHS_->values(), LHS_->stride(), &INFO_);
00609 if (INFO_!=0) return(INFO_);
00610 solved_ = true;
00611
00612 }
00613 int ierr1=0;
00614 if (refineSolution_ && !inverted()) ierr1 = applyRefinement();
00615 if (ierr1!=0)
00616 return(ierr1);
00617 else
00618 return(ierr);
00619
00620 if (equilibrate_) ierr1 = unequilibrateLHS();
00621 return(ierr1);
00622 }
00623
00624
00625 template<typename OrdinalType, typename ScalarType>
00626 int SerialDenseSolver<OrdinalType,ScalarType>::applyRefinement()
00627 {
00628 TEST_FOR_EXCEPTION(!solved(), std::logic_error,
00629 "SerialDenseSolver<T>::applyRefinement: Must have an existing solution!");
00630 TEST_FOR_EXCEPTION(A_==AF_, std::logic_error,
00631 "SerialDenseSolver<T>::applyRefinement: Cannot apply refinement if no original copy of A!");
00632
00633 OrdinalType NRHS = RHS_->numCols();
00634 FERR_.resize( NRHS );
00635 BERR_.resize( NRHS );
00636 allocateWORK();
00637 allocateIWORK();
00638
00639 INFO_ = 0;
00640 this->GERFS(ETranspChar[TRANS_], N_, NRHS, A_, LDA_, AF_, LDAF_, &IPIV_[0],
00641 RHS_->values(), RHS_->stride(), LHS_->values(), LHS_->stride(),
00642 &FERR_[0], &BERR_[0], &WORK_[0], &IWORK_[0], &INFO_);
00643
00644 solutionErrorsEstimated_ = true;
00645 reciprocalConditionEstimated_ = true;
00646 solutionRefined_ = true;
00647
00648 return(INFO_);
00649
00650 }
00651
00652
00653
00654 template<typename OrdinalType, typename ScalarType>
00655 int SerialDenseSolver<OrdinalType,ScalarType>::computeEquilibrateScaling()
00656 {
00657 if (R_.size()!=0) return(0);
00658
00659 R_.resize( M_ );
00660 C_.resize( N_ );
00661
00662 INFO_ = 0;
00663 this->GEEQU (M_, N_, AF_, LDAF_, &R_[0], &C_[0], &ROWCND_, &COLCND_, &AMAX_, &INFO_);
00664
00665 if (COLCND_<0.1*ScalarTraits<MagnitudeType>::one() ||
00666 ROWCND_<0.1*ScalarTraits<MagnitudeType>::one() ||
00667 AMAX_ < ScalarTraits<ScalarType>::rmin() || AMAX_ > ScalarTraits<ScalarType>::rmax() )
00668 shouldEquilibrate_ = true;
00669
00670 return(INFO_);
00671 }
00672
00673
00674
00675 template<typename OrdinalType, typename ScalarType>
00676 int SerialDenseSolver<OrdinalType,ScalarType>::equilibrateMatrix()
00677 {
00678 OrdinalType i, j;
00679 int ierr = 0;
00680
00681 if (equilibratedA_) return(0);
00682 if (R_.size()==0) ierr = computeEquilibrateScaling();
00683 if (ierr!=0) return(ierr);
00684 if (A_==AF_) {
00685 ScalarType * ptr;
00686 for (j=0; j<N_; j++) {
00687 ptr = A_ + j*LDA_;
00688 ScalarType s1 = C_[j];
00689 for (i=0; i<M_; i++) {
00690 *ptr = *ptr*s1*R_[i];
00691 ptr++;
00692 }
00693 }
00694 }
00695 else {
00696 ScalarType * ptr;
00697 ScalarType * ptr1;
00698 for (j=0; j<N_; j++) {
00699 ptr = A_ + j*LDA_;
00700 ptr1 = AF_ + j*LDAF_;
00701 ScalarType s1 = C_[j];
00702 for (i=0; i<M_; i++) {
00703 *ptr = *ptr*s1*R_[i];
00704 ptr++;
00705 *ptr1 = *ptr1*s1*R_[i];
00706 ptr1++;
00707 }
00708 }
00709 }
00710
00711 equilibratedA_ = true;
00712
00713 return(0);
00714 }
00715
00716
00717
00718 template<typename OrdinalType, typename ScalarType>
00719 int SerialDenseSolver<OrdinalType,ScalarType>::equilibrateRHS()
00720 {
00721 OrdinalType i, j;
00722 int ierr = 0;
00723
00724 if (equilibratedB_) return(0);
00725 if (R_.size()==0) ierr = computeEquilibrateScaling();
00726 if (ierr!=0) return(ierr);
00727
00728 ScalarType * R_tmp = &R_[0];
00729 if (transpose_) R_tmp = &C_[0];
00730
00731 OrdinalType LDB = RHS_->stride(), NRHS = RHS_->numCols();
00732 ScalarType * B = RHS_->values();
00733 ScalarType * ptr;
00734 for (j=0; j<NRHS; j++) {
00735 ptr = B + j*LDB;
00736 for (i=0; i<M_; i++) {
00737 *ptr = *ptr*R_tmp[i];
00738 ptr++;
00739 }
00740 }
00741
00742 equilibratedB_ = true;
00743
00744 return(0);
00745 }
00746
00747
00748
00749 template<typename OrdinalType, typename ScalarType>
00750 int SerialDenseSolver<OrdinalType,ScalarType>::unequilibrateLHS()
00751 {
00752 OrdinalType i, j;
00753
00754 if (!equilibratedB_) return(0);
00755
00756 ScalarType * C_tmp = &C_[0];
00757 if (transpose_) C_tmp = &R_[0];
00758
00759 OrdinalType LDX = RHS_->stride(), NRHS = RHS_->numCols();
00760 ScalarType * X = RHS_->values();
00761 ScalarType * ptr;
00762 for (j=0; j<NRHS; j++) {
00763 ptr = X + j*LDX;
00764 for (i=0; i<N_; i++) {
00765 *ptr = *ptr*C_tmp[i];
00766 ptr++;
00767 }
00768 }
00769
00770 return(0);
00771 }
00772
00773
00774
00775 template<typename OrdinalType, typename ScalarType>
00776 int SerialDenseSolver<OrdinalType,ScalarType>::invert()
00777 {
00778 if (!factored()) factor();
00779
00780
00781
00782
00783
00784
00785
00786
00787
00788
00789
00790
00791
00792 allocateWORK();
00793
00794 INFO_ = 0;
00795 this->GETRI( N_, AF_, LDAF_, &IPIV_[0], &WORK_[0], LWORK_, &INFO_);
00796
00797 inverted_ = true;
00798 factored_ = false;
00799
00800 return(INFO_);
00801 }
00802
00803
00804
00805 template<typename OrdinalType, typename ScalarType>
00806 int SerialDenseSolver<OrdinalType,ScalarType>::reciprocalConditionEstimate(MagnitudeType & Value)
00807 {
00808 if (reciprocalConditionEstimated()) {
00809 Value = RCOND_;
00810 return(0);
00811 }
00812
00813 if ( ANORM_<ScalarTraits<MagnitudeType>::zero() ) ANORM_ = Matrix_->normOne();
00814
00815 int ierr = 0;
00816 if (!factored()) ierr = factor();
00817 if (ierr!=0) return(ierr);
00818
00819 allocateWORK();
00820 allocateIWORK();
00821
00822
00823 INFO_ = 0;
00824 this->GECON( '1', N_, AF_, LDAF_, ANORM_, &RCOND_, &WORK_[0], &IWORK_[0], &INFO_);
00825 reciprocalConditionEstimated_ = true;
00826 Value = RCOND_;
00827
00828 return(INFO_);
00829 }
00830
00831
00832 template<typename OrdinalType, typename ScalarType>
00833 void SerialDenseSolver<OrdinalType,ScalarType>::Print(std::ostream& os) const {
00834
00835 if (Matrix_!=Teuchos::null) os << "Solver Matrix" << std::endl << *Matrix_ << std::endl;
00836 if (Factor_!=Teuchos::null) os << "Solver Factored Matrix" << std::endl << *Factor_ << std::endl;
00837 if (LHS_ !=Teuchos::null) os << "Solver LHS" << std::endl << *LHS_ << std::endl;
00838 if (RHS_ !=Teuchos::null) os << "Solver RHS" << std::endl << *RHS_ << std::endl;
00839
00840 }
00841
00842 }
00843
00844 #endif