I am working on a complex Finite Difference code and testing it's results vs. a standard (published) analytical solution. The code is from a reputable source (US government agency) and has been 'widely' used. When compiling the code with the Salford v254, I get the correct result for the test problem. However with the latest Silverfrost v8 32b compiler I get a significantly different and incorrect result. Could anyone advise me if there are some initialization or precision differences between these 2 compiler versions and/or if there are any compile options I could use to address this.

The results should be the same if both are compiled with FTN95.

Are you using the same FTN95 command line options in both cases?

If you are using Ver 8 32-bit, then both tests should be using x87 computation. If they are different, then there may be different optimisation in place. (what compile options are you using) You need to determine if the differences are due to round-off variability or are significantly different. For significantly different, you need to identify where the differences come from. You could try FTN95 /64, as this is a different instruction set, that doesn't use x87 instructions.

I would suggest you use /checkmate to utilise the run time checking available in FTN95.

Unfortunately 'complex Finite Difference' suggests a difficult task

It might not be the last time that old legacy codes show up long hidden bugs.

As far as OI am aware I am using teh same options - here is the detail:

Here is the Salford 254 command line to compile

C:\Progra~2\salford\FTN95 mod_array.f95 /LIST mod_array.f95.lis

and than SLINK simply load all obj files and create exe file

For the Silverfrost v8 I used Plato w project properties as following ticked / added

Miscellaneous

• extra compiler options - /FIXED_FORMAT
• output file name - HST2D.exe
• output file type - EXE
• pre-process source files is ticked

Switches /FPP /CFPP /FIXED_FORMAT

all other compiler/linker options are NOT ticked

Let me know if there is a difference after all.

I did use the 32 b and 64 b compilation for Silverfros v8 and I also used release and checkmate.

All the above gave the same result ie. different from Salford 254.

I noticed the name 'HST2D' in your post. It so happens that I worked with the successor USGS code HST3D 2.2.16 in 2016. I found several bugs (at least 15 of them) in the source code, using FTN95 and /checkmate. There were variables that needed initialisation, arrays that were too small or were referenced with indices outside their bounds, and variables that needed to be SAVEd.

You may see some posts in the thread https://forums.silverfrost.com/Forum/Topic/2878 that touch HST3D, although in that thread we were more concerned with the (then) new 64 bit FTN95 compiler and speed.

Assuming that HST3D is a corrected and augmented version of HST2D, I suspect that HST2D may contain similar errors as well. I cannot locate the source of HST2D, but HST3D is available at https://wwwbrr.cr.usgs.gov/projects/GW_Solute/hst/index.shtml .

Such bugs are not unusual, and finding them need not cause loss of respect for the authors of the software. Note that the Silverfrost compilers may also have had bugs that affected the results, as you may see from the list at https://www.silverfrost.com/default.aspx?id=19 . A bug in the compiler can hide a bug in the software being compiled!

If you can provide a link to the HST2D source (the original from USGS or your modified version) along with instructions to reproduce the behaviour that you noticed, I may be able to help.

Quoted from John-Silver ... which if true is worrying. Optimisation should only touch the efficiency of execution not accuracy of the results surely. Or am I missing something.

For floating point calculations, Optimisation can change the instruction set or change the order of calculations. Especially for FTN95 /32, x87 register use can also change what is stored in registers and when values are transferred (truncated) to memory. This can result in different round-off errors, which then need to be assessed for their significance. The worst round-off might not occur at the end of the calculation, resulting in greater than expected errors. ( eg subtracting two very similar values with differing round-off errors ) For iterative solutions, this can change the number of cycles to convergence, which can induce a significant change in the calculated values due to an extra iteration. Optimisation is many faceted, definitely not just optimum accuracy.

Then if there is a bug in the code, eg uninitialised values, optimisation can make a meal of this.

Quoted from John-Silver Optimisation should only touch the efficiency of execution not accuracy of the results surely.

Only if exact arithmetic is used (i.e., integer, character and boolean variables) and no integer overflow occurs. If finite precision floating point calculations are performed, the usual rules of algebra, e.g., a+b-a = b, are not always obeyed. Here is an example program to illustrate this.

program inexact
implicit none
real a,b,s
integer i
a=1.0
b=2e-7
do i=1,5
   s=a+b
   print '(1x,i2,2ES15.7)',i,b,s-a
   b=0.5*b
end do
end program

FTN95, and Intel Ifort /Od (optimisations disabled) give

  1  2.0000000E-07  2.3841858E-07
  2  1.0000000E-07  1.1920929E-07
  3  5.0000001E-08  0.0000000E+00
  4  2.5000000E-08  0.0000000E+00
  5  1.2500000E-08  0.0000000E+00

Note that the second and third columns, which show b and (a+b)-a, do not agree at all.

If it does touch accuracy then it should be in the positive sense i.e. results should be 'better'.

Usually, the opposite is true. There is a trade-off between at least three factors: (i) speed of compilation and linking, (ii) speed of resulting program and (iii) accuracy of results. Many compilers provide 'aggresive optimisations' or 'unsafe optimisations'.

Sometimes, you may get more accuracy than you expect, but for the wrong reasons. For the above program, if optimisations are allowed, the results from Ifort are:

  1  2.0000000E-07  2.0000000E-07
  2  1.0000000E-07  1.0000000E-07
  3  5.0000001E-08  5.0000001E-08
  4  2.5000000E-08  2.5000000E-08
  5  1.2500000E-08  1.2500000E-08

Looks great, does it not? The compiler realised that the results can be precomputed at compile time, and produced a program that simply prints the precomputed results. The Fortran standard allows the compile time arithmetic to be completely different from the runtime arithmetic; in this case, the compile time arithmetic was probably done in double precision.

It is interesting what we consider to be a compiler bug.

This accusation is often applied to a new compiler, because the previous compiler tolerated, what we accepted to be my 'old' Fortran. A number of these include:

Un-initialised variables, because the old compiler set memory to zero.

Local variables going out of scope, because the old compiler defaulted all local variables to SAVE. (this is a big problem for multi-thread codes where SAVE should not be used)

Those who have used static allocation compilers or linkers may disagree, but to me these are two serious coding bugs that are repeatedly appearing in complaints about most modern Fortran compilers.

I am not immune to this accusation, as I struggle with Intel iFort's management of array sections, by using non-contiguous memory for array arguments. I am probably now wrong, but I have always assumed that arrays are contiguous, and the array argument transfers the start memory address of the array. (for me to break type mixing rules)

The point is that just because the old compiler got what you now accept as the correct answer, doesn't mean that there are no bugs in the code. The bugs are still there, but did not cause a problem with the old compiler.

John_C,

I'm something of a fan of static allocation, but SAVE in all its guises is an abomination, and it is the job of the programmer, not the compiler to initialise variables before first use. It's a lucky accident, no more, if variables are initialised automatically at run time.

Those coding bugs may be present in any code, and I'll give you a win at the argument because they surface in old code when it is being resurrected after a decade or more when nobody looked at it.

There's also a situation where the fixing of a compiler bug invalidates something that was accepted in the past, and now declares it to be an error. The nightmare situation is when a really old code has been 'bodged' ('botched' in some versions of English) and then becomes rather unreadable.

Eddie

Thanks for all comments so far.

As per my initial post, there must be an algorithm in the code that gets compiled by Salford v254 as intended, but with Silverfrost v8 somewhat differently making it working incorrectly.

Please review below document (link below) to see what the difference is. The code calculates the effect of an underground heat exchange cell where a negative heat influx reduces the temperature originally at 10.5 degrees C (Graph 1 - heat exchange cell). The second graph shows the temperature at an observation point at some distance.

https://1drv.ms/w/s!AuTT_gAwgmEIh4YoSh_vjzS_bwq7NA

As graph 2 (observation well) shows the calculation of the SF254 compiled code follows a published analytical solution very closely and the calculation of Silverfrost v8 shows a significant difference.

I looked at your results, which are (disturbingly) only slightly different.

Did you use ftn95 /64 or ftn95 32-bit ? 32-bit will use x87 80-bit calculations, while I think 64-bit will use 64-bit calculations, which have a lower accuracy.

Apart from that, I would be compiling with /checkmate to check for out of range addressing and possible undefined values.

I used for Silverfrost ftn95 v8 both 64 and 32 bit and obtained exactly the same (incorrect) result.

I also did compile Silverfrost ftn95 v8 both 64 and 32 bit checkmate This produced a number of warnings re un-initialized variables. After fixing this by initializing variables properly I still got the same incorrect result.

Note that the Salford ftn v254 compiled version, which produced the correct results, was NOT compiled with checkmate option. I am currently preparing a version that can be compiled w checkmate (needs the same fix to initialise variables as per above)

I finally managed to track down the issue after some further discussion with mecej4, who kindly recompiled the code.

It turns out the the version of Salford v2.54 which was passed on to me as a straight folder under c:\Program Files (x86) (ie not an install) compiled by default with double precision reals. I am not sure whether this is standard or influenced by the original setup. I am also not sure if the previous developed was aware of this. I had inherited BAT compile files from the previous developer which did not stipulate 'all double precision' compilation flags.

Silverfrost v8.3 compiles by default with single precision reals.

For most of the test problems I had at hand this did not cause any issue other than totally insignificant differences.

However, when I further analysed the test problem I posted (and that created the problem) did change the basic input via a unit conversion and single precision truncation significantly reduced the heat influx.

I am still not sure if the code require requires a blanket double precision for all other calcs, but obviously as a matter of safety it can't hurt. Most codes of this style I have worked with before had the double precision declarations explicitly inside the code, obviously to avoid the problem I ran into.

Lessons learned (for the umpteenth time): never make assumptions about defaults, and understanding the physics (which was pointed out to me by mecej4) is very valuable when tracking down bug.

All comments from various contributors to this thread are much appreciated.

These things are easy to miss. The original compilations may have used /config to configure the compiler to use /dreal by default.

Is there a standard fortran statement that will declare all reals as kind=2 (double precision) ?

The safest way may be to call SELECTED_REAL_KIND to get the KIND value for the compiler.

Quoted from jcherw Is there a standard fortran statement that will declare all reals as kind=2 (double precision) ?

First of all, note that kind numbers are NON-PORTABLE.

For real variables, the default kind numbers in FTN95 are 1 (32-bit real) and 2 (64-bit real). For FTN95, there is also the X87 10-byte real, for KIND=3. FTN95 provides the /ALT_KINDS option to enable the use of code that assumes kind numbers of 4 and 8 for 32-bit and 64-bit reals (with some limitation). Details may be seen in the FTN95 /help output.

As Paul wrote, in portable code you would use SELECTED_INT_KIND. For a project made up of a large number of source files, a somewhat shorter (but still portable) way is to put the following code into a separate module source file:

module floats
integer, parameter :: sp = kind(0.0),  dp = kind(0.0d0)
end module floats

In each subprogram of your project, you may add USE FLOATS and then change variable declarations from REAL to REAL(sp), if you choose single precision, or to REAL(dp) if you choose double precision.

If you choose this approach, where you have real constants in the code you may need to append '_sp' or '_dp'. You have the value of Pi specified in various places in the code, but you do not have enough digits to provide double precision for Pi.

Just for completeness, here is a small test program that clarifies the issue.

  Program testPrec
  

! Test precision of time conversion ! Heat flux (utu) unit coversion (/day to /sec) to comply with SI units ! which is used in remainder of program (solver etc)

  implicit none

  real*4 uqh4, utu4, res4 ! heat flux, conversion & result single precision
  real*8 uqh8, utu8, res8 ! heat flux, conversion & result double precision
  real*8 res48, dif8, tot8
  
  utu4 = 86400.
  utu8 = 86400.
  
  uqh4 = -.1036E+09
  uqh8 = -.1036E+09
  
  res4 = uqh4/utu4
  res8 = uqh8/utu8

  write (*,'(A25,F15.5, E15.5, F15.5)')  &
    ' Single utu4 uqh4 res4: ', utu4, uqh4, res4
  write (*,'(A25,F15.5, E15.5, F15.5)')  &
   ' Double utu8 uqh8 res8: ', utu8, uqh8, res8
  
  res48 = res4
  dif8 = res48-res8
  
  write (*,'(A24,1X,3F15.5)') ' res48 res8 dif8:', res48, res8, dif8
  

! Accumulated heat flux difference after 1000 days (run of simulation)

  tot8 = 1.D03*utu8*dif8
  write (*,'(1X,A45,9X,F15.5)') &
   ' Accumulated heat flux difference @ 1000 days:', tot8
  
  end

and the results

Single utu4 uqh4 res4: 86400.00000 -0.10360E+09 -1199.07410 Double utu8 uqh8 res8: 86400.00000 -0.10360E+09 -1199.07407 res48 res8 dif8: -1199.07410 -1199.07407 -0.00002 Accumulated heat flux difference @ 1000 days -1953.12499

The expected precision of 32-bit reals is seven significant decimal digits. You are taking the difference between two floating point numbers that differ only in their eighth digits; that difference in Q (between res48 and res8), as reported, is 2 or 3 X 10-6, in whatever physical units were used. You are then multiplying that by a large number (1000 days X 86400 seconds/day). The result of doing so, which is about 2,000, may appear large, but that is simply because of the units used. The fractional error between res48 and res8 is still about 2 X 10-8, and should not cause any worries.