RealComplexFFTLinearOp.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_REAL_COMPLEX_FFT_LINEAR_OP_HPP
00030 #define THYRA_REAL_COMPLEX_FFT_LINEAR_OP_HPP
00031 
00032 #include "ComplexFFTLinearOp.hpp"
00033 #include "Thyra_DefaultSpmdVectorSpace.hpp"
00034 #include "Thyra_DetachedVectorView.hpp"
00035 #include "Thyra_DefaultSerialVectorSpaceConverter.hpp"
00036 #include "serial_1D_FFT.hpp"
00037 
00045 template<class RealScalar>
00046 class RealComplexFFTLinearOp : virtual public Thyra::LinearOpWithSolveBase< std::complex<RealScalar>, RealScalar >
00047 {
00048 public:
00049 
00051   typedef std::complex<RealScalar>                                                 RangeScalar;
00053   typedef RealScalar                                                               DomainScalar;
00055   typedef typename Teuchos::PromotionTraits<RangeScalar,DomainScalar>::promote     Scalar;
00056 
00058   RealComplexFFTLinearOp( const int N );
00059 
00062 
00064   Teuchos::RefCountPtr< const Thyra::VectorSpaceBase<RangeScalar> > range() const;
00066   Teuchos::RefCountPtr< const Thyra::VectorSpaceBase<DomainScalar> > domain() const;
00068   bool applySupports( const Thyra::EConj conj ) const;
00070   bool applyTransposeSupports( const Thyra::EConj conj ) const;
00072   void apply(
00073     const Thyra::EConj                            conj
00074     ,const Thyra::MultiVectorBase<DomainScalar>   &X
00075     ,Thyra::MultiVectorBase<RangeScalar>          *Y
00076     ,const RangeScalar                            alpha
00077     ,const RangeScalar                           beta
00078     ) const;
00080   void applyTranspose(
00081     const Thyra::EConj                            conj
00082     ,const Thyra::MultiVectorBase<RangeScalar>    &X
00083     ,Thyra::MultiVectorBase<DomainScalar>         *Y
00084     ,const DomainScalar                           alpha
00085     ,const DomainScalar                           beta
00086     ) const;
00087 
00089 
00092 
00094   bool solveSupportsConj(Thyra::EConj conj) const;
00096   bool solveTransposeSupportsConj(Thyra::EConj conj) const;
00098   bool solveSupportsSolveMeasureType(Thyra::EConj conj, const Thyra::SolveMeasureType solveMeasureType) const;
00100   bool solveTransposeSupportsSolveMeasureType(Thyra::EConj conj, const Thyra::SolveMeasureType solveMeasureType) const;
00102   void solve(
00103     const Thyra::EConj                           conj
00104     ,const Thyra::MultiVectorBase<RangeScalar>   &B
00105     ,Thyra::MultiVectorBase<DomainScalar>        *X
00106     ,const int                                   numBlocks
00107     ,const Thyra::BlockSolveCriteria<Scalar>     blockSolveCriteria[]
00108     ,Thyra::SolveStatus<Scalar>                  blockSolveStatus[]
00109     ) const;
00111   void solveTranspose(
00112     const Thyra::EConj                           conj
00113     ,const Thyra::MultiVectorBase<DomainScalar>  &B
00114     ,Thyra::MultiVectorBase<RangeScalar>         *X
00115     ,const int                                   numBlocks
00116     ,const Thyra::BlockSolveCriteria<Scalar>     blockSolveCriteria[]
00117     ,Thyra::SolveStatus<Scalar>                  blockSolveStatus[]
00118     ) const;
00119 
00121 
00122 private:
00123 
00124   Thyra::DefaultSerialVectorSpaceConverter<DomainScalar,RangeScalar>      realToComplexConverter_;
00125   ComplexFFTLinearOp<RealScalar>                                          complexFFTOp_;
00126   Teuchos::RefCountPtr< const Thyra::VectorSpaceBase< RealScalar > >      domain_;
00127   
00129   static void updateReal(
00130     const Thyra::MultiVectorBase<RangeScalar>      &V_complex
00131     ,const DomainScalar                            &beta
00132     ,Thyra::MultiVectorBase<DomainScalar>          *V
00133     );
00134 
00135 };
00136 
00137 // /////////////////////////
00138 // Definitions
00139 
00140 template<class RealScalar>
00141 RealComplexFFTLinearOp<RealScalar>::RealComplexFFTLinearOp( const int N )
00142   :complexFFTOp_(N), domain_(realToComplexConverter_.createVectorSpaceFrom(*complexFFTOp_.domain()))
00143 {}
00144 
00145 // Overridden from LinearOpBase
00146 
00147 template<class RealScalar>
00148 Teuchos::RefCountPtr< const Thyra::VectorSpaceBase<typename RealComplexFFTLinearOp<RealScalar>::RangeScalar > >
00149 RealComplexFFTLinearOp<RealScalar>::range() const
00150 {
00151   return complexFFTOp_.range();
00152 }
00153 
00154 template<class RealScalar>
00155 Teuchos::RefCountPtr< const Thyra::VectorSpaceBase<typename RealComplexFFTLinearOp<RealScalar>::DomainScalar > >
00156 RealComplexFFTLinearOp<RealScalar>::domain() const
00157 {
00158   return domain_;
00159 }
00160 
00161 template<class RealScalar>
00162 bool RealComplexFFTLinearOp<RealScalar>::applySupports( const Thyra::EConj conj ) const
00163 {
00164   return (conj == Thyra::NONCONJ_ELE);
00165 }
00166 
00167 template<class RealScalar>
00168 bool RealComplexFFTLinearOp<RealScalar>::applyTransposeSupports( const Thyra::EConj conj ) const
00169 {
00170   return (conj == Thyra::CONJ_ELE);
00171 }
00172 
00173 template<class RealScalar>
00174 void RealComplexFFTLinearOp<RealScalar>::apply(
00175   const Thyra::EConj                            conj
00176   ,const Thyra::MultiVectorBase<DomainScalar>   &X
00177   ,Thyra::MultiVectorBase<RangeScalar>          *Y
00178   ,const RangeScalar                            alpha
00179   ,const RangeScalar                            beta
00180   ) const
00181 {
00182   Teuchos::RefCountPtr< Thyra::MultiVectorBase<RangeScalar> >
00183     X_complex = Thyra::createMembers(complexFFTOp_.domain(),X.domain()->dim());
00184   realToComplexConverter_.convert( X, &*X_complex );
00185   Thyra::apply(complexFFTOp_,conj,*X_complex,Y,alpha,beta);
00186 }
00187 
00188 template<class RealScalar>
00189 void RealComplexFFTLinearOp<RealScalar>::applyTranspose(
00190   const Thyra::EConj                            conj
00191   ,const Thyra::MultiVectorBase<RangeScalar>    &X
00192   ,Thyra::MultiVectorBase<DomainScalar>         *Y
00193   ,const DomainScalar                           alpha
00194   ,const DomainScalar                           beta
00195   ) const
00196 {
00197   Teuchos::RefCountPtr< Thyra::MultiVectorBase<RangeScalar> >
00198     Y_complex = Thyra::createMembers(complexFFTOp_.range(),X.domain()->dim());
00199   Thyra::applyTranspose(complexFFTOp_,conj,X,&*Y_complex,RangeScalar(alpha));
00200   updateReal( *Y_complex, beta, Y );
00201 }
00202 
00203 // Overridden from LinearOpWithSolveBase
00204 
00205 template<class RealScalar>
00206 bool RealComplexFFTLinearOp<RealScalar>::solveSupportsConj(Thyra::EConj conj) const
00207 {
00208   return (conj == Thyra::NONCONJ_ELE);
00209 }
00210 
00211 template<class RealScalar>
00212 bool RealComplexFFTLinearOp<RealScalar>::solveTransposeSupportsConj(Thyra::EConj conj) const
00213 {
00214   return (conj == Thyra::CONJ_ELE);
00215 }
00216 
00217 template<class RealScalar>
00218 bool RealComplexFFTLinearOp<RealScalar>::solveSupportsSolveMeasureType(Thyra::EConj conj, const Thyra::SolveMeasureType solveMeasureType) const
00219 {
00220   return (conj == Thyra::NONCONJ_ELE);
00221 }
00222 
00223 template<class RealScalar>
00224 bool RealComplexFFTLinearOp<RealScalar>::solveTransposeSupportsSolveMeasureType(Thyra::EConj conj, const Thyra::SolveMeasureType solveMeasureType) const
00225 {
00226   return (conj == Thyra::CONJ_ELE);
00227 }
00228 
00229 template<class RealScalar>
00230 void RealComplexFFTLinearOp<RealScalar>::solve(
00231   const Thyra::EConj                           conj
00232   ,const Thyra::MultiVectorBase<RangeScalar>   &B
00233   ,Thyra::MultiVectorBase<DomainScalar>        *X
00234   ,const int                                   numBlocks
00235   ,const Thyra::BlockSolveCriteria<Scalar>     blockSolveCriteria[]
00236   ,Thyra::SolveStatus<Scalar>                  blockSolveStatus[]
00237   ) const
00238 {
00239   Teuchos::RefCountPtr< Thyra::MultiVectorBase<RangeScalar> >
00240     X_complex = Thyra::createMembers(complexFFTOp_.domain(),B.domain()->dim());
00241   Thyra::solve(complexFFTOp_,conj,B,&*X_complex,numBlocks,blockSolveCriteria,blockSolveStatus);
00242   updateReal( *X_complex, Teuchos::ScalarTraits<DomainScalar>::zero(), X );
00243 }
00244 
00245 template<class RealScalar>
00246 void RealComplexFFTLinearOp<RealScalar>::solveTranspose(
00247   const Thyra::EConj                           conj
00248   ,const Thyra::MultiVectorBase<DomainScalar>  &B
00249   ,Thyra::MultiVectorBase<RangeScalar>         *X
00250   ,const int                                   numBlocks
00251   ,const Thyra::BlockSolveCriteria<Scalar>     blockSolveCriteria[]
00252   ,Thyra::SolveStatus<Scalar>                  blockSolveStatus[]
00253   ) const
00254 {
00255   Teuchos::RefCountPtr< Thyra::MultiVectorBase<RangeScalar> >
00256     B_complex = Thyra::createMembers(complexFFTOp_.domain(),B.domain()->dim());
00257   realToComplexConverter_.convert( B, &*B_complex );
00258   Thyra::solveTranspose(complexFFTOp_,conj,*B_complex,X,numBlocks,blockSolveCriteria,blockSolveStatus);
00259 }
00260 
00261 // private
00262 
00263 template<class RealScalar>
00264 void RealComplexFFTLinearOp<RealScalar>::updateReal(
00265   const Thyra::MultiVectorBase<RangeScalar>      &V_complex
00266   ,const DomainScalar                            &beta
00267   ,Thyra::MultiVectorBase<DomainScalar>          *V
00268   )
00269 {
00270 #ifdef TEUCHOS_DEBUG
00271   TEST_FOR_EXCEPT( V==NULL );
00272 #endif  
00273   typedef Teuchos::ScalarTraits<Scalar> ST;
00274   Thyra::ConstDetachedMultiVectorView<RangeScalar>           ev_V_complex(V_complex);
00275   Thyra::DetachedMultiVectorView<DomainScalar>   ev_V(*V);
00276 #ifdef TEUCHOS_DEBUG
00277   TEST_FOR_EXCEPT( ev_V_complex.subDim() != ev_V.subDim() || ev_V_complex.numSubCols() != ev_V.numSubCols() );
00278 #endif  
00279   if( beta == ST::zero() ) {
00280     for( Thyra::Index j = 0; j < ev_V.numSubCols(); ++j ) {
00281       for( Thyra::Index i = 0; i < ev_V.subDim(); ++i ) {
00282         ev_V(i,j) = ev_V_complex(i,j).real();
00283       }
00284     }
00285   }
00286   else if( beta == ST::one() ) {
00287     for( Thyra::Index j = 0; j < ev_V.numSubCols(); ++j ) {
00288       for( Thyra::Index i = 0; i < ev_V.subDim(); ++i ) {
00289         ev_V(i,j) += ev_V_complex(i,j).real();
00290       }
00291     }
00292   }
00293   else {
00294     for( Thyra::Index j = 0; j < ev_V.numSubCols(); ++j ) {
00295       for( Thyra::Index i = 0; i < ev_V.subDim(); ++i ) {
00296         ev_V(i,j) = beta*ev_V(i,j) + ev_V_complex(i,j).real();
00297       }
00298     }
00299   }
00300 }
00301 
00302 #endif  // THYRA_REAL_COMPLEX_FFT_LINEAR_OP_HPP

Generated on Thu Sep 18 12:33:01 2008 for Thyra Package Browser (Single Doxygen Collection) by doxygen 1.3.9.1