Welcome to our forums
Thanks for the feedback.
Regarding the subscript overrun in XSDRIV1: it is my mistake. Please change the allocation statement to:
allocate (work(lenw), y(n+1))
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.
! 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
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'.
mecej4, I had some fun yesterday evening examining IFLAG from FZERO and came up with another example for you to consider.
! 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
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:
Another implementation, in the CERN library:
https://cds.cern.ch/record/2050906/files/c345.pdf
I gave the online NAG program the input
Bessel Roots
0.0,10,1
but it would not run (error code 127).
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.
! 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
Kenneth,
The argument IPVT needs to be an array of dimension (3).
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!
Kenneth,
If it is any consolation, 7/22 is a good approximation to 1 / Pi.
The bug relating to test04 above has now been fixed for the next release of FTN95.
Thank you, Paul.
mecej4,
Below is a link to another example (a well known problem in the electrical world) which uses the some of the SLATEC eigenvalue and eigenvector routines.
https://www.dropbox.com/s/l8mgxhbeiyjdix5/ieee1.f95?dl=0
I resisted the temptation to plot the mode shapes with %pl. The results all look to be correct to me.
Here is a better, tidied up version:
Kenneth,
I think that a reference/link to the source of the 'magic' constants in arrays HC and TS should be added. How about https://link.springer.com/content/pdf/bbm:978-1-4615-5633-6/1.pdf ? Table A3 there lists the constants.
The comments in the code refer to routines SPPDI and SDEEV, but there are no calls to those in that code.
The code would be a bit easier to read and understand if array operations were used. For example:
! Do this to avoid very small negative numbers which cause problems with SQRT later
where (abs(hm) < eps) hm = 0.0
where (abs(km) < eps) km = 0.0
mecej4,
A revised version can be found here:
https://www.dropbox.com/s/9s6jenwi8qvc2xm/ieee3.f95?dl=0
At the end of the revised code, a clearwin+ window is opened which outputs all the results in graphical format, which is appropriate for a dll intended for use with FTN95
I have added a description of the general method for solving this problem, and more comments in the code, which is now separated into a number of distinct steps.
In this particular problem, the parameters for the rotating mass inertia and shaft stiffness (the “magic” numbers) are not in standard units, but expressed in the per unit system used by electrical engineers for dynamic studies. The definition of these units is explained, as is the necessary departures from the “general” method. As a consequence of this, the intermediate values of the eigenvalues reported by the program do not correspond to omega**2. The system base frequency is correctly factored in at the end when reporting the actual torsional frequencies in Hz, i.e. this takes us back to the results that would be obtained in we started with rotating mass inertia J and shaft torsional stiffness K in their SI units.
The parameters in the benchmark model correspond to the large turbine generators at Mohave, Modes 2, 4, and 5 can be easily excited if their complementary electrical frequencies are present in the connected power system, due to the high participation of the generator mass in these modes as show by the mode shapes (eigenvalues). Two shaft failures occurred at Mohave in the early 1970s, which led to the development of the benchmark model for analysis purposes. If the electrical system produces frequencies complementary to one or more the torsional modes you can very rapidly fatigue the mechanical string.
Upon reflection, this was not perhaps a good example, as somebody blindly reusing the code for a torsional problem with J in kg,^2 and K in Nm/rad will run into problems – system electrical frequency is not relevant with this approach. Nevertheless it does demonstrate the basic methodology, the use of some of the SLATEC routines, and was an interesting problem for the “wee small hours” of the last few days. A refresher – I’ve not looked at this problem in earnest for at least a decade. I’m a bit fatigued after this one.
PS I can see a way of fixing all of these issues. But it will be some time before I can implement this, as there a a few other urgent things on the to-do list at the moment.
Kenneth, thanks for writing up a detailed explanation. I found a set of slides for a presentation on the Mohave turbine failure: https://www.geenergyconsulting.com/content/dam/Energy_Consulting/global/en_US/pdfs/GE_-_SSO_Risk_Analysis_Protection_and_Mitigation_Techniques.pdf
One of the figures in the presentation shows mode shapes, and it struck me that it was very similar to the graph that your program produces. Neat.
A more detailed discussion of SSR (Sub-Synchronous Resonance):
Since the installation where the failure occurred is in the US, recasting the description into SI units could be considered sacrilegious.
Mecej4,
Thanks for sharing both of these links.
The Energiforsk document appears to be mostly drawn from an earlier EPRI report which is on my desk, updated to reflect the more recent advances in SSR protection methods.
The GE document which you shared confirms something I have long expected. The use of series filters to block troublesome frequencies takes up a lot of physical space. The photograph of the filters at Navajo (plant now decommissioned) confirms this. A very good example of a purely analogue (R/L/C) solution to a complex problem.
Some pictures of the more conventional shunt filters on the IFA1 HVDC link on the UK side at Sellindge which also have a large footprint. The rack mounted capacitors and the cylindrical air cored reactors can be clearly seen.
https://www.geograph.org.uk/photo/6964452
My first attempts to understand the torsional problem were considered sacrilegious by the mechanical engineering team; I converted the mechanical system to its electrical equivalent circuit (series L to represent the shafts and C to ground to represent the inertias). The torsional frequencies appeared as local maxima/minima in the impedance vs. frequency scans at each of the nodes in the equivalent circuit, and from the off diagonal elements of the full impedance matrix at these specific frequencies the mode shapes were derived. Even although we agreed that the system differential equations were essentially the same, they remained sceptical for a long time!
There are a lot of fast Fourier transform subroutines in this library. Can anyone share experience and demo examples with these routines or any other libraries?
Need the ones which will deal with
Dan, it is nice to see you expressing interest in the SLATEC DLLs for FTN95.
Of the 54 tests that come with the package, Test51 is devoted to the FFT routines. The comments in the code indicate that the FFT routines originated in NCAR (U.S. Natl. Center for Atmospheric Research), and are dated 1985 or earlier. From the comment lines in the pertinent source file, FFTQX, we are able to note:
! This Program Tests The Package Of Fast Fourier
! Transforms For Both Complex And Real Periodic Sequences And
! Certain Other Symmetric Sequences That Are Listed Below.
!
! 1. Rffti Initialize Rfftf And Rfftb
! 2. Rfftf Forward Transform Of A Real Periodic Sequence
! 3. Rfftb Backward Transform Of A Real Coefficient Array
!
! 4. Ezffti Initialize Ezfftf And Ezfftb
! 5. Ezfftf A Simplified Real Periodic Forward Transform
! 6. Ezfftb A Simplified Real Periodic Backward Transform
!
! 7. Sinti Initialize Sint
! 8. Sint Sine Transform Of A Real Odd Sequence
!
! 9. Costi Initialize Cost
! 10. Cost Cosine Transform Of A Real Even Sequence
!
! 11. Sinqi Initialize Sinqf And Sinqb
! 12. Sinqf Forward Sine Transform With Odd Wave Numbers
! 13. Sinqb Unnormalized Inverse Of Sinqf
!
! 14. Cosqi Initialize Cosqf And Cosqb
! 15. Cosqf Forward Cosine Transform With Odd Wave Numbers
! 16. Cosqb Unnormalized Inverse Of Cosqf
!
! 17. Cffti Initialize Cfftf And Cfftb
! 18. Cfftf Forward Transform Of A Complex Periodic Sequence
! 19. Cfftb Unnormalized Inverse Of Cfftf
I do not think that SLATEC should be regarded as a replacement for specialized or more recent packages for solving problems in any of the numerous areas that it covers. On the other hand, it is free, high quality but old, and the availability of DLLs may make it easier for someone such as yourself to try it out. The download and documentation links that I provided above should get you started.
As for demos, other than the 54 quick checks included in the package, we have a handful, including those that Kenneth has contributed. If you can prepare one for FFTs, go for it! As Kenneth has demonstrated, less than 10 lines of Fortran code will do to produce very good graphs using FTN95. Such demos may persuade prospective users to continue with Fortran instead of switching to Python, Julia, etc., for the graphics.
mecej4,
Final version of the eigenvalue example can be found here.
https://www.dropbox.com/s/lggp8hvaufa7hm5/ieee4.f95?dl=0
Incorrect labels (mode numbers) in the graphics output corrected and the description now contains a physical interpretation of the mode shapes.
Thanks, Kenneth.
I am waiting for Paul's fix for the complex arithmetic bug to be released. When that occurs, I can rebuild the SLATEC DLLs and make them available, along with the modified sources.
In the meantime, I am replacing a large number of code segments in the SLATEC sources similar to
Ws(iw+1) = zero
call scopy(n,Ws(iw+1),0,Ws(iw+1),1)
by the modern Fortran equivalent
Ws(iw+1:iw+n) = Zero