forums.silverfrost.com Forum Index forums.silverfrost.com
Welcome to the Silverfrost forums
 
 FAQFAQ   SearchSearch   MemberlistMemberlist   UsergroupsUsergroups   RegisterRegister 
 ProfileProfile   Log in to check your private messagesLog in to check your private messages   Log inLog in 

Porting SLATEC Library -- Undocumented routines in SALFLIBC
Goto page Previous  1, 2, 3, 4  Next
 
Post new topic   Reply to topic    forums.silverfrost.com Forum Index -> Support
View previous topic :: View next topic  
Author Message
Kenneth_Smith



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

PostPosted: Wed Mar 08, 2023 9:50 am    Post subject: Reply with quote

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
View user's profile Send private message Visit poster's website
mecej4



Joined: 31 Oct 2006
Posts: 1885

PostPosted: Wed Mar 08, 2023 12:10 pm    Post subject: Reply with quote

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 FTN95-compiled 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
View user's profile Send private message
Kenneth_Smith



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

PostPosted: Wed Mar 08, 2023 8:17 pm    Post subject: Reply with quote

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
View user's profile Send private message Visit poster's website
mecej4



Joined: 31 Oct 2006
Posts: 1885

PostPosted: Thu Mar 09, 2023 12:29 am    Post subject: SLATEC DLLs and demo problem drivers Reply with quote

Here is a link to a Zip file containing 32- and 64-bit 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/1-EGPTavu0RMBuu8ZHY2aTzeGn7vcBJj0/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
View user's profile Send private message
Kenneth_Smith



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

PostPosted: Thu Mar 09, 2023 2:12 pm    Post subject: Reply with quote

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
View user's profile Send private message Visit poster's website
mecej4



Joined: 31 Oct 2006
Posts: 1885

PostPosted: Thu Mar 09, 2023 3:02 pm    Post subject: Reply with quote

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
View user's profile Send private message
Kenneth_Smith



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

PostPosted: Thu Mar 09, 2023 5:29 pm    Post subject: Reply with quote

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**2-2.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**2-2.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
View user's profile Send private message Visit poster's website
mecej4



Joined: 31 Oct 2006
Posts: 1885

PostPosted: Fri Mar 10, 2023 12:18 am    Post subject: Example of Gaussian Integration Reply with quote

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
View user's profile Send private message
Kenneth_Smith



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

PostPosted: Sat Mar 11, 2023 12:32 pm    Post subject: Reply with quote

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
View user's profile Send private message Visit poster's website
mecej4



Joined: 31 Oct 2006
Posts: 1885

PostPosted: Sat Mar 11, 2023 3:33 pm    Post subject: Reply with quote

Nice example, added to collection.

The SLATEC FZERO routine does not require an initial bracket, and your trial values lead to success.

In general, the root found is not necessarily the one that is closest to the trial value provided, but in this case that is not a problem.

It so happens that the NAG library has a routine to do exactly this calculation:

http://monet.nag.co.uk/nagexample/?nagroutine=s17alf+%3A+nagf_specfun_bessel_zeros&loaddata=Load+default+data+file&inputdata=S17ALF+Example+Program+Data%0D%0A+0.0++5++1+%3A+Values+of+A%2C+N+and+MODE%0D%0A

Another implementation, in the CERN library:

https://cds.cern.ch/record/2050906/files/c345.pdf

I gave the online NAG program the input
Code:
Bessel Roots
0.0,10,1

but it would not run (error code 127).
Back to top
View user's profile Send private message
Kenneth_Smith



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

PostPosted: Sun Mar 12, 2023 8:08 pm    Post subject: Reply with quote

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
View user's profile Send private message Visit poster's website
mecej4



Joined: 31 Oct 2006
Posts: 1885

PostPosted: Sun Mar 12, 2023 9:00 pm    Post subject: Reply with quote

Kenneth,

The argument IPVT needs to be an array of dimension (3).
Back to top
View user's profile Send private message
Kenneth_Smith



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

PostPosted: Sun Mar 12, 2023 9:25 pm    Post subject: Reply with quote

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
View user's profile Send private message Visit poster's website
mecej4



Joined: 31 Oct 2006
Posts: 1885

PostPosted: Mon Mar 13, 2023 1:13 am    Post subject: Reply with quote

Kenneth,

If it is any consolation, 7/22 is a good approximation to 1 / Pi.
Back to top
View user's profile Send private message
PaulLaidler
Site Admin


Joined: 21 Feb 2005
Posts: 7916
Location: Salford, UK

PostPosted: Mon Mar 13, 2023 10:10 am    Post subject: Reply with quote

The bug relating to test04 above has now been fixed for the next release of FTN95.
Back to top
View user's profile Send private message AIM Address
Display posts from previous:   
Post new topic   Reply to topic    forums.silverfrost.com Forum Index -> Support All times are GMT + 1 Hour
Goto page Previous  1, 2, 3, 4  Next
Page 2 of 4

 
Jump to:  
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