
forums.silverfrost.com Welcome to the Silverfrost forums

View previous topic :: View next topic 
Author 
Message 
Kenneth_Smith
Joined: 18 May 2012 Posts: 682 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: 7770 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: 682 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
