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
00030
00031
00032 #ifndef TEUCHOS_SERIALSPDDENSESOLVER_H
00033 #define TEUCHOS_SERIALSPDDENSESOLVER_H
00034
00039 #include "Teuchos_CompObject.hpp"
00040 #include "Teuchos_BLAS.hpp"
00041 #include "Teuchos_LAPACK.hpp"
00042 #include "Teuchos_RCP.hpp"
00043 #include "Teuchos_ConfigDefs.hpp"
00044 #include "Teuchos_SerialSymDenseMatrix.hpp"
00045 #include "Teuchos_SerialDenseMatrix.hpp"
00046
00114 namespace Teuchos {
00115
00116 template<typename OrdinalType, typename ScalarType>
00117 class SerialSpdDenseSolver : public CompObject, public Object, public BLAS<OrdinalType, ScalarType>,
00118 public LAPACK<OrdinalType, ScalarType>
00119 {
00120 public:
00121
00122 typedef typename Teuchos::ScalarTraits<ScalarType>::magnitudeType MagnitudeType;
00123
00125
00126
00127 SerialSpdDenseSolver();
00128
00129
00131 virtual ~SerialSpdDenseSolver();
00133
00135
00136
00138 int setMatrix(const RCP<SerialSymDenseMatrix<OrdinalType,ScalarType> >& A_in);
00139
00141
00144 int setVectors(const RCP<SerialDenseMatrix<OrdinalType, ScalarType> >& X,
00145 const RCP<SerialDenseMatrix<OrdinalType, ScalarType> >& B);
00147
00149
00150
00152
00154 void factorWithEquilibration(bool flag) {equilibrate_ = flag; return;}
00155
00157 void solveToRefinedSolution(bool flag) {refineSolution_ = flag; return;}
00158
00160
00163 void estimateSolutionErrors(bool flag);
00165
00167
00168
00170
00173 int factor();
00174
00176
00179 int solve();
00180
00182
00187 int invert();
00188
00190
00193 int computeEquilibrateScaling();
00194
00196
00199 int equilibrateMatrix();
00200
00202
00205 int equilibrateRHS();
00206
00207
00209
00212 int applyRefinement();
00213
00215
00218 int unequilibrateLHS();
00219
00221
00227 int reciprocalConditionEstimate(MagnitudeType & Value);
00229
00231
00232
00234 bool transpose() {return(transpose_);}
00235
00237 bool factored() {return(factored_);}
00238
00240 bool equilibratedA() {return(equilibratedA_);}
00241
00243 bool equilibratedB() {return(equilibratedB_);}
00244
00246 bool shouldEquilibrate() {computeEquilibrateScaling(); return(shouldEquilibrate_);}
00247
00249 bool solutionErrorsEstimated() {return(solutionErrorsEstimated_);}
00250
00252 bool inverted() {return(inverted_);}
00253
00255 bool reciprocalConditionEstimated() {return(reciprocalConditionEstimated_);}
00256
00258 bool solved() {return(solved_);}
00259
00261 bool solutionRefined() {return(solutionRefined_);}
00263
00265
00266
00268 RCP<SerialSymDenseMatrix<OrdinalType, ScalarType> > getMatrix() const {return(Matrix_);}
00269
00271 RCP<SerialSymDenseMatrix<OrdinalType, ScalarType> > getFactoredMatrix() const {return(Factor_);}
00272
00274 RCP<SerialDenseMatrix<OrdinalType, ScalarType> > getLHS() const {return(LHS_);}
00275
00277 RCP<SerialDenseMatrix<OrdinalType, ScalarType> > getRHS() const {return(RHS_);}
00278
00280 OrdinalType numRows() const {return(numRowCols_);}
00281
00283 OrdinalType numCols() const {return(numRowCols_);}
00284
00286 MagnitudeType ANORM() const {return(ANORM_);}
00287
00289 MagnitudeType RCOND() const {return(RCOND_);}
00290
00292
00294 MagnitudeType SCOND() {return(SCOND_);};
00295
00297 MagnitudeType AMAX() const {return(AMAX_);}
00298
00300 std::vector<ScalarType> FERR() const {return(FERR_);}
00301
00303 std::vector<ScalarType> BERR() const {return(BERR_);}
00304
00306 std::vector<ScalarType> R() const {return(R_);}
00307
00309
00311
00312
00313 void Print(std::ostream& os) const;
00315
00316 protected:
00317
00318 void allocateWORK() { LWORK_ = 4*numRowCols_; WORK_.resize( LWORK_ ); return;}
00319 void allocateIWORK() { IWORK_.resize( numRowCols_ ); return;}
00320 void resetMatrix();
00321 void resetVectors();
00322
00323 bool equilibrate_;
00324 bool shouldEquilibrate_;
00325 bool equilibratedA_;
00326 bool equilibratedB_;
00327 bool transpose_;
00328 bool factored_;
00329 bool estimateSolutionErrors_;
00330 bool solutionErrorsEstimated_;
00331 bool solved_;
00332 bool inverted_;
00333 bool reciprocalConditionEstimated_;
00334 bool refineSolution_;
00335 bool solutionRefined_;
00336
00337 OrdinalType numRowCols_;
00338 OrdinalType LDA_;
00339 OrdinalType LDAF_;
00340 OrdinalType INFO_;
00341 OrdinalType LWORK_;
00342
00343 std::vector<int> IWORK_;
00344
00345 MagnitudeType ANORM_;
00346 MagnitudeType RCOND_;
00347 MagnitudeType SCOND_;
00348 MagnitudeType AMAX_;
00349
00350 RCP<SerialSymDenseMatrix<OrdinalType, ScalarType> > Matrix_;
00351 RCP<SerialDenseMatrix<OrdinalType, ScalarType> > LHS_;
00352 RCP<SerialDenseMatrix<OrdinalType, ScalarType> > RHS_;
00353 RCP<SerialSymDenseMatrix<OrdinalType, ScalarType> > Factor_;
00354
00355 ScalarType * A_;
00356 ScalarType * AF_;
00357 std::vector<ScalarType> FERR_;
00358 std::vector<ScalarType> BERR_;
00359 std::vector<ScalarType> WORK_;
00360 std::vector<ScalarType> R_;
00361
00362 private:
00363
00364
00365 SerialSpdDenseSolver(const SerialSpdDenseSolver<OrdinalType, ScalarType>& Source);
00366 SerialSpdDenseSolver & operator=(const SerialSpdDenseSolver<OrdinalType, ScalarType>& Source);
00367
00368 };
00369
00370
00371
00372 template<typename OrdinalType, typename ScalarType>
00373 SerialSpdDenseSolver<OrdinalType,ScalarType>::SerialSpdDenseSolver()
00374 : CompObject(),
00375 Object("Teuchos::SerialSpdDenseSolver"),
00376 equilibrate_(false),
00377 shouldEquilibrate_(false),
00378 equilibratedA_(false),
00379 equilibratedB_(false),
00380 transpose_(false),
00381 factored_(false),
00382 estimateSolutionErrors_(false),
00383 solutionErrorsEstimated_(false),
00384 solved_(false),
00385 inverted_(false),
00386 reciprocalConditionEstimated_(false),
00387 refineSolution_(false),
00388 solutionRefined_(false),
00389 numRowCols_(0),
00390 LDA_(0),
00391 LDAF_(0),
00392 INFO_(0),
00393 LWORK_(0),
00394 ANORM_(ScalarTraits<ScalarType>::zero()),
00395 RCOND_(ScalarTraits<ScalarType>::zero()),
00396 SCOND_(ScalarTraits<ScalarType>::zero()),
00397 AMAX_(ScalarTraits<ScalarType>::zero()),
00398 A_(0),
00399 AF_(0)
00400 {
00401 resetMatrix();
00402 }
00403
00404
00405 template<typename OrdinalType, typename ScalarType>
00406 SerialSpdDenseSolver<OrdinalType,ScalarType>::~SerialSpdDenseSolver()
00407 {}
00408
00409
00410
00411 template<typename OrdinalType, typename ScalarType>
00412 void SerialSpdDenseSolver<OrdinalType,ScalarType>::resetVectors()
00413 {
00414 LHS_ = Teuchos::null;
00415 RHS_ = Teuchos::null;
00416 reciprocalConditionEstimated_ = false;
00417 solutionRefined_ = false;
00418 solved_ = false;
00419 solutionErrorsEstimated_ = false;
00420 equilibratedB_ = false;
00421 }
00422
00423
00424 template<typename OrdinalType, typename ScalarType>
00425 void SerialSpdDenseSolver<OrdinalType,ScalarType>::resetMatrix()
00426 {
00427 resetVectors();
00428 equilibratedA_ = false;
00429 factored_ = false;
00430 inverted_ = false;
00431 numRowCols_ = 0;
00432 LDA_ = 0;
00433 LDAF_ = 0;
00434 ANORM_ = -ScalarTraits<MagnitudeType>::one();
00435 RCOND_ = -ScalarTraits<MagnitudeType>::one();
00436 SCOND_ = -ScalarTraits<MagnitudeType>::one();
00437 AMAX_ = -ScalarTraits<MagnitudeType>::one();
00438 A_ = 0;
00439 AF_ = 0;
00440 INFO_ = 0;
00441 LWORK_ = 0;
00442 R_.resize(0);
00443 }
00444
00445
00446 template<typename OrdinalType, typename ScalarType>
00447 int SerialSpdDenseSolver<OrdinalType,ScalarType>::setMatrix(const RCP<SerialSymDenseMatrix<OrdinalType,ScalarType> >& A)
00448 {
00449 resetMatrix();
00450 Matrix_ = A;
00451 Factor_ = A;
00452 numRowCols_ = A->numRows();
00453 LDA_ = A->stride();
00454 LDAF_ = LDA_;
00455 A_ = A->values();
00456 AF_ = A->values();
00457 return(0);
00458 }
00459
00460
00461 template<typename OrdinalType, typename ScalarType>
00462 int SerialSpdDenseSolver<OrdinalType,ScalarType>::setVectors(const RCP<SerialDenseMatrix<OrdinalType,ScalarType> >& X,
00463 const RCP<SerialDenseMatrix<OrdinalType,ScalarType> >& B)
00464 {
00465
00466 TEST_FOR_EXCEPTION(B->numRows()!=X->numRows() || B->numCols() != X->numCols(), std::invalid_argument,
00467 "SerialSpdDenseSolver<T>::setVectors: X and B are not the same size!");
00468 TEST_FOR_EXCEPTION(B->values()==0, std::invalid_argument,
00469 "SerialSpdDenseSolver<T>::setVectors: B is an empty SerialDenseMatrix<T>!");
00470 TEST_FOR_EXCEPTION(X->values()==0, std::invalid_argument,
00471 "SerialSpdDenseSolver<T>::setVectors: X is an empty SerialDenseMatrix<T>!");
00472 TEST_FOR_EXCEPTION(B->stride()<1, std::invalid_argument,
00473 "SerialSpdDenseSolver<T>::setVectors: B has an invalid stride!");
00474 TEST_FOR_EXCEPTION(X->stride()<1, std::invalid_argument,
00475 "SerialSpdDenseSolver<T>::setVectors: X has an invalid stride!");
00476
00477 resetVectors();
00478 LHS_ = X;
00479 RHS_ = B;
00480 return(0);
00481 }
00482
00483
00484 template<typename OrdinalType, typename ScalarType>
00485 void SerialSpdDenseSolver<OrdinalType,ScalarType>::estimateSolutionErrors(bool flag)
00486 {
00487 estimateSolutionErrors_ = flag;
00488
00489
00490 refineSolution_ = refineSolution_ || flag;
00491 }
00492
00493
00494 template<typename OrdinalType, typename ScalarType>
00495 int SerialSpdDenseSolver<OrdinalType,ScalarType>::factor() {
00496
00497 if (factored()) return(0);
00498
00499 TEST_FOR_EXCEPTION(inverted(), std::logic_error,
00500 "SerialSpdDenseSolver<T>::factor: Cannot factor an inverted matrix!");
00501
00502 ANORM_ = Matrix_->normOne();
00503
00504
00505
00506
00507
00508 if (A_ == AF_)
00509 if (refineSolution_ ) {
00510 Factor_ = rcp( new SerialSymDenseMatrix<OrdinalType,ScalarType>(*Matrix_) );
00511 AF_ = Factor_->values();
00512 LDAF_ = Factor_->stride();
00513 }
00514
00515 int ierr = 0;
00516 if (equilibrate_) ierr = equilibrateMatrix();
00517
00518 if (ierr!=0) return(ierr);
00519
00520 INFO_ = 0;
00521 this->POTRF(Matrix_->UPLO(), numRowCols_, AF_, LDAF_, &INFO_);
00522 factored_ = true;
00523
00524 return(INFO_);
00525 }
00526
00527
00528
00529 template<typename OrdinalType, typename ScalarType>
00530 int SerialSpdDenseSolver<OrdinalType,ScalarType>::solve() {
00531
00532
00533
00534
00535
00536
00537
00538
00539
00540
00541 int ierr = 0;
00542 if (equilibrate_) {
00543 ierr = equilibrateRHS();
00544 equilibratedB_ = true;
00545 }
00546 if (ierr != 0) return(ierr);
00547
00548 TEST_FOR_EXCEPTION( (equilibratedA_ && !equilibratedB_) || (!equilibratedA_ && equilibratedB_) ,
00549 std::logic_error, "SerialSpdDenseSolver<T>::solve: Matrix and vectors must be similarly scaled!");
00550 TEST_FOR_EXCEPTION( RHS_==Teuchos::null, std::invalid_argument,
00551 "SerialSpdDenseSolver<T>::solve: No right-hand side vector (RHS) has been set for the linear system!");
00552 TEST_FOR_EXCEPTION( LHS_==Teuchos::null, std::invalid_argument,
00553 "SerialSpdDenseSolver<T>::solve: No solution vector (LHS) has been set for the linear system!");
00554
00555 if (shouldEquilibrate() && !equilibratedA_)
00556 std::cout << "WARNING! SerialSpdDenseSolver<T>::solve: System should be equilibrated!" << std::endl;
00557
00558 if (inverted()) {
00559
00560 TEST_FOR_EXCEPTION( RHS_->values() == LHS_->values(), std::invalid_argument,
00561 "SerialSpdDenseSolver<T>::solve: X and B must be different vectors if matrix is inverted.");
00562
00563 INFO_ = 0;
00564 this->GEMM(Teuchos::NO_TRANS, Teuchos::NO_TRANS, numRowCols_, RHS_->numCols(),
00565 numRowCols_, 1.0, AF_, LDAF_, RHS_->values(), RHS_->stride(), 0.0,
00566 LHS_->values(), LHS_->stride());
00567 if (INFO_!=0) return(INFO_);
00568 solved_ = true;
00569 }
00570 else {
00571
00572 if (!factored()) factor();
00573
00574 if (RHS_->values()!=LHS_->values()) {
00575 (*LHS_) = (*RHS_);
00576 }
00577 INFO_ = 0;
00578 this->POTRS(Matrix_->UPLO(), numRowCols_, RHS_->numCols(), AF_, LDAF_, LHS_->values(), LHS_->stride(), &INFO_);
00579 if (INFO_!=0) return(INFO_);
00580 solved_ = true;
00581
00582 }
00583 int ierr1=0;
00584 if (refineSolution_ && !inverted()) ierr1 = applyRefinement();
00585 if (ierr1!=0)
00586 return(ierr1);
00587 else
00588 return(ierr);
00589
00590 if (equilibrate_) ierr1 = unequilibrateLHS();
00591 return(ierr1);
00592 }
00593
00594
00595 template<typename OrdinalType, typename ScalarType>
00596 int SerialSpdDenseSolver<OrdinalType,ScalarType>::applyRefinement()
00597 {
00598 TEST_FOR_EXCEPTION(!solved(), std::logic_error,
00599 "SerialSpdDenseSolver<T>::applyRefinement: Must have an existing solution!");
00600 TEST_FOR_EXCEPTION(A_==AF_, std::logic_error,
00601 "SerialSpdDenseSolver<T>::applyRefinement: Cannot apply refinement if no original copy of A!");
00602
00603 OrdinalType NRHS = RHS_->numCols();
00604 FERR_.resize( NRHS );
00605 BERR_.resize( NRHS );
00606 allocateWORK();
00607 allocateIWORK();
00608
00609 INFO_ = 0;
00610 this->PORFS(Matrix_->UPLO(), numRowCols_, NRHS, A_, LDA_, AF_, LDAF_,
00611 RHS_->values(), RHS_->stride(), LHS_->values(), LHS_->stride(),
00612 &FERR_[0], &BERR_[0], &WORK_[0], &IWORK_[0], &INFO_);
00613
00614 solutionErrorsEstimated_ = true;
00615 reciprocalConditionEstimated_ = true;
00616 solutionRefined_ = true;
00617
00618 return(INFO_);
00619
00620 }
00621
00622
00623
00624 template<typename OrdinalType, typename ScalarType>
00625 int SerialSpdDenseSolver<OrdinalType,ScalarType>::computeEquilibrateScaling()
00626 {
00627 if (R_.size()!=0) return(0);
00628
00629 R_.resize( numRowCols_ );
00630
00631 INFO_ = 0;
00632 this->POEQU (numRowCols_, AF_, LDAF_, &R_[0], &SCOND_, &AMAX_, &INFO_);
00633 if (SCOND_<0.1*ScalarTraits<MagnitudeType>::one() ||
00634 AMAX_ < ScalarTraits<ScalarType>::rmin() ||
00635 AMAX_ > ScalarTraits<ScalarType>::rmax())
00636 shouldEquilibrate_ = true;
00637
00638 return(INFO_);
00639 }
00640
00641
00642
00643 template<typename OrdinalType, typename ScalarType>
00644 int SerialSpdDenseSolver<OrdinalType,ScalarType>::equilibrateMatrix()
00645 {
00646 OrdinalType i, j;
00647 int ierr = 0;
00648
00649 if (equilibratedA_) return(0);
00650 if (R_.size()==0) ierr = computeEquilibrateScaling();
00651 if (ierr!=0) return(ierr);
00652 if (Matrix_->upper()) {
00653 if (A_==AF_) {
00654 ScalarType * ptr;
00655 for (j=0; j<numRowCols_; j++) {
00656 ptr = A_ + j*LDA_;
00657 ScalarType s1 = R_[j];
00658 for (i=0; i<=j; i++) {
00659 *ptr = *ptr*s1*R_[i];
00660 ptr++;
00661 }
00662 }
00663 }
00664 else {
00665 ScalarType * ptr;
00666 ScalarType * ptr1;
00667 for (j=0; j<numRowCols_; j++) {
00668 ptr = A_ + j*LDA_;
00669 ptr1 = AF_ + j*LDAF_;
00670 ScalarType s1 = R_[j];
00671 for (i=0; i<=j; i++) {
00672 *ptr = *ptr*s1*R_[i];
00673 ptr++;
00674 *ptr1 = *ptr1*s1*R_[i];
00675 ptr1++;
00676 }
00677 }
00678 }
00679 }
00680 else {
00681 if (A_==AF_) {
00682 ScalarType * ptr;
00683 for (j=0; j<numRowCols_; j++) {
00684 ptr = A_ + j + j*LDA_;
00685 ScalarType s1 = R_[j];
00686 for (i=j; i<numRowCols_; i++) {
00687 *ptr = *ptr*s1*R_[i];
00688 ptr++;
00689 }
00690 }
00691 }
00692 else {
00693 ScalarType * ptr;
00694 ScalarType * ptr1;
00695 for (j=0; j<numRowCols_; j++) {
00696 ptr = A_ + j + j*LDA_;
00697 ptr1 = AF_ + j + j*LDAF_;
00698 ScalarType s1 = R_[j];
00699 for (i=j; i<numRowCols_; i++) {
00700 *ptr = *ptr*s1*R_[i];
00701 ptr++;
00702 *ptr1 = *ptr1*s1*R_[i];
00703 ptr1++;
00704 }
00705 }
00706 }
00707 }
00708
00709 equilibratedA_ = true;
00710
00711 return(0);
00712 }
00713
00714
00715
00716 template<typename OrdinalType, typename ScalarType>
00717 int SerialSpdDenseSolver<OrdinalType,ScalarType>::equilibrateRHS()
00718 {
00719 OrdinalType i, j;
00720 int ierr = 0;
00721
00722 if (equilibratedB_) return(0);
00723 if (R_.size()==0) ierr = computeEquilibrateScaling();
00724 if (ierr!=0) return(ierr);
00725
00726 OrdinalType LDB = RHS_->stride(), NRHS = RHS_->numCols();
00727 ScalarType * B = RHS_->values();
00728 ScalarType * ptr;
00729 for (j=0; j<NRHS; j++) {
00730 ptr = B + j*LDB;
00731 for (i=0; i<numRowCols_; i++) {
00732 *ptr = *ptr*R_[i];
00733 ptr++;
00734 }
00735 }
00736
00737 equilibratedB_ = true;
00738
00739 return(0);
00740 }
00741
00742
00743
00744 template<typename OrdinalType, typename ScalarType>
00745 int SerialSpdDenseSolver<OrdinalType,ScalarType>::unequilibrateLHS()
00746 {
00747 OrdinalType i, j;
00748
00749 if (!equilibratedB_) return(0);
00750
00751 OrdinalType LDX = RHS_->stride(), NRHS = RHS_->numCols();
00752 ScalarType * X = RHS_->values();
00753 ScalarType * ptr;
00754 for (j=0; j<NRHS; j++) {
00755 ptr = X + j*LDX;
00756 for (i=0; i<numRowCols_; i++) {
00757 *ptr = *ptr*R_[i];
00758 ptr++;
00759 }
00760 }
00761
00762 return(0);
00763 }
00764
00765
00766
00767 template<typename OrdinalType, typename ScalarType>
00768 int SerialSpdDenseSolver<OrdinalType,ScalarType>::invert()
00769 {
00770 if (!factored()) factor();
00771
00772 INFO_ = 0;
00773 this->POTRI( Matrix_->UPLO(), numRowCols_, AF_, LDAF_, &INFO_);
00774
00775
00776 if (Matrix_->upper())
00777 Matrix_->setLower();
00778 else
00779 Matrix_->setUpper();
00780
00781 inverted_ = true;
00782 factored_ = false;
00783
00784 return(INFO_);
00785 }
00786
00787
00788
00789 template<typename OrdinalType, typename ScalarType>
00790 int SerialSpdDenseSolver<OrdinalType,ScalarType>::reciprocalConditionEstimate(MagnitudeType & Value)
00791 {
00792 if (reciprocalConditionEstimated()) {
00793 Value = RCOND_;
00794 return(0);
00795 }
00796
00797 if ( ANORM_<ScalarTraits<MagnitudeType>::zero() ) ANORM_ = Matrix_->normOne();
00798
00799 int ierr = 0;
00800 if (!factored()) ierr = factor();
00801 if (ierr!=0) return(ierr);
00802
00803 allocateWORK();
00804 allocateIWORK();
00805
00806
00807 INFO_ = 0;
00808 this->POCON( Matrix_->UPLO(), numRowCols_, AF_, LDAF_, ANORM_, &RCOND_, &WORK_[0], &IWORK_[0], &INFO_);
00809 reciprocalConditionEstimated_ = true;
00810 Value = RCOND_;
00811
00812 return(INFO_);
00813 }
00814
00815
00816 template<typename OrdinalType, typename ScalarType>
00817 void SerialSpdDenseSolver<OrdinalType,ScalarType>::Print(std::ostream& os) const {
00818
00819 if (Matrix_!=Teuchos::null) os << "Solver Matrix" << std::endl << *Matrix_ << std::endl;
00820 if (Factor_!=Teuchos::null) os << "Solver Factored Matrix" << std::endl << *Factor_ << std::endl;
00821 if (LHS_ !=Teuchos::null) os << "Solver LHS" << std::endl << *LHS_ << std::endl;
00822 if (RHS_ !=Teuchos::null) os << "Solver RHS" << std::endl << *RHS_ << std::endl;
00823
00824 }
00825
00826 }
00827
00828 #endif