 |
forums.silverfrost.com
Welcome to the Silverfrost forums
|
| View previous topic :: View next topic |
| Author |
Message |
Kenneth_Smith
Joined: 18 May 2012
Posts: 726
Location: Hamilton, Lanarkshire, Scotland.
|
 Posted: Sun Apr 09, 2023 9:11 pm Post subject: Complex issue with kind=1 |
 |
|
The following program returns the correct results when compiled with Win32.
When compiled with X64 the results are incorrect. The error occurs at line 8, in which there is the product of two complex numbers, the latter defined by two nested function calls.
| Code: |
module problem_mod
integer, parameter :: wp = kind(1.0)
contains
! Rotate phasor by specified angle (deg)
complex(kind=wp) function crotphasordeg(phasor,angle)
complex(kind=wp), intent(in) :: phasor
real(kind=wp), intent(in) :: angle
crotphasordeg = phasor*cang2phasor(sdeg2rad(angle)) ! ERROR OCCURS IN THIS EXPRESSION
! crotphasordeg = cang2phasor(sdeg2rad(angle)) ! This alternative form gives the
! crotphasordeg = phasor*crotphasordeg ! correct result.
end function crotphasordeg
! Define the unit phasor at the specified angle
complex(kind=wp) function cang2phasor(angle)
real(kind=wp), intent(in) :: angle
cang2phasor = exp(cmplx(0.0,1.0,kind=wp)*angle) ! e^(j*angle_rad)
end function cang2phasor
! Convert supplied angle in deg to radians
real(kind=wp) function sdeg2rad(angle)
real(kind=wp), intent(in) :: angle
sdeg2rad = angle*(4.0*atan(1.0))/180.0
end function sdeg2rad
end module problem_mod
program p
use problem_mod
complex(kind=wp) z, z1
real(kind=wp) angle
z = cmplx(1.0,0.0,kind=wp)
angle = 30.0
z1 = crotphasordeg(z,angle)
print*, 'The phasor ', z, 'rotated by ', angle, 'degrees is', z1
angle = -30.0
z1 = crotphasordeg(z,angle)
print*, 'The phasor ', z, 'rotated by ', angle, 'degrees is', z1
print*
end program p
|
If the module parameter wp is changed to kind(1.d0) everything is as expected.
Evaluating the required values in a single line of code (kind=1), i.e. no function calls returns the correct results for both win32 and x64.
| Code: | print*, '+30', cmplx(1.0,0.0)*exp(cmplx(0.0,1.0)*(+30.0*(4.0*atan(1.0))/180.0))
print*, '-30', cmplx(1.0,0.0)*exp(cmplx(0.0,1.0)*(-30.0*(4.0*atan(1.0))/180.0))
end |
|
|
| Back to top |
|
 |
PaulLaidler
Site Admin
Joined: 21 Feb 2005
Posts: 8037
Location: Salford, UK
|
| Posted: Mon Apr 10, 2023 7:13 am Post subject: |
|
|
Ken
Thank you for the bug report. It appears to have been already fixed and I am assuming that it is the same bug as that recently reported by mecej4 elsewhere on the Forum.
This fix will be in the next release of FTN95. |
|
| Back to top |
|
|
Kenneth_Smith
Joined: 18 May 2012
Posts: 726
Location: Hamilton, Lanarkshire, Scotland.
|
| Posted: Mon Apr 10, 2023 10:38 am Post subject: |
|
|
Paul, that�s good news, i.e. you don�t see the real part of the result going off to an very large number.
I didn�t make the connection with mecej4�s post in the SLATEC thread, but looking at it this morning I can see some similarities. In particular, the error I was wrestling with yesterday is sensitive to the requested rotation angle.
I�ll park this one for now, pending the next FTN95 release.
Ken |
|
| Back to top |
|
|
|
|
You cannot post new topics in this forum
You cannot reply to topics in this forum
You cannot edit your posts in this forum
You cannot delete your posts in this forum
You cannot vote in polls in this forum
|
Powered by phpBB © 2001, 2005 phpBB Group
|
FAQ
Search
Memberlist
Usergroups
Register
Profile
Log in to check your private messages
Log in
