Принятый выше ответ является разумным решением, если вам нужно только что-то близкое к двойной точности (Real (8)), поскольку вы явно определили Pi для числа цифр, близких к Real (8). Фактическое значение числа Пи до полного реального (8) равно
3.1415926535897931
ср. ваш параметр
3.14159265358979
Если вы хотите иметь более общую полную точность, соответствующую результату _cp, вы можете использовать что-то вроде
Real(cP), Parameter :: Pi = 4_cp*ATan(1_cp)
***Редактировать: спасибо, francescalus, за обнаружение опечаток, это должен быть ATan(1._cp), хотя на самом деле, как объяснено ниже, следует использовать параметризованные числа и избегать использования _cp в аргументах и т. д.
Теперь, если _cp равно, скажем, 16, ваш Pi будет автоматически:
3.14159265358979323846264338327950280
Теперь всякий раз, когда вы изменяете _cp, ваш код будет автоматически иметь полную точность, соответствующую _cp.
Кстати, вы также можете параметризовать некоторые базовые числа, например:
Real(cp), Parameter :: Zero = 0_cp
Real(cp), Parameter :: One = 1_cp
:
... так далее
*** Редактировать: спасибо, francescalus, за обнаружение опечаток, но в этих конкретных выражениях, хотя 0.0_cp и 1.0_cp могут быть лучше, это не имеет значения, так как декларация Params позаботится об этом ›
Затем в своем коде вы можете написать, например:
...
Real(cp), Parameter :: Pi = Four*ATan(One)
...
S = Cmplx(Zero, Zero)*Exp(-Two) ! etc etc
и не нужно беспокоиться о добавлении _cp и т. д. повсюду, и его намного легче читать/отлаживать и т. д.
У этой стратегии есть и другие преимущества, особенно в унаследованном коде, такие как
If( a == 1.0 ) Then
vs.
If( a == One ) Then
... Но это уже другая история.
Между прочим, и в некоторой степени вопрос стиля, в Фортране арифметика по умолчанию использует наибольшую точность/тип в выражении, поэтому
S = S + x(k)*exp(real(-2.0,cp)*PI*j*real(k-1,cp)*real(i-1,cp)/real(N,cp))
должно быть эквивалентно
S = S + x(k)*exp(-2*PI*j*(k-1)*(i-1)/N)
или еще лучше
S = S + x(k)*exp(-Two*PI*j*(k-1)*(i-1)/N)
Лично я считаю, что чем легче читать, тем легче понять правильно.
... хотя, как указано выше, если ваша арифметика содержит только целые числа, могут применяться дополнительные соображения. Однако, если вы Parameter'ize только Real, то у вас никогда не должно быть риска (например, div/0) использования 2 там, где действительно требуется 2.0, если вы используете Two.
Наконец, поскольку у вас есть Parameter'ized Pi, почему бы не использовать также 2*Pi, например:
...
Real(cp), Parameter :: TwoPi = Two*Pi
теперь выражение
S = S + x(k)*exp(-Two*PI*j*(k-1)*(i-1)/N)
возможно
S = S + x(k)*exp(-TwoPI*j*(k-1)*(i-1)/N)
что немного экономит процессор. Фактически, расширение этой логики может еще больше повысить производительность вашего кода.
... в совершенно другом месте, знаете ли вы, что ваш var Temp() может быть автоматическим, скажем,
complex(cp) :: temp(Size(x, Dim = 1))
тогда нет необходимости в Allocatable и т.д. Хотя с того места, где я сижу, кажется, что можно написать свой код полностью без S и Temp(:), обеспечивая гораздо более простой и быстрый код.
***** Приложение: после некоторых комментариев к моему ответу я подумал, что это может помочь показать возможные изменения в коде OP для повышения производительности и, в некоторой степени, дальнейшего повышения точности.
Однако до этого в ОП не указывалось, почему требуется определенная степень точности. Это может быть важным вопросом во многих промышленных/реальных условиях. Например, в некоторых случаях скорость гораздо важнее точности, если точность достаточно хороша, чтобы обеспечить контекстно-зависимую надежность/принятие решений. Например, возможно, что в некоторых случаях все вычисления должны быть в Real(4). Этот вопрос требует отдельного и нетривиального обсуждения.
Как бы то ни было, в то время как ответ Росса исправил проблему точности с точки зрения представления в памяти, а мой первоначальный ответ улучшил точность с точки зрения того, что вам действительно нужно предоставить достаточное количество sig figs для начала (даже если объявления верны), уменьшая количество FLOPS не только повышает скорость, но также может помочь с точностью, поскольку каждый FLOP вводит возможность усечения/округления, особенно для длинной последовательности FLOP в одном выражении.
Кроме того, эту проблему OP можно было бы хорошо решить с помощью формулировки идентичности (4 * ATan (1)). Во многих случаях это либо невозможно, либо практично, и тогда требуются менее элегантные методы, но я оставляю их на другой день.
Ниже приведены только одна/две возможные альтернативы.
Я не проверял это.
поэтому может потребоваться настройка; пожалуйста, не стесняйтесь предлагать улучшения.
Если повезет, это может уменьшить FLOP примерно на 80%.
... Я был бы признателен ОП за сравнение различных реализаций и отчет о результатах, если это удобно. Вы также можете провести тест для _cp = 4, 8, 16, просто чтобы продемонстрировать компромисс между скоростью и точностью.
Эта альтернативная версия потребует очевидного обновления вызывающего s/r.
module fft_mod
!
public :: fft1D
!
Integer,parameter :: cp = selected_real_kind(14) ! just out of curiosity, why "14", see my comments on "reality vs. precision"
!>>real(cp),parameter :: PI = real(3.14159265358979,cp)
!
Real(cp), Parameter :: Zero = 0_cp ! these types of declarations can be created in a separate Module,
! say "MyDeclarationModule", and then just access with "Use"
!
!XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX
! NOTE: With respect to francescalus comments on issues with "1_cp" syntax, the usage here works fine as shown in the previous result,
! though francescalus' comments are fair in suggesting that other approaches make for better coding.
! However, and to be clear, I don't actually use "_xx" style myslef. I have used it here ONLY to be consistent with the
! the OP's and earlier answers. The usage of it here is NOT a recommendation. A discussion of the alternatives
! is much too long and strays from the immediate issues too far.
! ... Perhaps francescalus will take up the mantle and re-write the OP's code for alternatives
!XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX
!
Real(cp), Parameter :: One = 1_cp
Real(cp), Parameter :: Two = 2_cp
Real(cp), Parameter :: Four = 4_cp
!
! ... after this point, there is no need to see "_cp" again, in this example
!
Real(cp), Parameter :: Pi = Four*ATan(One) ! this guarrantees maximal precision for Pi, up to "_cp"
!
Real(cp), Parameter :: TwoPi = Two*Pi ! Vladimir: now, and only now (that I've talem most everything out of the loop),
! this declaration has less value ... unlike previously, when it was
! still in the inner NxN, and when it saved approx 15 - 20% of FLOPs
! Crucially, if you do a lot of computational mathematics, TwoPi, PiBy2, Root2Pi, etc etc
! arise with considerable frequency, and creating these as Parameters makes for much
! improvement, and I would think it a practice to be encouraged.
!
Complex(cp), Parameter :: SqrtM1 = Cmplx(Zero, One) ! sqrt(-1), this is your "j",
! sorry but "j" sounds just too much like an "integer" to me
! (for us "old timers"), so I changed the name to something more meaningful to me
!
!
Contains
!
Pure Subroutine fft1D(x, n, Temp) ! OPTIONAL, return the results explicitly to save vector assignment, and preserve original data
! OPTIONAL, send the size of x() explicitly as n, much cleaner
! you could leave it as s/r fft1D(x), and most of the savings below would still apply
! ... also, my practice is to make everything Pure/Elemental where reasonable
!
Integer, Intent(In) :: n ! Optional, but cleaner
! complex(cp), Intent(InOut) :: x(n) ! consider as just Intent(In), and return Temp() as Intent(Out)
complex(cp), Intent(In) :: x(n) ! consider as just Intent(In), and return Temp() as Intent(Out)
!
! Locals
! ------
!
! Complex(cp), Allocatable :: temp(:) ! no need for Allocatable
! Complex(cp) :: temp(Size(x), Dim = 1) ! I prefer to pass n explicitly
Complex(cp), Intent(Out) :: temp(n) ! make Intent(Out) to preserve data and save a vector assignment
!
! complex(cp) :: S
! integer :: N
! complex(cp) :: j ! sqrt(-1)
!
Integer :: i, k
!
Complex(cp) :: TPSdivN ! new reduce 4 nxn mults to 1
!
Complex(cp) :: TPSdivNiM1 ! new, to reduce 5 nxn mults to 1 nx1
!
Real(cp) :: rKM1(n) ! new, to obviate nxn k-1's in inner loop
!
! N=size(x) ! my preference is to pass N explicitly,
! but can be done this way if prefered (need to make appropriate changes above)
!
! allocate(temp(N)) ! Temp() can be either automatic or an Arg, no need for Allocate
! j = cmplx(0.0_cp,1.0_cp,cp) ! make Parameter, and rename to SqrtM1 (personal pref)
! S = cmplx(0.0_cp,0.0_cp,cp) ! Not required with "improved" inner loop
!
!
temp(1:n) = Cmplx(Zero, Zero) ! S is not needed, just start with Temp() directly
!
TPSdivN = -TwoPi*SqrtM1/n ! new, move everything out all loops that can be
!
ForAll( k=1:n) rKM1(k) = k - 1 ! new, this allows the elimination of NxN "k-1's" in the inner
! my preference is for ForAll's, but you could use Do,
!
do i = 1,N
!
TPSdivNiM1 = TPSdivN*(i-1) ! new. everything out of inner loop that can be
!
! improved, inner do, but could use "no-Do" alternative as illustrated below
!
! do k = 1,N
!>> S = S + x(k)*exp(real(-2.0,cp)*PI*j*real(k-1,cp)*real(i-1,cp)/real(N,cp)) ! original, many unneccessary nxn FLOPs, type conversions, etc
! temp(i) = temp(i) + x(k)*Exp(TPSdivNiM1*rKM1(k))) ! use array of k-1's, then the inner k-loop is no longer required,
! and can use Product()/Sum() or Dot_Product() intrinsics, see example below
! End Do
!
! alternative "direct array" approach to inner loop, or "no-DO" version ... there are other possibilities.
!
Temp(i) = Dot_Product( x(1:n), Exp(TPSdivNiM1*rKM1(1:n)) ) ! this approach can have certain advantages depending on circumstances
!
! temp(i) = S ! not needed
! S = cmplx(0.0_cp,0.0_cp,cp) ! not needed
!
End Do
!
! x(1:n) = temp(1:n) ! not need if Temp() is Intent(Out) that save this vector assignment, and the original data
! deallocate(temp) ! not needed
!
End Subroutine fft1D
!
! end module
End Module fft_mod
person
DrOli
schedule
03.09.2015
