serial_1D_FFT.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_SERIAL_1D_FFT_HPP
00030 #define THYRA_SERIAL_1D_FFT_HPP
00031 
00032 #include "Thyra_LinearOpWithSolveBaseDecl.hpp"
00033 #include "Thyra_SingleRhsLinearOpWithSolveBaseDecl.hpp"
00034 
00044 template<class RealScalar>
00045 void serial_1D_FFT( RealScalar data[], unsigned long nn, int isign )
00046 {
00047   unsigned long n,mmax,m,j,istep,i;
00048   RealScalar wtemp,wr,wpr,wpi,wi,theta;
00049   RealScalar tempr,tempi;
00050   n=nn << 1;
00051   j=1;
00052   for (i=1;i<n;i+=2) { // This is the bit-reversal section of the routine.
00053     if (j > i) {
00054       std::swap(data[j],data[i]);
00055       std::swap(data[j+1],data[i+1]);
00056     }
00057     m=nn;
00058     while (m >= 2 && j > m) {
00059       j -= m;
00060       m >>= 1;
00061     }
00062     j += m;
00063   }
00064   // Here begins the Danielson-Lanczos section of the routine.
00065   mmax=2;
00066   while (n > mmax) { // Outer loop executed log2(nn) times.
00067     istep=mmax << 1;
00068     theta=isign*(6.28318530717959/mmax); // Initialize the trigonometric recurrence.
00069     wtemp=sin(0.5*theta);
00070     wpr = -2.0*wtemp*wtemp;
00071     wpi=sin(theta);
00072     wr=1.0;
00073     wi=0.0;
00074     for (m=1;m<mmax;m+=2) {
00075       for (i=m;i<=n;i+=istep) {
00076         j=i+mmax;
00077         tempr=wr*data[j]-wi*data[j+1];
00078         tempi=wr*data[j+1]+wi*data[j];
00079         data[j]=data[i]-tempr;
00080         data[j+1]=data[i+1]-tempi;
00081         data[i] += tempr;
00082         data[i+1] += tempi;
00083       }
00084       wr=(wtemp=wr)*wpr-wi*wpi+wr; // Trigonometric recurrence.
00085       wi=wi*wpr+wtemp*wpi+wi;
00086     }
00087     mmax=istep;
00088   }
00089 }
00090 
00091 #endif  // THYRA_SERIAL_1D_FFT_HPP

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