Teuchos - Trilinos Tools Package Version of the Day
Loading...
Searching...
No Matches
Teuchos_SerialSpdDenseSolver.hpp
Go to the documentation of this file.
1
2// @HEADER
3// *****************************************************************************
4// Teuchos: Common Tools Package
5//
6// Copyright 2004 NTESS and the Teuchos contributors.
7// SPDX-License-Identifier: BSD-3-Clause
8// *****************************************************************************
9// @HEADER
10
11#ifndef TEUCHOS_SERIALSPDDENSESOLVER_H
12#define TEUCHOS_SERIALSPDDENSESOLVER_H
13
18#include "Teuchos_CompObject.hpp"
19#include "Teuchos_BLAS.hpp"
20#include "Teuchos_LAPACK.hpp"
21#include "Teuchos_RCP.hpp"
22#include "Teuchos_ConfigDefs.hpp"
23#include "Teuchos_SerialSymDenseMatrix.hpp"
24#include "Teuchos_SerialDenseSolver.hpp"
25#include "Teuchos_SerialDenseMatrix.hpp"
26
27#ifdef HAVE_TEUCHOSNUMERICS_EIGEN
28#include "Eigen/Dense"
29#endif
30
98namespace Teuchos {
99
100 template<typename OrdinalType, typename ScalarType>
101 class SerialSpdDenseSolver : public CompObject, public Object, public BLAS<OrdinalType, ScalarType>,
102 public LAPACK<OrdinalType, ScalarType>
103 {
104 public:
105
106 typedef typename Teuchos::ScalarTraits<ScalarType>::magnitudeType MagnitudeType;
107#ifdef HAVE_TEUCHOSNUMERICS_EIGEN
108 typedef typename Eigen::Matrix<ScalarType,Eigen::Dynamic,Eigen::Dynamic,Eigen::ColMajor> EigenMatrix;
109 typedef typename Eigen::Matrix<ScalarType,Eigen::Dynamic,1> EigenVector;
110 typedef typename Eigen::InnerStride<Eigen::Dynamic> EigenInnerStride;
111 typedef typename Eigen::OuterStride<Eigen::Dynamic> EigenOuterStride;
112 typedef typename Eigen::Map<EigenVector,0,EigenInnerStride > EigenVectorMap;
113 typedef typename Eigen::Map<const EigenVector,0,EigenInnerStride > EigenConstVectorMap;
114 typedef typename Eigen::Map<EigenMatrix,0,EigenOuterStride > EigenMatrixMap;
115 typedef typename Eigen::Map<const EigenMatrix,0,EigenOuterStride > EigenConstMatrixMap;
116 typedef typename Eigen::PermutationMatrix<Eigen::Dynamic,Eigen::Dynamic,OrdinalType> EigenPermutationMatrix;
117 typedef typename Eigen::Array<OrdinalType,Eigen::Dynamic,1> EigenOrdinalArray;
118 typedef typename Eigen::Map<EigenOrdinalArray> EigenOrdinalArrayMap;
119 typedef typename Eigen::Array<ScalarType,Eigen::Dynamic,1> EigenScalarArray;
120 typedef typename Eigen::Map<EigenScalarArray> EigenScalarArrayMap;
121#endif
122
124
125
126 SerialSpdDenseSolver();
127
128
130 virtual ~SerialSpdDenseSolver();
132
134
135
137 int setMatrix(const RCP<SerialSymDenseMatrix<OrdinalType,ScalarType> >& A_in);
138
140
143 int setVectors(const RCP<SerialDenseMatrix<OrdinalType, ScalarType> >& X,
144 const RCP<SerialDenseMatrix<OrdinalType, ScalarType> >& B);
146
148
149
151
153 void factorWithEquilibration(bool flag) {equilibrate_ = flag; return;}
154
156 void solveToRefinedSolution(bool flag) {refineSolution_ = flag; return;}
157
159
162 void estimateSolutionErrors(bool flag);
164
166
167
169
172 int factor();
173
175
178 int solve();
179
181
186 int invert();
187
189
192 int computeEquilibrateScaling();
193
195
198 int equilibrateMatrix();
199
201
204 int equilibrateRHS();
205
207
210 int applyRefinement();
211
213
216 int unequilibrateLHS();
217
219
225 int reciprocalConditionEstimate(MagnitudeType & Value);
227
229
230
232 bool transpose() {return(transpose_);}
233
235 bool factored() {return(factored_);}
236
238 bool equilibratedA() {return(equilibratedA_);}
239
241 bool equilibratedB() {return(equilibratedB_);}
242
244 bool shouldEquilibrate() {computeEquilibrateScaling(); return(shouldEquilibrate_);}
245
247 bool solutionErrorsEstimated() {return(solutionErrorsEstimated_);}
248
250 bool inverted() {return(inverted_);}
251
253 bool reciprocalConditionEstimated() {return(reciprocalConditionEstimated_);}
254
256 bool solved() {return(solved_);}
257
259 bool solutionRefined() {return(solutionRefined_);}
261
263
264
266 RCP<SerialSymDenseMatrix<OrdinalType, ScalarType> > getMatrix() const {return(Matrix_);}
267
269 RCP<SerialSymDenseMatrix<OrdinalType, ScalarType> > getFactoredMatrix() const {return(Factor_);}
270
272 RCP<SerialDenseMatrix<OrdinalType, ScalarType> > getLHS() const {return(LHS_);}
273
275 RCP<SerialDenseMatrix<OrdinalType, ScalarType> > getRHS() const {return(RHS_);}
276
278 OrdinalType numRows() const {return(numRowCols_);}
279
281 OrdinalType numCols() const {return(numRowCols_);}
282
284 MagnitudeType ANORM() const {return(ANORM_);}
285
287 MagnitudeType RCOND() const {return(RCOND_);}
288
290
292 MagnitudeType SCOND() {return(SCOND_);};
293
295 MagnitudeType AMAX() const {return(AMAX_);}
296
298 std::vector<MagnitudeType> FERR() const {return(FERR_);}
299
301 std::vector<MagnitudeType> BERR() const {return(BERR_);}
302
304 std::vector<MagnitudeType> R() const {return(R_);}
305
307
309
310
311 void Print(std::ostream& os) const;
313
314 protected:
315
316 void allocateWORK() { LWORK_ = 4*numRowCols_; WORK_.resize( LWORK_ ); return;}
317 void allocateIWORK() { IWORK_.resize( numRowCols_ ); return;}
318 void resetMatrix();
319 void resetVectors();
320
321 bool equilibrate_;
322 bool shouldEquilibrate_;
323 bool equilibratedA_;
324 bool equilibratedB_;
325 bool transpose_;
326 bool factored_;
327 bool estimateSolutionErrors_;
328 bool solutionErrorsEstimated_;
329 bool solved_;
330 bool inverted_;
331 bool reciprocalConditionEstimated_;
332 bool refineSolution_;
333 bool solutionRefined_;
334
335 OrdinalType numRowCols_;
336 OrdinalType LDA_;
337 OrdinalType LDAF_;
338 OrdinalType INFO_;
339 OrdinalType LWORK_;
340
341 std::vector<int> IWORK_;
342
343 MagnitudeType ANORM_;
344 MagnitudeType RCOND_;
345 MagnitudeType SCOND_;
346 MagnitudeType AMAX_;
347
348 RCP<SerialSymDenseMatrix<OrdinalType, ScalarType> > Matrix_;
349 RCP<SerialDenseMatrix<OrdinalType, ScalarType> > LHS_;
350 RCP<SerialDenseMatrix<OrdinalType, ScalarType> > RHS_;
351 RCP<SerialSymDenseMatrix<OrdinalType, ScalarType> > Factor_;
352
353 ScalarType * A_;
354 ScalarType * AF_;
355 std::vector<MagnitudeType> FERR_;
356 std::vector<MagnitudeType> BERR_;
357 std::vector<ScalarType> WORK_;
358 std::vector<MagnitudeType> R_;
359#ifdef HAVE_TEUCHOSNUMERICS_EIGEN
360 Eigen::LLT<EigenMatrix,Eigen::Upper> lltu_;
361 Eigen::LLT<EigenMatrix,Eigen::Lower> lltl_;
362#endif
363
364 private:
365 // SerialSpdDenseSolver copy constructor (put here because we don't want user access)
366
367 SerialSpdDenseSolver(const SerialSpdDenseSolver<OrdinalType, ScalarType>& Source);
368 SerialSpdDenseSolver & operator=(const SerialSpdDenseSolver<OrdinalType, ScalarType>& Source);
369
370 };
371
372 // Helper traits to distinguish work arrays for real and complex-valued datatypes.
373 using namespace details;
374
375//=============================================================================
376
377template<typename OrdinalType, typename ScalarType>
378SerialSpdDenseSolver<OrdinalType,ScalarType>::SerialSpdDenseSolver()
379 : CompObject(),
380 equilibrate_(false),
381 shouldEquilibrate_(false),
382 equilibratedA_(false),
383 equilibratedB_(false),
384 transpose_(false),
385 factored_(false),
386 estimateSolutionErrors_(false),
387 solutionErrorsEstimated_(false),
388 solved_(false),
389 inverted_(false),
390 reciprocalConditionEstimated_(false),
391 refineSolution_(false),
392 solutionRefined_(false),
393 numRowCols_(0),
394 LDA_(0),
395 LDAF_(0),
396 INFO_(0),
397 LWORK_(0),
398 ANORM_(ScalarTraits<MagnitudeType>::zero()),
399 RCOND_(ScalarTraits<MagnitudeType>::zero()),
400 SCOND_(ScalarTraits<MagnitudeType>::zero()),
401 AMAX_(ScalarTraits<MagnitudeType>::zero()),
402 A_(0),
403 AF_(0)
404{
405 resetMatrix();
406}
407//=============================================================================
408
409template<typename OrdinalType, typename ScalarType>
412
413//=============================================================================
414
415template<typename OrdinalType, typename ScalarType>
416void SerialSpdDenseSolver<OrdinalType,ScalarType>::resetVectors()
417{
418 LHS_ = Teuchos::null;
419 RHS_ = Teuchos::null;
420 reciprocalConditionEstimated_ = false;
421 solutionRefined_ = false;
422 solved_ = false;
423 solutionErrorsEstimated_ = false;
424 equilibratedB_ = false;
425}
426//=============================================================================
427
428template<typename OrdinalType, typename ScalarType>
429void SerialSpdDenseSolver<OrdinalType,ScalarType>::resetMatrix()
430{
431 resetVectors();
432 equilibratedA_ = false;
433 factored_ = false;
434 inverted_ = false;
435 numRowCols_ = 0;
436 LDA_ = 0;
437 LDAF_ = 0;
438 ANORM_ = -ScalarTraits<MagnitudeType>::one();
439 RCOND_ = -ScalarTraits<MagnitudeType>::one();
440 SCOND_ = -ScalarTraits<MagnitudeType>::one();
441 AMAX_ = -ScalarTraits<MagnitudeType>::one();
442 A_ = 0;
443 AF_ = 0;
444 INFO_ = 0;
445 LWORK_ = 0;
446 R_.resize(0);
447}
448//=============================================================================
449
450template<typename OrdinalType, typename ScalarType>
451int SerialSpdDenseSolver<OrdinalType,ScalarType>::setMatrix(const RCP<SerialSymDenseMatrix<OrdinalType,ScalarType> >& A)
452{
453 resetMatrix();
454 Matrix_ = A;
455 Factor_ = A;
456 numRowCols_ = A->numRows();
457 LDA_ = A->stride();
458 LDAF_ = LDA_;
459 A_ = A->values();
460 AF_ = A->values();
461 return(0);
462}
463//=============================================================================
464
465template<typename OrdinalType, typename ScalarType>
466int SerialSpdDenseSolver<OrdinalType,ScalarType>::setVectors(const RCP<SerialDenseMatrix<OrdinalType,ScalarType> >& X,
467 const RCP<SerialDenseMatrix<OrdinalType,ScalarType> >& B)
468{
469 // Check that these new vectors are consistent.
470 TEUCHOS_TEST_FOR_EXCEPTION(B->numRows()!=X->numRows() || B->numCols() != X->numCols(), std::invalid_argument,
471 "SerialSpdDenseSolver<T>::setVectors: X and B are not the same size!");
472 TEUCHOS_TEST_FOR_EXCEPTION(B->values()==0, std::invalid_argument,
473 "SerialSpdDenseSolver<T>::setVectors: B is an empty SerialDenseMatrix<T>!");
474 TEUCHOS_TEST_FOR_EXCEPTION(X->values()==0, std::invalid_argument,
475 "SerialSpdDenseSolver<T>::setVectors: X is an empty SerialDenseMatrix<T>!");
476 TEUCHOS_TEST_FOR_EXCEPTION(B->stride()<1, std::invalid_argument,
477 "SerialSpdDenseSolver<T>::setVectors: B has an invalid stride!");
478 TEUCHOS_TEST_FOR_EXCEPTION(X->stride()<1, std::invalid_argument,
479 "SerialSpdDenseSolver<T>::setVectors: X has an invalid stride!");
480
481 resetVectors();
482 LHS_ = X;
483 RHS_ = B;
484 return(0);
485}
486//=============================================================================
487
488template<typename OrdinalType, typename ScalarType>
489void SerialSpdDenseSolver<OrdinalType,ScalarType>::estimateSolutionErrors(bool flag)
490{
491 estimateSolutionErrors_ = flag;
492
493 // If the errors are estimated, this implies that the solution must be refined
494 refineSolution_ = refineSolution_ || flag;
495}
496//=============================================================================
497
498template<typename OrdinalType, typename ScalarType>
499int SerialSpdDenseSolver<OrdinalType,ScalarType>::factor() {
500
501 if (factored()) return(0); // Already factored
502
503 TEUCHOS_TEST_FOR_EXCEPTION(inverted(), std::logic_error,
504 "SerialSpdDenseSolver<T>::factor: Cannot factor an inverted matrix!");
505
506 ANORM_ = Matrix_->normOne(); // Compute 1-Norm of A
507
508
509 // If we want to refine the solution, then the factor must
510 // be stored separatedly from the original matrix
511
512 if (A_ == AF_)
513 if (refineSolution_ ) {
514 Factor_ = rcp( new SerialSymDenseMatrix<OrdinalType,ScalarType>(*Matrix_) );
515 AF_ = Factor_->values();
516 LDAF_ = Factor_->stride();
517 }
518
519 int ierr = 0;
520 if (equilibrate_) ierr = equilibrateMatrix();
521
522 if (ierr!=0) return(ierr);
523
524 INFO_ = 0;
525
526#ifdef HAVE_TEUCHOSNUMERICS_EIGEN
527 EigenMatrixMap aMap( AF_, numRowCols_, numRowCols_, EigenOuterStride(LDAF_) );
528
529 if (Matrix_->upper())
530 lltu_.compute(aMap);
531 else
532 lltl_.compute(aMap);
533
534#else
535 this->POTRF(Matrix_->UPLO(), numRowCols_, AF_, LDAF_, &INFO_);
536#endif
537
538 factored_ = true;
539
540 return(INFO_);
541}
542
543//=============================================================================
544
545template<typename OrdinalType, typename ScalarType>
546int SerialSpdDenseSolver<OrdinalType,ScalarType>::solve() {
547
548 // We will call one of four routines depending on what services the user wants and
549 // whether or not the matrix has been inverted or factored already.
550 //
551 // If the matrix has been inverted, use DGEMM to compute solution.
552 // Otherwise, if the user want the matrix to be equilibrated or wants a refined solution, we will
553 // call the X interface.
554 // Otherwise, if the matrix is already factored we will call the TRS interface.
555 // Otherwise, if the matrix is unfactored we will call the SV interface.
556
557 int ierr = 0;
558 if (equilibrate_) {
559 ierr = equilibrateRHS();
560 }
561 if (ierr != 0) return(ierr); // Can't equilibrate B, so return.
562
563 TEUCHOS_TEST_FOR_EXCEPTION( RHS_==Teuchos::null, std::invalid_argument,
564 "SerialSpdDenseSolver<T>::solve: No right-hand side vector (RHS) has been set for the linear system!");
565 TEUCHOS_TEST_FOR_EXCEPTION( LHS_==Teuchos::null, std::invalid_argument,
566 "SerialSpdDenseSolver<T>::solve: No solution vector (LHS) has been set for the linear system!");
567
568 if (inverted()) {
569
570 TEUCHOS_TEST_FOR_EXCEPTION( RHS_->values() == LHS_->values(), std::invalid_argument,
571 "SerialSpdDenseSolver<T>::solve: X and B must be different vectors if matrix is inverted.");
572 TEUCHOS_TEST_FOR_EXCEPTION( (equilibratedA_ && !equilibratedB_) || (!equilibratedA_ && equilibratedB_) ,
573 std::logic_error, "SerialSpdDenseSolver<T>::solve: Matrix and vectors must be similarly scaled!");
574
575 INFO_ = 0;
576 this->GEMM(Teuchos::NO_TRANS, Teuchos::NO_TRANS, numRowCols_, RHS_->numCols(),
577 numRowCols_, 1.0, AF_, LDAF_, RHS_->values(), RHS_->stride(), 0.0,
578 LHS_->values(), LHS_->stride());
579 if (INFO_!=0) return(INFO_);
580 solved_ = true;
581 }
582 else {
583
584 if (!factored()) factor(); // Matrix must be factored
585
586 TEUCHOS_TEST_FOR_EXCEPTION( (equilibratedA_ && !equilibratedB_) || (!equilibratedA_ && equilibratedB_) ,
587 std::logic_error, "SerialSpdDenseSolver<T>::solve: Matrix and vectors must be similarly scaled!");
588
589 if (RHS_->values()!=LHS_->values()) {
590 (*LHS_) = (*RHS_); // Copy B to X if needed
591 }
592 INFO_ = 0;
593
594#ifdef HAVE_TEUCHOSNUMERICS_EIGEN
595 EigenMatrixMap rhsMap( RHS_->values(), RHS_->numRows(), RHS_->numCols(), EigenOuterStride(RHS_->stride()) );
596 EigenMatrixMap lhsMap( LHS_->values(), LHS_->numRows(), LHS_->numCols(), EigenOuterStride(LHS_->stride()) );
597
598 if (Matrix_->upper())
599 lhsMap = lltu_.solve(rhsMap);
600 else
601 lhsMap = lltl_.solve(rhsMap);
602
603#else
604 this->POTRS(Matrix_->UPLO(), numRowCols_, RHS_->numCols(), AF_, LDAF_, LHS_->values(), LHS_->stride(), &INFO_);
605#endif
606
607 if (INFO_!=0) return(INFO_);
608 solved_ = true;
609
610 }
611
612 // Warn users that the system should be equilibrated, but it wasn't.
613 if (shouldEquilibrate() && !equilibratedA_)
614 std::cout << "WARNING! SerialSpdDenseSolver<T>::solve: System should be equilibrated!" << std::endl;
615
616 int ierr1=0;
617 if (refineSolution_ && !inverted()) ierr1 = applyRefinement();
618 if (ierr1!=0)
619 return(ierr1);
620
621 if (equilibrate_) ierr1 = unequilibrateLHS();
622 return(ierr1);
623}
624//=============================================================================
625
626template<typename OrdinalType, typename ScalarType>
627int SerialSpdDenseSolver<OrdinalType,ScalarType>::applyRefinement()
628{
629 TEUCHOS_TEST_FOR_EXCEPTION(!solved(), std::logic_error,
630 "SerialSpdDenseSolver<T>::applyRefinement: Must have an existing solution!");
631 TEUCHOS_TEST_FOR_EXCEPTION(A_==AF_, std::logic_error,
632 "SerialSpdDenseSolver<T>::applyRefinement: Cannot apply refinement if no original copy of A!");
633
634
635#ifdef HAVE_TEUCHOSNUMERICS_EIGEN
636 // Implement templated PORFS or use Eigen.
637 return (-1);
638#else
639 OrdinalType NRHS = RHS_->numCols();
640 FERR_.resize( NRHS );
641 BERR_.resize( NRHS );
642 allocateWORK();
643 allocateIWORK();
644
645 INFO_ = 0;
646 std::vector<typename details::lapack_traits<ScalarType>::iwork_type> PORFS_WORK( numRowCols_ );
647 this->PORFS(Matrix_->UPLO(), numRowCols_, NRHS, A_, LDA_, AF_, LDAF_,
648 RHS_->values(), RHS_->stride(), LHS_->values(), LHS_->stride(),
649 &FERR_[0], &BERR_[0], &WORK_[0], &PORFS_WORK[0], &INFO_);
650
651 solutionErrorsEstimated_ = true;
652 reciprocalConditionEstimated_ = true;
653 solutionRefined_ = true;
654
655 return(INFO_);
656#endif
657}
658
659//=============================================================================
660
661template<typename OrdinalType, typename ScalarType>
662int SerialSpdDenseSolver<OrdinalType,ScalarType>::computeEquilibrateScaling()
663{
664 if (R_.size()!=0) return(0); // Already computed
665
666 R_.resize( numRowCols_ );
667
668 INFO_ = 0;
669 this->POEQU (numRowCols_, AF_, LDAF_, &R_[0], &SCOND_, &AMAX_, &INFO_);
670 if (SCOND_<0.1*ScalarTraits<MagnitudeType>::one() ||
671 AMAX_ < ScalarTraits<ScalarType>::rmin() ||
672 AMAX_ > ScalarTraits<ScalarType>::rmax())
673 shouldEquilibrate_ = true;
674
675 return(INFO_);
676}
677
678//=============================================================================
679
680template<typename OrdinalType, typename ScalarType>
681int SerialSpdDenseSolver<OrdinalType,ScalarType>::equilibrateMatrix()
682{
683 OrdinalType i, j;
684 int ierr = 0;
685
686 if (equilibratedA_) return(0); // Already done.
687 if (R_.size()==0) ierr = computeEquilibrateScaling(); // Compute R if needed.
688 if (ierr!=0) return(ierr); // If return value is not zero, then can't equilibrate.
689 if (Matrix_->upper()) {
690 if (A_==AF_) {
691 ScalarType * ptr;
692 for (j=0; j<numRowCols_; j++) {
693 ptr = A_ + j*LDA_;
694 ScalarType s1 = R_[j];
695 for (i=0; i<=j; i++) {
696 *ptr = *ptr*s1*R_[i];
697 ptr++;
698 }
699 }
700 }
701 else {
702 ScalarType * ptr;
703 ScalarType * ptr1;
704 for (j=0; j<numRowCols_; j++) {
705 ptr = A_ + j*LDA_;
706 ptr1 = AF_ + j*LDAF_;
707 ScalarType s1 = R_[j];
708 for (i=0; i<=j; i++) {
709 *ptr = *ptr*s1*R_[i];
710 ptr++;
711 *ptr1 = *ptr1*s1*R_[i];
712 ptr1++;
713 }
714 }
715 }
716 }
717 else {
718 if (A_==AF_) {
719 ScalarType * ptr;
720 for (j=0; j<numRowCols_; j++) {
721 ptr = A_ + j + j*LDA_;
722 ScalarType s1 = R_[j];
723 for (i=j; i<numRowCols_; i++) {
724 *ptr = *ptr*s1*R_[i];
725 ptr++;
726 }
727 }
728 }
729 else {
730 ScalarType * ptr;
731 ScalarType * ptr1;
732 for (j=0; j<numRowCols_; j++) {
733 ptr = A_ + j + j*LDA_;
734 ptr1 = AF_ + j + j*LDAF_;
735 ScalarType s1 = R_[j];
736 for (i=j; i<numRowCols_; i++) {
737 *ptr = *ptr*s1*R_[i];
738 ptr++;
739 *ptr1 = *ptr1*s1*R_[i];
740 ptr1++;
741 }
742 }
743 }
744 }
745
746 equilibratedA_ = true;
747
748 return(0);
749}
750
751//=============================================================================
752
753template<typename OrdinalType, typename ScalarType>
754int SerialSpdDenseSolver<OrdinalType,ScalarType>::equilibrateRHS()
755{
756 OrdinalType i, j;
757 int ierr = 0;
758
759 if (equilibratedB_) return(0); // Already done.
760 if (R_.size()==0) ierr = computeEquilibrateScaling(); // Compute R if needed.
761 if (ierr!=0) return(ierr); // Can't count on R being computed.
762
763 OrdinalType LDB = RHS_->stride(), NRHS = RHS_->numCols();
764 ScalarType * B = RHS_->values();
765 ScalarType * ptr;
766 for (j=0; j<NRHS; j++) {
767 ptr = B + j*LDB;
768 for (i=0; i<numRowCols_; i++) {
769 *ptr = *ptr*R_[i];
770 ptr++;
771 }
772 }
773
774 equilibratedB_ = true;
775
776 return(0);
777}
778
779//=============================================================================
780
781template<typename OrdinalType, typename ScalarType>
782int SerialSpdDenseSolver<OrdinalType,ScalarType>::unequilibrateLHS()
783{
784 OrdinalType i, j;
785
786 if (!equilibratedB_) return(0); // Nothing to do
787
788 OrdinalType LDX = LHS_->stride(), NLHS = LHS_->numCols();
789 ScalarType * X = LHS_->values();
790 ScalarType * ptr;
791 for (j=0; j<NLHS; j++) {
792 ptr = X + j*LDX;
793 for (i=0; i<numRowCols_; i++) {
794 *ptr = *ptr*R_[i];
795 ptr++;
796 }
797 }
798
799 return(0);
800}
801
802//=============================================================================
803
804template<typename OrdinalType, typename ScalarType>
805int SerialSpdDenseSolver<OrdinalType,ScalarType>::invert()
806{
807#ifdef HAVE_TEUCHOSNUMERICS_EIGEN
808 // Implement templated inverse or use Eigen.
809 return (-1);
810#else
811 if (!factored()) factor(); // Need matrix factored.
812
813 INFO_ = 0;
814 this->POTRI( Matrix_->UPLO(), numRowCols_, AF_, LDAF_, &INFO_);
815
816 // Copy lower/upper triangle to upper/lower triangle to make full inverse.
817 if (Matrix_->upper())
818 Matrix_->setLower();
819 else
820 Matrix_->setUpper();
821
822 inverted_ = true;
823 factored_ = false;
824
825 return(INFO_);
826#endif
827}
828
829//=============================================================================
830
831template<typename OrdinalType, typename ScalarType>
832int SerialSpdDenseSolver<OrdinalType,ScalarType>::reciprocalConditionEstimate(MagnitudeType & Value)
833{
834#ifdef HAVE_TEUCHOSNUMERICS_EIGEN
835 // Implement templated GECON or use Eigen. Eigen currently doesn't have a condition estimation function.
836 return (-1);
837#else
838 if (reciprocalConditionEstimated()) {
839 Value = RCOND_;
840 return(0); // Already computed, just return it.
841 }
842
843 if ( ANORM_<ScalarTraits<MagnitudeType>::zero() ) ANORM_ = Matrix_->normOne();
844
845 int ierr = 0;
846 if (!factored()) ierr = factor(); // Need matrix factored.
847 if (ierr!=0) return(ierr);
848
849 allocateWORK();
850 allocateIWORK();
851
852 // We will assume a one-norm condition number
853 INFO_ = 0;
854 std::vector<typename details::lapack_traits<ScalarType>::iwork_type> POCON_WORK( numRowCols_ );
855 this->POCON( Matrix_->UPLO(), numRowCols_, AF_, LDAF_, ANORM_, &RCOND_, &WORK_[0], &POCON_WORK[0], &INFO_);
856 reciprocalConditionEstimated_ = true;
857 Value = RCOND_;
858
859 return(INFO_);
860#endif
861}
862//=============================================================================
863
864template<typename OrdinalType, typename ScalarType>
865void SerialSpdDenseSolver<OrdinalType,ScalarType>::Print(std::ostream& os) const {
866
867 if (Matrix_!=Teuchos::null) os << "Solver Matrix" << std::endl << *Matrix_ << std::endl;
868 if (Factor_!=Teuchos::null) os << "Solver Factored Matrix" << std::endl << *Factor_ << std::endl;
869 if (LHS_ !=Teuchos::null) os << "Solver LHS" << std::endl << *LHS_ << std::endl;
870 if (RHS_ !=Teuchos::null) os << "Solver RHS" << std::endl << *RHS_ << std::endl;
871
872}
873
874} // namespace Teuchos
875
876#endif /* TEUCHOS_SERIALSPDDENSESOLVER_H */
Teuchos_BLAS.hpp
Templated interface class to BLAS routines.
Teuchos_CompObject.hpp
Object for storing data and providing functionality that is common to all computational classes.
Teuchos_ConfigDefs.hpp
Teuchos header file which uses auto-configuration information to include necessary C++ headers.
Teuchos_LAPACK.hpp
Templated interface class to LAPACK routines.
Teuchos_RCP.hpp
Reference-counted pointer class and non-member templated function implementations.
Teuchos_SerialDenseMatrix.hpp
Templated serial dense matrix class.
Teuchos_SerialDenseSolver.hpp
Templated class for solving dense linear problems.
Teuchos_SerialSymDenseMatrix.hpp
Templated serial, dense, symmetric matrix class.
Teuchos::BLAS
Templated BLAS wrapper.
Definition Teuchos_BLAS.hpp:213
Teuchos::CompObject
Functionality and data that is common to all computational classes.
Definition Teuchos_CompObject.hpp:34
Teuchos::LAPACK
The Templated LAPACK Wrapper Class.
Definition Teuchos_LAPACK.hpp:65
Teuchos::Object
The base Teuchos class.
Definition Teuchos_Object.hpp:36
Teuchos::RCP
Smart reference counting pointer class for automatic garbage collection.
Definition Teuchos_RCPDecl.hpp:397
Teuchos::RCP::rcp
RCP< T > rcp(const boost::shared_ptr< T > &sptr)
Conversion function that takes in a boost::shared_ptr object and spits out a Teuchos::RCP object.
Teuchos::RCP::ptr
Ptr< T > ptr() const
Get a safer wrapper raw C++ pointer to the underlying object.
Definition Teuchos_RCP.hpp:382
Teuchos::SerialSpdDenseSolver
A class for constructing and using Hermitian positive definite dense matrices.
Definition Teuchos_SerialSpdDenseSolver.hpp:103
Teuchos::SerialSpdDenseSolver::SCOND
MagnitudeType SCOND()
Ratio of smallest to largest equilibration scale factors for the this matrix (returns -1 if not yet c...
Definition Teuchos_SerialSpdDenseSolver.hpp:292
Teuchos::SerialSpdDenseSolver::setMatrix
int setMatrix(const RCP< SerialSymDenseMatrix< OrdinalType, ScalarType > > &A_in)
Sets the pointers for coefficient matrix.
Definition Teuchos_SerialSpdDenseSolver.hpp:451
Teuchos::SerialSpdDenseSolver::numCols
OrdinalType numCols() const
Returns column dimension of system.
Definition Teuchos_SerialSpdDenseSolver.hpp:281
Teuchos::SerialSpdDenseSolver::FERR
std::vector< MagnitudeType > FERR() const
Returns a pointer to the forward error estimates computed by LAPACK.
Definition Teuchos_SerialSpdDenseSolver.hpp:298
Teuchos::SerialSpdDenseSolver::applyRefinement
int applyRefinement()
Apply Iterative Refinement.
Definition Teuchos_SerialSpdDenseSolver.hpp:627
Teuchos::SerialSpdDenseSolver::transpose
bool transpose()
Returns true if transpose of this matrix has and will be used.
Definition Teuchos_SerialSpdDenseSolver.hpp:232
Teuchos::SerialSpdDenseSolver::solve
int solve()
Computes the solution X to AX = B for the this matrix and the B provided to SetVectors()....
Definition Teuchos_SerialSpdDenseSolver.hpp:546
Teuchos::SerialSpdDenseSolver::R
std::vector< MagnitudeType > R() const
Returns a pointer to the row scaling vector used for equilibration.
Definition Teuchos_SerialSpdDenseSolver.hpp:304
Teuchos::SerialSpdDenseSolver::setVectors
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).
Definition Teuchos_SerialSpdDenseSolver.hpp:466
Teuchos::SerialSpdDenseSolver::factor
int factor()
Computes the in-place Cholesky factorization of the matrix using the LAPACK routine DPOTRF.
Definition Teuchos_SerialSpdDenseSolver.hpp:499
Teuchos::SerialSpdDenseSolver::getRHS
RCP< SerialDenseMatrix< OrdinalType, ScalarType > > getRHS() const
Returns pointer to current RHS.
Definition Teuchos_SerialSpdDenseSolver.hpp:275
Teuchos::SerialSpdDenseSolver::estimateSolutionErrors
void estimateSolutionErrors(bool flag)
Causes all solves to estimate the forward and backward solution error.
Definition Teuchos_SerialSpdDenseSolver.hpp:489
Teuchos::SerialSpdDenseSolver::solutionRefined
bool solutionRefined()
Returns true if the current set of vectors has been refined.
Definition Teuchos_SerialSpdDenseSolver.hpp:259
Teuchos::SerialSpdDenseSolver::~SerialSpdDenseSolver
virtual ~SerialSpdDenseSolver()
SerialSpdDenseSolver destructor.
Definition Teuchos_SerialSpdDenseSolver.hpp:410
Teuchos::SerialSpdDenseSolver::getMatrix
RCP< SerialSymDenseMatrix< OrdinalType, ScalarType > > getMatrix() const
Returns pointer to current matrix.
Definition Teuchos_SerialSpdDenseSolver.hpp:266
Teuchos::SerialSpdDenseSolver::solved
bool solved()
Returns true if the current set of vectors has been solved.
Definition Teuchos_SerialSpdDenseSolver.hpp:256
Teuchos::SerialSpdDenseSolver::reciprocalConditionEstimated
bool reciprocalConditionEstimated()
Returns true if the condition number of the this matrix has been computed (value available via Recipr...
Definition Teuchos_SerialSpdDenseSolver.hpp:253
Teuchos::SerialSpdDenseSolver::getFactoredMatrix
RCP< SerialSymDenseMatrix< OrdinalType, ScalarType > > getFactoredMatrix() const
Returns pointer to factored matrix (assuming factorization has been performed).
Definition Teuchos_SerialSpdDenseSolver.hpp:269
Teuchos::SerialSpdDenseSolver::solveToRefinedSolution
void solveToRefinedSolution(bool flag)
Causes all solves to compute solution to best ability using iterative refinement.
Definition Teuchos_SerialSpdDenseSolver.hpp:156
Teuchos::SerialSpdDenseSolver::numRows
OrdinalType numRows() const
Returns row dimension of system.
Definition Teuchos_SerialSpdDenseSolver.hpp:278
Teuchos::SerialSpdDenseSolver::ANORM
MagnitudeType ANORM() const
Returns the 1-Norm of the this matrix (returns -1 if not yet computed).
Definition Teuchos_SerialSpdDenseSolver.hpp:284
Teuchos::SerialSpdDenseSolver::solutionErrorsEstimated
bool solutionErrorsEstimated()
Returns true if forward and backward error estimated have been computed (available via FERR() and BER...
Definition Teuchos_SerialSpdDenseSolver.hpp:247
Teuchos::SerialSpdDenseSolver::AMAX
MagnitudeType AMAX() const
Returns the absolute value of the largest entry of the this matrix (returns -1 if not yet computed).
Definition Teuchos_SerialSpdDenseSolver.hpp:295
Teuchos::SerialSpdDenseSolver::factored
bool factored()
Returns true if matrix is factored (factor available via AF() and LDAF()).
Definition Teuchos_SerialSpdDenseSolver.hpp:235
Teuchos::SerialSpdDenseSolver::Print
void Print(std::ostream &os) const
Print service methods; defines behavior of ostream << operator.
Definition Teuchos_SerialSpdDenseSolver.hpp:865
Teuchos::SerialSpdDenseSolver::SerialSpdDenseSolver
SerialSpdDenseSolver()
Default constructor; matrix should be set using setMatrix(), LHS and RHS set with setVectors().
Definition Teuchos_SerialSpdDenseSolver.hpp:378
Teuchos::SerialSpdDenseSolver::BERR
std::vector< MagnitudeType > BERR() const
Returns a pointer to the backward error estimates computed by LAPACK.
Definition Teuchos_SerialSpdDenseSolver.hpp:301
Teuchos::SerialSpdDenseSolver::unequilibrateLHS
int unequilibrateLHS()
Unscales the solution vectors if equilibration was used to solve the system.
Definition Teuchos_SerialSpdDenseSolver.hpp:782
Teuchos::SerialSpdDenseSolver::factorWithEquilibration
void factorWithEquilibration(bool flag)
Causes equilibration to be called just before the matrix factorization as part of the call to factor.
Definition Teuchos_SerialSpdDenseSolver.hpp:153
Teuchos::SerialSpdDenseSolver::equilibrateMatrix
int equilibrateMatrix()
Equilibrates the this matrix.
Definition Teuchos_SerialSpdDenseSolver.hpp:681
Teuchos::SerialSpdDenseSolver::reciprocalConditionEstimate
int reciprocalConditionEstimate(MagnitudeType &Value)
Returns the reciprocal of the 1-norm condition number of the this matrix.
Definition Teuchos_SerialSpdDenseSolver.hpp:832
Teuchos::SerialSpdDenseSolver::equilibratedB
bool equilibratedB()
Returns true if RHS is equilibrated (RHS available via B() and LDB()).
Definition Teuchos_SerialSpdDenseSolver.hpp:241
Teuchos::SerialSpdDenseSolver::getLHS
RCP< SerialDenseMatrix< OrdinalType, ScalarType > > getLHS() const
Returns pointer to current LHS.
Definition Teuchos_SerialSpdDenseSolver.hpp:272
Teuchos::SerialSpdDenseSolver::inverted
bool inverted()
Returns true if matrix inverse has been computed (inverse available via AF() and LDAF()).
Definition Teuchos_SerialSpdDenseSolver.hpp:250
Teuchos::SerialSpdDenseSolver::computeEquilibrateScaling
int computeEquilibrateScaling()
Computes the scaling vector S(i) = 1/sqrt(A(i,i) of the this matrix.
Definition Teuchos_SerialSpdDenseSolver.hpp:662
Teuchos::SerialSpdDenseSolver::RCOND
MagnitudeType RCOND() const
Returns the reciprocal of the condition number of the this matrix (returns -1 if not yet computed).
Definition Teuchos_SerialSpdDenseSolver.hpp:287
Teuchos::SerialSpdDenseSolver::equilibrateRHS
int equilibrateRHS()
Equilibrates the current RHS.
Definition Teuchos_SerialSpdDenseSolver.hpp:754
Teuchos::SerialSpdDenseSolver::shouldEquilibrate
bool shouldEquilibrate()
Returns true if the LAPACK general rules for equilibration suggest you should equilibrate the system.
Definition Teuchos_SerialSpdDenseSolver.hpp:244
Teuchos::SerialSpdDenseSolver::equilibratedA
bool equilibratedA()
Returns true if factor is equilibrated (factor available via AF() and LDAF()).
Definition Teuchos_SerialSpdDenseSolver.hpp:238
TEUCHOS_TEST_FOR_EXCEPTION
#define TEUCHOS_TEST_FOR_EXCEPTION(throw_exception_test, Exception, msg)
Macro for throwing an exception with breakpointing to ease debugging.
Definition Teuchos_TestForException.hpp:146
Teuchos
The Teuchos namespace contains all of the classes, structs and enums used by Teuchos,...
details
Teuchos implementation details.
Teuchos::ScalarTraits
This structure defines some basic traits for a scalar field type.
Definition Teuchos_ScalarTraitsDecl.hpp:59
Teuchos::ScalarTraits::one
static T one()
Returns representation of one for this scalar type.
Definition Teuchos_ScalarTraitsDecl.hpp:104
Teuchos::ScalarTraits::magnitudeType
T magnitudeType
Mandatory typedef for result of magnitude.
Definition Teuchos_ScalarTraitsDecl.hpp:61
Generated by doxygen 1.9.8