Thyra Package Browser (Single Doxygen Collection) Version of the Day
Thyra_DefaultMultiVectorLinearOpWithSolve_def.hpp
Go to the documentation of this file.
00001 // @HEADER
00002 // ***********************************************************************
00003 // 
00004 //    Thyra: Interfaces and Support for Abstract Numerical Algorithms
00005 //                 Copyright (2004) Sandia Corporation
00006 // 
00007 // Under terms of Contract DE-AC04-94AL85000, there is a non-exclusive
00008 // license for use of this work by or on behalf of the U.S. Government.
00009 // 
00010 // This library is free software; you can redistribute it and/or modify
00011 // it under the terms of the GNU Lesser General Public License as
00012 // published by the Free Software Foundation; either version 2.1 of the
00013 // License, or (at your option) any later version.
00014 //  
00015 // This library is distributed in the hope that it will be useful, but
00016 // WITHOUT ANY WARRANTY; without even the implied warranty of
00017 // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
00018 // Lesser General Public License for more details.
00019 //  
00020 // You should have received a copy of the GNU Lesser General Public
00021 // License along with this library; if not, write to the Free Software
00022 // Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307
00023 // USA
00024 // Questions? Contact Michael A. Heroux (maherou@sandia.gov) 
00025 // 
00026 // ***********************************************************************
00027 // @HEADER
00028 
00029 #ifndef THYRA_MULTI_VECTOR_LINEAR_OP_WITH_SOLVE_HPP
00030 #define THYRA_MULTI_VECTOR_LINEAR_OP_WITH_SOLVE_HPP
00031 
00032 #include "Thyra_DefaultMultiVectorLinearOpWithSolve_decl.hpp"
00033 #include "Thyra_DefaultDiagonalLinearOp.hpp"
00034 #include "Thyra_LinearOpWithSolveBase.hpp"
00035 #include "Thyra_DefaultMultiVectorProductVectorSpace.hpp"
00036 #include "Thyra_DefaultMultiVectorProductVector.hpp"
00037 #include "Thyra_AssertOp.hpp"
00038 #include "Teuchos_dyn_cast.hpp"
00039 
00040 
00041 namespace Thyra {
00042 
00043 
00044 // Constructors/initializers/accessors
00045 
00046 
00047 template<class Scalar>
00048 DefaultMultiVectorLinearOpWithSolve<Scalar>::DefaultMultiVectorLinearOpWithSolve()
00049 {}
00050 
00051 
00052 template<class Scalar>
00053 void DefaultMultiVectorLinearOpWithSolve<Scalar>::nonconstInitialize(
00054   const RCP<LinearOpWithSolveBase<Scalar> > &lows,
00055   const RCP<const DefaultMultiVectorProductVectorSpace<Scalar> > &multiVecRange,
00056   const RCP<const DefaultMultiVectorProductVectorSpace<Scalar> > &multiVecDomain
00057   )
00058 {
00059   validateInitialize(lows,multiVecRange,multiVecDomain);
00060   lows_ = lows;
00061   multiVecRange_ = multiVecRange;
00062   multiVecDomain_ = multiVecDomain;
00063 }
00064 
00065 
00066 template<class Scalar>
00067 void DefaultMultiVectorLinearOpWithSolve<Scalar>::initialize(
00068   const RCP<const LinearOpWithSolveBase<Scalar> > &lows,
00069   const RCP<const DefaultMultiVectorProductVectorSpace<Scalar> > &multiVecRange,
00070   const RCP<const DefaultMultiVectorProductVectorSpace<Scalar> > &multiVecDomain
00071   )
00072 {
00073   validateInitialize(lows,multiVecRange,multiVecDomain);
00074   lows_ = lows;
00075   multiVecRange_ = multiVecRange;
00076   multiVecDomain_ = multiVecDomain;
00077 }
00078 
00079 
00080 template<class Scalar>
00081 RCP<LinearOpWithSolveBase<Scalar> >
00082 DefaultMultiVectorLinearOpWithSolve<Scalar>::getNonconstLinearOpWithSolve()
00083 {
00084   return lows_.getNonconstObj();
00085 }
00086 
00087 
00088 template<class Scalar>
00089 RCP<const LinearOpWithSolveBase<Scalar> >
00090 DefaultMultiVectorLinearOpWithSolve<Scalar>::getLinearOpWithSolve() const
00091 {
00092   return lows_.getConstObj();
00093 }
00094 
00095 
00096 template<class Scalar>
00097 void DefaultMultiVectorLinearOpWithSolve<Scalar>::uninitialize()
00098 {
00099   lows_.uninitialize();
00100   multiVecRange_ = Teuchos::null;
00101   multiVecDomain_ = Teuchos::null;
00102 }
00103 
00104 
00105 // Overridden from LinearOpBase
00106 
00107 
00108 template<class Scalar>
00109 RCP< const VectorSpaceBase<Scalar> >
00110 DefaultMultiVectorLinearOpWithSolve<Scalar>::range() const
00111 {
00112   return multiVecRange_;
00113 }
00114 
00115 
00116 template<class Scalar>
00117 RCP< const VectorSpaceBase<Scalar> >
00118 DefaultMultiVectorLinearOpWithSolve<Scalar>::domain() const
00119 {
00120   return multiVecDomain_;
00121 }
00122 
00123 
00124 template<class Scalar>
00125 RCP<const LinearOpBase<Scalar> >
00126 DefaultMultiVectorLinearOpWithSolve<Scalar>::clone() const
00127 {
00128   return Teuchos::null; // ToDo: Implement if needed ???
00129 }
00130 
00131 
00132 // protected
00133 
00134 
00135 // Overridden from LinearOpBase
00136 
00137 
00138 template<class Scalar>
00139 bool DefaultMultiVectorLinearOpWithSolve<Scalar>::opSupportedImpl(
00140   EOpTransp M_trans
00141   ) const
00142 {
00143   return Thyra::opSupported(*lows_.getConstObj(),M_trans);
00144 }
00145 
00146 
00147 template<class Scalar>
00148 void DefaultMultiVectorLinearOpWithSolve<Scalar>::applyImpl(
00149   const EOpTransp M_trans,
00150   const MultiVectorBase<Scalar> &XX,
00151   const Ptr<MultiVectorBase<Scalar> > &YY,
00152   const Scalar alpha,
00153   const Scalar beta
00154   ) const
00155 {
00156 
00157   using Teuchos::dyn_cast;
00158   typedef DefaultMultiVectorProductVector<Scalar> MVPV;
00159 
00160   const Ordinal numCols = XX.domain()->dim();
00161 
00162   for (Ordinal col_j = 0; col_j < numCols; ++col_j) {
00163 
00164     const RCP<const VectorBase<Scalar> > x = XX.col(col_j);
00165     const RCP<VectorBase<Scalar> > y = YY->col(col_j);
00166  
00167     RCP<const MultiVectorBase<Scalar> >
00168       X = dyn_cast<const MVPV>(*x).getMultiVector().assert_not_null();
00169     RCP<MultiVectorBase<Scalar> >
00170       Y = dyn_cast<MVPV>(*y).getNonconstMultiVector().assert_not_null();
00171     
00172     Thyra::apply( *lows_.getConstObj(), M_trans, *X, Y.ptr(), alpha, beta );
00173 
00174   }
00175 
00176 }
00177 
00178 
00179 // Overridden from LinearOpWithSolveBase
00180 
00181 
00182 template<class Scalar>
00183 bool
00184 DefaultMultiVectorLinearOpWithSolve<Scalar>::solveSupportsImpl(
00185   EOpTransp M_trans
00186   ) const
00187 {
00188   return Thyra::solveSupports(*lows_.getConstObj(),M_trans);
00189 }
00190 
00191 
00192 template<class Scalar>
00193 bool
00194 DefaultMultiVectorLinearOpWithSolve<Scalar>::solveSupportsSolveMeasureTypeImpl(
00195   EOpTransp M_trans, const SolveMeasureType& solveMeasureType
00196   ) const
00197 {
00198   return Thyra::solveSupportsSolveMeasureType(
00199     *lows_.getConstObj(),M_trans,solveMeasureType);
00200 }
00201 
00202 
00203 template<class Scalar>
00204 SolveStatus<Scalar>
00205 DefaultMultiVectorLinearOpWithSolve<Scalar>::solveImpl(
00206   const EOpTransp transp,
00207   const MultiVectorBase<Scalar> &BB,
00208   const Ptr<MultiVectorBase<Scalar> > &XX,
00209   const Ptr<const SolveCriteria<Scalar> > solveCriteria
00210   ) const
00211 {
00212 
00213   using Teuchos::dyn_cast;
00214   using Teuchos::outArg;
00215   using Teuchos::inOutArg;
00216   typedef DefaultMultiVectorProductVector<Scalar> MVPV;
00217 
00218   const Ordinal numCols = BB.domain()->dim();
00219 
00220   SolveStatus<Scalar> overallSolveStatus;
00221   accumulateSolveStatusInit(outArg(overallSolveStatus));
00222   
00223   for (Ordinal col_j = 0; col_j < numCols; ++col_j) {
00224 
00225     const RCP<const VectorBase<Scalar> > b = BB.col(col_j);
00226     const RCP<VectorBase<Scalar> > x = XX->col(col_j);
00227 
00228     RCP<const MultiVectorBase<Scalar> >
00229       B = dyn_cast<const MVPV>(*b).getMultiVector().assert_not_null();
00230     RCP<MultiVectorBase<Scalar> >
00231       X = dyn_cast<MVPV>(*x).getNonconstMultiVector().assert_not_null();
00232 
00233     const SolveStatus<Scalar> solveStatus =
00234       Thyra::solve(*lows_.getConstObj(), transp, *B, X.ptr(), solveCriteria);
00235 
00236     accumulateSolveStatus(
00237       SolveCriteria<Scalar>(), // Never used
00238       solveStatus, inOutArg(overallSolveStatus) );
00239 
00240   }
00241   
00242   return overallSolveStatus;
00243 
00244 }
00245 
00246 
00247 // private
00248 
00249 
00250 template<class Scalar>
00251 void DefaultMultiVectorLinearOpWithSolve<Scalar>::validateInitialize(
00252   const RCP<const LinearOpWithSolveBase<Scalar> > &lows,
00253   const RCP<const DefaultMultiVectorProductVectorSpace<Scalar> > &multiVecRange,
00254   const RCP<const DefaultMultiVectorProductVectorSpace<Scalar> > &multiVecDomain
00255   )
00256 {
00257 #ifdef TEUCHOS_DEBUG
00258   TEST_FOR_EXCEPT(is_null(lows));
00259   TEST_FOR_EXCEPT(is_null(multiVecRange));
00260   TEST_FOR_EXCEPT(is_null(multiVecDomain));
00261   TEST_FOR_EXCEPT( multiVecRange->numBlocks() != multiVecDomain->numBlocks() );
00262   THYRA_ASSERT_VEC_SPACES(
00263     "DefaultMultiVectorLinearOpWithSolve<Scalar>::initialize(lows,multiVecRange,multiVecDomain)",
00264     *lows->range(), *multiVecRange->getBlock(0) );
00265   THYRA_ASSERT_VEC_SPACES(
00266     "DefaultMultiVectorLinearOpWithSolve<Scalar>::initialize(lows,multiVecRange,multiVecDomain)",
00267     *lows->domain(), *multiVecDomain->getBlock(0) );
00268 #endif
00269 }
00270 
00271 
00272 } // end namespace Thyra
00273 
00274 
00275 #endif  // THYRA_MULTI_VECTOR_LINEAR_OP_WITH_SOLVE_HPP
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Defines