- Mark as New
- Bookmark
- Subscribe
- Mute
- Subscribe to RSS Feed
- Permalink
- Report Inappropriate Content
Hi there
I think I found a bug in Visual Fortran Compiler 10.1.025. It is about a subroutine used to solve nonlinear equations with MFP approach (see below). The brief code is to apply equation Y = X^2 - 25 just for checking purpose.
The problem is that when I add two printing commands on lines 10 and 20, the code works perfectly, but when I omit them it does not! In fact, it can be checked by inputting 0.001, 100, 1 and 10 respectively to get X= 5.
Could someone let me know if it is a bug or I made a mistake? I compiled the subroutine with another compiler and it worked even after omitting the two printing commands.
Many thanks
Hamid
1 READ*, ERS, IMAX, EPSDPER, EPSDMAX
EPSL = EPSDPER; EPSU = EPSDMAX
CALL MFP(EPSL, EPSU, EPSD, ERS, ERA, ITER, IMAX)
GO TO 1
END
FUNCTION F(X)
F = X**2-25
RETURN
END
SUBROUTINE MFP(EPSL, EPSU, EPSD, ERS, ERA, ITER, IMAX)
ITER = 0
FL = F(EPSL)
FU = F(EPSU)
10 PRINT*, 'FL=', FL, 'FU=', FU
DO
EPSDOLD = EPSD
EPSD = EPSU - FU * (EPSL - EPSU) / (FL - FU)
20 PRINT*, 'EPSL=', EPSL, 'EPSU=', EPSU
FR = F(EPSD)
ITER = ITER + 1
C print*, 'iter=', iter
C100 PRINT*, 'EPSD=', EPSD
IF (EPSD.NE.0) ERA = ABS((EPSD - EPSDOLD) / EPSD) * 100
PRINT*, 'ERA=', ERA
TEST = FL * FR
IF (TEST.LT.0) THEN
EPSU = EPSD
FU = F(EPSU)
IU = 0
IL = IL + 1
IF (IL.GE.2) FL = FL / 2
ELSE IF (TEST.GT.0) THEN
EPSL = EPSD
FL = F(EPSL)
IL = 0
IU = IU + 1
IF (IU.GE.2) FU = FU / 2
ELSE
ERA = 0
END IF
print*, 'ERA=', ERA
IF (ERA.LT.ERS .OR. ITER.GT.IMAX) EXIT
END DO
AMFP1 = EPSD
PRINT*, 'X=', AMFP1
RETURN
END
- Mark as New
- Bookmark
- Subscribe
- Mute
- Subscribe to RSS Feed
- Permalink
- Report Inappropriate Content
Hi Rich
Appreciated very much for your kind consideration of the problem.
The posted program gives the same wrong answer (X= 3.181818) after commenting the printing statements. However, as you annotated EPSDOLD is not initialized (Jim also mentioned the same about EPSD - thank you Jim). So, when I initialize EPSD or EPSDOLD to zero, the program works perfectly even if the floating point numbers are not double precision.
Interestingly, when I do not initialize them to zero and print them, they are zero! But sometimes the program does not work correctly depending on the printing statements. For example, in the following code, when lines 31, 32 and 33 are stated, the result is right:
.001
100
1
10
epsd 0.0000000E+00
epsd 0.0000000E+00
epsd 0.0000000E+00
epsd 3.181818
epsd 4.310345
epsd 5.142132
epsd 4.989630
epsd 4.999855
epsd 5.000141
X= 5.000000
but, when just one of the printing lines is commented, the result is wrong:
.001
100
1
10
epsd 0.0000000E+00
epsd 0.0000000E+00
X= 3.181818
If you could forward it to the developers, it may be interesting for them. The compiler I am using is: Intel Visual Fortran Compiler Professional for applications running on IA-32, Version 11.0 Build 20080930 Package ID: w_cprof_p_11.0.061.
FUNCTION F(X)
F = X**2-25
RETURN
END
SUBROUTINE MFP(EPSL, EPSU, EPSD, ERS, ERA, ITER, IMAX)
ITER = 0; IL = 0; IU = 0; !EPSD = 0
31 print*, 'epsd',epsd
FL = F(EPSL)
FU = F(EPSU)
32 print*, 'epsd',epsd
DO
33 print*, 'epsd',epsd
EPSDOLD = EPSD
EPSD = EPSU - FU * (EPSL - EPSU) / (FL - FU)
FR = F(EPSD)
ITER = ITER + 1
IF (EPSD.NE.0) ERA = ABS((EPSD - EPSDOLD) / EPSD) * 100
TEST = FL * FR
IF (TEST.LT.0) THEN
EPSU = EPSD
FU = F(EPSU)
IU = 0
IL = IL + 1
IF (IL.GE.2) FL = FL / 2
ELSE IF (TEST.GT.0) THEN
EPSL = EPSD
FL = F(EPSL)
IL = 0
IU = IU + 1
IF (IU.GE.2) FU = FU / 2
ELSE
ERA = 0
END IF
IF (ERA.LT.ERS .OR. ITER.GT.IMAX) EXIT
END DO
AMFP1 = EPSD
PRINT*, 'X=', AMFP1
RETURN
END
1 READ*, ERS, IMAX, EPSDPER, EPSDMAX
EPSL = EPSDPER; EPSU = EPSDMAX
CALL MFP(EPSL, EPSU, EPSD, ERS, ERA, ITER, IMAX)
GO TO 1
END
Many thanks
Hamid
Hi Hamid,
It looks like EPSD is not being updated as it's supposed to be (2 lines after the line numbered 33). If it's staying at zero (the default value) this means that your convergence condition is not going to be properly evaluated:
[cpp]IF(EPSD.NE.0) THEN
ERA = ABS((EPSD - EPSDOLD) / EPSD) * 100
END IF
[/cpp]
ERA is not evaluated and so stays at its default (uninitialized value) of 0 (which is less than ERS). Which would kick you out after the first pass with the current value of EPSD. But we know that the first evaluation is 3.18181...
It looks like a compiler bug. Can you run the debugger and track EPSD? Does this happen when you run the debugger (different compiler options - less optimization)?
Rich
Link Copied
- Mark as New
- Bookmark
- Subscribe
- Mute
- Subscribe to RSS Feed
- Permalink
- Report Inappropriate Content
I'm not sure what I should be seeing, as no matter how I try it, it prints X - 5. at the end. However, I do see a bug in your code. There's an IF-THEN-ELSE in the subroutine which sets variables IL and IU. If TEST is not zaero, then there will be a reference to an uninitialized variable IL or IU, depending on the sign of TEST, and this will affect the value of FL or FU. I'm not certain that this is relevant to the problem you're seeing, but it is suspicious.
- Mark as New
- Bookmark
- Subscribe
- Mute
- Subscribe to RSS Feed
- Permalink
- Report Inappropriate Content
I'm not sure what I should be seeing, as no matter how I try it, it prints X - 5. at the end. However, I do see a bug in your code. There's an IF-THEN-ELSE in the subroutine which sets variables IL and IU. If TEST is not zaero, then there will be a reference to an uninitialized variable IL or IU, depending on the sign of TEST, and this will affect the value of FL or FU. I'm not certain that this is relevant to the problem you're seeing, but it is suspicious.
Many thanks for this. But, there is still the same problem.
Could you check the following subroutine with a Visual Fortran Compiler. When it works with another compiler (like the second compiler I used), it seems that there is something wrong with the visual compiler. I do not know how the printing commands can affect the result.
1 READ*, ERS, IMAX, EPSDPER, EPSDMAX
EPSL = EPSDPER; EPSU = EPSDMAX
CALL MFP(EPSL, EPSU, EPSD, ERS, ERA, ITER, IMAX)
GO TO 1
END
FUNCTION F(X)
F = X**2-25
RETURN
END
SUBROUTINE MFP(EPSL, EPSU, EPSD, ERS, ERA, ITER, IMAX)
ITER = 0; IL = 0; IU = 0
FL = F(EPSL)
FU = F(EPSU)
C10 PRINT*, 'FL=', FL, 'FU=', FU
DO
EPSDOLD = EPSD
EPSD = EPSU - FU * (EPSL - EPSU) / (FL - FU)
C20 PRINT*, 'EPSL=', EPSL, 'EPSU=', EPSU
FR = F(EPSD)
ITER = ITER + 1
IF (EPSD.NE.0) ERA = ABS((EPSD - EPSDOLD) / EPSD) * 100
TEST = FL * FR
IF (TEST.LT.0) THEN
EPSU = EPSD
FU = F(EPSU)
IU = 0
IL = IL + 1
IF (IL.GE.2) FL = FL / 2
ELSE
IF (TEST.GT.0) THEN
EPSL = EPSD
FL = F(EPSL)
IL = 0
IU = IU + 1
IF (IU.GE.2) FU = FU / 2
ELSE
ERA = 0
END IF
END IF
IF (ERA.LT.ERS .OR. ITER.GT.IMAX) EXIT
END DO
AMFP1 = EPSD
PRINT*, 'X=', AMFP1
RETURN
END
- Mark as New
- Bookmark
- Subscribe
- Mute
- Subscribe to RSS Feed
- Permalink
- Report Inappropriate Content
What are the expected and bad results? I tried various combinations but did not see a problem. Please also tell me what compile options you use (if you're using Visual Studio, right click on the project, select Properties > Fortran > Command line and copy all of the All options text.)
- Mark as New
- Bookmark
- Subscribe
- Mute
- Subscribe to RSS Feed
- Permalink
- Report Inappropriate Content
What are the expected and bad results? I tried various combinations but did not see a problem. Please also tell me what compile options you use (if you're using Visual Studio, right click on the project, select Properties > Fortran > Command line and copy all of the All options text.)
Thank you for your help. By inputting 0.001, 100, 1 and 10 respectively , the correct answer must be X= 5 but when lines 10 and 20 are omitted, the answer is 3.181818.
I could not find the compile options yet and sadly I have to leave my office now.
Ill come back to find the compile options as soon as possible.
- Mark as New
- Bookmark
- Subscribe
- Mute
- Subscribe to RSS Feed
- Permalink
- Report Inappropriate Content
I can reproduce it now. I'll let you know what I find.
- Mark as New
- Bookmark
- Subscribe
- Mute
- Subscribe to RSS Feed
- Permalink
- Report Inappropriate Content
Thank you for your help. By inputting 0.001, 100, 1 and 10 respectively , the correct answer must be X= 5 but when lines 10 and 20 are omitted, the answer is 3.181818.
I could not find the compile options yet and sadly I have to leave my office know.
Ill come back to find the compile options as soon as possible.
Just a thought,
I've seen something similar before. Depending upon what compiler flags you have set, the optimizer may be playing some tricks with your code to keep your computations in a double length register for subsequent computations. The print statements clears the registers and forces the subsequent computations to be performed using stored values, hence losing the extra precision and changing results due to differing round-off errors.
I ran the code and didn't reproduce the results you cited (I didn't play with any of the compilation flags though). It worked fine both ways for me. The only thing that struck me was the use of implicit typing which means that your floating point numbers are defaulting to single precision which can make a big difference, particularly in algorithms like your Regula Falsi where your subtracting numbers that are getting close to one another as your interval converges (catastrophic cancellation) and then dividing by the result which tends to amplify errors. I'd start with declaring your variables and making sure that your floating point numbers are double precision.
Cheers
- Mark as New
- Bookmark
- Subscribe
- Mute
- Subscribe to RSS Feed
- Permalink
- Report Inappropriate Content
I can reproduce it now. I'll let you know what I find.
Hello again
Please note that the program was written with Visual Studio, but to compile it I used Visual Fortran Compiler applying ifort command. It means the only files I haveare mfp02.f, mfp02.obj and mfp02.exeand by right clicking cannot obtain the compile options. But, from Visual Studio (Project/Properties /Fortran/Command line) I could find the options:
/nologo /debug:full /Od /gen-interfaces /warn:interfaces /module:"Debug" /object:"Debug" /traceback /check:bounds /libs:static /threads /dbglibs /c
Regarding ender01's advice, I tried to define all the numbers as double procession with the following command:
implicit double precision (a-z)
However, the result was: X= NaN.
I would very much appreciate your consideration.
- Mark as New
- Bookmark
- Subscribe
- Mute
- Subscribe to RSS Feed
- Permalink
- Report Inappropriate Content
- Mark as New
- Bookmark
- Subscribe
- Mute
- Subscribe to RSS Feed
- Permalink
- Report Inappropriate Content
Thank you for checking. I got correct result when I useed another compiler, but the problem is with Visual Fortran Compiler 10.1.025 . I have also checked version 11.0.061 getting wrong answers after omitting lines 10 and 20.
- Mark as New
- Bookmark
- Subscribe
- Mute
- Subscribe to RSS Feed
- Permalink
- Report Inappropriate Content
Thank you for checking. I got correct result when I useed another compiler, but the problem is with Visual Fortran Compiler 10.1.025 . I have also checked version 11.0.061 getting wrong answers after omitting lines 10 and 20.
I wouldn't make everything double precision, just the floating point numbers. If you convert everything, including integers, you replace any integer math with real math, which, depending upon the equations, can be quite different (1./2 = 0.5 but 1/2 = 0). Also, implicit typing is a bad idea in general, pick the wrong variable names and you're doing integer math when you weren't planning on it (omega = sqrt(k/mass), for instance, will cause this with implicit typing. I would really start by making sure that all of your variables are appropriately typed, not just setting everything to double precision.
Cheers,
Rich
- Mark as New
- Bookmark
- Subscribe
- Mute
- Subscribe to RSS Feed
- Permalink
- Report Inappropriate Content
I wouldn't make everything double precision, just the floating point numbers. If you convert everything, including integers, you replace any integer math with real math, which, depending upon the equations, can be quite different (1./2 = 0.5 but 1/2 = 0). Also, implicit typing is a bad idea in general, pick the wrong variable names and you're doing integer math when you weren't planning on it (omega = sqrt(k/mass), for instance, will cause this with implicit typing. I would really start by making sure that all of your variables are appropriately typed, not just setting everything to double precision.
Cheers,
Rich
Dear Rich
Much appreciated your help.
I have just double checked the floating point numbers and everything seems to be correct. But there is still the same problem (X= NaN) even if only the floating point numers are made double precision.
Could you compile the following program and check it by inputting 0.001, 100, 1 and 10 respectively to get the correct answer X= 5.
Many thanks
Hamid
FUNCTION F(X)
F = X**2-25
RETURN
END
SUBROUTINE MFP(EPSL, EPSU, EPSD, ERS, ERA, ITER, IMAX)
ITER = 0; IL = 0; IU = 0
FL = F(EPSL)
FU = F(EPSU)
10 PRINT*, 'FL=', FL, 'FU=', FU
DO
EPSDOLD = EPSD
EPSD = EPSU - FU * (EPSL - EPSU) / (FL - FU)
20 PRINT*, 'EPSL=', EPSL, 'EPSU=', EPSU
FR = F(EPSD)
ITER = ITER + 1
IF (EPSD.NE.0) ERA = ABS((EPSD - EPSDOLD) / EPSD) * 100
PRINT*, 'ERA=', ERA
TEST = FL * FR
IF (TEST.LT.0) THEN
EPSU = EPSD
FU = F(EPSU)
IU = 0
IL = IL + 1
IF (IL.GE.2) FL = FL / 2
ELSE IF (TEST.GT.0) THEN
EPSL = EPSD
FL = F(EPSL)
IL = 0
IU = IU + 1
IF (IU.GE.2) FU = FU / 2
ELSE
ERA = 0
END IF
print*, 'ERA=', ERA
IF (ERA.LT.ERS .OR. ITER.GT.IMAX) EXIT
END DO
AMFP1 = EPSD
PRINT*, 'X=', AMFP1
RETURN
END
DOUBLE PRECISION EPSDPER, EPSDMAX, EPSL, EPSU, EPSD, EPSDOLD
DOUBLE PRECISION X, ERS, ERA, FL, FR, FU, AMFP1
! implicit double precision (a-z)
1 READ*, ERS, IMAX, EPSDPER, EPSDMAX
EPSL = EPSDPER; EPSU = EPSDMAX
CALL MFP(EPSL, EPSU, EPSD, ERS, ERA, ITER, IMAX)
GO TO 1
END
- Mark as New
- Bookmark
- Subscribe
- Mute
- Subscribe to RSS Feed
- Permalink
- Report Inappropriate Content
Hi there
I think I found a bug in Fortran Visual Compiler 10.1.025. It is about a subroutine used to solve nonlinear equations with MFP approach (see below). The brief code is to apply equation Y = X^2 - 25 just for checking purpose.
The problem is that when I add two printing commands on lines 10 and 20, the code works perfectly, but when I omit them it does not! In fact, it can be checked by inputting 0.001, 100, 1 and 10 respectively to get X= 5.
Could someone let me know if it is a bug or I made a mistake? I compiled the subroutine with another compiler and it worked even after omitting the two printing commands.
Many thanks
Hamid
1 READ*, ERS, IMAX, EPSDPER, EPSDMAX
EPSL = EPSDPER; EPSU = EPSDMAX
CALL MFP(EPSL, EPSU, EPSD, ERS, ERA, ITER, IMAX)
GO TO 1
END
FUNCTION F(X)
F = X**2-25
RETURN
END
SUBROUTINE MFP(EPSL, EPSU, EPSD, ERS, ERA, ITER, IMAX)
ITER = 0
FL = F(EPSL)
FU = F(EPSU)
10 PRINT*, 'FL=', FL, 'FU=', FU
DO
EPSDOLD = EPSD
EPSD = EPSU - FU * (EPSL - EPSU) / (FL - FU)
20 PRINT*, 'EPSL=', EPSL, 'EPSU=', EPSU
FR = F(EPSD)
ITER = ITER + 1
C print*, 'iter=', iter
C100 PRINT*, 'EPSD=', EPSD
IF (EPSD.NE.0) ERA = ABS((EPSD - EPSDOLD) / EPSD) * 100
PRINT*, 'ERA=', ERA
TEST = FL * FR
IF (TEST.LT.0) THEN
EPSU = EPSD
FU = F(EPSU)
IU = 0
IL = IL + 1
IF (IL.GE.2) FL = FL / 2
ELSE IF (TEST.GT.0) THEN
EPSL = EPSD
FL = F(EPSL)
IL = 0
IU = IU + 1
IF (IU.GE.2) FU = FU / 2
ELSE
ERA = 0
END IF
print*, 'ERA=', ERA
IF (ERA.LT.ERS .OR. ITER.GT.IMAX) EXIT
END DO
AMFP1 = EPSD
PRINT*, 'X=', AMFP1
RETURN
END
Hamid
On call to MFP the argument EPSD is uninitialized and used withing MFPD.
Do you expect uninitialized variables to be 0?
Although some compilers initialize uninitialized variables to 0 there is no Fortran requirement to do so.
Try setting EPSD = 0.0 or whatever you wish for initialization.
Jim Dempsey
- Mark as New
- Bookmark
- Subscribe
- Mute
- Subscribe to RSS Feed
- Permalink
- Report Inappropriate Content
Dear Rich
Much appreciated your help.
I have just double checked the floating point numbers and everything seems to be correct. But there is still the same problem (X= NaN) even if only the floating point numers are made double precision.
Could you compile the following program and check it by inputting 0.001, 100, 1 and 10 respectively to get the correct answer X= 5.
Many thanks
Hamid
FUNCTION F(X)
F = X**2-25
RETURN
END
SUBROUTINE MFP(EPSL, EPSU, EPSD, ERS, ERA, ITER, IMAX)
ITER = 0; IL = 0; IU = 0
FL = F(EPSL)
FU = F(EPSU)
10 PRINT*, 'FL=', FL, 'FU=', FU
DO
EPSDOLD = EPSD
EPSD = EPSU - FU * (EPSL - EPSU) / (FL - FU)
20 PRINT*, 'EPSL=', EPSL, 'EPSU=', EPSU
FR = F(EPSD)
ITER = ITER + 1
IF (EPSD.NE.0) ERA = ABS((EPSD - EPSDOLD) / EPSD) * 100
PRINT*, 'ERA=', ERA
TEST = FL * FR
IF (TEST.LT.0) THEN
EPSU = EPSD
FU = F(EPSU)
IU = 0
IL = IL + 1
IF (IL.GE.2) FL = FL / 2
ELSE IF (TEST.GT.0) THEN
EPSL = EPSD
FL = F(EPSL)
IL = 0
IU = IU + 1
IF (IU.GE.2) FU = FU / 2
ELSE
ERA = 0
END IF
print*, 'ERA=', ERA
IF (ERA.LT.ERS .OR. ITER.GT.IMAX) EXIT
END DO
AMFP1 = EPSD
PRINT*, 'X=', AMFP1
RETURN
END
DOUBLE PRECISION EPSDPER, EPSDMAX, EPSL, EPSU, EPSD, EPSDOLD
DOUBLE PRECISION X, ERS, ERA, FL, FR, FU, AMFP1
! implicit double precision (a-z)
1 READ*, ERS, IMAX, EPSDPER, EPSDMAX
EPSL = EPSDPER; EPSU = EPSDMAX
CALL MFP(EPSL, EPSU, EPSD, ERS, ERA, ITER, IMAX)
GO TO 1
END
Hi Hamid,
It won't compile as written, I get type mismatch errors on the call to MFP (error #6633: The type of the actual argument differs from the type of the dummy argument. [EPSU], [EPSD], [ERS], and [ERA]). Did it compile this way for you? If I add in:
DOUBLE PRECISION EPSU, EPSD, ERS, ERA
to the MFP subroutine and
DOUBLE PRECISION X
to the function F(X), then I start getting type mismatch errors on calls such as FL = F(EPSL). When I remove all of the type declaration, I get the correct answer of 5.000000. Also note, where you explicitly typed the variable is outside of the scope of the subroutine and function, all of those local variables are still implicitly typed.
What you have to remember is that modern optimizing compilers are pretty clever at fixing some programming errors. If you are using an older compiler, it may just be that your compiler isn't catching it and the newer ones are, hence other compilers aren't seeing the problem.
BTW, this code looks pretty old, where did you get it? You might be better off rewriting it in more modern Fortran. I'll bang on it for a little while and see if I can get you something that will work.
Cheers,
Rich
- Mark as New
- Bookmark
- Subscribe
- Mute
- Subscribe to RSS Feed
- Permalink
- Report Inappropriate Content
Dear Rich
Much appreciated your help.
I have just double checked the floating point numbers and everything seems to be correct. But there is still the same problem (X= NaN) even if only the floating point numers are made double precision.
Could you compile the following program and check it by inputting 0.001, 100, 1 and 10 respectively to get the correct answer X= 5.
Many thanks
Hamid
FUNCTION F(X)
F = X**2-25
RETURN
END
SUBROUTINE MFP(EPSL, EPSU, EPSD, ERS, ERA, ITER, IMAX)
ITER = 0; IL = 0; IU = 0
FL = F(EPSL)
FU = F(EPSU)
10 PRINT*, 'FL=', FL, 'FU=', FU
DO
EPSDOLD = EPSD
EPSD = EPSU - FU * (EPSL - EPSU) / (FL - FU)
20 PRINT*, 'EPSL=', EPSL, 'EPSU=', EPSU
FR = F(EPSD)
ITER = ITER + 1
IF (EPSD.NE.0) ERA = ABS((EPSD - EPSDOLD) / EPSD) * 100
PRINT*, 'ERA=', ERA
TEST = FL * FR
IF (TEST.LT.0) THEN
EPSU = EPSD
FU = F(EPSU)
IU = 0
IL = IL + 1
IF (IL.GE.2) FL = FL / 2
ELSE IF (TEST.GT.0) THEN
EPSL = EPSD
FL = F(EPSL)
IL = 0
IU = IU + 1
IF (IU.GE.2) FU = FU / 2
ELSE
ERA = 0
END IF
print*, 'ERA=', ERA
IF (ERA.LT.ERS .OR. ITER.GT.IMAX) EXIT
END DO
AMFP1 = EPSD
PRINT*, 'X=', AMFP1
RETURN
END
DOUBLE PRECISION EPSDPER, EPSDMAX, EPSL, EPSU, EPSD, EPSDOLD
DOUBLE PRECISION X, ERS, ERA, FL, FR, FU, AMFP1
! implicit double precision (a-z)
1 READ*, ERS, IMAX, EPSDPER, EPSDMAX
EPSL = EPSDPER; EPSU = EPSDMAX
CALL MFP(EPSL, EPSU, EPSD, ERS, ERA, ITER, IMAX)
GO TO 1
END
Hi Rich
I made mistakes in making double precisions that is the reason for the NoN result. Sorry about that. But in the corrected program (see below) has the first problem. X=3.18... instead of X=5.
I will check the uninitialized variables as Jim said. But I am almost sure the compiler I am using, makes them zero. I will check it tomorrow and tell you the results.
Many thanks
Hamid
DOUBLE PRECISION ERS,ERA,EPSDPER,EPSDMAX,EPSL,EPSU,EPSD,EPSDOLD
DOUBLE PRECISION F,X,FL,FU,FR,TEST,AMFP1
1 READ*, ERS, IMAX, EPSDPER, EPSDMAX
EPSL = EPSDPER; EPSU = EPSDMAX
CALL MFP(EPSL, EPSU, EPSD, ERS, ERA, ITER, IMAX)
GO TO 1
END
DOUBLE PRECISION FUNCTION F(X)
DOUBLE PRECISION X
F = X**2-25
RETURN
END
SUBROUTINE MFP(EPSL,EPSU,EPSD,ERS,ERA,ITER,IMAX)
DOUBLE PRECISION EPSL,EPSU,EPSD,ERS,ERA,FL,FU,FR,EPSDOLD
DOUBLE PRECISION TEST,AMFP1
ITER = 0; IL = 0; IU = 0
FL = F(EPSL)
FU = F(EPSU)
!10 PRINT*, 'FL=', FL, 'FU=', FU
DO
EPSDOLD = EPSD
EPSD = EPSU - FU * (EPSL - EPSU) / (FL - FU)
!20 PRINT*, 'EPSL=', EPSL, 'EPSU=', EPSU
FR = F(EPSD)
ITER = ITER + 1
IF (EPSD.NE.0) ERA = ABS((EPSD - EPSDOLD) / EPSD) * 100
PRINT*, 'ERA=', ERA
TEST = FL * FR
IF (TEST.LT.0) THEN
EPSU = EPSD
FU = F(EPSU)
IU = 0
IL = IL + 1
IF (IL.GE.2) FL = FL / 2
ELSE IF (TEST.GT.0) THEN
EPSL = EPSD
FL = F(EPSL)
IL = 0
IU = IU + 1
IF (IU.GE.2) FU = FU / 2
ELSE
ERA = 0
END IF
print*, 'ERA=', ERA
IF (ERA.LT.ERS .OR. ITER.GT.IMAX) EXIT
END DO
AMFP1 = EPSD
PRINT*, 'X=', AMFP1
RETURN
END
- Mark as New
- Bookmark
- Subscribe
- Mute
- Subscribe to RSS Feed
- Permalink
- Report Inappropriate Content
Hi Hamid,
It won't compile as written, I get type mismatch errors on the call to MFP (error #6633: The type of the actual argument differs from the type of the dummy argument. [EPSU], [EPSD], [ERS], and [ERA]). Did it compile this way for you? If I add in:
DOUBLE PRECISION EPSU, EPSD, ERS, ERA
to the MFP subroutine and
DOUBLE PRECISION X
to the function F(X), then I start getting type mismatch errors on calls such as FL = F(EPSL). When I remove all of the type declaration, I get the correct answer of 5.000000. Also note, where you explicitly typed the variable is outside of the scope of the subroutine and function, all of those local variables are still implicitly typed.
What you have to remember is that modern optimizing compilers are pretty clever at fixing some programming errors. If you are using an older compiler, it may just be that your compiler isn't catching it and the newer ones are, hence other compilers aren't seeing the problem.
BTW, this code looks pretty old, where did you get it? You might be better off rewriting it in more modern Fortran. I'll bang on it for a little while and see if I can get you something that will work.
Cheers,
Rich
Hi Hamid,
Here's a reworking of the code:
[cpp]! HamidCode.f90
!
! FUNCTIONS:
! HamidCode - Entry point of console application.
!
!****************************************************************************
!
! PROGRAM: HamidCode
!
! PURPOSE: Entry point for the console application.
!
!****************************************************************************
program HamidCode
implicit none
! Variables
! Precision constants
integer, parameter:: HI=selected_real_kind(12)
integer, parameter:: LO=selected_real_kind(5)
integer, parameter:: I_HI=selected_int_kind(9)
integer, parameter:: I_LO=selected_int_kind(4)
integer(kind = I_LO) :: ans
real(kind=HI) :: EPSDPER, EPSDMAX, EPSL, EPSU, EPSD, ERS, ERA
integer(kind=I_LO) :: IMAX, ITER
! Read in data
do
write(*,*) "Please enter the convergence tolerance, maximum number of &
&iterations, and the lower and upper bounds of the initial &
&interval, respectively"
read(*,*) ERS, IMAX, EPSDPER, EPSDMAX
EPSL = EPSDPER; EPSU = EPSDMAX
CALL MFP(EPSL, EPSU, EPSD, ERS, ERA, ITER, IMAX)
write(*,*) "Again? Enter 0 to continue"
read(*,*) ans
if (ans /= 0) exit
end do
contains
function F(X)
implicit none
real(kind=HI), intent(in) :: X
real(kind=HI) :: F
F = X**2-25._HI
return
end function F
subroutine MFP(EPSL, EPSU, EPSD, ERS, ERA, ITER, IMAX)
implicit none
real(kind=HI), intent(inout):: EPSL, EPSU
real(kind=HI), intent(out) :: EPSD, ERS, ERA
integer(kind=I_LO), intent(in) :: IMAX
integer(kind=I_LO), intent(out) :: ITER
integer(kind=I_LO) :: IL, IU
real(kind=HI) :: FL, FU, FR, TEST
real(kind=HI) :: EPSDOLD, AMFP1
! Initialize iteration variables
ITER = 0 ! Iteration count
IL = 0 ! Interval (bracket) lower bound
IU = 0 ! Interval (bracket) upper bound
! Evaluate F(X) at bracket bounds
FL = F(EPSL) ! F(X_lower)
FU = F(EPSU) ! F(X_upper)
! Echo function evaluation on interval bounds
10 write(*,*) 'FL=', FL, 'FU=', FU
! Iterate using Regula Falsi
do
! Store old value of zero crossing
! NOTE: This variable is not initialized on call to MFP
EPSDOLD = EPSD
! Compute updated estimate of zero crossing
EPSD = EPSU - FU * (EPSL - EPSU) / (FL - FU)
! Echo current bounds of interval
20 write(*,*) 'EPSL=', EPSL, 'EPSU=', EPSU
! Evaluate function at updated zero crossing estimate
FR = F(EPSD)
! Increment interation counter
ITER = ITER + 1
! Check for convergence
! Compute relative change in zero crossing estimate
if(EPSD.NE.0) then
! Avoid divide by zero
! NOTE: EPSDOLD not initialized on first pass
! Shouldn't matter in general but there are
! pathological cases
ERA = ABS((EPSD - EPSDOLD) / EPSD) * 100
end if
! Echo relative change in zero estimate
write(*,*) 'ERA=', ERA
! Use sign test to determine which end of interval to
! discard (want to bracket zero crossing, need one bound
! to evaluate to a positive number and the other bound
! to a negative number => if the function has a zero
! can usually expect a sign change).
TEST = FL * FR
if (TEST.LT.0) then
! If product is negative, new value has oposite sign
! of lower bound, update upper bound
EPSU = EPSD
! Re-evaluate function at new upper bound
! NOTE: This is superfluous, FR still holds this value
FU = F(EPSU)
! Oh, we have a modified Regula-Falsi of sorts...
IU = 0 ! Upper bound replaced, zero persistence counter
IL = IL + 1 ! Lower bound persisting, increment counter
! Most implementations do this every step rather than
! every other. I think that this is intended to stabilize
! the method but as long as you are bracketing a zero
! crossing (continuous or not) the method is stable (unlike
! Newton-Raphson methods where its usually a good idea to
! stabilize if you don't know the function)
! If lower bound is persisting, skew next estimate closer to
! lower bound by using half the function evaluation at the
! lower bound to compute next estimate
if (IL.GE.2) FL = FL / 2
else if (TEST.GT.0) then
! If product is positive, new value has same sign as lower
! lower bound, update lower bound
EPSL = EPSD
! Re-evaluate function at new lower bound
! NOTE: This is superfluous, FR still holds this value
FL = F(EPSL)
! Oh, we have a modified Regula-Falsi of sorts...
IL = 0 ! Lower bound replaced, zero persistence counter
IU = IU + 1 ! Upper bound persisting, increment counter
! If upper bound is persisting, skew next estimate closer to
! upper bound by using half the function evaluation at the
! upper bound to compute next estimate
if (IU.GE.2) FU = FU / 2
else
! We appear to have landed right on the zero!
ERA = 0._HI ! Percent change should be zero from here on out!
end if
! Echo change in zero crossing estimate
print*, 'ERA=', ERA
! Check convergence
! NOTE: In this method you should also check the size of the bracket
! (very steep functions will cause excessive iterations) as
! well as magnitude of the functions that are flat over part
! of the bracketed domain can cause convergence issues)
if (ERA.LT.ERS .OR. ITER.GT.IMAX) exit
end do
! Store final estimate of zero crossing
AMFP1 = EPSD
! Echo final estimate of zero crossing
write(*,*) 'X=', AMFP1
RETURN
end subroutine MFP
end program HamidCode[/cpp]
I've annotated it as well. There is, in addition, one variable which is used in a computation that is not initialized. Not all compilers handle uninitialized variables in the same way. I've called it out in the annotation. Let me know if that works. Also, you can see the correct typing for your variables should you want to go back and type them in the original code. This one works, with or without the print statements. If you still can't get it to work, send a step by step report of the values you have for each step of the code for an iteration or two and I'll match it up with my version and see if we can't narrow down where the misbehavior occurs.
Cheers,
Rich
- Mark as New
- Bookmark
- Subscribe
- Mute
- Subscribe to RSS Feed
- Permalink
- Report Inappropriate Content
EPSL= 1.000000 EPSU= 10.00000
EPSL= 3.181818 EPSU= 10.00000
EPSL= 4.310345 EPSU= 10.00000
EPSL= 4.310345 EPSU= 5.142132
EPSL= 4.989630 EPSU= 5.142132
EPSL= 4.999855 EPSU= 5.142132
EPSL= 4.999855 EPSU= 5.000141
X= 5.000000
Press any key to continue
So I can see where you are getting your answer 3.181818 from - the loop is exiting at the start of the second pass.
I believe that is where you should expend your debugging effort. It should be pretty straightforward.
- Mark as New
- Bookmark
- Subscribe
- Mute
- Subscribe to RSS Feed
- Permalink
- Report Inappropriate Content
Hi Hamid,
Here's a reworking of the code:
! HamidCode.f90 ! ! FUNCTIONS: ! HamidCode - Entry point of console application. ! !**************************************************************************** ! ! PROGRAM: HamidCode ! ! PURPOSE: Entry point for the console application. ! !**************************************************************************** program HamidCode implicit none ! Variables ! Precision constants integer, parameter:: HI=selected_real_kind(12) integer, parameter:: LO=selected_real_kind(5) integer, parameter:: I_HI=selected_int_kind(9) integer, parameter:: I_LO=selected_int_kind(4) integer(kind = I_LO) :: ans real(kind=HI) :: EPSDPER, EPSDMAX, EPSL, EPSU, EPSD, ERS, ERA integer(kind=I_LO) :: IMAX, ITER ! Read in data do write(*,*) "Please enter the convergence tolerance, maximum number of & &iterations, and the lower and upper bounds of the initial & &interval, respectively" read(*,*) ERS, IMAX, EPSDPER, EPSDMAX EPSL = EPSDPER; EPSU = EPSDMAX CALL MFP(EPSL, EPSU, EPSD, ERS, ERA, ITER, IMAX) write(*,*) "Again? Enter 0 to continue" read(*,*) ans if (ans /= 0) exit end do contains function F(X) implicit none real(kind=HI), intent(in) :: X real(kind=HI) :: F F = X**2-25._HI return end function F subroutine MFP(EPSL, EPSU, EPSD, ERS, ERA, ITER, IMAX) implicit none real(kind=HI), intent(inout):: EPSL, EPSU real(kind=HI), intent(out) :: EPSD, ERS, ERA integer(kind=I_LO), intent(in) :: IMAX integer(kind=I_LO), intent(out) :: ITER integer(kind=I_LO) :: IL, IU real(kind=HI) :: FL, FU, FR, TEST real(kind=HI) :: EPSDOLD, AMFP1 ! Initialize iteration variables ITER = 0 ! Iteration count IL = 0 ! Interval (bracket) lower bound IU = 0 ! Interval (bracket) upper bound ! Evaluate F(X) at bracket bounds FL = F(EPSL) ! F(X_lower) FU = F(EPSU) ! F(X_upper) ! Echo function evaluation on interval bounds 10 write(*,*) 'FL=', FL, 'FU=', FU ! Iterate using Regula Falsi do ! Store old value of zero crossing ! NOTE: This variable is not initialized on call to MFP EPSDOLD = EPSD ! Compute updated estimate of zero crossing EPSD = EPSU - FU * (EPSL - EPSU) / (FL - FU) ! Echo current bounds of interval 20 write(*,*) 'EPSL=', EPSL, 'EPSU=', EPSU ! Evaluate function at updated zero crossing estimate FR = F(EPSD) ! Increment interation counter ITER = ITER + 1 ! Check for convergence ! Compute relative change in zero crossing estimate if(EPSD.NE.0) then ! Avoid divide by zero ! NOTE: EPSDOLD not initialized on first pass ! Shouldn't matter in general but there are ! pathological cases ERA = ABS((EPSD - EPSDOLD) / EPSD) * 100 end if ! Echo relative change in zero estimate write(*,*) 'ERA=', ERA ! Use sign test to determine which end of interval to ! discard (want to bracket zero crossing, need one bound ! to evaluate to a positive number and the other bound ! to a negative number => if the function has a zero ! can usually expect a sign change). TEST = FL * FR if (TEST.LT.0) then ! If product is negative, new value has oposite sign ! of lower bound, update upper bound EPSU = EPSD ! Re-evaluate function at new upper bound ! NOTE: This is superfluous, FR still holds this value FU = F(EPSU) ! Oh, we have a modified Regula-Falsi of sorts... IU = 0 ! Upper bound replaced, zero persistence counter IL = IL + 1 ! Lower bound persisting, increment counter ! Most implementations do this every step rather than ! every other. I think that this is intended to stabilize ! the method but as long as you are bracketing a zero ! crossing (continuous or not) the method is stable (unlike ! Newton-Raphson methods where its usually a good idea to ! stabilize if you don't know the function) ! If lower bound is persisting, skew next estimate closer to ! lower bound by using half the function evaluation at the ! lower bound to compute next estimate if (IL.GE.2) FL = FL / 2 else if (TEST.GT.0) then ! If product is positive, new value has same sign as lower ! lower bound, update lower bound EPSL = EPSD ! Re-evaluate function at new lower bound ! NOTE: This is superfluous, FR still holds this value FL = F(EPSL) ! Oh, we have a modified Regula-Falsi of sorts... IL = 0 ! Lower bound replaced, zero persistence counter IU = IU + 1 ! Upper bound persisting, increment counter ! If upper bound is persisting, skew next estimate closer to ! upper bound by using half the function evaluation at the ! upper bound to compute next estimate if (IU.GE.2) FU = FU / 2 else ! We appear to have landed right on the zero! ERA = 0._HI ! Percent change should be zero from here on out! end if ! Echo change in zero crossing estimate print*, 'ERA=', ERA ! Check convergence ! NOTE: In this method you should also check the size of the bracket ! (very steep functions will cause excessive iterations) as ! well as magnitude of the functions that are flat over part ! of the bracketed domain can cause convergence issues) if (ERA.LT.ERS .OR. ITER.GT.IMAX) exit end do ! Store final estimate of zero crossing AMFP1 = EPSD ! Echo final estimate of zero crossing write(*,*) 'X=', AMFP1 RETURN end subroutine MFP end program HamidCode
I've annotated it as well. There is, in addition, one variable which is used in a computation that is not initialized. Not all compilers handle uninitialized variables in the same way. I've called it out in the annotation. Let me know if that works. Also, you can see the correct typing for your variables should you want to go back and type them in the original code. This one works, with or without the print statements. If you still can't get it to work, send a step by step report of the values you have for each step of the code for an iteration or two and I'll match it up with my version and see if we can't narrow down where the misbehavior occurs.
Cheers,
Rich
Hi Rich
Appreciated very much for your kind consideration of the problem.
The posted program gives the same wrong answer (X= 3.181818) after commenting the printing statements. However, as you annotated EPSDOLD is not initialized (Jim also mentioned the same about EPSD - thank you Jim). So, when I initialize EPSD or EPSDOLD to zero, the program works perfectly even if the floating point numbers are not double precision.
Interestingly, when I do not initialize them to zero and print them, they are zero! But sometimes the program does not work correctly depending on the printing statements. For example, in the following code, when lines 31, 32 and 33 are stated, the result is right:
.001
100
1
10
epsd 0.0000000E+00
epsd 0.0000000E+00
epsd 0.0000000E+00
epsd 3.181818
epsd 4.310345
epsd 5.142132
epsd 4.989630
epsd 4.999855
epsd 5.000141
X= 5.000000
but, when just one of the printing lines is commented, the result is wrong:
.001
100
1
10
epsd 0.0000000E+00
epsd 0.0000000E+00
X= 3.181818
If you could forward it to the developers, it may be interesting for them. The compiler I am using is: Intel Visual Fortran Compiler Professional for applications running on IA-32, Version 11.0 Build 20080930 Package ID: w_cprof_p_11.0.061.
FUNCTION F(X)
F = X**2-25
RETURN
END
SUBROUTINE MFP(EPSL, EPSU, EPSD, ERS, ERA, ITER, IMAX)
ITER = 0; IL = 0; IU = 0; !EPSD = 0
31 print*, 'epsd',epsd
FL = F(EPSL)
FU = F(EPSU)
32 print*, 'epsd',epsd
DO
33 print*, 'epsd',epsd
EPSDOLD = EPSD
EPSD = EPSU - FU * (EPSL - EPSU) / (FL - FU)
FR = F(EPSD)
ITER = ITER + 1
IF (EPSD.NE.0) ERA = ABS((EPSD - EPSDOLD) / EPSD) * 100
TEST = FL * FR
IF (TEST.LT.0) THEN
EPSU = EPSD
FU = F(EPSU)
IU = 0
IL = IL + 1
IF (IL.GE.2) FL = FL / 2
ELSE IF (TEST.GT.0) THEN
EPSL = EPSD
FL = F(EPSL)
IL = 0
IU = IU + 1
IF (IU.GE.2) FU = FU / 2
ELSE
ERA = 0
END IF
IF (ERA.LT.ERS .OR. ITER.GT.IMAX) EXIT
END DO
AMFP1 = EPSD
PRINT*, 'X=', AMFP1
RETURN
END
1 READ*, ERS, IMAX, EPSDPER, EPSDMAX
EPSL = EPSDPER; EPSU = EPSDMAX
CALL MFP(EPSL, EPSU, EPSD, ERS, ERA, ITER, IMAX)
GO TO 1
END
Many thanks
Hamid
- Mark as New
- Bookmark
- Subscribe
- Mute
- Subscribe to RSS Feed
- Permalink
- Report Inappropriate Content
Hi Rich
Appreciated very much for your kind consideration of the problem.
The posted program gives the same wrong answer (X= 3.181818) after commenting the printing statements. However, as you annotated EPSDOLD is not initialized (Jim also mentioned the same about EPSD - thank you Jim). So, when I initialize EPSD or EPSDOLD to zero, the program works perfectly even if the floating point numbers are not double precision.
Interestingly, when I do not initialize them to zero and print them, they are zero! But sometimes the program does not work correctly depending on the printing statements. For example, in the following code, when lines 31, 32 and 33 are stated, the result is right:
.001
100
1
10
epsd 0.0000000E+00
epsd 0.0000000E+00
epsd 0.0000000E+00
epsd 3.181818
epsd 4.310345
epsd 5.142132
epsd 4.989630
epsd 4.999855
epsd 5.000141
X= 5.000000
but, when just one of the printing lines is commented, the result is wrong:
.001
100
1
10
epsd 0.0000000E+00
epsd 0.0000000E+00
X= 3.181818
If you could forward it to the developers, it may be interesting for them. The compiler I am using is: Intel Visual Fortran Compiler Professional for applications running on IA-32, Version 11.0 Build 20080930 Package ID: w_cprof_p_11.0.061.
FUNCTION F(X)
F = X**2-25
RETURN
END
SUBROUTINE MFP(EPSL, EPSU, EPSD, ERS, ERA, ITER, IMAX)
ITER = 0; IL = 0; IU = 0; !EPSD = 0
31 print*, 'epsd',epsd
FL = F(EPSL)
FU = F(EPSU)
32 print*, 'epsd',epsd
DO
33 print*, 'epsd',epsd
EPSDOLD = EPSD
EPSD = EPSU - FU * (EPSL - EPSU) / (FL - FU)
FR = F(EPSD)
ITER = ITER + 1
IF (EPSD.NE.0) ERA = ABS((EPSD - EPSDOLD) / EPSD) * 100
TEST = FL * FR
IF (TEST.LT.0) THEN
EPSU = EPSD
FU = F(EPSU)
IU = 0
IL = IL + 1
IF (IL.GE.2) FL = FL / 2
ELSE IF (TEST.GT.0) THEN
EPSL = EPSD
FL = F(EPSL)
IL = 0
IU = IU + 1
IF (IU.GE.2) FU = FU / 2
ELSE
ERA = 0
END IF
IF (ERA.LT.ERS .OR. ITER.GT.IMAX) EXIT
END DO
AMFP1 = EPSD
PRINT*, 'X=', AMFP1
RETURN
END
1 READ*, ERS, IMAX, EPSDPER, EPSDMAX
EPSL = EPSDPER; EPSU = EPSDMAX
CALL MFP(EPSL, EPSU, EPSD, ERS, ERA, ITER, IMAX)
GO TO 1
END
Many thanks
Hamid
Hi Hamid,
It looks like EPSD is not being updated as it's supposed to be (2 lines after the line numbered 33). If it's staying at zero (the default value) this means that your convergence condition is not going to be properly evaluated:
[cpp]IF(EPSD.NE.0) THEN
ERA = ABS((EPSD - EPSDOLD) / EPSD) * 100
END IF
[/cpp]
ERA is not evaluated and so stays at its default (uninitialized value) of 0 (which is less than ERS). Which would kick you out after the first pass with the current value of EPSD. But we know that the first evaluation is 3.18181...
It looks like a compiler bug. Can you run the debugger and track EPSD? Does this happen when you run the debugger (different compiler options - less optimization)?
Rich
- Mark as New
- Bookmark
- Subscribe
- Mute
- Subscribe to RSS Feed
- Permalink
- Report Inappropriate Content
Hi Hamid,
It looks like EPSD is not being updated as it's supposed to be (2 lines after the line numbered 33). If it's staying at zero (the default value) this means that your convergence condition is not going to be properly evaluated:
IF(EPSD.NE.0) THEN ERA = ABS((EPSD - EPSDOLD) / EPSD) * 100 END IF
ERA is not evaluated and so stays at its default (uninitialized value) of 0 (which is less than ERS). Which would kick you out after the first pass with the current value of EPSD. But we know that the first evaluation is 3.18181...
It looks like a compiler bug. Can you run the debugger and track EPSD? Does this happen when you run the debugger (different compiler options - less optimization)?
Rich
The comparison between EPSD and 0 ismixed mode - what happens if you use 0. or 0D0? does EPSD get updated in this case?
David
- Mark as New
- Bookmark
- Subscribe
- Mute
- Subscribe to RSS Feed
- Permalink
- Report Inappropriate Content
The comparison between EPSD and 0 ismixed mode - what happens if you use 0. or 0D0? does EPSD get updated in this case?
David
Hi David
Thank you for this. Using 0. Or 0D0 did not make any difference.
Hamid
- Subscribe to RSS Feed
- Mark Topic as New
- Mark Topic as Read
- Float this Topic for Current User
- Bookmark
- Subscribe
- Printer Friendly Page