Belos Version of the Day
Loading...
Searching...
No Matches
BelosLinearProblem.hpp
Go to the documentation of this file.
1
2//@HEADER
3// ************************************************************************
4//
5// Belos: Block Linear Solvers Package
6// Copyright 2004 Sandia Corporation
7//
8// Under the terms of Contract DE-AC04-94AL85000 with Sandia Corporation,
9// the U.S. Government retains certain rights in this software.
10//
11// Redistribution and use in source and binary forms, with or without
12// modification, are permitted provided that the following conditions are
13// met:
14//
15// 1. Redistributions of source code must retain the above copyright
16// notice, this list of conditions and the following disclaimer.
17//
18// 2. Redistributions in binary form must reproduce the above copyright
19// notice, this list of conditions and the following disclaimer in the
20// documentation and/or other materials provided with the distribution.
21//
22// 3. Neither the name of the Corporation nor the names of the
23// contributors may be used to endorse or promote products derived from
24// this software without specific prior written permission.
25//
26// THIS SOFTWARE IS PROVIDED BY SANDIA CORPORATION "AS IS" AND ANY
27// EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
28// IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR
29// PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL SANDIA CORPORATION OR THE
30// CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL,
31// EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO,
32// PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR
33// PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF
34// LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING
35// NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
36// SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
37//
38// Questions? Contact Michael A. Heroux (maherou@sandia.gov)
39//
40// ************************************************************************
41//@HEADER
42
43#ifndef BELOS_LINEAR_PROBLEM_HPP
44#define BELOS_LINEAR_PROBLEM_HPP
45
51#include "Teuchos_ParameterList.hpp"
52#include "Teuchos_TimeMonitor.hpp"
53
54namespace Belos {
55
57
58
62 public:
63 LinearProblemError (const std::string& what_arg) :
64 BelosError(what_arg) {}
65 };
66
68
81 template <class ScalarType, class MV, class OP>
83 public:
84
86
87
94 LinearProblem (void);
95
102 LinearProblem (const Teuchos::RCP<const OP> &A,
103 const Teuchos::RCP<MV> &X,
104 const Teuchos::RCP<const MV> &B);
105
110
112 virtual ~LinearProblem (void) {}
113
115
117
118
122 void setOperator (const Teuchos::RCP<const OP> &A) {
123 A_ = A;
124 isSet_=false;
125 }
126
133 void setLHS (const Teuchos::RCP<MV> &X) {
134 X_ = X;
135 isSet_=false;
136 }
137
142 void setRHS (const Teuchos::RCP<const MV> &B) {
143 B_ = B;
144 isSet_=false;
145 }
146
152 void setInitResVec(const Teuchos::RCP<const MV> &R0) {
153 R0_user_ = R0;
154 isSet_=false;
155 }
156
162 void setInitPrecResVec(const Teuchos::RCP<const MV> &PR0) {
163 PR0_user_ = PR0;
164 isSet_=false;
165 }
166
170 void setLeftPrec(const Teuchos::RCP<const OP> &LP) { LP_ = LP; }
171
175 void setRightPrec(const Teuchos::RCP<const OP> &RP) { RP_ = RP; }
176
185 virtual void setCurrLS ();
186
196 virtual void setLSIndex (const std::vector<int>& index);
197
208 void setHermitian( bool isSym = true ) { isHermitian_ = isSym; }
209
216 void setLabel (const std::string& label);
217
254 virtual
255 Teuchos::RCP<MV>
256 updateSolution (const Teuchos::RCP<MV>& update = Teuchos::null,
257 bool updateLP = false,
258 ScalarType scale = Teuchos::ScalarTraits<ScalarType>::one());
259
277 virtual
278 Teuchos::RCP<MV>
279 updateSolution( const Teuchos::RCP<MV>& update = Teuchos::null,
280 ScalarType scale = Teuchos::ScalarTraits<ScalarType>::one() ) const
281 { return const_cast<LinearProblem<ScalarType,MV,OP> *>(this)->updateSolution( update, false, scale ); }
282
284
286
287
313 virtual
314 bool
315 setProblem (const Teuchos::RCP<MV> &newX = Teuchos::null,
316 const Teuchos::RCP<const MV> &newB = Teuchos::null);
317
319
321
322
324 Teuchos::RCP<const OP> getOperator() const { return(A_); }
325
327 Teuchos::RCP<MV> getLHS() const { return(X_); }
328
330 Teuchos::RCP<const MV> getRHS() const { return(B_); }
331
333 Teuchos::RCP<const MV> getInitResVec() const;
334
339 Teuchos::RCP<const MV> getInitPrecResVec() const;
340
355 Teuchos::RCP<MV> getCurrLHSVec();
356
371 Teuchos::RCP<const MV> getCurrRHSVec();
372
374 Teuchos::RCP<const OP> getLeftPrec() const { return(LP_); };
375
377 Teuchos::RCP<const OP> getRightPrec() const { return(RP_); };
378
400 const std::vector<int>& getLSIndex() const { return(rhsIndex_); }
401
408 int getLSNumber() const { return(lsNum_); }
409
416 Teuchos::Array<Teuchos::RCP<Teuchos::Time> > getTimers() const {
417 return Teuchos::tuple(timerOp_,timerPrec_);
418 }
419
420
422
424
425
433 bool isSolutionUpdated() const { return(solutionUpdated_); }
434
436 bool isProblemSet() const { return(isSet_); }
437
443 bool isHermitian() const { return(isHermitian_); }
444
446 bool isLeftPrec() const { return(LP_!=Teuchos::null); }
447
449 bool isRightPrec() const { return(RP_!=Teuchos::null); }
450
452
454
455
457
464 virtual void apply( const MV& x, MV& y ) const;
465
473 virtual void applyOp( const MV& x, MV& y ) const;
474
481 virtual void applyLeftPrec( const MV& x, MV& y ) const;
482
489 virtual void applyRightPrec( const MV& x, MV& y ) const;
490
492
496 virtual void computeCurrResVec( MV* R , const MV* X = 0, const MV* B = 0 ) const;
497
499
503 virtual void computeCurrPrecResVec( MV* R, const MV* X = 0, const MV* B = 0 ) const;
504
506
507 protected:
508
510 Teuchos::RCP<const OP> A_;
511
513 Teuchos::RCP<MV> X_;
514
516 Teuchos::RCP<MV> curX_;
517
519 Teuchos::RCP<const MV> B_;
520
522 Teuchos::RCP<const MV> curB_;
523
525 Teuchos::RCP<MV> R0_;
526
528 Teuchos::RCP<MV> PR0_;
529
531 Teuchos::RCP<const MV> R0_user_;
532
534 Teuchos::RCP<const MV> PR0_user_;
535
537 Teuchos::RCP<const OP> LP_;
538
540 Teuchos::RCP<const OP> RP_;
541
543 mutable Teuchos::RCP<Teuchos::Time> timerOp_, timerPrec_;
544
547
550
552 std::vector<int> rhsIndex_;
553
556
558
559
561 bool isSet_;
562
566
569
571
573 std::string label_;
574
575 private:
576
579 };
580
581 //--------------------------------------------
582 // Constructor Implementations
583 //--------------------------------------------
584
585 template <class ScalarType, class MV, class OP>
587 blocksize_(0),
588 num2Solve_(0),
589 rhsIndex_(0),
590 lsNum_(0),
591 isSet_(false),
592 isHermitian_(false),
593 solutionUpdated_(false),
594 label_("Belos")
595 {
596 }
597
598 template <class ScalarType, class MV, class OP>
599 LinearProblem<ScalarType,MV,OP>::LinearProblem(const Teuchos::RCP<const OP> &A,
600 const Teuchos::RCP<MV> &X,
601 const Teuchos::RCP<const MV> &B
602 ) :
603 A_(A),
604 X_(X),
605 B_(B),
606 blocksize_(0),
607 num2Solve_(0),
608 rhsIndex_(0),
609 lsNum_(0),
610 isSet_(false),
611 isHermitian_(false),
612 solutionUpdated_(false),
613 label_("Belos")
614 {
615 }
616
617 template <class ScalarType, class MV, class OP>
619 A_(Problem.A_),
620 X_(Problem.X_),
621 curX_(Problem.curX_),
622 B_(Problem.B_),
623 curB_(Problem.curB_),
624 R0_(Problem.R0_),
625 PR0_(Problem.PR0_),
626 R0_user_(Problem.R0_user_),
627 PR0_user_(Problem.PR0_user_),
628 LP_(Problem.LP_),
629 RP_(Problem.RP_),
630 timerOp_(Problem.timerOp_),
631 timerPrec_(Problem.timerPrec_),
632 blocksize_(Problem.blocksize_),
633 num2Solve_(Problem.num2Solve_),
634 rhsIndex_(Problem.rhsIndex_),
635 lsNum_(Problem.lsNum_),
636 isSet_(Problem.isSet_),
637 isHermitian_(Problem.isHermitian_),
638 solutionUpdated_(Problem.solutionUpdated_),
639 label_(Problem.label_)
640 {
641 }
642
643 template <class ScalarType, class MV, class OP>
644 void LinearProblem<ScalarType,MV,OP>::setLSIndex(const std::vector<int>& index)
645 {
646 // Set new linear systems using the indices in index.
647 rhsIndex_ = index;
648
649 // Compute the new block linear system.
650 // ( first clean up old linear system )
651 curB_ = Teuchos::null;
652 curX_ = Teuchos::null;
653
654 // Create indices for the new linear system.
655 int validIdx = 0, ivalidIdx = 0;
656 blocksize_ = rhsIndex_.size();
657 std::vector<int> vldIndex( blocksize_ );
658 std::vector<int> newIndex( blocksize_ );
659 std::vector<int> iIndex( blocksize_ );
660 for (int i=0; i<blocksize_; ++i) {
661 if (rhsIndex_[i] > -1) {
662 vldIndex[validIdx] = rhsIndex_[i];
663 newIndex[validIdx] = i;
664 validIdx++;
665 }
666 else {
667 iIndex[ivalidIdx] = i;
668 ivalidIdx++;
669 }
670 }
671 vldIndex.resize(validIdx);
672 newIndex.resize(validIdx);
673 iIndex.resize(ivalidIdx);
674 num2Solve_ = validIdx;
675
676 // Create the new linear system using index
677 if (num2Solve_ != blocksize_) {
678 newIndex.resize(num2Solve_);
679 vldIndex.resize(num2Solve_);
680 //
681 // First create multivectors of blocksize.
682 // Fill the RHS with random vectors LHS with zero vectors.
683 curX_ = MVT::Clone( *X_, blocksize_ );
684 MVT::MvInit(*curX_);
685 Teuchos::RCP<MV> tmpCurB = MVT::Clone( *B_, blocksize_ );
686 MVT::MvRandom(*tmpCurB);
687 //
688 // Now put in the part of B into curB
689 Teuchos::RCP<const MV> tptr = MVT::CloneView( *B_, vldIndex );
690 MVT::SetBlock( *tptr, newIndex, *tmpCurB );
691 curB_ = tmpCurB;
692 //
693 // Now put in the part of X into curX
694 tptr = MVT::CloneView( *X_, vldIndex );
695 MVT::SetBlock( *tptr, newIndex, *curX_ );
696 //
697 solutionUpdated_ = false;
698 }
699 else {
700 curX_ = MVT::CloneViewNonConst( *X_, rhsIndex_ );
701 curB_ = MVT::CloneView( *B_, rhsIndex_ );
702 }
703 //
704 // Increment the number of linear systems that have been loaded into this object.
705 //
706 lsNum_++;
707 }
708
709
710 template <class ScalarType, class MV, class OP>
712 {
713 //
714 // We only need to copy the solutions back if the linear systems of
715 // interest are less than the block size.
716 //
717 if (num2Solve_ < blocksize_) {
718 //
719 // Get a view of the current solutions and correction std::vector.
720 //
721 int validIdx = 0;
722 std::vector<int> newIndex( num2Solve_ );
723 std::vector<int> vldIndex( num2Solve_ );
724 for (int i=0; i<blocksize_; ++i) {
725 if ( rhsIndex_[i] > -1 ) {
726 vldIndex[validIdx] = rhsIndex_[i];
727 newIndex[validIdx] = i;
728 validIdx++;
729 }
730 }
731 Teuchos::RCP<MV> tptr = MVT::CloneViewNonConst( *curX_, newIndex );
732 MVT::SetBlock( *tptr, vldIndex, *X_ );
733 }
734 //
735 // Clear the current vectors of this linear system so that any future calls
736 // to get the vectors for this system return null pointers.
737 //
738 curX_ = Teuchos::null;
739 curB_ = Teuchos::null;
740 rhsIndex_.resize(0);
741 }
742
743
744 template <class ScalarType, class MV, class OP>
745 Teuchos::RCP<MV>
747 updateSolution (const Teuchos::RCP<MV>& update,
748 bool updateLP,
749 ScalarType scale)
750 {
751 using Teuchos::RCP;
752 using Teuchos::null;
753
754 RCP<MV> newSoln;
755 if (update.is_null())
756 { // The caller didn't supply an update vector, so we assume
757 // that the current solution curX_ is unchanged, and return a
758 // pointer to it.
759 newSoln = curX_;
760 }
761 else // the update vector is NOT null.
762 {
763 if (updateLP) // The caller wants us to update the linear problem.
764 {
765 if (RP_.is_null())
766 { // There is no right preconditioner.
767 // curX_ := curX_ + scale * update.
768 MVT::MvAddMv( 1.0, *curX_, scale, *update, *curX_ );
769 }
770 else
771 { // There is indeed a right preconditioner, so apply it
772 // before computing the new solution.
773 RCP<MV> rightPrecUpdate =
774 MVT::Clone (*update, MVT::GetNumberVecs (*update));
775 applyRightPrec( *update, *rightPrecUpdate );
776 // curX_ := curX_ + scale * rightPrecUpdate.
777 MVT::MvAddMv( 1.0, *curX_, scale, *rightPrecUpdate, *curX_ );
778 }
779 solutionUpdated_ = true;
780 newSoln = curX_;
781 }
782 else
783 { // Rather than updating our stored current solution curX_,
784 // we make a copy and compute the new solution in the
785 // copy, without modifying curX_.
786 newSoln = MVT::Clone (*update, MVT::GetNumberVecs (*update));
787 if (RP_.is_null())
788 { // There is no right preconditioner.
789 // newSoln := curX_ + scale * update.
790 MVT::MvAddMv( 1.0, *curX_, scale, *update, *newSoln );
791 }
792 else
793 { // There is indeed a right preconditioner, so apply it
794 // before computing the new solution.
795 RCP<MV> rightPrecUpdate =
796 MVT::Clone (*update, MVT::GetNumberVecs (*update));
797 applyRightPrec( *update, *rightPrecUpdate );
798 // newSoln := curX_ + scale * rightPrecUpdate.
799 MVT::MvAddMv( 1.0, *curX_, scale, *rightPrecUpdate, *newSoln );
800 }
801 }
802 }
803 return newSoln;
804 }
805
806 template <class ScalarType, class MV, class OP>
807 void LinearProblem<ScalarType,MV,OP>::setLabel(const std::string& label)
808 {
809 if (label != label_) {
810 label_ = label;
811 // Create new timers if they have already been created.
812 if (timerOp_ != Teuchos::null) {
813 std::string opLabel = label_ + ": Operation Op*x";
814#ifdef BELOS_TEUCHOS_TIME_MONITOR
815 timerOp_ = Teuchos::TimeMonitor::getNewCounter( opLabel );
816#endif
817 }
818 if (timerPrec_ != Teuchos::null) {
819 std::string precLabel = label_ + ": Operation Prec*x";
820#ifdef BELOS_TEUCHOS_TIME_MONITOR
821 timerPrec_ = Teuchos::TimeMonitor::getNewCounter( precLabel );
822#endif
823 }
824 }
825 }
826
827 template <class ScalarType, class MV, class OP>
828 bool
830 setProblem (const Teuchos::RCP<MV> &newX,
831 const Teuchos::RCP<const MV> &newB)
832 {
833 // Create timers if the haven't been created yet.
834 if (timerOp_ == Teuchos::null) {
835 std::string opLabel = label_ + ": Operation Op*x";
836#ifdef BELOS_TEUCHOS_TIME_MONITOR
837 timerOp_ = Teuchos::TimeMonitor::getNewCounter( opLabel );
838#endif
839 }
840 if (timerPrec_ == Teuchos::null) {
841 std::string precLabel = label_ + ": Operation Prec*x";
842#ifdef BELOS_TEUCHOS_TIME_MONITOR
843 timerPrec_ = Teuchos::TimeMonitor::getNewCounter( precLabel );
844#endif
845 }
846
847 // Set the linear system using the arguments newX and newB
848 if (newX != Teuchos::null)
849 X_ = newX;
850 if (newB != Teuchos::null)
851 B_ = newB;
852
853 // Invalidate the current linear system indices and multivectors.
854 rhsIndex_.resize(0);
855 curX_ = Teuchos::null;
856 curB_ = Teuchos::null;
857
858 // If we didn't set a matrix A, a left-hand side X, or a
859 // right-hand side B, then we didn't set the problem.
860 if (A_ == Teuchos::null || X_ == Teuchos::null || B_ == Teuchos::null) {
861 isSet_ = false;
862 return isSet_;
863 }
864
865 // Reset whether the solution has been updated. (We're just
866 // setting the problem now, so of course we haven't updated the
867 // solution yet.)
868 solutionUpdated_ = false;
869
870 // Compute the initial residuals.
871 if(Teuchos::is_null(R0_user_)) {
872 if (R0_==Teuchos::null || MVT::GetNumberVecs( *R0_ )!=MVT::GetNumberVecs( *B_ )) {
873 R0_ = MVT::Clone( *B_, MVT::GetNumberVecs( *B_ ) );
874 }
875 computeCurrResVec( &*R0_, &*X_, &*B_ );
876
877 if (LP_!=Teuchos::null) {
878 if (PR0_==Teuchos::null || (PR0_==R0_) || (MVT::GetNumberVecs(*PR0_)!=MVT::GetNumberVecs(*B_))) {
879 PR0_ = MVT::Clone( *B_, MVT::GetNumberVecs( *B_ ) );
880 }
881 applyLeftPrec( *R0_, *PR0_ );
882 }
883 else {
884 PR0_ = R0_;
885 }
886 }
887 else { // User specified the residuals
888 // If the user did not specify the right sized residual, create one and set R0_user_ to null.
889 if (MVT::GetNumberVecs( *R0_user_ )!=MVT::GetNumberVecs( *B_ )) {
890 Teuchos::RCP<MV> helper = MVT::Clone( *B_, MVT::GetNumberVecs( *B_ ) );
891 computeCurrResVec( &*helper, &*X_, &*B_ );
892 R0_user_ = Teuchos::null;
893 R0_ = helper;
894 }
895
896 if (LP_!=Teuchos::null) {
897 // If the user provided preconditioned residual is the wrong size or pointing at
898 // the wrong object, create one and set the PR0_user_ to null.
899 if (PR0_user_==Teuchos::null || (PR0_user_==R0_) || (PR0_user_==R0_user_)
900 || (MVT::GetNumberVecs(*PR0_user_)!=MVT::GetNumberVecs(*B_))) {
901 Teuchos::RCP<MV> helper = MVT::Clone( *B_, MVT::GetNumberVecs( *B_ ) );
902
903 // Get the initial residual from getInitResVec because R0_user_ may be null from above.
904 applyLeftPrec( *getInitResVec(), *helper );
905 PR0_user_ = Teuchos::null;
906 PR0_ = helper;
907 }
908 }
909 else {
910 // The preconditioned initial residual vector is the residual vector.
911 // NOTE: R0_user_ could be null if incompatible.
912 if (R0_user_!=Teuchos::null)
913 {
914 PR0_user_ = R0_user_;
915 }
916 else
917 {
918 PR0_user_ = Teuchos::null;
919 PR0_ = R0_;
920 }
921 }
922 }
923
924 // The problem has been set and is ready for use.
925 isSet_ = true;
926
927 // Return isSet.
928 return isSet_;
929 }
930
931 template <class ScalarType, class MV, class OP>
933 {
934 if(Teuchos::nonnull(R0_user_)) {
935 return R0_user_;
936 }
937 return(R0_);
938 }
939
940 template <class ScalarType, class MV, class OP>
942 {
943 if(Teuchos::nonnull(PR0_user_)) {
944 return PR0_user_;
945 }
946 return(PR0_);
947 }
948
949 template <class ScalarType, class MV, class OP>
951 {
952 if (isSet_) {
953 return curX_;
954 }
955 else {
956 return Teuchos::null;
957 }
958 }
959
960 template <class ScalarType, class MV, class OP>
962 {
963 if (isSet_) {
964 return curB_;
965 }
966 else {
967 return Teuchos::null;
968 }
969 }
970
971 template <class ScalarType, class MV, class OP>
972 void LinearProblem<ScalarType,MV,OP>::apply( const MV& x, MV& y ) const
973 {
974 using Teuchos::null;
975 using Teuchos::RCP;
976
977 const bool leftPrec = LP_ != null;
978 const bool rightPrec = RP_ != null;
979
980 // We only need a temporary vector for intermediate results if
981 // there is a left or right preconditioner. We really should just
982 // keep this temporary vector around instead of allocating it each
983 // time.
984 RCP<MV> ytemp = (leftPrec || rightPrec) ? MVT::Clone (y, MVT::GetNumberVecs (y)) : null;
985
986 //
987 // No preconditioning.
988 //
989 if (!leftPrec && !rightPrec) {
990 applyOp( x, y );
991 }
992 //
993 // Preconditioning is being done on both sides
994 //
995 else if( leftPrec && rightPrec )
996 {
997 applyRightPrec( x, y );
998 applyOp( y, *ytemp );
999 applyLeftPrec( *ytemp, y );
1000 }
1001 //
1002 // Preconditioning is only being done on the left side
1003 //
1004 else if( leftPrec )
1005 {
1006 applyOp( x, *ytemp );
1007 applyLeftPrec( *ytemp, y );
1008 }
1009 //
1010 // Preconditioning is only being done on the right side
1011 //
1012 else
1013 {
1014 applyRightPrec( x, *ytemp );
1015 applyOp( *ytemp, y );
1016 }
1017 }
1018
1019 template <class ScalarType, class MV, class OP>
1020 void LinearProblem<ScalarType,MV,OP>::applyOp( const MV& x, MV& y ) const {
1021 if (A_.get()) {
1022#ifdef BELOS_TEUCHOS_TIME_MONITOR
1023 Teuchos::TimeMonitor OpTimer(*timerOp_);
1024#endif
1025 OPT::Apply( *A_,x, y);
1026 }
1027 else {
1028 MVT::MvAddMv( Teuchos::ScalarTraits<ScalarType>::one(), x,
1029 Teuchos::ScalarTraits<ScalarType>::zero(), x, y );
1030 }
1031 }
1032
1033 template <class ScalarType, class MV, class OP>
1034 void LinearProblem<ScalarType,MV,OP>::applyLeftPrec( const MV& x, MV& y ) const {
1035 if (LP_!=Teuchos::null) {
1036#ifdef BELOS_TEUCHOS_TIME_MONITOR
1037 Teuchos::TimeMonitor PrecTimer(*timerPrec_);
1038#endif
1039 OPT::Apply( *LP_,x, y);
1040 }
1041 else {
1042 MVT::MvAddMv( Teuchos::ScalarTraits<ScalarType>::one(), x,
1043 Teuchos::ScalarTraits<ScalarType>::zero(), x, y );
1044 }
1045 }
1046
1047 template <class ScalarType, class MV, class OP>
1048 void LinearProblem<ScalarType,MV,OP>::applyRightPrec( const MV& x, MV& y ) const {
1049 if (RP_!=Teuchos::null) {
1050#ifdef BELOS_TEUCHOS_TIME_MONITOR
1051 Teuchos::TimeMonitor PrecTimer(*timerPrec_);
1052#endif
1053 OPT::Apply( *RP_,x, y);
1054 }
1055 else {
1056 MVT::MvAddMv( Teuchos::ScalarTraits<ScalarType>::one(), x,
1057 Teuchos::ScalarTraits<ScalarType>::zero(), x, y );
1058 }
1059 }
1060
1061 template <class ScalarType, class MV, class OP>
1062 void LinearProblem<ScalarType,MV,OP>::computeCurrPrecResVec( MV* R, const MV* X, const MV* B ) const {
1063
1064 if (R) {
1065 if (X && B) // The entries are specified, so compute the residual of Op(A)X = B
1066 {
1067 if (LP_!=Teuchos::null)
1068 {
1069 Teuchos::RCP<MV> R_temp = MVT::Clone( *B, MVT::GetNumberVecs( *B ) );
1070 applyOp( *X, *R_temp );
1071 MVT::MvAddMv( -1.0, *R_temp, 1.0, *B, *R_temp );
1072 applyLeftPrec( *R_temp, *R );
1073 }
1074 else
1075 {
1076 applyOp( *X, *R );
1077 MVT::MvAddMv( -1.0, *R, 1.0, *B, *R );
1078 }
1079 }
1080 else {
1081 // The solution and right-hand side may not be specified, check and use which ones exist.
1082 Teuchos::RCP<const MV> localB, localX;
1083 if (B)
1084 localB = Teuchos::rcp( B, false );
1085 else
1086 localB = curB_;
1087
1088 if (X)
1089 localX = Teuchos::rcp( X, false );
1090 else
1091 localX = curX_;
1092
1093 if (LP_!=Teuchos::null)
1094 {
1095 Teuchos::RCP<MV> R_temp = MVT::Clone( *localB, MVT::GetNumberVecs( *localB ) );
1096 applyOp( *localX, *R_temp );
1097 MVT::MvAddMv( -1.0, *R_temp, 1.0, *localB, *R_temp );
1098 applyLeftPrec( *R_temp, *R );
1099 }
1100 else
1101 {
1102 applyOp( *localX, *R );
1103 MVT::MvAddMv( -1.0, *R, 1.0, *localB, *R );
1104 }
1105 }
1106 }
1107 }
1108
1109
1110 template <class ScalarType, class MV, class OP>
1111 void LinearProblem<ScalarType,MV,OP>::computeCurrResVec( MV* R, const MV* X, const MV* B ) const {
1112
1113 if (R) {
1114 if (X && B) // The entries are specified, so compute the residual of Op(A)X = B
1115 {
1116 applyOp( *X, *R );
1117 MVT::MvAddMv( -1.0, *R, 1.0, *B, *R );
1118 }
1119 else {
1120 // The solution and right-hand side may not be specified, check and use which ones exist.
1121 Teuchos::RCP<const MV> localB, localX;
1122 if (B)
1123 localB = Teuchos::rcp( B, false );
1124 else
1125 localB = curB_;
1126
1127 if (X)
1128 localX = Teuchos::rcp( X, false );
1129 else
1130 localX = curX_;
1131
1132 applyOp( *localX, *R );
1133 MVT::MvAddMv( -1.0, *R, 1.0, *localB, *R );
1134 }
1135 }
1136 }
1137
1138} // end Belos namespace
1139
1140#endif /* BELOS_LINEAR_PROBLEM_HPP */
1141
1142
Declaration of basic traits for the multivector type.
Class which defines basic traits for the operator type.
Parent class to all Belos exceptions.
Definition: BelosTypes.hpp:60
Exception thrown to signal error with the Belos::LinearProblem object.
LinearProblemError(const std::string &what_arg)
A linear system to solve, and its associated information.
Teuchos::RCP< const MV > PR0_user_
User-defined preconditioned initial residual of the linear system.
std::vector< int > rhsIndex_
Indices of current linear systems being solver for.
bool isSet_
Has the linear problem to solve been set?
void setLabel(const std::string &label)
Set the label prefix used by the timers in this object.
bool isRightPrec() const
Whether the linear system is being preconditioned on the right.
void setInitPrecResVec(const Teuchos::RCP< const MV > &PR0)
Set the user-defined preconditioned residual of linear problem .
int blocksize_
Current block size of linear system.
virtual void applyLeftPrec(const MV &x, MV &y) const
Apply ONLY the left preconditioner to x, returning y.
bool isSolutionUpdated() const
Has the current approximate solution been updated?
std::string label_
Linear problem label that prefixes the timer labels.
Teuchos::RCP< const MV > curB_
Current right-hand side of the linear system.
void setLeftPrec(const Teuchos::RCP< const OP > &LP)
Set left preconditioner (LP) of linear problem .
Teuchos::RCP< const MV > B_
Right-hand side of linear system.
virtual void computeCurrPrecResVec(MV *R, const MV *X=0, const MV *B=0) const
Compute a residual R for this operator given a solution X, and right-hand side B.
virtual ~LinearProblem(void)
Destructor (declared virtual for memory safety of derived classes).
virtual void apply(const MV &x, MV &y) const
Apply the composite operator of this linear problem to x, returning y.
void setLHS(const Teuchos::RCP< MV > &X)
Set left-hand-side X of linear problem .
virtual void setLSIndex(const std::vector< int > &index)
Tell the linear problem which linear system(s) need to be solved next.
bool solutionUpdated_
Has the current approximate solution been updated?
Teuchos::RCP< MV > PR0_
Preconditioned initial residual of the linear system.
virtual void applyRightPrec(const MV &x, MV &y) const
Apply ONLY the right preconditioner to x, returning y.
Teuchos::RCP< const OP > LP_
Left preconditioning operator of linear system.
Teuchos::RCP< const OP > getLeftPrec() const
Get a pointer to the left preconditioner.
Teuchos::RCP< const MV > R0_user_
User-defined initial residual of the linear system.
Teuchos::RCP< MV > curX_
Current solution vector of the linear system.
Teuchos::RCP< MV > R0_
Initial residual of the linear system.
void setRightPrec(const Teuchos::RCP< const OP > &RP)
Set right preconditioner (RP) of linear problem .
virtual void setCurrLS()
Tell the linear problem that the solver is finished with the current linear system.
int getLSNumber() const
The number of linear systems that have been set.
void setHermitian(bool isSym=true)
Tell the linear problem that the (preconditioned) operator is Hermitian.
void setOperator(const Teuchos::RCP< const OP > &A)
Set the operator A of the linear problem .
bool isProblemSet() const
Whether the problem has been set.
int num2Solve_
Number of linear systems that are currently being solver for ( <= blocksize_ )
Teuchos::RCP< Teuchos::Time > timerPrec_
Teuchos::RCP< const MV > getInitPrecResVec() const
A pointer to the preconditioned initial residual vector.
Teuchos::RCP< const OP > getRightPrec() const
Get a pointer to the right preconditioner.
Teuchos::RCP< const OP > getOperator() const
A pointer to the (unpreconditioned) operator A.
Teuchos::RCP< Teuchos::Time > timerOp_
Timers.
const std::vector< int > & getLSIndex() const
(Zero-based) indices of the linear system(s) currently being solved.
Teuchos::Array< Teuchos::RCP< Teuchos::Time > > getTimers() const
The timers for this object.
bool isHermitian_
Whether the operator A is symmetric (in real arithmetic, or Hermitian in complex arithmetic).
LinearProblem(void)
Default constructor.
void setRHS(const Teuchos::RCP< const MV > &B)
Set right-hand-side B of linear problem .
Teuchos::RCP< const MV > getRHS() const
A pointer to the right-hand side B.
bool isHermitian() const
Whether the (preconditioned) operator is Hermitian.
Teuchos::RCP< MV > getLHS() const
A pointer to the left-hand side X.
virtual void applyOp(const MV &x, MV &y) const
Apply ONLY the operator to x, returning y.
Teuchos::RCP< const MV > getInitResVec() const
A pointer to the initial unpreconditioned residual vector.
Teuchos::RCP< MV > getCurrLHSVec()
Get a pointer to the current left-hand side (solution) of the linear system.
bool isLeftPrec() const
Whether the linear system is being preconditioned on the left.
Teuchos::RCP< const MV > getCurrRHSVec()
Get a pointer to the current right-hand side of the linear system.
virtual void computeCurrResVec(MV *R, const MV *X=0, const MV *B=0) const
Compute a residual R for this operator given a solution X, and right-hand side B.
int lsNum_
Number of linear systems that have been loaded in this linear problem object.
Teuchos::RCP< MV > X_
Solution vector of linear system.
Teuchos::RCP< const OP > A_
Operator of linear system.
virtual Teuchos::RCP< MV > updateSolution(const Teuchos::RCP< MV > &update=Teuchos::null, bool updateLP=false, ScalarType scale=Teuchos::ScalarTraits< ScalarType >::one())
Compute the new solution to the linear system using the given update vector.
void setInitResVec(const Teuchos::RCP< const MV > &R0)
Set the user-defined residual of linear problem .
virtual Teuchos::RCP< MV > updateSolution(const Teuchos::RCP< MV > &update=Teuchos::null, ScalarType scale=Teuchos::ScalarTraits< ScalarType >::one()) const
Compute the new solution to the linear system using the given update vector.
virtual bool setProblem(const Teuchos::RCP< MV > &newX=Teuchos::null, const Teuchos::RCP< const MV > &newB=Teuchos::null)
Set up the linear problem manager.
Teuchos::RCP< const OP > RP_
Right preconditioning operator of linear system.
Traits class which defines basic operations on multivectors.
Class which defines basic traits for the operator type.

Generated for Belos by doxygen 1.9.6