Teuchos - Trilinos Tools Package  Version of the Day
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_
44 
48 #include "Teuchos_CompObject.hpp"
49 #include "Teuchos_BLAS.hpp"
50 #include "Teuchos_LAPACK.hpp"
51 #include "Teuchos_RCP.hpp"
52 #include "Teuchos_ConfigDefs.hpp"
54 #include "Teuchos_ScalarTraits.hpp"
55 
56 #ifdef HAVE_TEUCHOSNUMERICS_EIGEN
57 #include "Eigen/Dense"
58 #endif
59 
134 namespace 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 
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; if (trans != Teuchos::NO_TRANS) { transpose_ = true; } }
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  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 
455 template<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 
491 template<typename OrdinalType, typename ScalarType>
493 {}
494 
495 //=============================================================================
496 
497 template<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 
510 template<typename OrdinalType, typename ScalarType>
511 void 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 
536 template<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 
553 template<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 
576 template<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 
586 template<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_ ) {
602  Factor_ = rcp( new SerialDenseMatrix<OrdinalType,ScalarType>(*Matrix_) );
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_) );
618  EigenPermutationMatrix p;
619  EigenOrdinalArray ind;
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 
646 template<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  equilibratedB_ = true;
662  }
663  if (ierr != 0) return(ierr); // Can't equilibrate B, so return.
664 
665  TEUCHOS_TEST_FOR_EXCEPTION( (equilibratedA_ && !equilibratedB_) || (!equilibratedA_ && equilibratedB_) ,
666  std::logic_error, "SerialDenseSolver<T>::solve: Matrix and vectors must be similarly scaled!");
667  TEUCHOS_TEST_FOR_EXCEPTION( RHS_==Teuchos::null, std::invalid_argument,
668  "SerialDenseSolver<T>::solve: No right-hand side vector (RHS) has been set for the linear system!");
669  TEUCHOS_TEST_FOR_EXCEPTION( LHS_==Teuchos::null, std::invalid_argument,
670  "SerialDenseSolver<T>::solve: No solution vector (LHS) has been set for the linear system!");
671 
672  if (shouldEquilibrate() && !equilibratedA_)
673  std::cout << "WARNING! SerialDenseSolver<T>::solve: System should be equilibrated!" << std::endl;
674 
675  if (inverted()) {
676 
677  TEUCHOS_TEST_FOR_EXCEPTION( RHS_->values() == LHS_->values(), std::invalid_argument,
678  "SerialDenseSolver<T>::solve: X and B must be different vectors if matrix is inverted.");
679 
680  INFO_ = 0;
681  this->GEMM(TRANS_, Teuchos::NO_TRANS, N_, RHS_->numCols(), N_, 1.0, AF_, LDAF_,
682  RHS_->values(), RHS_->stride(), 0.0, LHS_->values(), LHS_->stride());
683  if (INFO_!=0) return(INFO_);
684  solved_ = true;
685  }
686  else {
687 
688  if (!factored()) factor(); // Matrix must be factored
689 
690  if (RHS_->values()!=LHS_->values()) {
691  (*LHS_) = (*RHS_); // Copy B to X if needed
692  }
693  INFO_ = 0;
694 
695 #ifdef HAVE_TEUCHOSNUMERICS_EIGEN
696  EigenMatrixMap rhsMap( RHS_->values(), RHS_->numRows(), RHS_->numCols(), EigenOuterStride(RHS_->stride()) );
697  EigenMatrixMap lhsMap( LHS_->values(), LHS_->numRows(), LHS_->numCols(), EigenOuterStride(LHS_->stride()) );
698  if ( TRANS_ == Teuchos::NO_TRANS ) {
699  lhsMap = lu_.solve(rhsMap);
700  } else if ( TRANS_ == Teuchos::TRANS ) {
701  lu_.matrixLU().template triangularView<Eigen::Upper>().transpose().solveInPlace(lhsMap);
702  lu_.matrixLU().template triangularView<Eigen::UnitLower>().transpose().solveInPlace(lhsMap);
703  EigenMatrix x = lu_.permutationP().transpose() * lhsMap;
704  lhsMap = x;
705  } else {
706  lu_.matrixLU().template triangularView<Eigen::Upper>().adjoint().solveInPlace(lhsMap);
707  lu_.matrixLU().template triangularView<Eigen::UnitLower>().adjoint().solveInPlace(lhsMap);
708  EigenMatrix x = lu_.permutationP().transpose() * lhsMap;
709  lhsMap = x;
710  }
711 #else
712  this->GETRS(ETranspChar[TRANS_], N_, RHS_->numCols(), AF_, LDAF_, &IPIV_[0], LHS_->values(), LHS_->stride(), &INFO_);
713 #endif
714 
715  if (INFO_!=0) return(INFO_);
716  solved_ = true;
717 
718  }
719  int ierr1=0;
720  if (refineSolution_ && !inverted()) ierr1 = applyRefinement();
721  if (ierr1!=0)
722  return(ierr1);
723 
724  if (equilibrate_) ierr1 = unequilibrateLHS();
725  return(ierr1);
726 }
727 //=============================================================================
728 
729 template<typename OrdinalType, typename ScalarType>
731 {
732  TEUCHOS_TEST_FOR_EXCEPTION(!solved(), std::logic_error,
733  "SerialDenseSolver<T>::applyRefinement: Must have an existing solution!");
734  TEUCHOS_TEST_FOR_EXCEPTION(A_==AF_, std::logic_error,
735  "SerialDenseSolver<T>::applyRefinement: Cannot apply refinement if no original copy of A!");
736 
737 #ifdef HAVE_TEUCHOSNUMERICS_EIGEN
738  // Implement templated GERFS or use Eigen.
739  return (-1);
740 #else
741 
742  OrdinalType NRHS = RHS_->numCols();
743  FERR_.resize( NRHS );
744  BERR_.resize( NRHS );
745  allocateWORK();
746 
747  INFO_ = 0;
748  std::vector<typename details::lapack_traits<ScalarType>::iwork_type> GERFS_WORK( N_ );
749  this->GERFS(ETranspChar[TRANS_], N_, NRHS, A_, LDA_, AF_, LDAF_, &IPIV_[0],
750  RHS_->values(), RHS_->stride(), LHS_->values(), LHS_->stride(),
751  &FERR_[0], &BERR_[0], &WORK_[0], &GERFS_WORK[0], &INFO_);
752 
753  solutionErrorsEstimated_ = true;
754  reciprocalConditionEstimated_ = true;
755  solutionRefined_ = true;
756 
757  return(INFO_);
758 #endif
759 }
760 
761 //=============================================================================
762 
763 template<typename OrdinalType, typename ScalarType>
765 {
766  if (R_.size()!=0) return(0); // Already computed
767 
768  R_.resize( M_ );
769  C_.resize( N_ );
770 
771  INFO_ = 0;
772  this->GEEQU (M_, N_, AF_, LDAF_, &R_[0], &C_[0], &ROWCND_, &COLCND_, &AMAX_, &INFO_);
773 
774  if (COLCND_<0.1*ScalarTraits<MagnitudeType>::one() ||
775  ROWCND_<0.1*ScalarTraits<MagnitudeType>::one() ||
777  shouldEquilibrate_ = true;
778 
779  return(INFO_);
780 }
781 
782 //=============================================================================
783 
784 template<typename OrdinalType, typename ScalarType>
786 {
787  OrdinalType i, j;
788  int ierr = 0;
789 
790  if (equilibratedA_) return(0); // Already done.
791  if (R_.size()==0) ierr = computeEquilibrateScaling(); // Compute R and C if needed.
792  if (ierr!=0) return(ierr); // If return value is not zero, then can't equilibrate.
793  if (A_==AF_) {
794  ScalarType * ptr;
795  for (j=0; j<N_; j++) {
796  ptr = A_ + j*LDA_;
797  ScalarType s1 = C_[j];
798  for (i=0; i<M_; i++) {
799  *ptr = *ptr*s1*R_[i];
800  ptr++;
801  }
802  }
803  }
804  else {
805  ScalarType * ptr;
806  ScalarType * ptr1;
807  for (j=0; j<N_; j++) {
808  ptr = A_ + j*LDA_;
809  ptr1 = AF_ + j*LDAF_;
810  ScalarType s1 = C_[j];
811  for (i=0; i<M_; i++) {
812  *ptr = *ptr*s1*R_[i];
813  ptr++;
814  *ptr1 = *ptr1*s1*R_[i];
815  ptr1++;
816  }
817  }
818  }
819 
820  equilibratedA_ = true;
821 
822  return(0);
823 }
824 
825 //=============================================================================
826 
827 template<typename OrdinalType, typename ScalarType>
829 {
830  OrdinalType i, j;
831  int ierr = 0;
832 
833  if (equilibratedB_) return(0); // Already done.
834  if (R_.size()==0) ierr = computeEquilibrateScaling(); // Compute R and C if needed.
835  if (ierr!=0) return(ierr); // Can't count on R and C being computed.
836 
837  MagnitudeType * R_tmp = &R_[0];
838  if (transpose_) R_tmp = &C_[0];
839 
840  OrdinalType LDB = RHS_->stride(), NRHS = RHS_->numCols();
841  ScalarType * B = RHS_->values();
842  ScalarType * ptr;
843  for (j=0; j<NRHS; j++) {
844  ptr = B + j*LDB;
845  for (i=0; i<M_; i++) {
846  *ptr = *ptr*R_tmp[i];
847  ptr++;
848  }
849  }
850 
851  equilibratedB_ = true;
852 
853  return(0);
854 }
855 
856 //=============================================================================
857 
858 template<typename OrdinalType, typename ScalarType>
860 {
861  OrdinalType i, j;
862 
863  if (!equilibratedB_) return(0); // Nothing to do
864 
865  MagnitudeType * C_tmp = &C_[0];
866  if (transpose_) C_tmp = &R_[0];
867 
868  OrdinalType LDX = LHS_->stride(), NLHS = LHS_->numCols();
869  ScalarType * X = LHS_->values();
870  ScalarType * ptr;
871  for (j=0; j<NLHS; j++) {
872  ptr = X + j*LDX;
873  for (i=0; i<N_; i++) {
874  *ptr = *ptr*C_tmp[i];
875  ptr++;
876  }
877  }
878 
879  return(0);
880 }
881 
882 //=============================================================================
883 
884 template<typename OrdinalType, typename ScalarType>
886 {
887 
888  if (!factored()) factor(); // Need matrix factored.
889 
890 #ifdef HAVE_TEUCHOSNUMERICS_EIGEN
891  EigenMatrixMap invMap( AF_, M_, N_, EigenOuterStride(LDAF_) );
892  EigenMatrix invMat = lu_.inverse();
893  for (OrdinalType j=0; j<invMap.outerSize(); j++) {
894  for (OrdinalType i=0; i<invMap.innerSize(); i++) {
895  invMap(i,j) = invMat(i,j);
896  }
897  }
898 #else
899 
900  /* This section work with LAPACK Version 3.0 only
901  // Setting LWORK = -1 and calling GETRI will return optimal work space size in WORK_TMP
902  OrdinalType LWORK_TMP = -1;
903  ScalarType WORK_TMP;
904  GETRI( N_, AF_, LDAF_, &IPIV_[0], &WORK_TMP, &LWORK_TMP, &INFO_);
905  LWORK_TMP = WORK_TMP; // Convert to integer
906  if (LWORK_TMP>LWORK_) {
907  LWORK_ = LWORK_TMP;
908  WORK_.resize( LWORK_ );
909  }
910  */
911  // This section will work with any version of LAPACK
912  allocateWORK();
913 
914  INFO_ = 0;
915  this->GETRI( N_, AF_, LDAF_, &IPIV_[0], &WORK_[0], LWORK_, &INFO_);
916 #endif
917 
918  inverted_ = true;
919  factored_ = false;
920 
921  return(INFO_);
922 
923 }
924 
925 //=============================================================================
926 
927 template<typename OrdinalType, typename ScalarType>
929 {
930 #ifdef HAVE_TEUCHOSNUMERICS_EIGEN
931  // Implement templated GECON or use Eigen. Eigen currently doesn't have a condition estimation function.
932  return (-1);
933 #else
934  if (reciprocalConditionEstimated()) {
935  Value = RCOND_;
936  return(0); // Already computed, just return it.
937  }
938 
939  if ( ANORM_<ScalarTraits<MagnitudeType>::zero() ) ANORM_ = Matrix_->normOne();
940 
941  int ierr = 0;
942  if (!factored()) ierr = factor(); // Need matrix factored.
943  if (ierr!=0) return(ierr);
944 
945  allocateWORK();
946 
947  // We will assume a one-norm condition number
948  INFO_ = 0;
949  std::vector<typename details::lapack_traits<ScalarType>::iwork_type> GECON_WORK( 2*N_ );
950  this->GECON( '1', N_, AF_, LDAF_, ANORM_, &RCOND_, &WORK_[0], &GECON_WORK[0], &INFO_);
951 
952  reciprocalConditionEstimated_ = true;
953  Value = RCOND_;
954 
955  return(INFO_);
956 #endif
957 }
958 //=============================================================================
959 
960 template<typename OrdinalType, typename ScalarType>
962 
963  if (Matrix_!=Teuchos::null) os << "Solver Matrix" << std::endl << *Matrix_ << std::endl;
964  if (Factor_!=Teuchos::null) os << "Solver Factored Matrix" << std::endl << *Factor_ << std::endl;
965  if (LHS_ !=Teuchos::null) os << "Solver LHS" << std::endl << *LHS_ << std::endl;
966  if (RHS_ !=Teuchos::null) os << "Solver RHS" << std::endl << *RHS_ << std::endl;
967 
968 }
969 
970 } // namespace Teuchos
971 
972 #endif /* _TEUCHOS_SERIALDENSESOLVER_HPP_ */
std::vector< MagnitudeType > C() const
Returns a pointer to the column scale vector used for equilibration.
int equilibrateMatrix()
Equilibrates the this matrix.
T magnitudeType
Mandatory typedef for result of magnitude.
Templated serial dense matrix class.
void solveWithTranspose(bool flag)
If flag is true, causes all subsequent function calls to work with the transpose of this matrix...
Templated interface class to BLAS routines.
int solve()
Computes the solution X to AX = B for the this matrix and the B provided to SetVectors()..
std::vector< MagnitudeType > FERR() const
Returns a pointer to the forward error estimates computed by LAPACK.
MagnitudeType AMAX() const
Returns the absolute value of the largest entry of the this matrix (returns -1 if not yet computed)...
RCP< SerialDenseMatrix< OrdinalType, ScalarType > > getMatrix() const
Returns pointer to current matrix.
#define TEUCHOS_TEST_FOR_EXCEPTION(throw_exception_test, Exception, msg)
Macro for throwing an exception with breakpointing to ease debugging.
bool solved()
Returns true if the current set of vectors has been solved.
bool solutionRefined()
Returns true if the current set of vectors has been refined.
Teuchos header file which uses auto-configuration information to include necessary C++ headers...
Templated BLAS wrapper.
RCP< SerialDenseMatrix< OrdinalType, ScalarType > > getLHS() const
Returns pointer to current LHS.
bool reciprocalConditionEstimated()
Returns true if the condition number of the this matrix has been computed (value available via Recipr...
Object for storing data and providing functionality that is common to all computational classes...
OrdinalType numRows() const
Returns row dimension of system.
This structure defines some basic traits for a scalar field type.
int unequilibrateLHS()
Unscales the solution vectors if equilibration was used to solve the system.
MagnitudeType COLCND() const
Ratio of smallest to largest column scale factors for 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. ...
RCP< SerialDenseMatrix< OrdinalType, ScalarType > > getRHS() const
Returns pointer to current RHS.
bool factored()
Returns true if matrix is factored (factor available via getFactoredMatrix()).
Functionality and data that is common to all computational classes.
Templated interface class to LAPACK routines.
TEUCHOS_DEPRECATED RCP< T > rcp(T *p, Dealloc_T dealloc, bool owns_mem)
Deprecated.
Teuchos implementation details.
bool shouldEquilibrate()
Returns true if the LAPACK general rules for equilibration suggest you should equilibrate the system...
bool inverted()
Returns true if matrix inverse has been computed (inverse available via getFactoredMatrix()).
int reciprocalConditionEstimate(MagnitudeType &Value)
Returns the reciprocal of the 1-norm condition number of the this matrix.
The base Teuchos class.
void factorWithEquilibration(bool flag)
Causes equilibration to be called just before the matrix factorization as part of the call to factor...
std::vector< OrdinalType > IPIV() const
Returns pointer to pivot vector (if factorization has been computed), zero otherwise.
void solveWithTransposeFlag(Teuchos::ETransp trans)
All subsequent function calls will work with the transpose-type set by this method (Teuchos::NO_TRANS...
int invert()
Inverts the this matrix.
bool transpose()
Returns true if transpose of this matrix has and will be used.
int applyRefinement()
Apply Iterative Refinement.
int computeEquilibrateScaling()
Computes the scaling vector S(i) = 1/sqrt(A(i,i)) of the this matrix.
void estimateSolutionErrors(bool flag)
Causes all solves to estimate the forward and backward solution error.
bool equilibratedB()
Returns true if RHS is equilibrated (RHS available via getRHS()).
SerialDenseSolver()
Default constructor; matrix should be set using setMatrix(), LHS and RHS set with setVectors()...
bool solutionErrorsEstimated()
Returns true if forward and backward error estimated have been computed (available via FERR() and BER...
The Templated LAPACK Wrapper Class.
MagnitudeType ROWCND() const
Ratio of smallest to largest row scale factors for the this matrix (returns -1 if not yet computed)...
void Print(std::ostream &os) const
Print service methods; defines behavior of ostream << operator.
virtual ~SerialDenseSolver()
SerialDenseSolver destructor.
MagnitudeType ANORM() const
Returns the 1-Norm of the this matrix (returns -1 if not yet computed).
int equilibrateRHS()
Equilibrates the current RHS.
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.
The Teuchos namespace contains all of the classes, structs and enums used by Teuchos, as well as a number of utility routines.
std::vector< MagnitudeType > BERR() const
Returns a pointer to the backward error estimates computed by LAPACK.
MagnitudeType RCOND() const
Returns the reciprocal of the condition number of the this matrix (returns -1 if not yet computed)...
Defines basic traits for the scalar field type.
std::vector< MagnitudeType > R() const
Returns a pointer to the row scaling vector used for equilibration.
Smart reference counting pointer class for automatic garbage collection.
int factor()
Computes the in-place LU factorization of the matrix using the LAPACK routine _GETRF.
bool equilibratedA()
Returns true if factor is equilibrated (factor available via getFactoredMatrix()).
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 > > getFactoredMatrix() const
Returns pointer to factored matrix (assuming factorization has been performed).
A class for solving dense linear problems.
int setMatrix(const RCP< SerialDenseMatrix< OrdinalType, ScalarType > > &A)
Sets the pointers for coefficient matrix.
This class creates and provides basic support for dense rectangular matrix of templated type...