Teuchos - Trilinos Tools Package  Version of the Day
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
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; transpose_ = (trans != Teuchos::NO_TRANS) ? true : false; }
200 
202 
204  void solveToRefinedSolution(bool flag) {refineSolution_ = flag; return;}
205 
207 
211  void estimateSolutionErrors(bool flag);
213 
215 
216 
218 
221  int factor();
222 
224 
227  int solve();
228 
230 
233  int invert();
234 
236 
240 
242 
246  int equilibrateMatrix();
247 
249 
253  int equilibrateRHS();
254 
256 
260  int applyRefinement();
261 
263 
267  int unequilibrateLHS();
268 
270 
276  int reciprocalConditionEstimate(MagnitudeType & Value);
278 
280 
281 
283  bool transpose() {return(transpose_);}
284 
286  bool factored() {return(factored_);}
287 
289  bool equilibratedA() {return(equilibratedA_);}
290 
292  bool equilibratedB() {return(equilibratedB_);}
293 
295  bool shouldEquilibrate() {computeEquilibrateScaling(); return(shouldEquilibrate_);}
296 
298  bool solutionErrorsEstimated() {return(solutionErrorsEstimated_);}
299 
301  bool inverted() {return(inverted_);}
302 
304  bool reciprocalConditionEstimated() {return(reciprocalConditionEstimated_);}
305 
307  bool solved() {return(solved_);}
308 
310  bool solutionRefined() {return(solutionRefined_);}
312 
314 
315 
318 
321 
324 
327 
329  OrdinalType numRows() const {return(M_);}
330 
332  OrdinalType numCols() const {return(N_);}
333 
335  std::vector<OrdinalType> IPIV() const {return(IPIV_);}
336 
338  MagnitudeType ANORM() const {return(ANORM_);}
339 
341  MagnitudeType RCOND() const {return(RCOND_);}
342 
344 
346  MagnitudeType ROWCND() const {return(ROWCND_);}
347 
349 
351  MagnitudeType COLCND() const {return(COLCND_);}
352 
354  MagnitudeType AMAX() const {return(AMAX_);}
355 
357  std::vector<MagnitudeType> FERR() const {return(FERR_);}
358 
360  std::vector<MagnitudeType> BERR() const {return(BERR_);}
361 
363  std::vector<MagnitudeType> R() const {return(R_);}
364 
366  std::vector<MagnitudeType> C() const {return(C_);}
368 
370 
371  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>(Teuchos::Copy, *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  }
662  if (ierr != 0) return(ierr); // Can't equilibrate B, so return.
663 
664  TEUCHOS_TEST_FOR_EXCEPTION( RHS_==Teuchos::null, std::invalid_argument,
665  "SerialDenseSolver<T>::solve: No right-hand side vector (RHS) has been set for the linear system!");
666  TEUCHOS_TEST_FOR_EXCEPTION( LHS_==Teuchos::null, std::invalid_argument,
667  "SerialDenseSolver<T>::solve: No solution vector (LHS) has been set for the linear system!");
668 
669  if (inverted()) {
670 
671  TEUCHOS_TEST_FOR_EXCEPTION( RHS_->values() == LHS_->values(), std::invalid_argument,
672  "SerialDenseSolver<T>::solve: X and B must be different vectors if matrix is inverted.");
673  TEUCHOS_TEST_FOR_EXCEPTION( (equilibratedA_ && !equilibratedB_) || (!equilibratedA_ && equilibratedB_) ,
674  std::logic_error, "SerialDenseSolver<T>::solve: Matrix and vectors must be similarly scaled!");
675 
676  INFO_ = 0;
677  this->GEMM(TRANS_, Teuchos::NO_TRANS, N_, RHS_->numCols(), N_, 1.0, AF_, LDAF_,
678  RHS_->values(), RHS_->stride(), 0.0, LHS_->values(), LHS_->stride());
679  if (INFO_!=0) return(INFO_);
680  solved_ = true;
681  }
682  else {
683 
684  if (!factored()) factor(); // Matrix must be factored
685 
686  TEUCHOS_TEST_FOR_EXCEPTION( (equilibratedA_ && !equilibratedB_) || (!equilibratedA_ && equilibratedB_) ,
687  std::logic_error, "SerialDenseSolver<T>::solve: Matrix and vectors must be similarly scaled!");
688 
689  if (RHS_->values()!=LHS_->values()) {
690  (*LHS_) = (*RHS_); // Copy B to X if needed
691  }
692  INFO_ = 0;
693 
694 #ifdef HAVE_TEUCHOSNUMERICS_EIGEN
695  EigenMatrixMap rhsMap( RHS_->values(), RHS_->numRows(), RHS_->numCols(), EigenOuterStride(RHS_->stride()) );
696  EigenMatrixMap lhsMap( LHS_->values(), LHS_->numRows(), LHS_->numCols(), EigenOuterStride(LHS_->stride()) );
697  if ( TRANS_ == Teuchos::NO_TRANS ) {
698  lhsMap = lu_.solve(rhsMap);
699  } else if ( TRANS_ == Teuchos::TRANS ) {
700  lu_.matrixLU().template triangularView<Eigen::Upper>().transpose().solveInPlace(lhsMap);
701  lu_.matrixLU().template triangularView<Eigen::UnitLower>().transpose().solveInPlace(lhsMap);
702  EigenMatrix x = lu_.permutationP().transpose() * lhsMap;
703  lhsMap = x;
704  } else {
705  lu_.matrixLU().template triangularView<Eigen::Upper>().adjoint().solveInPlace(lhsMap);
706  lu_.matrixLU().template triangularView<Eigen::UnitLower>().adjoint().solveInPlace(lhsMap);
707  EigenMatrix x = lu_.permutationP().transpose() * lhsMap;
708  lhsMap = x;
709  }
710 #else
711  this->GETRS(ETranspChar[TRANS_], N_, RHS_->numCols(), AF_, LDAF_, &IPIV_[0], LHS_->values(), LHS_->stride(), &INFO_);
712 #endif
713 
714  if (INFO_!=0) return(INFO_);
715  solved_ = true;
716 
717  }
718 
719  // Warn the user that the matrix should be equilibrated, if it isn't being done already.
720  if (shouldEquilibrate() && !equilibratedA_)
721  std::cout << "WARNING! SerialDenseSolver<T>::solve: System should be equilibrated!" << std::endl;
722 
723  int ierr1=0;
724  if (refineSolution_ && !inverted()) ierr1 = applyRefinement();
725  if (ierr1!=0)
726  return(ierr1);
727 
728  if (equilibrate_) ierr1 = unequilibrateLHS();
729  return(ierr1);
730 }
731 //=============================================================================
732 
733 template<typename OrdinalType, typename ScalarType>
735 {
736  TEUCHOS_TEST_FOR_EXCEPTION(!solved(), std::logic_error,
737  "SerialDenseSolver<T>::applyRefinement: Must have an existing solution!");
738  TEUCHOS_TEST_FOR_EXCEPTION(A_==AF_, std::logic_error,
739  "SerialDenseSolver<T>::applyRefinement: Cannot apply refinement if no original copy of A!");
740 
741 #ifdef HAVE_TEUCHOSNUMERICS_EIGEN
742  // Implement templated GERFS or use Eigen.
743  return (-1);
744 #else
745 
746  OrdinalType NRHS = RHS_->numCols();
747  FERR_.resize( NRHS );
748  BERR_.resize( NRHS );
749  allocateWORK();
750 
751  INFO_ = 0;
752  std::vector<typename details::lapack_traits<ScalarType>::iwork_type> GERFS_WORK( N_ );
753  this->GERFS(ETranspChar[TRANS_], N_, NRHS, A_, LDA_, AF_, LDAF_, &IPIV_[0],
754  RHS_->values(), RHS_->stride(), LHS_->values(), LHS_->stride(),
755  &FERR_[0], &BERR_[0], &WORK_[0], &GERFS_WORK[0], &INFO_);
756 
757  solutionErrorsEstimated_ = true;
758  reciprocalConditionEstimated_ = true;
759  solutionRefined_ = true;
760 
761  return(INFO_);
762 #endif
763 }
764 
765 //=============================================================================
766 
767 template<typename OrdinalType, typename ScalarType>
769 {
770  if (R_.size()!=0) return(0); // Already computed
771 
772  R_.resize( M_ );
773  C_.resize( N_ );
774 
775  INFO_ = 0;
776  this->GEEQU (M_, N_, AF_, LDAF_, &R_[0], &C_[0], &ROWCND_, &COLCND_, &AMAX_, &INFO_);
777 
778  if (COLCND_<0.1*ScalarTraits<MagnitudeType>::one() ||
779  ROWCND_<0.1*ScalarTraits<MagnitudeType>::one() ||
781  shouldEquilibrate_ = true;
782 
783  return(INFO_);
784 }
785 
786 //=============================================================================
787 
788 template<typename OrdinalType, typename ScalarType>
790 {
791  OrdinalType i, j;
792  int ierr = 0;
793 
794  if (equilibratedA_) return(0); // Already done.
795  if (R_.size()==0) ierr = computeEquilibrateScaling(); // Compute R and C if needed.
796  if (ierr!=0) return(ierr); // If return value is not zero, then can't equilibrate.
797  if (A_==AF_) {
798  ScalarType * ptr;
799  for (j=0; j<N_; j++) {
800  ptr = A_ + j*LDA_;
801  ScalarType s1 = C_[j];
802  for (i=0; i<M_; i++) {
803  *ptr = *ptr*s1*R_[i];
804  ptr++;
805  }
806  }
807  }
808  else {
809  ScalarType * ptr;
810  ScalarType * ptr1;
811  for (j=0; j<N_; j++) {
812  ptr = A_ + j*LDA_;
813  ptr1 = AF_ + j*LDAF_;
814  ScalarType s1 = C_[j];
815  for (i=0; i<M_; i++) {
816  *ptr = *ptr*s1*R_[i];
817  ptr++;
818  *ptr1 = *ptr1*s1*R_[i];
819  ptr1++;
820  }
821  }
822  }
823 
824  equilibratedA_ = true;
825 
826  return(0);
827 }
828 
829 //=============================================================================
830 
831 template<typename OrdinalType, typename ScalarType>
833 {
834  OrdinalType i, j;
835  int ierr = 0;
836 
837  if (equilibratedB_) return(0); // Already done.
838  if (R_.size()==0) ierr = computeEquilibrateScaling(); // Compute R and C if needed.
839  if (ierr!=0) return(ierr); // Can't count on R and C being computed.
840 
841  MagnitudeType * R_tmp = &R_[0];
842  if (transpose_) R_tmp = &C_[0];
843 
844  OrdinalType LDB = RHS_->stride(), NRHS = RHS_->numCols();
845  ScalarType * B = RHS_->values();
846  ScalarType * ptr;
847  for (j=0; j<NRHS; j++) {
848  ptr = B + j*LDB;
849  for (i=0; i<M_; i++) {
850  *ptr = *ptr*R_tmp[i];
851  ptr++;
852  }
853  }
854 
855  equilibratedB_ = true;
856 
857  return(0);
858 }
859 
860 //=============================================================================
861 
862 template<typename OrdinalType, typename ScalarType>
864 {
865  OrdinalType i, j;
866 
867  if (!equilibratedB_) return(0); // Nothing to do
868 
869  MagnitudeType * C_tmp = &C_[0];
870  if (transpose_) C_tmp = &R_[0];
871 
872  OrdinalType LDX = LHS_->stride(), NLHS = LHS_->numCols();
873  ScalarType * X = LHS_->values();
874  ScalarType * ptr;
875  for (j=0; j<NLHS; j++) {
876  ptr = X + j*LDX;
877  for (i=0; i<N_; i++) {
878  *ptr = *ptr*C_tmp[i];
879  ptr++;
880  }
881  }
882 
883  return(0);
884 }
885 
886 //=============================================================================
887 
888 template<typename OrdinalType, typename ScalarType>
890 {
891 
892  if (!factored()) factor(); // Need matrix factored.
893 
894 #ifdef HAVE_TEUCHOSNUMERICS_EIGEN
895  EigenMatrixMap invMap( AF_, M_, N_, EigenOuterStride(LDAF_) );
896  EigenMatrix invMat = lu_.inverse();
897  for (OrdinalType j=0; j<invMap.outerSize(); j++) {
898  for (OrdinalType i=0; i<invMap.innerSize(); i++) {
899  invMap(i,j) = invMat(i,j);
900  }
901  }
902 #else
903 
904  /* This section work with LAPACK Version 3.0 only
905  // Setting LWORK = -1 and calling GETRI will return optimal work space size in WORK_TMP
906  OrdinalType LWORK_TMP = -1;
907  ScalarType WORK_TMP;
908  GETRI( N_, AF_, LDAF_, &IPIV_[0], &WORK_TMP, &LWORK_TMP, &INFO_);
909  LWORK_TMP = WORK_TMP; // Convert to integer
910  if (LWORK_TMP>LWORK_) {
911  LWORK_ = LWORK_TMP;
912  WORK_.resize( LWORK_ );
913  }
914  */
915  // This section will work with any version of LAPACK
916  allocateWORK();
917 
918  INFO_ = 0;
919  this->GETRI( N_, AF_, LDAF_, &IPIV_[0], &WORK_[0], LWORK_, &INFO_);
920 #endif
921 
922  inverted_ = true;
923  factored_ = false;
924 
925  return(INFO_);
926 
927 }
928 
929 //=============================================================================
930 
931 template<typename OrdinalType, typename ScalarType>
933 {
934 #ifdef HAVE_TEUCHOSNUMERICS_EIGEN
935  // Implement templated GECON or use Eigen. Eigen currently doesn't have a condition estimation function.
936  return (-1);
937 #else
938  if (reciprocalConditionEstimated()) {
939  Value = RCOND_;
940  return(0); // Already computed, just return it.
941  }
942 
943  if ( ANORM_<ScalarTraits<MagnitudeType>::zero() ) ANORM_ = Matrix_->normOne();
944 
945  int ierr = 0;
946  if (!factored()) ierr = factor(); // Need matrix factored.
947  if (ierr!=0) return(ierr);
948 
949  allocateWORK();
950 
951  // We will assume a one-norm condition number
952  INFO_ = 0;
953  std::vector<typename details::lapack_traits<ScalarType>::iwork_type> GECON_WORK( 2*N_ );
954  this->GECON( '1', N_, AF_, LDAF_, ANORM_, &RCOND_, &WORK_[0], &GECON_WORK[0], &INFO_);
955 
956  reciprocalConditionEstimated_ = true;
957  Value = RCOND_;
958 
959  return(INFO_);
960 #endif
961 }
962 //=============================================================================
963 
964 template<typename OrdinalType, typename ScalarType>
966 
967  if (Matrix_!=Teuchos::null) os << "Solver Matrix" << std::endl << *Matrix_ << std::endl;
968  if (Factor_!=Teuchos::null) os << "Solver Factored Matrix" << std::endl << *Factor_ << std::endl;
969  if (LHS_ !=Teuchos::null) os << "Solver LHS" << std::endl << *LHS_ << std::endl;
970  if (RHS_ !=Teuchos::null) os << "Solver RHS" << std::endl << *RHS_ << std::endl;
971 
972 }
973 
974 } // namespace Teuchos
975 
976 #endif /* _TEUCHOS_SERIALDENSESOLVER_HPP_ */
int equilibrateMatrix()
Equilibrates the this matrix.
void Print(std::ostream &os) const
Print service methods; defines behavior of ostream &lt;&lt; operator.
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...
std::vector< MagnitudeType > R() const
Returns a pointer to the row scaling vector used for equilibration.
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()..
#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.
bool reciprocalConditionEstimated()
Returns true if the condition number of the this matrix has been computed (value available via Recipr...
RCP< SerialDenseMatrix< OrdinalType, ScalarType > > getFactoredMatrix() const
Returns pointer to factored matrix (assuming factorization has been performed).
Object for storing data and providing functionality that is common to all computational classes...
RCP< SerialDenseMatrix< OrdinalType, ScalarType > > getRHS() const
Returns pointer to current RHS.
std::vector< MagnitudeType > FERR() const
Returns a pointer to the forward error estimates computed by LAPACK.
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.
void solveToRefinedSolution(bool flag)
Causes all solves to compute solution to best ability using iterative refinement. ...
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.
RCP< SerialDenseMatrix< OrdinalType, ScalarType > > getLHS() const
Returns pointer to current LHS.
TEUCHOS_DEPRECATED RCP< T > rcp(T *p, Dealloc_T dealloc, bool owns_mem)
Deprecated.
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...
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.
MagnitudeType ANORM() const
Returns the 1-Norm of the this matrix (returns -1 if not yet computed).
bool equilibratedB()
Returns true if RHS is equilibrated (RHS available via getRHS()).
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.
virtual ~SerialDenseSolver()
SerialDenseSolver destructor.
std::vector< OrdinalType > IPIV() const
Returns pointer to pivot vector (if factorization has been computed), zero otherwise.
OrdinalType numRows() const
Returns row dimension of system.
std::vector< MagnitudeType > C() const
Returns a pointer to the column scale vector used for equilibration.
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).
MagnitudeType COLCND() const
Ratio of smallest to largest column scale factors for the this matrix (returns -1 if not yet computed...
MagnitudeType ROWCND() const
Ratio of smallest to largest row scale factors for the this matrix (returns -1 if not yet computed)...
std::vector< MagnitudeType > BERR() const
Returns a pointer to the backward 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)...
Defines basic traits for the scalar field type.
OrdinalType numCols() const
Returns column dimension of system.
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()).
RCP< SerialDenseMatrix< OrdinalType, ScalarType > > getMatrix() const
Returns pointer to current matrix.
Reference-counted pointer class and non-member templated function implementations.
static T one()
Returns representation of one for this scalar type.
A class for solving dense linear problems.
MagnitudeType RCOND() const
Returns the reciprocal of the condition number of the this matrix (returns -1 if not yet computed)...
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...