Teuchos - Trilinos Tools Package  Version of the Day
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
Teuchos_SerialBandDenseSolver.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_SERIALBANDDENSESOLVER_HPP_
43 #define _TEUCHOS_SERIALBANDDENSESOLVER_HPP_
44 
49 #include "Teuchos_CompObject.hpp"
50 #include "Teuchos_BLAS.hpp"
51 #include "Teuchos_LAPACK.hpp"
52 #include "Teuchos_RCP.hpp"
53 #include "Teuchos_ConfigDefs.hpp"
57 #include "Teuchos_ScalarTraits.hpp"
58 
59 #ifdef HAVE_TEUCHOSNUMERICS_EIGEN
60 #include "Eigen/Dense"
61 #endif
62 
63 namespace Teuchos {
64 
165  template<typename OrdinalType, typename ScalarType>
166  class SerialBandDenseSolver : public CompObject, public Object, public BLAS<OrdinalType, ScalarType>,
167  public LAPACK<OrdinalType, ScalarType>
168  {
169 
170  public:
171 
172  typedef typename Teuchos::ScalarTraits<ScalarType>::magnitudeType MagnitudeType;
173 #ifdef HAVE_TEUCHOSNUMERICS_EIGEN
174  typedef typename Eigen::Matrix<ScalarType,Eigen::Dynamic,Eigen::Dynamic,Eigen::ColMajor> EigenMatrix;
175  typedef typename Eigen::internal::BandMatrix<ScalarType,Eigen::Dynamic,Eigen::Dynamic,Eigen::ColMajor> EigenBandMatrix;
176  typedef typename Eigen::Matrix<ScalarType,Eigen::Dynamic,1> EigenVector;
177  typedef typename Eigen::InnerStride<Eigen::Dynamic> EigenInnerStride;
178  typedef typename Eigen::OuterStride<Eigen::Dynamic> EigenOuterStride;
179  typedef typename Eigen::Map<EigenVector,0,EigenInnerStride > EigenVectorMap;
180  typedef typename Eigen::Map<const EigenVector,0,EigenInnerStride > EigenConstVectorMap;
181  typedef typename Eigen::Map<EigenMatrix,0,EigenOuterStride > EigenMatrixMap;
182  typedef typename Eigen::Map<const EigenMatrix,0,EigenOuterStride > EigenConstMatrixMap;
183  typedef typename Eigen::PermutationMatrix<Eigen::Dynamic,Eigen::Dynamic,OrdinalType> EigenPermutationMatrix;
184  typedef typename Eigen::Array<OrdinalType,Eigen::Dynamic,1> EigenOrdinalArray;
185  typedef typename Eigen::Map<EigenOrdinalArray> EigenOrdinalArrayMap;
186  typedef typename Eigen::Array<ScalarType,Eigen::Dynamic,1> EigenScalarArray;
187  typedef typename Eigen::Map<EigenScalarArray> EigenScalarArrayMap;
188 #endif
189 
191 
192 
195 
197  virtual ~SerialBandDenseSolver();
198 
200 
202 
205 
207 
212 
214 
216 
219  void factorWithEquilibration(bool flag) {equilibrate_ = flag; return;}
220 
224  void solveWithTranspose(bool flag) {transpose_ = flag; if (flag) TRANS_ = Teuchos::TRANS; else TRANS_ = Teuchos::NO_TRANS; return;}
225 
228  void solveWithTransposeFlag(Teuchos::ETransp trans) {TRANS_ = trans; transpose_ = (trans != Teuchos::NO_TRANS) ? true : false;}
229 
231  void solveToRefinedSolution(bool flag) {refineSolution_ = flag; return;}
232 
234 
237  void estimateSolutionErrors(bool flag);
239 
241 
242 
244 
247  int factor();
248 
250 
253  int solve();
254 
256 
260 
262 
265  int equilibrateMatrix();
266 
268 
271  int equilibrateRHS();
272 
274 
277  int applyRefinement();
278 
280 
283  int unequilibrateLHS();
284 
286 
292  int reciprocalConditionEstimate(MagnitudeType & Value);
294 
296 
297 
299  bool transpose() {return(transpose_);}
300 
302  bool factored() {return(factored_);}
303 
305  bool equilibratedA() {return(equilibratedA_);}
306 
308  bool equilibratedB() {return(equilibratedB_);}
309 
311  bool shouldEquilibrate() {computeEquilibrateScaling(); return(shouldEquilibrate_);}
312 
314  bool solutionErrorsEstimated() {return(solutionErrorsEstimated_);}
315 
317  bool reciprocalConditionEstimated() {return(reciprocalConditionEstimated_);}
318 
320  bool solved() {return(solved_);}
321 
323  bool solutionRefined() {return(solutionRefined_);}
325 
327 
328 
331 
334 
337 
340 
342  OrdinalType numRows() const {return(M_);}
343 
345  OrdinalType numCols() const {return(N_);}
346 
348  OrdinalType numLower() const {return(KL_);}
349 
351  OrdinalType numUpper() const {return(KU_);}
352 
354  std::vector<OrdinalType> IPIV() const {return(IPIV_);}
355 
357  MagnitudeType ANORM() const {return(ANORM_);}
358 
360  MagnitudeType RCOND() const {return(RCOND_);}
361 
363 
365  MagnitudeType ROWCND() const {return(ROWCND_);}
366 
368 
370  MagnitudeType COLCND() const {return(COLCND_);}
371 
373  MagnitudeType AMAX() const {return(AMAX_);}
374 
376  std::vector<MagnitudeType> FERR() const {return(FERR_);}
377 
379  std::vector<MagnitudeType> BERR() const {return(BERR_);}
380 
382  std::vector<MagnitudeType> R() const {return(R_);}
383 
385  std::vector<MagnitudeType> C() const {return(C_);}
387 
389 
390  void Print(std::ostream& os) const;
393  protected:
394 
395  void allocateWORK() { LWORK_ = 3*N_; WORK_.resize( LWORK_ ); return;}
396  void resetMatrix();
397  void resetVectors();
398 
399 
400  bool equilibrate_;
401  bool shouldEquilibrate_;
402  bool equilibratedA_;
403  bool equilibratedB_;
404  bool transpose_;
405  bool factored_;
406  bool estimateSolutionErrors_;
407  bool solutionErrorsEstimated_;
408  bool solved_;
409  bool reciprocalConditionEstimated_;
410  bool refineSolution_;
411  bool solutionRefined_;
412 
413  Teuchos::ETransp TRANS_;
414 
415  OrdinalType M_;
416  OrdinalType N_;
417  OrdinalType KL_;
418  OrdinalType KU_;
419  OrdinalType Min_MN_;
420  OrdinalType LDA_;
421  OrdinalType LDAF_;
422  OrdinalType INFO_;
423  OrdinalType LWORK_;
424 
425  std::vector<OrdinalType> IPIV_;
426  std::vector<int> IWORK_;
427 
428  MagnitudeType ANORM_;
429  MagnitudeType RCOND_;
430  MagnitudeType ROWCND_;
431  MagnitudeType COLCND_;
432  MagnitudeType AMAX_;
433 
434  RCP<SerialBandDenseMatrix<OrdinalType, ScalarType> > Matrix_;
435  RCP<SerialDenseMatrix<OrdinalType, ScalarType> > LHS_;
436  RCP<SerialDenseMatrix<OrdinalType, ScalarType> > RHS_;
437  RCP<SerialBandDenseMatrix<OrdinalType, ScalarType> > Factor_;
438 
439  ScalarType * A_;
440  ScalarType * AF_;
441  std::vector<MagnitudeType> FERR_;
442  std::vector<MagnitudeType> BERR_;
443  std::vector<ScalarType> WORK_;
444  std::vector<MagnitudeType> R_;
445  std::vector<MagnitudeType> C_;
446 #ifdef HAVE_TEUCHOSNUMERICS_EIGEN
447  Eigen::PartialPivLU<EigenMatrix> lu_;
448 #endif
449 
450  private:
451 
452  // SerialBandDenseSolver copy constructor (put here because we don't want user access)
453  SerialBandDenseSolver(const SerialBandDenseSolver<OrdinalType, ScalarType>& Source);
454  SerialBandDenseSolver & operator=(const SerialBandDenseSolver<OrdinalType, ScalarType>& Source);
455 
456  };
457 
458  // Helper traits to distinguish work arrays for real and complex-valued datatypes.
459  using namespace details;
460 
461 //=============================================================================
462 
463 template<typename OrdinalType, typename ScalarType>
465  : CompObject(),
466  equilibrate_(false),
467  shouldEquilibrate_(false),
468  equilibratedA_(false),
469  equilibratedB_(false),
470  transpose_(false),
471  factored_(false),
472  estimateSolutionErrors_(false),
473  solutionErrorsEstimated_(false),
474  solved_(false),
475  reciprocalConditionEstimated_(false),
476  refineSolution_(false),
477  solutionRefined_(false),
478  TRANS_(Teuchos::NO_TRANS),
479  M_(0),
480  N_(0),
481  KL_(0),
482  KU_(0),
483  Min_MN_(0),
484  LDA_(0),
485  LDAF_(0),
486  INFO_(0),
487  LWORK_(0),
488  ANORM_(0.0),
489  RCOND_(ScalarTraits<MagnitudeType>::zero()),
490  ROWCND_(ScalarTraits<MagnitudeType>::zero()),
491  COLCND_(ScalarTraits<MagnitudeType>::zero()),
492  AMAX_(ScalarTraits<MagnitudeType>::zero()),
493  A_(NULL),
494  AF_(NULL)
495 {
496  resetMatrix();
497 }
498 //=============================================================================
499 
500 template<typename OrdinalType, typename ScalarType>
502 {}
503 
504 //=============================================================================
505 
506 template<typename OrdinalType, typename ScalarType>
508 {
509  LHS_ = Teuchos::null;
510  RHS_ = Teuchos::null;
511  reciprocalConditionEstimated_ = false;
512  solutionRefined_ = false;
513  solved_ = false;
514  solutionErrorsEstimated_ = false;
515  equilibratedB_ = false;
516 }
517 //=============================================================================
518 
519 template<typename OrdinalType, typename ScalarType>
520 void SerialBandDenseSolver<OrdinalType,ScalarType>::resetMatrix()
521 {
522  resetVectors();
523  equilibratedA_ = false;
524  factored_ = false;
525  M_ = 0;
526  N_ = 0;
527  KL_ = 0;
528  KU_ = 0;
529  Min_MN_ = 0;
530  LDA_ = 0;
531  LDAF_ = 0;
536  A_ = 0;
537  AF_ = 0;
538  INFO_ = 0;
539  LWORK_ = 0;
540  R_.resize(0);
541  C_.resize(0);
542 }
543 //=============================================================================
544 
545 template<typename OrdinalType, typename ScalarType>
547 {
548 
549  OrdinalType m = AB->numRows();
550  OrdinalType n = AB->numCols();
551  OrdinalType kl = AB->lowerBandwidth();
552  OrdinalType ku = AB->upperBandwidth();
553 
554  // Check that the new matrix is consistent.
555  TEUCHOS_TEST_FOR_EXCEPTION(AB->values()==0, std::invalid_argument,
556  "SerialBandDenseSolver<T>::setMatrix: A is an empty SerialBandDenseMatrix<T>!");
557 
558  resetMatrix();
559  Matrix_ = AB;
560  Factor_ = AB;
561  M_ = m;
562  N_ = n;
563  KL_ = kl;
564  KU_ = ku-kl;
565  Min_MN_ = TEUCHOS_MIN(M_,N_);
566  LDA_ = AB->stride();
567  LDAF_ = LDA_;
568  A_ = AB->values();
569  AF_ = AB->values();
570 
571  return(0);
572 }
573 //=============================================================================
574 
575 template<typename OrdinalType, typename ScalarType>
578 {
579  // Check that these new vectors are consistent.
580  TEUCHOS_TEST_FOR_EXCEPTION(B->numRows()!=X->numRows() || B->numCols() != X->numCols(), std::invalid_argument,
581  "SerialBandDenseSolver<T>::setVectors: X and B are not the same size!");
582  TEUCHOS_TEST_FOR_EXCEPTION(B->values()==0, std::invalid_argument,
583  "SerialBandDenseSolver<T>::setVectors: B is an empty SerialDenseMatrix<T>!");
584  TEUCHOS_TEST_FOR_EXCEPTION(X->values()==0, std::invalid_argument,
585  "SerialBandDenseSolver<T>::setVectors: X is an empty SerialDenseMatrix<T>!");
586  TEUCHOS_TEST_FOR_EXCEPTION(B->stride()<1, std::invalid_argument,
587  "SerialBandDenseSolver<T>::setVectors: B has an invalid stride!");
588  TEUCHOS_TEST_FOR_EXCEPTION(X->stride()<1, std::invalid_argument,
589  "SerialBandDenseSolver<T>::setVectors: X has an invalid stride!");
590 
591  resetVectors();
592  LHS_ = X;
593  RHS_ = B;
594  return(0);
595 }
596 //=============================================================================
597 
598 template<typename OrdinalType, typename ScalarType>
600 {
601  estimateSolutionErrors_ = flag;
602 
603  // If the errors are estimated, this implies that the solution must be refined
604  refineSolution_ = refineSolution_ || flag;
605 }
606 //=============================================================================
607 
608 template<typename OrdinalType, typename ScalarType>
610 
611  if (factored()) return(0); // Already factored
612 
613  ANORM_ = Matrix_->normOne(); // Compute 1-Norm of A
614 
615  // If we want to refine the solution, then the factor must
616  // be stored separatedly from the original matrix
617 
618  if (A_ == AF_)
619  if (refineSolution_ ) {
620  Factor_ = rcp( new SerialBandDenseMatrix<OrdinalType,ScalarType>(*Matrix_) );
621  AF_ = Factor_->values();
622  LDAF_ = Factor_->stride();
623  }
624 
625  int ierr = 0;
626  if (equilibrate_) ierr = equilibrateMatrix();
627 
628  if (ierr!=0) return(ierr);
629 
630  if (IPIV_.size()==0) IPIV_.resize( N_ ); // Allocated Pivot vector if not already done.
631 
632  INFO_ = 0;
633 
634 #ifdef HAVE_TEUCHOSNUMERICS_EIGEN
635  EigenMatrixMap aMap( AF_+KL_, KL_+KU_+1, N_, EigenOuterStride(LDAF_) );
636  EigenMatrix fullMatrix(M_, N_);
637  for (OrdinalType j=0; j<N_; j++) {
638  for (OrdinalType i=TEUCHOS_MAX(0,j-KU_); i<=TEUCHOS_MIN(M_-1, j+KL_); i++) {
639  fullMatrix(i,j) = aMap(KU_+i-j, j);
640  }
641  }
642 
643  EigenPermutationMatrix p;
644  EigenOrdinalArray ind;
645  EigenOrdinalArrayMap ipivMap( &IPIV_[0], Min_MN_ );
646 
647  lu_.compute(fullMatrix);
648  p = lu_.permutationP();
649  ind = p.indices();
650 
651  for (OrdinalType i=0; i<ipivMap.innerSize(); i++) {
652  ipivMap(i) = ind(i);
653  }
654 
655 #else
656  this->GBTRF(M_, N_, KL_, KU_, AF_, LDAF_, &IPIV_[0], &INFO_);
657 #endif
658 
659  factored_ = true;
660 
661  return(INFO_);
662 }
663 
664 //=============================================================================
665 
666 template<typename OrdinalType, typename ScalarType>
668 
669  // If the user want the matrix to be equilibrated or wants a refined solution, we will
670  // call the X interface.
671  // Otherwise, if the matrix is already factored we will call the TRS interface.
672  // Otherwise, if the matrix is unfactored we will call the SV interface.
673 
674  int ierr = 0;
675  if (equilibrate_) {
676  ierr = equilibrateRHS();
677  }
678  if (ierr != 0) return(ierr); // Can't equilibrate B, so return.
679 
680  TEUCHOS_TEST_FOR_EXCEPTION( RHS_==Teuchos::null, std::invalid_argument,
681  "SerialBandDenseSolver<T>::solve: No right-hand side vector (RHS) has been set for the linear system!");
682  TEUCHOS_TEST_FOR_EXCEPTION( LHS_==Teuchos::null, std::invalid_argument,
683  "SerialBandDenseSolver<T>::solve: No solution vector (LHS) has been set for the linear system!");
684 
685  if (!factored()) factor(); // Matrix must be factored
686 
687  TEUCHOS_TEST_FOR_EXCEPTION( (equilibratedA_ && !equilibratedB_) || (!equilibratedA_ && equilibratedB_) ,
688  std::logic_error, "SerialBandDenseSolver<T>::solve: Matrix and vectors must be similarly scaled!");
689 
690  if (shouldEquilibrate() && !equilibratedA_)
691  std::cout << "WARNING! SerialBandDenseSolver<T>::solve: System should be equilibrated!" << std::endl;
692 
693  if (RHS_->values()!=LHS_->values()) {
694  (*LHS_) = (*RHS_); // Copy B to X if needed
695  }
696  INFO_ = 0;
697 
698 #ifdef HAVE_TEUCHOSNUMERICS_EIGEN
699  EigenMatrixMap rhsMap( RHS_->values(), RHS_->numRows(), RHS_->numCols(), EigenOuterStride(RHS_->stride()) );
700  EigenMatrixMap lhsMap( LHS_->values(), LHS_->numRows(), LHS_->numCols(), EigenOuterStride(LHS_->stride()) );
701  if ( TRANS_ == Teuchos::NO_TRANS ) {
702  lhsMap = lu_.solve(rhsMap);
703  } else if ( TRANS_ == Teuchos::TRANS ) {
704  lu_.matrixLU().template triangularView<Eigen::Upper>().transpose().solveInPlace(lhsMap);
705  lu_.matrixLU().template triangularView<Eigen::UnitLower>().transpose().solveInPlace(lhsMap);
706  EigenMatrix x = lu_.permutationP().transpose() * lhsMap;
707  lhsMap = x;
708  } else {
709  lu_.matrixLU().template triangularView<Eigen::Upper>().adjoint().solveInPlace(lhsMap);
710  lu_.matrixLU().template triangularView<Eigen::UnitLower>().adjoint().solveInPlace(lhsMap);
711  EigenMatrix x = lu_.permutationP().transpose() * lhsMap;
712  lhsMap = x;
713  }
714 #else
715  this->GBTRS(ETranspChar[TRANS_], N_, KL_, KU_, RHS_->numCols(), AF_, LDAF_, &IPIV_[0], LHS_->values(), LHS_->stride(), &INFO_);
716 #endif
717 
718  if (INFO_!=0) return(INFO_);
719  solved_ = true;
720 
721  int ierr1=0;
722  if (refineSolution_) ierr1 = applyRefinement();
723  if (ierr1!=0)
724  return(ierr1);
725 
726  if (equilibrate_) ierr1 = unequilibrateLHS();
727 
728  return(ierr1);
729 }
730 //=============================================================================
731 
732 template<typename OrdinalType, typename ScalarType>
734 {
735  TEUCHOS_TEST_FOR_EXCEPTION(!solved(), std::logic_error,
736  "SerialBandDenseSolver<T>::applyRefinement: Must have an existing solution!");
737  TEUCHOS_TEST_FOR_EXCEPTION(A_==AF_, std::logic_error,
738  "SerialBandDenseSolver<T>::applyRefinement: Cannot apply refinement if no original copy of A!");
739 
740 #ifdef HAVE_TEUCHOSNUMERICS_EIGEN
741  // Implement templated GERFS or use Eigen.
742  return (-1);
743 #else
744 
745  OrdinalType NRHS = RHS_->numCols();
746  FERR_.resize( NRHS );
747  BERR_.resize( NRHS );
748  allocateWORK();
749 
750  INFO_ = 0;
751  std::vector<typename details::lapack_traits<ScalarType>::iwork_type> GBRFS_WORK( N_ );
752  this->GBRFS(ETranspChar[TRANS_], N_, KL_, KU_, NRHS, A_+KL_, LDA_, AF_, LDAF_, &IPIV_[0],
753  RHS_->values(), RHS_->stride(), LHS_->values(), LHS_->stride(),
754  &FERR_[0], &BERR_[0], &WORK_[0], &GBRFS_WORK[0], &INFO_);
755 
756  solutionErrorsEstimated_ = true;
757  reciprocalConditionEstimated_ = true;
758  solutionRefined_ = true;
759 
760  return(INFO_);
761 #endif
762 }
763 
764 //=============================================================================
765 
766 template<typename OrdinalType, typename ScalarType>
768 {
769  if (R_.size()!=0) return(0); // Already computed
770 
771  R_.resize( M_ );
772  C_.resize( N_ );
773 
774  INFO_ = 0;
775  this->GBEQU (M_, N_, KL_, KU_, AF_+KL_, LDAF_, &R_[0], &C_[0], &ROWCND_, &COLCND_, &AMAX_, &INFO_);
776 
777  if (COLCND_<0.1*ScalarTraits<MagnitudeType>::one() ||
778  ROWCND_<0.1*ScalarTraits<MagnitudeType>::one() ||
780  shouldEquilibrate_ = true;
781 
782  return(INFO_);
783 }
784 
785 //=============================================================================
786 
787 template<typename OrdinalType, typename ScalarType>
789 {
790  OrdinalType i, j;
791  int ierr = 0;
792 
793  if (equilibratedA_) return(0); // Already done.
794  if (R_.size()==0) ierr = computeEquilibrateScaling(); // Compute R and C if needed.
795  if (ierr!=0) return(ierr); // If return value is not zero, then can't equilibrate.
796 
797  if (A_==AF_) {
798 
799  ScalarType * ptr;
800  for (j=0; j<N_; j++) {
801  ptr = A_ + KL_ + j*LDA_ + TEUCHOS_MAX(KU_-j, 0);
802  ScalarType s1 = C_[j];
803  for (i=TEUCHOS_MAX(0,j-KU_); i<=TEUCHOS_MIN(M_-1,j+KL_); i++) {
804  *ptr = *ptr*s1*R_[i];
805  ptr++;
806  }
807  }
808  } else {
809 
810  ScalarType * ptr;
811  ScalarType * ptrL;
812  ScalarType * ptrU;
813  for (j=0; j<N_; j++) {
814  ptr = A_ + KL_ + j*LDA_ + TEUCHOS_MAX(KU_-j, 0);
815  ScalarType s1 = C_[j];
816  for (i=TEUCHOS_MAX(0,j-KU_); i<=TEUCHOS_MIN(M_-1,j+KL_); i++) {
817  *ptr = *ptr*s1*R_[i];
818  ptr++;
819  }
820  }
821  for (j=0; j<N_; j++) {
822  ptrU = AF_ + j*LDAF_ + TEUCHOS_MAX(KL_+KU_-j, 0);
823  ptrL = AF_ + KL_ + KU_ + 1 + j*LDAF_;
824  ScalarType s1 = C_[j];
825  for (i=TEUCHOS_MAX(0,j-(KL_+KU_)); i<=TEUCHOS_MIN(M_-1,j); i++) {
826  *ptrU = *ptrU*s1*R_[i];
827  ptrU++;
828  }
829  for (i=TEUCHOS_MAX(0,j); i<=TEUCHOS_MIN(M_-1,j+KL_); i++) {
830  *ptrL = *ptrL*s1*R_[i];
831  ptrL++;
832  }
833  }
834  }
835 
836  equilibratedA_ = true;
837 
838  return(0);
839 }
840 
841 //=============================================================================
842 
843 template<typename OrdinalType, typename ScalarType>
845 {
846  OrdinalType i, j;
847  int ierr = 0;
848 
849  if (equilibratedB_) return(0); // Already done.
850  if (R_.size()==0) ierr = computeEquilibrateScaling(); // Compute R and C if needed.
851  if (ierr!=0) return(ierr); // Can't count on R and C being computed.
852 
853  MagnitudeType * R_tmp = &R_[0];
854  if (transpose_) R_tmp = &C_[0];
855 
856  OrdinalType LDB = RHS_->stride(), NRHS = RHS_->numCols();
857  ScalarType * B = RHS_->values();
858  ScalarType * ptr;
859  for (j=0; j<NRHS; j++) {
860  ptr = B + j*LDB;
861  for (i=0; i<M_; i++) {
862  *ptr = *ptr*R_tmp[i];
863  ptr++;
864  }
865  }
866 
867  equilibratedB_ = true;
868 
869  return(0);
870 }
871 
872 
873 //=============================================================================
874 
875 template<typename OrdinalType, typename ScalarType>
877 {
878  OrdinalType i, j;
879 
880  if (!equilibratedB_) return(0); // Nothing to do
881 
882  MagnitudeType * C_tmp = &C_[0];
883  if (transpose_) C_tmp = &R_[0];
884 
885  OrdinalType LDX = LHS_->stride(), NLHS = LHS_->numCols();
886  ScalarType * X = LHS_->values();
887  ScalarType * ptr;
888  for (j=0; j<NLHS; j++) {
889  ptr = X + j*LDX;
890  for (i=0; i<N_; i++) {
891  *ptr = *ptr*C_tmp[i];
892  ptr++;
893  }
894  }
895 
896  return(0);
897 }
898 
899 //=============================================================================
900 
901 template<typename OrdinalType, typename ScalarType>
903 {
904 #ifdef HAVE_TEUCHOSNUMERICS_EIGEN
905  // Implement templated GECON or use Eigen. Eigen currently doesn't have a condition estimation function.
906  return (-1);
907 #else
908 
909  if (reciprocalConditionEstimated()) {
910  Value = RCOND_;
911  return(0); // Already computed, just return it.
912  }
913 
914  if ( ANORM_<ScalarTraits<MagnitudeType>::zero() ) ANORM_ = Matrix_->normOne();
915 
916  int ierr = 0;
917  if (!factored()) ierr = factor(); // Need matrix factored.
918  if (ierr!=0) return(ierr);
919 
920  allocateWORK();
921 
922  // We will assume a one-norm condition number
923  INFO_ = 0;
924  std::vector<typename details::lapack_traits<ScalarType>::iwork_type> GBCON_WORK( N_ );
925  this->GBCON( '1', N_, KL_, KU_, AF_, LDAF_, &IPIV_[0], ANORM_, &RCOND_, &WORK_[0], &GBCON_WORK[0], &INFO_);
926  reciprocalConditionEstimated_ = true;
927  Value = RCOND_;
928 
929  return(INFO_);
930 #endif
931 }
932 //=============================================================================
933 
934 template<typename OrdinalType, typename ScalarType>
936 
937  if (Matrix_!=Teuchos::null) os << "Solver Matrix" << std::endl << *Matrix_ << std::endl;
938  if (Factor_!=Teuchos::null) os << "Solver Factored Matrix" << std::endl << *Factor_ << std::endl;
939  if (LHS_ !=Teuchos::null) os << "Solver LHS" << std::endl << *LHS_ << std::endl;
940  if (RHS_ !=Teuchos::null) os << "Solver RHS" << std::endl << *RHS_ << std::endl;
941 
942 }
943 
944 //=============================================================================
945 
946 
947 //=============================================================================
948 
949 
950 } // namespace Teuchos
951 
952 #endif /* _TEUCHOS_SERIALBANDDENSESOLVER_HPP_ */
bool equilibratedA()
Returns true if factor is equilibrated (factor available via AF() and LDAF()).
MagnitudeType RCOND() const
Returns the reciprocal of the condition number of the this matrix (returns -1 if not yet computed)...
RCP< SerialBandDenseMatrix< 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.
T magnitudeType
Mandatory typedef for result of magnitude.
Templated serial dense matrix class.
int applyRefinement()
Apply Iterative Refinement.
void solveToRefinedSolution(bool flag)
Set whether or not to use iterative refinement to improve solutions to linear systems.
int setMatrix(const RCP< SerialBandDenseMatrix< OrdinalType, ScalarType > > &AB)
Sets the pointer for coefficient matrix.
Templated interface class to BLAS routines.
std::vector< MagnitudeType > R() const
Returns a pointer to the row scaling vector used for equilibration.
bool reciprocalConditionEstimated()
Returns true if the condition number of the this matrix has been computed (value available via Recipr...
MagnitudeType ROWCND() const
Ratio of smallest to largest row scale factors for the this matrix (returns -1 if not yet computed)...
int reciprocalConditionEstimate(MagnitudeType &Value)
Returns the reciprocal of the 1-norm condition number of the this matrix.
#define TEUCHOS_TEST_FOR_EXCEPTION(throw_exception_test, Exception, msg)
Macro for throwing an exception with breakpointing to ease debugging.
Teuchos header file which uses auto-configuration information to include necessary C++ headers...
Templated BLAS wrapper.
std::vector< OrdinalType > IPIV() const
Returns pointer to pivot vector (if factorization has been computed), zero otherwise.
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).
int unequilibrateLHS()
Unscales the solution vectors if equilibration was used to solve the system.
int equilibrateMatrix()
Equilibrates the this matrix.
Object for storing data and providing functionality that is common to all computational classes...
This structure defines some basic traits for a scalar field type.
int computeEquilibrateScaling()
Computes the scaling vector S(i) = 1/sqrt(A(i,i)) of the this matrix.
bool solutionRefined()
Returns true if the current set of vectors has been refined.
Templated class for solving dense linear problems.
Functionality and data that is common to all computational classes.
Templated interface class to LAPACK routines.
Templated serial dense matrix class.
TEUCHOS_DEPRECATED RCP< T > rcp(T *p, Dealloc_T dealloc, bool owns_mem)
Deprecated.
bool solutionErrorsEstimated()
Returns true if forward and backward error estimated have been computed (available via FERR() and BER...
OrdinalType numLower() const
Returns lower bandwidth of system.
The base Teuchos class.
std::vector< MagnitudeType > FERR() const
Returns a pointer to the forward error estimates computed by LAPACK.
This class creates and provides basic support for banded dense matrices of templated type...
bool equilibratedB()
Returns true if RHS is equilibrated (RHS available via B() and LDB()).
A class for representing and solving banded dense linear systems.
RCP< SerialDenseMatrix< OrdinalType, ScalarType > > getLHS() const
Returns pointer to current LHS.
The Templated LAPACK Wrapper Class.
void Print(std::ostream &os) const
Print service methods; defines behavior of ostream &lt;&lt; operator.
bool transpose()
Returns true if transpose of this matrix has and will be used.
OrdinalType numRows() const
Returns row dimension of system.
int solve()
Computes the solution X to AX = B for the this matrix and the B provided to SetVectors()..
MagnitudeType ANORM() const
Returns the 1-Norm of the this matrix (returns -1 if not yet computed).
bool factored()
Returns true if matrix is factored (factor available via AF() and LDAF()).
bool solved()
Returns true if the current set of vectors has been solved.
int factor()
Computes the in-place LU factorization of the matrix using the LAPACK routine _GBTRF.
bool shouldEquilibrate()
Returns true if the LAPACK general rules for equilibration suggest you should equilibrate the system...
MagnitudeType COLCND() const
Ratio of smallest to largest column scale factors for the this matrix (returns -1 if not yet computed...
virtual ~SerialBandDenseSolver()
SerialBandDenseSolver destructor.
SerialBandDenseSolver()
Default constructor; matrix should be set using setMatrix(), LHS and RHS set with setVectors()...
MagnitudeType AMAX() const
Returns the absolute value of the largest entry of the this matrix (returns -1 if not yet computed)...
Defines basic traits for the scalar field type.
void solveWithTransposeFlag(Teuchos::ETransp trans)
Smart reference counting pointer class for automatic garbage collection.
std::vector< MagnitudeType > BERR() const
Returns a pointer to the backward error estimates computed by LAPACK.
OrdinalType numUpper() const
Returns upper bandwidth of system.
std::vector< MagnitudeType > C() const
Returns a pointer to the column scale vector used for equilibration.
RCP< SerialBandDenseMatrix< OrdinalType, ScalarType > > getMatrix() const
Returns pointer to current matrix.
int equilibrateRHS()
Equilibrates the current RHS.
OrdinalType numCols() const
Returns column dimension of system.
Reference-counted pointer class and non-member templated function implementations.
static T one()
Returns representation of one for this scalar type.
RCP< SerialDenseMatrix< OrdinalType, ScalarType > > getRHS() const
Returns pointer to current RHS.
This class creates and provides basic support for dense rectangular matrix of templated type...