More info on file and errors on code ported from VBA/Excel

[John_B_8]

Attached is the latest version of the ported code from VBA/Excel and a listing of errors that the compiler is throwing up

[John_B_8]

Sorry, the original query is awaiting the forum administrators' approval before it is posted. I'm trying to port some working VBA/Excel code to fortran, but I am falling foul of some of the Fortran rules , I have read around but I still can not see why the code is throwing up these particular errors.

Thanks in advance

[John_B_8]

I thought I could 'end' the main code then have the functions appended but that did not seem to work, too many end statements, but taking out the first end statement then seems to run into the problem of having funtions after executable statemnets

[John_B_8]

I thought I could 'end' the main code then have the functions appended but that did not seem to work, too many end statements, but taking out the first end statement then seems to run into the problem of having funtions after executable statemnets

[Les_Neilson]

John,

In the help look at CONTAINS and try having something like :

program console2

..... your code

contains

function poisp

...

end function poisp

function poislow
end function poislow

...

end program console2

But also check that the functions have a value when the error conditions occur. (e.g. poisp = 0.0 at the start of the function poisp)

You should also get into the habit of specifying "implicit none" in your code. It will help identify typos for example.

HTH

Les

[John_B_8]

Thanks Les,

In the meantime (without trying contain)  the code is now building, but is not working correctly. I named  the main program Poisson, then ended it  with an

 end program poisson statement,

 then placed my user defined function after that end statement.

Its not giving the right answers though. It is as if it it not going to those user defined function to compute the value of the functions.

[Steven_L_Intel1]

There are several problems here.  The first, as Les suggested, is that you are not ending one program unit before starting another. Les' suggestion of the contains is a good one - the contains line would go in place of the "199 continue".

Les also is hinting that you have argument type mismatches due to the lack of declarations in the functions. You also have the functions returning single precision REAL but you're assigning values to double precision variables - you should be consistent throughout.

[Casey]

John,

I have attached a modified version of your fortran code that addresses many of the problems Steve and Les have pointed out with the goal of being helpful and illustrating what they are saying.  The code sample compiles and runs.

Changes I made include:

1) removing goto statements in favor of better use of if constructs and return statements.

2) using a parameter to declare the precision of your reals and using this throughout your code and on all your floating point constants.

3) Put the "contains" line where it needs to go and moved "end program" where it needs to go.

4) Re-worked the variable declarations on your internal functions.  You only need to declare the dummy arguments and local variables.  It is not neccesary to declare variables you are not using, nor should you declare the names of your internal functions.  A second note, your internal functions can access variables in their host (your program), though I do consider it good practice to pass those variables (as you are doing) rather than accessing the host scope.

5) Took the simplest approach to making your argument types consistent.  It would be more proper (in my opinion) to make x1, x2 integers to remove the need for checking that in your functions, but that would require reworking more code than I cared to spent time on.

Take a look at the code and verify it performs properly and feel free to ask questions about my changes.

[John_B_8]

Thanks Steve and Casey,

I was going to revert to old textbook example of the user defined function average to discover how to call such functions.

I'll check the type consistency event should be integer, although Poisson events are always integer, their mean can and often is real. Double Precision floating point is needed in most variable instances, otherwise the results will be unreliable.

Thanks again for the code, I'll read and digest and apply and see how it goes.

[John_B_8]

Casey,

I'm afraid my Fortran is a bit stuck in the past, normally I dip into Fortran 77 'Principles of programming' by Jerrold L. Wagner, then use Clive Page's, 'Professional Programmer's guide to Fortran 77', with another short document he did ,'Fortran90 for Fortran77 Programmers', and that has in section 2.3 'new forms for specification statements' info about new Type statements:-

The DOUBLE PRECISION data type is now just a special case of REAL so all facilities are identical (???); this means that double precison complex is fully standardised.

(In my old '77 manual it appears that there was a somewhat inconsistent prohibition of the use of  DOUBLE PRECISION in certain operations with complex numbers in F77.)

"new form with double-colon allows all attributes of variables to be specified at once."

There was no mention however of real(kind=wp).

In the attached document I have put a selection of  information that I found about this topic. It may be of use to others.

I found that 'wp' is short for working precision in the section that I lifted from a NAG document (attached), this gives an explanation, but I will have to read it very carefully in order to understand it, as it is not described as clealry as one would hope.

[John_B_8]

I have also found a document online by Bo Einarsson and Yuij Shokin, "Fortran 90 for the Fortan 77 Programmer",

(aside- copy and paste does not seem to work when posting to this forum. I get a pop-up window, but it is non-functional.)

Anyway, section 13 of the above document 'The new Precision Concept'   the statement REAL(KIND=WP) is explained quite nicely.

[John_B_8]

Arrgh, my earlier comment/post on REAL(KIND=WP) has been queued again by the admin procedures on this forum, so they will be out of order, sorry.

[John_B_8]

Casey,

There must be a mathematical error in part of my porting of the VBA code, which I'll try to fix, but the part that is working is giving a slightly different answer to the VBA/Excel code, which also using a similar 'R' statistical package code to compute the same value, seems to be the more correct answer, than mine. These values are notorious for the demands they place on the precision of numerical calculations, ref. L. Knusel, 'Computation of  Statistical Distributions - documentation of the Program ELV'. I see that there is now, sort of a double double precision data type in Fortran, and your wp 'trick', seems to be a very convenient way to set the precision. What are the limits on the acceptable values for wp ?

Thanks,

John

ps in the Einarsson document they give the following example

     INTEGER, PARAMETER       :: WP = SELECTED_REAL_KIND(14,300)

    REAL(KIND=WP)                  :: working_precision_variable

But they note "Regrettably now we have to give all floating point constants (with?) the additional _WP,

for example

            REAL(KIND=WP)           :: PI

            PI =3.141592653589793_WP

while sine the the intrinsic functions are generic, they will automatically use the corrcet data type and kind, which means that the argument determines which 'kind' the result should have (usually the same as the argument).

With this method you will in practice obtain double precision on systems based on IEEE 754 abnd single precision on coputers like Cray or compures based on the Digital Equipment Alpha-processor, which in all cases means a precision of about 15 significant digits."

[John_B_8]

Casey,

I sorted out the error in the PoisP routine where I should have had .GE. I had incorrectly written .GT. in the logic test. That got the probablities to match those of R and of some online probablity calculators.  The PoisHigh function was also working ok, but not PoisLow, on its first few circuiots round the loop it was reporting that x2 was not a whole number.

I looked through it and could not find any errors in comparison with the original VBA/excel code, but one thing I did note was that in PoisP I was using the same variable names as for the main program unit.

Now like subroutines in the function itself we can use new names and thus not worry about any local change to a variable affecting global or main program variables.

i.e. if we send PoisP(x,y,z) we can use in the function PoisP(a,b,c) the function just uses the position of the variable to translate between them so that

in this case global x = local a, global y = local b etc.

So I thought I would use this idea and simply rename the repeated names within the PoisP as x1 -> px1, x2 -> px1, z -> pz.

So now there should be no danger of 'crosstalk' causing weird (but actually logical) behaviour witin the code.

Problem was after that simple change, the code refused to compile giving a linking error cannot open file 'Debug/Casy.exe'.

So I'm stuck and I have no idea what is happening.

Is it something simple like I need to rename a file in the debug folder ?

ANSWER NO this latter problem is because I had an instance of the code running (command window open). When I exited that DOS/command box. The code has compiled, now to see what happens !

The offending code is attached,

ps I have now made the majority the function variables 'unique' so using the passing of dummy arguments to subroutines method.

[John_B_8]

PoisLow still not working, everything else looks ok, any ideas why it is failing mathematically, the porting of the functions looks sound, perhaps I'm getting code blindness though ?

NB confint.xls can be used to check the correct working of the fortran code.

[John_B_8]

Now this gets strange.

Instead of using 'function' I thought I woulod code up a version using a PoisP as a subroutine to PoisLow just to test the code in this way. Doing it this way seems to have worked ??!!!

How come ?

[Casey]

John,

There are quite a few comments to work through.  Just a note, my perspecive is a bit opposite to yours, and I primarily deal with fortran 90 and use features of 03 and 08 when I can and do not deal much with '77 code.  Books I keep on my shelf related to fortran are "Modern Fortran explained" by Metcalf and "Modern Fortran in Practice" by Markus.  

The (kind=wp) and the _wp attached to constants is just a convientent way to ensure everything is the same type and kind in a way that you can easily change the precision of all the variables in one place.  You'll see I set wp=kind(1.d0), which just makes wp mean double precision.  You can also select an arbitrary precision by using wp = selected_real_kind(sig_digits, exponent_rang), and let fortran worry about what type it will need to hold that number.  E.g. wp = selected_real_kind(10,40) will ensure your variables of that kind have enough precision to represent 10 digits and an exponenet range 10*(+/- 40).  The name "wp" is arbitrary, but as you found is a convention for "working precision".   Your book might call adding _wp to your constants regrettably, but knowing that all math, named variables and temporary variables will be of your desired precision means that you are not losing precision unkowningly somewhere.  

I'll also comment on scope.  In your internal functions (functions between "contains" and "end program") you have access to all variables in your host scope, but if you have dummy arguments or local variables of the same name, you will block access to those in the host scope and only have access to the local ones.  The dummy argument names should not be a problem as they were named (nor are they with the new names).

I don't have time now deal with the mathematical accuracy issues, but if I find some later I may take a shot at it.

[John_B_8]

Casey,

I found the problem (although I'm not sure why it causing the problem) sort of by accident.

This might sound a bit lame, but I use in conjunction with the microsoft intel compiler an old freebie Fortran compiler called Force. Although it has one foot stuck in the old punch card environament, it is very good when it comes to defining errors. Anyway it siad that the BigNum 'integer' !!?? was too large. This has the form of 1 followed by ten zeroes, now in the excel/VBA this causes no problems and no warnings at any stage when doing the lower limit with 3 as the number of events and 0.025 as probability limit. In Fortran this triggered the repeated warning "x2 not a whole number" (i.e. when the PoisLow function was in turn using the PoisP function) and then gave a number way too high. So in Force to get it to compile I had to drop a zero off the length of BigNum, when I did that the subroutine version ran without any warnings. Then re-loading the function version back into the Intel complier and knocking a zero off BigNum, got that to run with out the warnings and gave the correct answer.

I've no idea why that number is ok in excel/VBA but causes probelms in Fortran.

Now I can move on to getting the code to work with input and output files instead of interactive input.

There aren't going to be any open unit surprises, are there ?

Thanks for your help !

[Casey]

I'm not surprised by that actually.  I did not mention it (i forgot), but in the code I sent you I took the easy way out and declared BigNum to be a real(kind=wp) rather than an integer.   If you want BigNum to be an integer you can make it type integer(kind=selected_int_kind(11)) to ensure it is allocated enough storage to hold a number that big.  A default 4 byte integer can only hold numbers +/- 2**31 (~2 billion), which is not enough to hold your BigNum. 

[Casey]

John B. wrote:

Now I can move on to getting the code to work with input and output files instead of interactive input.

 There aren't going to be any open unit surprises, are there ?

 Thanks for your help !

The only potential surprise I can think of is regarding record lengths in unformatted files.  If you are specifying a record length (recl=) in bytes, you'll want to add the compile flag /assume:byterecl because the intel default is not bytes (as it is on many other fortran compilers).  

[mecej4]

In this thread much has been said about the various versions of Fortran code that is intended to calculate the  Poisson cumulative distribution function and its inverse. Another side of the question is whether the mathematics is correct. Unless I am mistaken, the calculations are missing an important element: the Poisson density function is e^-λλ^x/x! , but I see no instance of the EXP function in any of the codes. What is the explanation?

My second comment is concerning the use of bisection of the interval in solving for the inverse. It is simple to code, but very inefficient, taking roughly three calls to POISP() to obtain each decimal digit in the answer.

The third comment is regarding BigNum.  Because of the exponential function in the PDF, for values of λ from 0 to 20, using BigNum = 40 is a good approximation to the usage of BigNum = ∞. Therefore, there is little reason to worry about using 8-byte integers, etc. Even single byte integers are sufficient!