
forums.silverfrost.com Welcome to the Silverfrost forums

View previous topic :: View next topic 
Author 
Message 
Kenneth_Smith
Joined: 18 May 2012 Posts: 709 Location: Hamilton, Lanarkshire, Scotland.

Posted: Wed Mar 08, 2023 9:50 am Post subject: 


Mecej4,
Great work you are doing on this. I will certainly give the DLL a try when it becomes available. The functions for modified Bessel functions with complex arguments are of particular interest to myself (they are used extensively in Schelkunoff�s model of coaxial cables with cylindrical shields and/or armours). 

Back to top 


mecej4
Joined: 31 Oct 2006 Posts: 1896

Posted: Wed Mar 08, 2023 12:10 pm Post subject: 


Kenneth_Smith,
Thanks for your kind words.
Yesterday, I posted a note about my progress:
http://forums.silverfrost.com/viewtopic.php?t=4745
I am very pleased with the quality of the SLATEC software, as well as the ability of my software tools to help me detect and fix the small number of bugs in the package ( a score or two, but that is minute in a package with 1900 source files).
The FTN95compiled 64 bit DLL works very well, both in regard to accuracy as well as speed. There is one compiler bug (relating to complex variables, noted in this thread above) for Silverfrost to fix before the DLL can pass all the 54 SLATEC tests.
Do you need to evaluate Modified Bessel functions with double precision complex arguments? Or will single precision do? Note that there was no general support for double precision complex ("DCOMPLEX") in Fortran 77, so SLATEC may not provide the support that you may want. However, we can discuss the issues and consider whether it is worthwhile to supplement SLATEC by adding double precision complex functions using Fortran 200X intrinsics to the extent possible.
If you are willing to test the DLL and the few usage examples that I have prepared, notwithstanding the bug in FTN95, I can send you download links for the DLL and the example sources.
Yesterday I found this documentation page for SLATEC, and you may wish to browse through it.
http://www.cs.yorku.ca/~roumani/fortran/slatecAPI/keylist.index.html 

Back to top 


Kenneth_Smith
Joined: 18 May 2012 Posts: 709 Location: Hamilton, Lanarkshire, Scotland.

Posted: Wed Mar 08, 2023 8:17 pm Post subject: 


mecej4,
This is interesting. The first Fortran 77 I used in the mid 1980s supported COMPLEX*16, and I was under the mistaken impression that this was part of the Fortran 77 standard for nearly 40 years!
Almost certainly the application I have in mind for the modified Bessel functions would run into floating point issues if the calculations were performed in single precision.
I have found the Fortran 77 routines from �FORTRAN routines for computation of Special Functions� by Shanjie Zhang, and Jianming Jin, to be useful (these do use COMPLEX*16 � perhaps reinforcing my erroneous understanding of what the 77 standard specified) . The original link for the source code files where I downloaded the routines does not work anymore. These routines are not in the public domain.
They do appear to have been collated into a single file by John Burkardt (I have not checked so see if they reflect the information I have in the errata document in the download I obtained about 8 years ago).
https://people.sc.fsu.edu/~jburkardt/f77_src/special_functions/special_functions.f
Thanks for the link to the SLATEC documentation page, that is very useful � and easy to search.
I have some time at present to experiment with the DLL and I am more than happy to help with this where I can. 

Back to top 


mecej4
Joined: 31 Oct 2006 Posts: 1896

Posted: Thu Mar 09, 2023 12:29 am Post subject: SLATEC DLLs and demo problem drivers 


Here is a link to a Zip file containing 32 and 64bit SLATEC DLLs for use with FTN95. The Zip also contains a small number of test programs that exercise the DLLs.
https://drive.google.com/file/d/1EGPTavu0RMBuu8ZHY2aTzeGn7vcBJj0/view?usp=sharing
If interested, you may download the various testNN.f files from https://netlib.org/slatec/chk/ and run them with the DLLs. You DO NOT need any of the other source files in that directory, since those routines are already included in the DLL.
To build and run, for example, Test25, the commands are
Code:  ftn95 test25.f
slink test25.obj slatec32.dll
echo 3  test25 
The "echo 3" is needed since each of the SLATEC test programs reads one integer before running. That integer specifies the verbosity of the output, and allowable values are 1, 2 or 3.
Feedback is welcome. 

Back to top 


Kenneth_Smith
Joined: 18 May 2012 Posts: 709 Location: Hamilton, Lanarkshire, Scotland.

Posted: Thu Mar 09, 2023 2:12 pm Post subject: 


mecej4,
I have tested each of your examples and only have one comment:
xsdriv1.f90 Array subscripts out of bounds at line 22.
All the other examples run as expected.
I have also compared the SLATEC routines for J0, J1, Y0, and Y1 (real single and double precision arguments) with the corresponding intrinsic BESSEL_J0 etc and don�t see any unexpected results, with the 32 or 64 bit dlls.
It appears to be relatively straightforward to access the library functions.
I shall do some more testing of other functions later. I will also look at what happens when running from within a Plato project environment. I don�t use FTN95 from the command line that often. 

Back to top 


mecej4
Joined: 31 Oct 2006 Posts: 1896

Posted: Thu Mar 09, 2023 3:02 pm Post subject: 


Thanks for the feedback.
Regarding the subscript overrun in XSDRIV1: it is my mistake. Please change the allocation statement to:
Code:  allocate (work(lenw), y(n+1)) 


Back to top 


Kenneth_Smith
Joined: 18 May 2012 Posts: 709 Location: Hamilton, Lanarkshire, Scotland.

Posted: Thu Mar 09, 2023 5:29 pm Post subject: 


mecej4,
The integration subroutine DGAUS8 is also fine for both 32 and 64 bit dlls.
The following program (which you may wish to adapt as a further example) demonstrates DGAUS8 and was developed and run from within Plato.
Code:  ! EXAMPLE (Integrate a real function of one variable over a finite interval)
!
! The problem is to integrate the function
!
! f(x) = (x**2)*(x**22.d0)*sin(x)
!
! over the interval a to b
!
! Output from the DGAUS8 subroutine is compared with the exact solution.
!
! ##########
!
program test_dgaus8
implicit none
integer i, ierr
double precision a, b, err, ans, exact, delta
double precision, external :: fun
double precision funt, x
funt(x) = 4.d0*x*((x**2)7.d0)*sin(x)((x**4)14.d0*(x**2)+28.d0)*cos(x)
a = 0.d0
b = 0.d0
do i = 1, 12
err = 1.0d06
call DGAUS8 (FUN, A, B, ERR, ANS, IERR)
exact = funt(b)  funt(a)
delta = ans  exact
write(6,*)
write(6,'(A20,F16.12)') 'Lower limit ', A
write(6,'(A20,F16.12)') 'Upper limit ', B
write(6,'(A20,F16.12)') 'Solution from DGUAS8', ANS
write(6,'(A20,F16.12)') 'Exact solution ', exact
write(6,'(A20,F16.12)') 'Delta ', delta
b = b + 4.d0*atan(1.d0)/12.d0
end do
end program test_dgaus8
double precision function fun(x)
double precision x
fun = (x**2)*(x**22.d0)*sin(x)
end function fun

Last edited by Kenneth_Smith on Fri Mar 10, 2023 12:00 pm; edited 1 time in total 

Back to top 


mecej4
Joined: 31 Oct 2006 Posts: 1896

Posted: Fri Mar 10, 2023 12:18 am Post subject: Example of Gaussian Integration 


Excellent example of the power of integration by parts. I have added it to the examples collection.
Change the last line to "end" or "end function" or "end function fun". 

Back to top 


Kenneth_Smith
Joined: 18 May 2012 Posts: 709 Location: Hamilton, Lanarkshire, Scotland.

Posted: Sat Mar 11, 2023 12:32 pm Post subject: 


mecej4,
I had some fun yesterday evening examining IFLAG from FZERO and came up with another example for you to consider.
Code:  ! EXAMPLE (Use of Slatec routines BESJ0 and FZERO)
!
! The problem is to find the first ten zeros of:
!
! f(x) = BESJ0(x) for x > 0.0
!
! ##########
program xbesj0fzero
implicit none
integer iflag,nroots
real b,c,r,abserr,relerr,bstart,cstart,dx
real, external :: besj0
external :: fzero
relerr = 0.0001
abserr = 0.00005
dx = 0.1
bstart = 0.0
cstart = bstart + dx
nroots = 0
write(6,'(1x,a)') 'Roots of besj0(x):'
write(6,'(1x,t17,a,t27,a)') 'x','f(x)'
! Begin search
do while (nroots .lt. 10)
b = bstart
c = cstart
call FZERO(BESJ0,b,c,r,relerr,abserr,iflag)
! Exammine iflag to determine next step
if (iflag .eq. 4) then
! No change in sign of F(X) was found  expand search range
cstart = cstart + dx
cycle
else if (iflag .eq. 1) then
! B is within the requested tolerance of a zero.
nroots = nroots + 1
write(6,'(1x,a,i3,2(4x,2f10.4))') 'Root', nroots, b, BESJ0(b)
! Now set search range for next root beyond b
bstart = b + dx
cstart = bstart + dx
cycle
else if (iflag .eq. 5) then
STOP 'Too many iterations in FZERO'
end if
end do
end program xbesj0fzero



Back to top 


mecej4
Joined: 31 Oct 2006 Posts: 1896


Back to top 


Kenneth_Smith
Joined: 18 May 2012 Posts: 709 Location: Hamilton, Lanarkshire, Scotland.

Posted: Sun Mar 12, 2023 8:08 pm Post subject: 


mecej4,
My initial version was much longer and included a call to an initial routine to bracket the roots, rather than expanding the search space based on iflag returned from FZERO.
I was not ware of these other routines to do the same calculation. I chose BESJ0 since it was simply another function in the library which was well behaved for this purpose.
Here is another example which uses the routines CGECO and CGEDI. I did not expect this to run correctly due to the issues with complex arithmetic. What is interesting is that in this particular case, the results returned using the 64 bit compiler are correct while those with the 32 bit compiler are not correct. This is not what I was expecting (I thought it would be the other way round). Perhaps you can try this example?
Code corrected, as per next post. Now runs correctly.
Code:  ! EXAMPLE (Use of Slatec routines CGECO and CGEDI)
!
! The problem is to find the inverse of the complex 3x3 matrix
!
!  1.0 1.0 1.0 
! h =  1.0 a**2 a  where a = EXP(j*2*pi/3)
!  1.0 a a**2 
!
! (a is approximately equal to 0.500 + j*0.866).
!
! The solution to the problem is
!
!  1.0 1.0 1.0 
! hinv = (1/3) *  1.0 a a**2 
!  1.0 a**2 a 
!
! ##########
program xcgecocgedi
implicit none
integer ipvt(3),job
real rcond, pi
complex a,a2,h(3,3),c1,w(3,3),det(3,3)
external :: CGECO, CGEDI
! Set up the h matrix
pi = 4.0*atan(1.0)
a = EXP(cmplx(0.0,1.0)*2.0*pi/3.0)
a2 = a*a
c1 = cmplx(1.0,0.0)
write(6,'(a,1x,"(",f10.6,",",f10.6,")")') 'a = ', a
write(6,'(a,1x,"(",f10.6,",",f10.6,")")') 'a**2 = ', a2
write(6,'(a,1x,"(",f10.6,",",f10.6,")")') 'c1 = ', c1
write(6,*)
h = transpose(reshape([c1,c1,c1,c1,a2,a,c1,a,a2],shape(h)))
write(6,'(a)') 'H matrix'
call print_cmat_s(h,3)
! Now call CGECO and CGEDI to obtain inverse of h
call CGECO(h, 3, 3, IPVT, RCOND, w)
job = 01
call CGEDI(h, 3, 3, IPVT, DET, w, JOB)
! h now contains the inverse of the original h
write(6,*)
write(6,'(a)') 'Inverse of H matrix calculated using CGECO and CGEDI'
call print_cmat_s(h,3)
! Here we have the correct solution
h = (1.0/3.0)*transpose(reshape([c1,c1,c1,c1,a,a2,c1,a2,a],shape(h)))
write(6,*)
write(6,'(a)') 'Correct inverse of H matrix'
call print_cmat_s(h,3)
end program xcgecocgedi
subroutine print_cmat_s(a,n)
implicit none
integer, intent(in) :: n
complex, intent(in) :: a(n,n)
integer i, j
do i = 1, n
write(6,'(3(1x,"(",f10.6,",",f10.6,")"))') (a(i,j),j=1,n)
end do
end subroutine print_cmat_s 
Error in print_cmat_s also corrected J=1,n was previously j=1,3
Last edited by Kenneth_Smith on Mon Mar 13, 2023 11:58 am; edited 4 times in total 

Back to top 


mecej4
Joined: 31 Oct 2006 Posts: 1896

Posted: Sun Mar 12, 2023 9:00 pm Post subject: 


Kenneth,
The argument IPVT needs to be an array of dimension (3). 

Back to top 


Kenneth_Smith
Joined: 18 May 2012 Posts: 709 Location: Hamilton, Lanarkshire, Scotland.

Posted: Sun Mar 12, 2023 9:25 pm Post subject: 


mecej4,
Thanks I should have spotted that. So no problems with this example.
The Scotland vs. Ireland rugby match as on the television when I was looking at this, which may have distracted me. The final score was very upsetting to a Scotsman! 

Back to top 


mecej4
Joined: 31 Oct 2006 Posts: 1896

Posted: Mon Mar 13, 2023 1:13 am Post subject: 


Kenneth,
If it is any consolation, 7/22 is a good approximation to 1 / Pi. 

Back to top 


PaulLaidler Site Admin
Joined: 21 Feb 2005 Posts: 8011 Location: Salford, UK

Posted: Mon Mar 13, 2023 10:10 am Post subject: 


The bug relating to test04 above has now been fixed for the next release of FTN95. 

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
