ComplexFFTLinearOp.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_COMPLEX_FFT_LINEAR_OP_HPP
00030 #define THYRA_COMPLEX_FFT_LINEAR_OP_HPP
00031 
00032 #include "Thyra_LinearOpWithSolveBase.hpp"
00033 #include "Thyra_SingleRhsLinearOpWithSolveBase.hpp"
00034 #include "Thyra_DefaultSpmdVectorSpace.hpp"
00035 #include "Thyra_DetachedVectorView.hpp"
00036 #include "serial_1D_FFT.hpp"
00037 
00045 template<class RealScalar>
00046 class ComplexFFTLinearOp
00047   : virtual public Thyra::LinearOpWithSolveBase< std::complex<RealScalar> >
00048   , virtual protected Thyra::SingleRhsLinearOpWithSolveBase< std::complex<RealScalar> >
00049 {
00050 public:
00051 
00053   typedef std::complex<RealScalar> Scalar;
00055   ComplexFFTLinearOp( const int N );
00056 
00060   Teuchos::RefCountPtr< const Thyra::VectorSpaceBase< std::complex<RealScalar> > > range() const;
00062   Teuchos::RefCountPtr< const Thyra::VectorSpaceBase< std::complex<RealScalar> > > domain() const;
00064 
00065 protected:
00066 
00070   bool opSupported(Thyra::ETransp M_trans) const;
00072 
00076   void apply(
00077     const Thyra::ETransp                                    M_trans
00078     ,const Thyra::VectorBase< std::complex<RealScalar> >    &x
00079     ,Thyra::VectorBase< std::complex<RealScalar> >          *y
00080     ,const std::complex<RealScalar>                         alpha
00081     ,const std::complex<RealScalar>                         beta
00082     ) const;
00084 
00088   bool solveSupportsTrans(Thyra::ETransp M_trans) const;
00090   bool solveSupportsSolveMeasureType(
00091     Thyra::ETransp M_trans, const Thyra::SolveMeasureType& solveMeasureType
00092     ) const;
00094 
00098   Thyra::SolveStatus< std::complex<RealScalar> > solve(
00099     const Thyra::ETransp                                             M_trans
00100     ,const Thyra::VectorBase< std::complex<RealScalar> >             &b
00101     ,Thyra::VectorBase< std::complex<RealScalar> >                   *x
00102     ,const Thyra::SolveCriteria< std::complex<RealScalar> >          *solveCriteria
00103     ) const;
00105 
00106 private:
00107 
00108   Teuchos::RefCountPtr< const Thyra::VectorSpaceBase< std::complex<RealScalar> > > space_;
00109 
00110 };
00111 
00112 // /////////////////////////
00113 // Definitions
00114 
00115 template<class RealScalar>
00116 ComplexFFTLinearOp<RealScalar>::ComplexFFTLinearOp( const int N )
00117 {
00118   space_ = Teuchos::rcp(
00119     new Thyra::DefaultSpmdVectorSpace< std::complex<RealScalar> >(
00120       int(std::pow(double(2),N))
00121       )
00122     );
00123 }
00124 
00125 // Overridden from LinearOpBase
00126 
00127 template<class RealScalar>
00128 Teuchos::RefCountPtr< const Thyra::VectorSpaceBase< std::complex<RealScalar> > >
00129 ComplexFFTLinearOp<RealScalar>::range() const
00130 {
00131   return space_;
00132 }
00133 
00134 template<class RealScalar>
00135 Teuchos::RefCountPtr< const Thyra::VectorSpaceBase< std::complex<RealScalar> > >
00136 ComplexFFTLinearOp<RealScalar>::domain() const
00137 {
00138   return space_;
00139 }
00140 
00141 // Overridden from SingleScalarLinearOpBase
00142 
00143 template<class RealScalar>
00144 bool ComplexFFTLinearOp<RealScalar>::opSupported(Thyra::ETransp M_trans) const
00145 {
00146   return ( M_trans == Thyra::NOTRANS || M_trans == Thyra::CONJTRANS );
00147 }
00148 
00149 // Overridden from SingleRhsLinearOpBase
00150 
00151 template<class RealScalar>
00152 void ComplexFFTLinearOp<RealScalar>::apply(
00153   const Thyra::ETransp                                    M_trans
00154   ,const Thyra::VectorBase< std::complex<RealScalar> >    &x
00155   ,Thyra::VectorBase< std::complex<RealScalar> >          *y
00156   ,const std::complex<RealScalar>                         alpha
00157   ,const std::complex<RealScalar>                         beta
00158   ) const
00159 {
00160   typedef Teuchos::ScalarTraits< std::complex<RealScalar> > ST;
00161 #ifdef TEUCHOS_DEBUG
00162   TEST_FOR_EXCEPT( !( M_trans == Thyra::NOTRANS || M_trans == Thyra::CONJTRANS ) );
00163 #endif
00164   // Update y first
00165   Thyra::Vt_S( y, beta );
00166   // Translate from input x into one long array with data[] that will be
00167   // passed to and from the FFT routine
00168   const Thyra::ConstDetachedVectorView<Scalar>  x_ev(x);
00169   std::vector<RealScalar> data(2*x_ev.subDim());
00170   for( int k = 0; k < x_ev.subDim(); ++k ) {
00171     data[2*k]   = x_ev[k].real();
00172     data[2*k+1] = x_ev[k].imag();
00173   }
00174   // Call the FFT rountine
00175   serial_1D_FFT(
00176     &data[0]-1                                            // This function is 1-based of all things!
00177     ,x_ev.subDim()                                        // 1/2 length of data[]
00178     ,Thyra::real_trans(M_trans)==Thyra::NOTRANS ? +1 : -1 // +1 = fwd, -1 = adjoint
00179     );
00180   // Add the scaled result into y
00181   const Thyra::DetachedVectorView<Scalar>  y_ev(*y);
00182   const Scalar scalar = alpha * Scalar(1)/ST::squareroot(x_ev.subDim()); // needed to make adjoint == inverse!
00183   for( int k = 0; k < y_ev.subDim(); ++k ) {
00184     y_ev[k] += ( scalar * Scalar(data[2*k],data[2*k+1]) );
00185   }
00186 }
00187 
00188 // Overridden from SingleScalarLinearOpWithSolveBase
00189 
00190 template<class RealScalar>
00191 bool ComplexFFTLinearOp<RealScalar>::solveSupportsTrans(Thyra::ETransp M_trans) const
00192 {
00193   return ( M_trans == Thyra::NOTRANS || M_trans == Thyra::CONJTRANS );
00194 }
00195 
00196 template<class RealScalar>
00197 bool ComplexFFTLinearOp<RealScalar>::solveSupportsSolveMeasureType(
00198   Thyra::ETransp M_trans, const Thyra::SolveMeasureType& solveMeasureType
00199   ) const
00200 {
00201   return ( M_trans == Thyra::NOTRANS || M_trans == Thyra::CONJTRANS );
00202 }
00203 
00204 // Overridden from SingleRhsLinearOpWithSolveBase
00205 
00206 template<class RealScalar>
00207 Thyra::SolveStatus< std::complex<RealScalar> >
00208 ComplexFFTLinearOp<RealScalar>::solve(
00209   const Thyra::ETransp                                             M_trans
00210   ,const Thyra::VectorBase< std::complex<RealScalar> >             &b
00211   ,Thyra::VectorBase< std::complex<RealScalar> >                   *x
00212   ,const Thyra::SolveCriteria< std::complex<RealScalar> >          *solveCriteria
00213   ) const
00214 {
00215   typedef Teuchos::ScalarTraits< std::complex<RealScalar> > ST;
00216 #ifdef TEUCHOS_DEBUG
00217   TEST_FOR_EXCEPT( !( M_trans == Thyra::NOTRANS || M_trans == Thyra::CONJTRANS ) );
00218 #endif
00219   Thyra::apply( *this, M_trans==Thyra::NOTRANS?Thyra::CONJTRANS:Thyra::NOTRANS, b, x );
00220   typedef Thyra::SolveCriteria< std::complex<RealScalar> >  SC;
00221   typedef Thyra::SolveStatus< std::complex<RealScalar> >    SS;
00222   SS solveStatus;
00223   solveStatus.solveStatus = Thyra::SOLVE_STATUS_CONVERGED;
00224   solveStatus.achievedTol = SS::unknownTolerance();
00225   return solveStatus;
00226 }
00227 
00228 #endif  // THYRA_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