URGENT Help Needed with Some Questions on a FORTRAN Code

[Deb_Chatterjee]

All,

I need some help in understanding a why a variable is not givingits correct numerical value in my present FORTRANcode which I am developing for calculations related to EMwave propagation in layered media. I am just copying therelevant lines from the source code.

.
.
.

[bash]CKP=-JJ/CR0
CZ=CON+(CON/ERF1)
write(*,*)' CZ=',CZ !!! Debugging line !!!
write(*,*)' CKP=',CKP,' before multiplied by CZ' !!! Debugging line !!!
CKP=CZ*CKP
write(*,*)' CKP=',CKP,' after multiplied by CZ' !!! Debugging line !!![/bash]

.
.
.

When I now start to look at the numerical value of the CKP variable, I get some startling results. I have checked the subroutine and the declaration of the variables CZ, CKP, JJ, CR0, CON, ERF1 are all correct (i.e. declared as COMPLEX*16 - double precision complex).

Before the (complex) multiplication the values (approximately) are:

CZ=1.306+j3.06E-06

CKP=-1.345-j0.5742

(These are approximate, because I have rounded off while typing.)

The funny part is that just after multiplication, I am again typing out the value of CKP. This time the result (???) from is:

CKP=0.000000000-j0.75

The CORRECT result, if you multiply both CZ and CKP using the standard scientific calculator you get

CKP=-6.05E-08-j0.75

I am wondering why this discrepancy - the real part is typed out as zero instead of being negative real ?

Can anyone please point out what should I do next ? Please help !!!

- thanks,

Deb Chatterjee

[JohnNichols]

Deb:

Suggestions:
One
If you have Maple or Mathematica try in to see if you get the same answer.
Two
Scientific calculators can give different answers, I use this sometimes in my exams, as the HP 49g can give different answers at the eighth + place from the TI, and it drives the students mad as they have problems with precision and accuracy. never met a TI calculator I did not hate.
Three
More beer.
Four:
Put a bit more of the code into a F90 file and let someone try it.

Regards

John Nichols

[mecej4]

> I am just copying therelevant lines from the source code.

It is possible that there is more that is relevant. In cases where the code is not doing what you expect, not only are bugs in the code to be suspected, but also your judgement of what is relevant.

Before the (complex) multiplication the values (approximately) are:

CZ=1.306+j3.06E-06

CKP=-1.345-j0.5742

(These are approximate, because I have rounded off while typing.)

Please do not do such gratuitous rounding in cases where you are trying to pin down the value of a number that is very small compared to others in the same expression.

You probably made a mistake, since the product of these values of CZ and CKP should be

-1.756568 - j 0.7499093

The real part does not agree with either the value that you say you got with your calculator nor with what you say the program gave you. It is hard to believe that with list directed output the program gave you exactly zero.

[Deb_Chatterjee]

You are right the values of the numbers were as follows:

(before multiplication)

CZ=1.30581039449541+j0.305810394495413E-00004

CKP=-0.134475289582368E-00004-j0.574209490914052

After multiplication of the above two complex numbers I get (from the FORTRAN code):

CKP=CKP*CZ=0.000000000000000 - j0.749808722264728

I used the format G25.15E5,2X,G25.15E5 to print a complex number.

Now I (re)multiplied manually the two numbers on my scientific calculator (available on my desktop). This gives the following (using the preceding numerical values for CZ and CKP):

CKP=CKP*CZ=(-1.7559923093943645674499293088E-05+1.75599230939436508778691243476E-05)
+j(-0.7498087218540979054148423532-4.1123941357068860170861677984E-10)

CKP=5.2033698312596E-21 - j0.74980872226533731898553095490862

My concern is regarding the real part of CKP. I know it should be very small for the input geometry I entered. But it cannotbe exactly zero. So, is there any way I can print the very small (yet non-zero) real part of CKP ?

[Steven_L_Intel1]

Try using E24.16E3 insead of G25.15E5 in the format. When G format switches to F format, you can get more rounding.

[mecej4]

>But it cannotbe exactly zero.
It is not exactly zero unless you used an appropriate format. You may be misinterpreting the output obtained using fixed format. Hint: there is no E00 or E000 or E+00 after the mantissa. To fix the idea in your mind, try the following

real :: x,y
...
x = 0
y=1e-6
write(*,10)x,y
10 format(1x,2F10.3)
...

> So, is there any way I can print the very small (yet non-zero) real part of CKP ?

As Steve has already told you, do not use F or G format for such situations, but use E, ES or EN.

You have another problem in the code snippets (these may or may not be present in your real code) : the variables are declared COMPLEX*16; however, the constant values assigned to them are going to be trimmed to about seven digits since they are single-precision. Subsequently, subtractive cancellation destroys any pretense of accuracy in the result, giving 1.53E-11 for the real part of the product. Compared to what the calculator gave you, this is off by 10 orders of magnitude. If you really aim to use this number, you have to be more careful in your work.

[Steven_L_Intel1]

Now that I think about it some more, the value probably is exactly zero. With G format, if the value was very small, say 1E-100, it would be displayed in E format. Since it was displayed in F format, the value was probably zero. It may have been rounded to zero during some computation. Without seeing a buildable and runnable test case, I can't be sure.

[mecej4]

Until the original poster recognizes the importance of reporting the code in toto, we can never be sure. If he had, for example,

cz=cmplx(1e-20,0.5)*cmplx(1e-20,0.0)
write(*,20)cz
20 format(1x,(G15.3,2x,G15.3))

underflow would occur and a zero would be printed even with E or G format.

[Steven_L_Intel1]

Change the project property Fortran > Floating Point > Floating Point Model to either "Source" or "Strict". What is happening is that the compiler is optimizing the multiplication in a way that causes the very small real part to go to zero. I haven't studied the code enough to see which operation made this happen, but I assume it is some sort of cancellation error due to the order of operations.

When you do this, you'll get the results you want. Of course, a value in the 10E-21 range is still pretty darned small.

[Deb_Chatterjee]

Steve,

It WORKED !!!!! Thanks a zillion !!!!

- Deb

P.S.: Can you please let me know my knowledge with coding in FORTRAN ? I am assuming it is pretty bad. I had exposure to FORTRAN 77 and some F 90, and I am a EE with background in Numerical Methods with applications to EM Theory. I have almost no sophisticated knowledge of coding like EE folks of today who takemore sophisticatedcourses in Computer Science to sharpen their programming skills. So, given that and what you saw in my code, what improvements should I make in my code writing and how ?

[Steven_L_Intel1 · Accepted solution]

Deb, here are some suggestions. Your code doesn't look bad, but you're missing out on some modern usages that can help you along.

• Use free-form source instead of fixed-form. This means a file type of .f90, use of ! for comments and & at the end of lines for continuation. Fixed-form source can mask coding errors and is restrictive as to line length.
• Use contained procedures rather than separate program and subroutines. This will also allow you to use uplevel references to variables rather than using COMMON. Or, better, put all the subroutines in a MODULE and have module variables instead of COMMON. But if you must use COMMON, put the COMMON definitions in an INCLUDE file and INCLUDE them, rather than repeating them in each routine, as that is error-prone.
• Instead of using REAL*8 and COMPLEX*16 (which are extensions), use REAL(DP) and COMPLEX(DP) where DP is defined as "INTEGER, PARAMETER :: DP = SELECTED_REAL_KIND(15)"
• Use generic intrinsics such as SQRT and ATAN2 instead of specifics CDSQRT and DATAN2. Generics will assume the proper datatype based on the arguments. Using the specifics can hurt maintainability of the code.
• Don't use INTRINSIC and EXTERNAL statements unless required
• Add IMPLICIT NONE to all program units
• Use DO..END DO rather than labeled CONTINUE statements
• Indent your IF..THEN..ELSE..END IF blocks to make it clear what the nesting level is
• Add more comments!

[mecej4]

I posted a reply to your #8. However, since #8 was marked "private", my reply is also classified as such, with no sayso from me. If you consent, we can request Steve to change my reply to "public".

By looking at your expressions in the program, it was possible to establish that the quantity in question, CKP, must be purely imaginary when tan1=tan2, as is true for the input data that you gave. Please see my reply #14 for details.

[Deb_Chatterjee]

Mecej4:

Your comments are correct on Tan1=Tan2 and hence on CKP. The following is in order:

1. I cannot make the code listing public because this code (under development) is a part of US Navy work.

2.Steve can see our comments because I had also sent him the code listing. Sure Steve can see our communications. He made some suggestions which I did not know and I'll implement them this week.

This code and its future reincarnations shall run into thousands of lines and in such cases my codes (developed for research work) become unmanageable. I need (per Steve's suggestion see his reply # 13) to modualrize my code. There will be several subroutines and IMSL based integrations of rapidly oscillatory functions etc., and I have to plan on how best to modularize. (In fact, Iam not an expert programmer, but implementing my own formulation gives me some fun [??] even at this50 years of age.)

Regarding IMSL Bessel functions (for complexarguments and integerorders), I disagree with you. I am not saying that these Bessel functions are wrong, but there are better algorithms which havebeen developed by researchersand they need to be used. Also, I am not too fond of recursion methods because they result in round-off errors.

Anyway, yours (& Steve) inputs/feedback is most welcome. Let me implement them andrun my code.

One thing: how to convert a *.FOR code to a *.F90 one ? I really don't know. Any quick pointers ?

- Deb

[mecej4]

> I cannot make the code listing public

I did not ask for that. I want only my response, which has none of your code, to be made public.

> I am not saying that these Bessel functions are wrong, but there are better algorithms which havebeen developed by researchersand they need to be used.

You are not using one of those anonymous "better algorithms". You are using numerical quadrature, which can rarely match (either in speed or accuracy) a good polynomial or Pade' approximation.

Also, I am not too fond of recursion methods because they result in round-off errors.

Only if the recursion level is quite deep and the error amplification factor is >> 1. For the small orders of Bessel functions that you use, this consideration does not apply.

However, you have the right to reject any and all advice given here.

[Steven_L_Intel1]

I cannot make a private reply public. What you can do is reply to one of the public posts, and that will be public.

[mecej4]

> What you can do is reply to one of the public posts, and that will be public.

Feature noted; thanks, Steve. I'll leave things as they are as far as this thread is concerned -- we don't want Navy Seals showing up in our front yard tomorrow.

[Deb_Chatterjee]

I am open to all and any advice. Regarding Bessel functions, and this specific code uses low orders only, Iexperienced that judicious use of asymptotic approximations, mindful of their regions of applicability and validity, are very useful.I was using the IMSL integration routine here to see if the numericalvalues (for real arguments) agree with ones in Abramowtiz and Stegun. In the finalversion of the code - which will take time - these Bessel functions will be approximated by their large argument formsand hopefully closed form solutions will be obtained for much complicated situations. Thiscode is the precursor to thelatter one to be born.

One question: how do I convert a *.FOR code to the *.F90 code asSteve suggested ? (I really need to know.)

- Deb

[Steven_L_Intel1]

You don't have to convert to free-form source, but if you wanted to, change the file type to .f90. Replace C in column 1 comments with !. If you use continuation, replace the character in column 6 with an & and add an & at the end of the previous line. (The & in column 6 is optional, and doesn't have to be in 6, but any nonblank character in 6 needs to be removed otherwise.)

[mecej4]

> One question: how do I convert a *.FOR code to the *.F90 code asSteve suggested ? (I really need to know.)

There are commercial tools to do restructuring, adding interfaces and declarations, but it will be difficult for you to debug the converted program without knowing Fortran 95. The use of IMSL (for which you do not have source code and may not have Fortran 9x modules for the interfaces) is an additional complication.

Given the small size of your program, I suggest that you obtain a good Fortran 95 book, and with its help convert the program manually, one subprogram at a time. At the end you will have not only converted your program but will have learned Fortran 95 along the way.

[Deb_Chatterjee]

mecej4:

Will do as you have recommended.

Can you recommend me a good FORTRAN 95 book, now that you (and Steve) have a good idea of my FORTRAN background ?

- Deb

[mecej4]

A compact book that I have found suited to my inclinations:

Fortran 95/2003 Explained, Metcalf M., Reid J. and Cohen M., 2004, Oxford University Press.

There is a more recent edition , Modern Fortran Explained, also from Oxford University Press.

There are older editions that cover only Fortran 95/90.

You can see a list of other books at www.fortran.com .