| // add term to sum
y_re[k]+=term_re;
y_im[k]+=term_im;
}
}
// transform y to polar coordinates
for(k=0;k
{
// compute amplitude of y(k)
a[k]=sqrt(y_re[k]*y_re[k]+y_im[k]*y_im[k]);
// compute phase of y(k)
phi[k]=atan2(y_im[k],y_re[k]);
}
}
// initialize input data
void init(double a[],int *m)
{
int j,n=*m;
for(j=0;j
{
a[j]=sin(1./sqrt((double)j+1.0));
}
}
// Dummy function to use result, preventing compiler from
// optimizing away the computation.
void consume(double a[],double b[],double c[])
{
}
附录 3 – 驱动器程序
下面是一个用于计时 DFT 代码的 Fortran 驱动器程序。
program main
interface
subroutine dft(x,a,phi,n)
real*8 x(n),a(n),phi(n)
integer n
end subroutine
subroutine init(a,n)
real*8 a(n)
integer n
end subroutine
subroutine consume(a,b,c)
real*8 a(*),b(*),c(*)
end subroutine
end interface
! Parameters:
! nmax is the problem size.
! nrep is the number of repetitions of the
! problem. This should be chosen so that
! the elapsed time is long enough to give
! sufficient timing resolution.
! cyc is the clock frequency in Hz for the
! processor that the program is to be run on.
! (Can be found from AIX command pmcycles.)
integer, parameter :: nmax=1000
integer, parameter :: nrep=100
real*8, parameter :: cyc=4704000000.d0
real*8 x(nmax), a(nmax), phi(nmax)
real*8 tx, ty, accum, del(4)
intrinsic sin, sqrt
real*8 rtc
del(4)=0.d0
acc = 0.d0
do k=1,nrep
tx=rtc()
call init(x,nmax)
call dft(x,a,phi,nmax)
ty=rtc()
call consume(x,a,phi)
do j=1,nmax
acc = acc + a(j) + phi(j)
end do
del(4) = del(4) + (ty-tx)
end do
del(1) = del(4)/real(nmax,8)
del(2) = del(1)/real(nrep,8)
del(3) = cyc*del(2)
print *,'acc=',acc/real(nrep,8),' n=',nmax,
& ' r=',nrep, ' a=',del(1),' b=',del(2),
& ' c=',del(3),' w=',del(4)
end program
声明
获得的精确性可能会随着使用处理器的模型及其配置而发生变化，同时依赖的因素还有所用汇编器的版本。因此，在不同的操作条件下您将获得不同的使用经验。 |
