Difficulty with Relatively Simple IF Structure
monir
I've an integrable smooth funcG(ups), and I need to replace it with a
stepped function with 15 equal Gstep values, and find the
corresponding XR (related to ups by the ACOS function below).
Gstep of segment i between upsn and upsp is given by:
.....abs {integral [upsn to upsp] funcG} / (upsprev-upsnod) - (i-1)
*Gstep
or..abs(res/(upsprev-upsnod))-(i-1)*Gstep
XR of interest is 1.0 to 0.7043... where G has its abs max value.
Starting from XR = 1.0 and move towards 0.704....
The following (modified) simplified code (F 77) works fine for other
scenarios, but not for this one. The difficulty here appears to be in
the nested IF structure.
My idea is to systematically approach the correct XRnod for the
current segment
by retaining the latest good value of XRnod (before the result of the
test is reversed) and:
- move inward {if the current step > (i-1)*Gstep} or
- move outward {if the current step < (i-1)*Gstep}
The assignment of xrtemp is either incomplete or incorrect!!
C
****************************************************************
SUBROUTINE SHTSBD
C
****************************************************************
C declarations .....
C
pi = 3.141592654E0
XH = 0.20
XHM = 0.5*(1.-XH)
XHP = 0.5*(1.+XH)
tol = 1.E-7
XRPREV = 1.0
XRGMAX = 0.70432264
UPSPREV = pi
GPREV = 0.0
NFTIP = 15
GMAX = -0.033839315
Gstep = abs(GMAX/NFTIP)
DO 10 I=2, NFTIP + 1
XRNOD = XRGMAX
UPSNOD = ACOS(XHP/XHM-XRNOD/XHM)
11 CALL QROMO (FUNCG, upsnod, upsprev,res,MIDPNT)
if (abs(abs(res/(upsprev-upsnod))-(i-1)*Gstep) .gt. tol) then
if (abs(res/(upsprev-upsnod)) .gt. (i-1)*Gstep) then !
secant inward towards 1.0
xrtemp = xrnod
xrnod =(xrnod+xrprev)/2.
UPSNOD = ACOS(XHP/XHM-XRNOD/XHM)
goto 11
else !secant
outward
xrnod =(xrnod+xrtemp)/2.
UPSNOD = ACOS(XHP/XHM-XRNOD/XHM)
goto 11
end if
end if
GNOD = SIGN(real(res/(upsprev-upsnod)),1.)
GPREV = GNOD
XRPREV = XRNOD
10 CONTINUE
RETURN
END
Your expert help would be greatly appreciated.
Kind regards.
Monir
e p chandler
There is nothing wrong with a nested IF, IF you understand what it is
doing. Without being familiar with your problem and the subroutines
you are trying to use it is difficult to be of much help. Of course I
can not run your program fragment, so I have to guess... I suspect
that your program does not work correctly for other functions, even
though you think that it does.
0. You have not described or shown EXACTLY how your program fails.
What is the faulty output?
1. Have you used all of the debugging and checking options available
from your compiler?
1a. I highly recommend adding IMPLICIT NONE to all of your code AND
then declaring the types of ALL variables.
2. GNOD = SIGN(real(res/(upsprev-upsnod)),1.)
This statement DOES NOT do what you think it does. As you use it, it
returns ABS(real(res/(upsprev-upsnod)). SIGN(x,y) returns ABS(x) with
the sign of Y affixed. The sign of 1. is always +.
RES, UPSPREV and UPSNOD are all IMPLICITLY typed (A-H,O-Z) as REAL,
hence the REAL() is not needed. see 1a too.
3. I don't know exactly how your two-way branch works, but it should
be symmetric in some sense.
4. (i-1)*gstep is a CONSTANT for each value of I. You might see things
more clearly if you set a variable equal to its value, say BETA.
5. abs(res/(upsprev-upsnod)) is a COMMON SUB EXPRESSION in your nested
IFs. You might see more clearly if you set a variable equal to its
value, say ALPHA.
You may want code something like this (it's not completely Fortran):
DO 10 I=2, NFTIP + 1
XRNOD = XRGMAX
UPSNOD = ACOS(XHP/XHM-XRNOD/XHM)
ALPHA = (I-1)*Gstep !##
11 CALL QROMO (FUNCG, upsnod, upsprev,res,MIDPNT)
BETA = ABS(RES/(UPSPREV-UPSNOD)) !##
DELTA = ALPHA - BETA !##
! if (abs(abs(res/(upsprev-upsnod))-(i-1)*Gstep) .gt. tol) then
! if (abs(res/(upsprev-upsnod)) .gt. (i-1)*Gstep) then !
secant inward towards 1.0
if ( abs(DELTA) .gt. tol) then
if ( DELTA .GT. 0. ) then ! secant inward towards 1.0
xrtemp = xrnod
xrnod =(xrnod+xrprev)/2.
UPSNOD = ACOS(XHP/XHM-XRNOD/XHM)
goto 11
else IF (DELTA .LT. 0) THEN !secant outward
xrnod =(xrnod+xrtemp)/2.
UPSNOD = ACOS(XHP/XHM-XRNOD/XHM)
goto 11
else
WHAT DO YOU WANT HERE IF DELTA .EQ. 0?
end if
end if
!#### GNOD = SIGN(real(res/(upsprev-upsnod)),1.)
WHAT YOU REALLY WANT HERE
GPREV = GNOD
XRPREV = XRNOD
10 CONTINUE
See this diagram for conditions:
BETA-TOL BETA BETA+TOL
ALPHA: <-----| |----->
ALPHA: <-------------|------------>
HTH
- e
monir
>
> - Show quoted text -
Hi e p chandler;
Thank you for your comments, and I'll try to address your concerns as
follows:
1) I agree that there's nothing wrong with a nested IF, but I believe
its structure in the provided modified subroutine is questionable.
The problem as described in my OP is simplified further below so you
don't have to guess.
2) The program works perfectly for other scenarios where that
particular IF block was not invoked in that subroutine.
3) The program fails to produce the "expected" results of subdividing
a given smooth function into segments of equal step values over a
desired range of the function independent variable. Will be addressed
later with some numerical results.
4) I've used print statements in the subroutine for debugging
purposes, and not the full debugging and checking options available in
g95 compiler
5) I've all the correct declaration in the FULL version of the
subroutine. The posted (tested) version of the sub is a simplified
almost self-contained (apart from calling the integrator QROMO)
version. Had I posted the FULL version of the sub with the
declarations, data blocks, equivalence statement, etc., it would've
confused and distracted the reader.
6) Let us leave the SIGN() statement for now as is. Will look at it
later once the IF two-way branch problem is resolved.
7) Your coding suggestions were more or less already there in the FULL
version of the sub. Nevertheless, I've used your suggested additional
variable names in the latest abbreviated version below.
8) Here's the latest version of the subroutine:
C
****************************************************************
SUBROUTINE SHTSBD2
C
****************************************************************
C declarations .....res is DB
C
pi = 3.141592654E0
XH = 0.20
XHM = 0.5*(1.-XH)
XHP = 0.5*(1.+XH)
tol = 1.E-7
XRPREV = 1.0
XRGMAX = 0.70432264
UPSPREV = pi
GPREV = 0.0
NFTIP = 15
GMAX = -0.033839315
Gstep = abs(GMAX/NFTIP)
DO 10 I=2, NFTIP + 1
alpha = (i-1)*Gstep
xrtemp1=xrgmax
xrtemp2=xrprev
XRNOD = XRGMAX
XRNOD = XRGMAX
UPSNOD = ACOS(XHP/XHM-XRNOD/XHM)
! res is the result of integrating FUNCG from upsnod to upsprev
11 CALL QROMO (FUNCG, upsnod, upsprev,res,MIDPNT)
beta = ABS(res/(upsprev-upsnod))
delta = beta - alpha
if (abs(delta) .gt. tol) then
if(delta .gt. 0) then ! secant inward towards 1.0
xrtemp1 = xrnod
xrnod =(xrnod+xrtemp2)/2.
UPSNOD = ACOS(XHP/XHM-XRNOD/XHM)
goto 11
else !secant
outward
xrtemp2 = xrnod
xrnod =(xrnod+xrtemp1)/2.
UPSNOD = ACOS(XHP/XHM-XRNOD/XHM)
goto 11
end if
end if
GNOD = SIGN(real(res/(upsprev-upsnod)),1.)
GPREV = GNOD
XRPREV = XRNOD
10 CONTINUE
RETURN
END
Will post shortly the numerical results for a simple example where
calling QROMO is optional and the function integration can be
performed instead by its analytical close form. This way, the results
could be easily reproduced.
Regards.
Monir
monir
>
> - Show quoted text -
Hi;
Continuation to my earlier post:
8) NUMERICAL EXAMPLE.
Let us assume I've a smooth function represented by the Fourier sine
series:
FuncG = sum(m=1 to 9)GM(m)*SIN(m*ups)
where:
GM( 1) = -.329564E-01
GM( 2) = 0.457676E-02
GM( 3) = -.265454E-03
GM( 4) = 0.101691E-03
GM( 5) = 0.523956E-04
GM( 6) = 0.112262E-04
GM( 7) = -.596746E-06
GM( 8) = 0.863187E-05
GM( 9) = 0.398432E-05
The integral of FuncG from ups2 to ups1 is given by:
res = - sum(m=1 to 9) GM(m)*[cos(m*ups1)-cos(m*ups2)]/m
So one doesn't need QROMO for the integration.
9) The |max| of FuncG is 0.03383932
and it's at XR = 0.70432264
The relation between variable "ups" and "XR" is given by:
UPSNOD = ACOS(XHP/XHM-XR/XHM)
where:
XHM = 0.400
XHP = 0.600
10) I need to replace the smooth FuncG by a stepped function with 15
equal Gstep values, and find the
corresponding XR over the range XR=1.0 (ups=pi) to XR = 0.70432264
First step at XR = 1.0
From Item 9) above:
...Gstep = (0.03383932/15) = 0.002255954
11) Subroutine SHTSBD2 (F 77, posted earlier) and a different Excel
VBA procedure produce almost identical (but unexpected and strange)
results:
(XRNOD vs stepped '-G' below are ready for plotting as a stepped
function )
....I...XRNOD.....sgmt no..sgmt span...stepped '-G'.....Gstep
1.0
0 0.00225595
2 1.0 1 0.0021915 0.002256
2 0.9978085 1 0.002256
0.00225590
3 0.9978085 2 1.51E-05 0.0045119
3 0.9977934 2 0.0045119
0.00225598
4 0.9977934 3 0.006689 0.0067678
4 0.9911044 3 0.0067678
0.00225599
5 0.9911044 4 0.0001274 0.0090238
5 0.990977 4 0.0090238
0.00225597
6 0.990977 5 0.0115604 0.0112798
6 0.9794166 5 0.0112798
0.00225585
7 0.9794166 6 0.0004701 0.0135356
7 0.9789465 6 0.0135356
0.00225602
8 0.9789465 7 0.0171849 0.0157917
8 0.9617616 7 0.0157917
0.00225603
9 0.9617616 8 0.0012885 0.0180477
9 0.9604731 8 0.0180477
0.00225597
10 0.9604731 9 0.02424864 0.0203037
10 0.93622446 9 0.0203037
0.00225583
11 0.93622446 10 0.0031352 0.0225595
11 0.93308926 10 0.0225595 0.00225597
12 0.93308926 11 0.03440786 0.0248155
12 0.8986814 11 0.0248155
0.00225605
13 0.8986814 12 0.007779 0.0270715
13 0.8909024 12 0.0270715
0.00225598
14 0.8909024 13 0.05386484 0.0293275
14 0.83703756 13 0.0293275 0.00225584
15 0.83703756 14 0.02669046 0.0315833
15 0.8103471 14 0.0315833
0.00225600
16 not converging
should be:
16 0.8103471 15 0.03383932
16 0.70432264 15 0.03383932 0
12) Some observations:
a) one would expect the segments spans to increase from XR=1.0 moving
inward as the function flattens out towards the max G at XR = 0.704.
By examining the above numerical results, you'll find the rate of
increase is different for the odd and even no. segments. Why ??
b) The odd number segments 1, 3, 5, ... appear to have the correct
spans, while even no. segments 2, 4, 6, ... have very small spans.
For example, segment no. 2 has a span of 1.51E-05 surrounded by
segment no. 1 of span 0.0021915 and segment no. 3 of span 0.006689.
The function is a smooth and gradually increasing from 0.0 at XR=1.0
to max at the inner bound of XR range-of-interest. What is the cause
of the instability ?? Is the logic in the Subroutine incorrect ?? or
is it rather a mathematical problem ??
c) all the segments have the correct desired Gstep (within tol)
despite the non-homogeneous spans. Strange !!
d) Sub SHTSBD2 would not converge for the last value of Index I=NFTIP
+1 in DO 10.
I can't change the Index of the DO loop, since the program package
works perfectly for other scenarios where that particular IF block is
not invoked.
13) I would like to suggest, if it's not too much trouble, for someone
independently to try the simple problem outlined in Items 8, 9, 10
above (possibly using one of those system Algebraic packages) to
compare the results of XRNOD with those tabulated in Item 11) above.
Unfortunately, I don't have access to any of such packages.
Your help would be greatly appreciated.
Regards.
Monir
e p chandler
"monir" <mon...@mondenet.com> wrote in message
GNOD = SIGN(real(res/(upsprev-upsnod)),1.)
RETURN
END
--------------- end of quoted text ----------------
Well you've corrected the problem with xrtemp being un-initialized. Also you
have a solution that looks (anti) symmetric between the two cases. I don't
know if your algorithm is correct, it's not my field.
Really you only need a simple 3 way branch
1. if beta > alpha + tol
2. if beta < alpha - tol
3. otherwise.
-- elliot
e p chandler
[big snip]
Hi;
Regards.
Monir
------------- end of message -----------------
Now you've got a problem in numerical analysis. Your G function changes sign
(oscillates). Perhaps the secant method (whatever that is) is unstable under
those conditions.
For example, it's easy to cook up a root finder using Newton's method. With
a bad quess it might not converge. For some functions, it might not
converge, etc.
Perhaps you need the help of specialists in numerical analysis, even though
some of them do hang out here.
-- elliot
George
What the helliot, elliot.
Having immersed myself in fortran, in my stupid pain day, I'd like to see
if Neumann fails, an instance of which I have not reckoned.
The symbols seem willfully bad like "ups". What's more, you'd need Tony
Romo to dig out of this syntax deficit.
--
George
America is a friend to the people of Iraq. Our demands are directed only at
the regime that enslaves them and threatens us. When these demands are met,
the first and greatest benefit will come to Iraqi men, women and children.
George W. Bush
Picture of the Day http://apod.nasa.gov/apod/
monir
>
> - Show quoted text -
Hi Elliot;
Thank you for taking the time. The sample G function (Item 8 of my
post) DOES NOT change sign within the range-of-interest.
By the way, the table (Item 11) appears to have lost its formatting in
the process of posting. My apologies if it was unreadable!!
Will post the question in Google sci.math.num-analysis DG.
Regards.
Monir
robin
On Nov 18, 9:46 pm, e p chandler <e...@juno.com> wrote:
>4) I've used print statements in the subroutine for debugging
>purposes, and not the full debugging and checking options available in
>g95 compiler
Why not, after e. p. chandler suggested that you do?
>5) I've all the correct declaration in the FULL version of the
>subroutine.
How do you know?
Do you have IMPLICIT NONE in this subroutine?
It's not shown here.
monir
Thank you all for your interest, and I apologise for the delay in
responding.
1) As I indicated earlier, the program has been extensively tested and
works perfectly for other scenarios where that particular nested IF
block is not invoked in Subroutine SHTSBD2.
2) The numerical method works perfectly as well for all other tested
smooth, monotonic, integrable and well-defined 1D functions over the
XR-range-of-interest.
3) The difficulty appears to be "conceptual" rather than programming
or numerical. The current problem situation is somewhat analogous
(albeit a weak one!) to specifying X1 and X2 believed to bracket a
root of a function, and it turned out to be a false bracket. In such
situations, one can't blame the algorithm!!
4) I'm currently testing an alternative concept using Excel VBA. The
results so far are very promising.
5) Once the extensive numerical experimentation is done, will return
to the Fortran "modified" nested IF structure in Sub SHTSBD2 where
your help would be greatly appreciated. Will focus only on the IF
block, variables names, initial values, etc.
Regards.
Monir
robin
~On Nov 22, 7:11 pm, "robin" <robi...@bigpond.com> wrote:
~> "monir" <mon...@mondenet.com> wrote in message
~>
~> news:a1c21aaa-c40d-4737...@v4g2000yqa.googlegroups.com...
~> On Nov 18, 9:46 pm, e p chandler <e...@juno.com> wrote:
~> >4) I've used print statements in the subroutine for debugging
~> >purposes, and not the full debugging and checking options available in
~> >g95 compiler
~> Why not, after e. p. chandler suggested that you do?
~> >5) I've all the correct declaration in the FULL version of the
~> >subroutine.
~> How do you know?
~> Do you have IMPLICIT NONE in this subroutine?
~> It's not shown here.
~Thank you all for your interest, and I apologise for the delay in
~responding.
~1) As I indicated earlier, the program has been extensively tested and
~works perfectly for other scenarios where that particular nested IF
~block is not invoked in Subroutine SHTSBD2.
Again, you've not provided the evidence requested.
And just because your program has been "extensively tested"
doesn't mean that it has no bugs.
By your own admission, the program "works" EXCEPT
for that particular nested IF!
monir
1) Robin: "Again, you've not provided the evidence requested. And
just because your program has been "extensively tested" doesn't mean
that it has no bugs."
I believe I've all the correct declaration in the FULL version of the
subroutine. The posted (tested) version of Sub SHTSBD2b is a
simplified almost self-contained version (apart from calling the
integrator QROMO). Posting the FULL version of the Sub (together
with the other 50 or 60 subroutines of the program) with all the
declarations, data blocks, equivalence statements, 1000s of lines of
code, etc., would be confusing to the reader and impractical.
2) Robin: "By your own admission, the program "works" EXCEPT for that
particular nested IF!"
The program works perfectly for all other scenarios when that
particular recently-added nested IF block is not required or not
invoked in Subroutine SHTSBD2b.
3) I've successfully tested an alternative concept using Excel VBA.
So it's the time now to return to the Fortran nested IF structure in
Sub SHTSBD2b where your help would be greatly appreciated.
Let us focus only on the IF block, variables names, initial values,
etc.
4) Here's the logic behind the IF structure:
Let us say you have a straight line extending between 1.0 and 0.704...
and you need to divide it to NFTIP segments starting from 1.0
Each segment has to meet certain condition specified in the outer IF.
You start (i=2) with the bounds XRprev = 1.0 and XRgmax = 0.704...,
and use the bisection method in succession:
1st attempt: XRnod = (1.0 + 0.704...)/2 = 0.85...
If the outer IF condition is met, then:
End outer IF
and assign XRprev = 0.85... and go to the next i
If not and delta > tol, first branch of the inner IF, then move
towards 1.0 and try:
...2nd attempt: XRnod = (1.0+0.85...)/2
...go to 11 and test again
If not and delta < tol, second branch of the inner IF, then move
towards 0.704... and try:
...2nd attempt: XRnod = (0.704...+0.85...)/2
...go to 11 and test again
End inner IF
Repeat until successful for the current i=2
If successful with, say, XRnod 0.97..., then
End outer IF
and assign XRprev = 0.97... and go for the next i loop, with the XR
range available now becomes:
... the latest successful XRnod = 0.97... and XRgmax = 0.704...
And so on.
The index of DO 10 MUST remain unchanged ( 2 to NFTIP+1).
C ****************************************************************
SUBROUTINE SHTSBD2b
C ****************************************************************
C declarations .....res is DB
C
pi = 3.141592654E0
XH = 0.20
XHM = 0.5*(1.-XH)
XHP = 0.5*(1.+XH)
tol = 1.E-7
XRPREV = 1.0
XRGMAX = 0.70432264
UPSPREV = pi
NFTIP = 15
GMAX = -0.033839315
Gstep = abs(0.995*GMAX/(NFTIP-1))
DO 10 I=2, NFTIP + 1
if (i .eq. nftip+1) then
alpha = (nftip-1)*Gstep
else
alpha = (i-1.5)*Gstep
end if
xrtemp1=xrgmax
xrtemp2=xrprev
XRNOD = XRGMAX
UPSNOD = ACOS(XHP/XHM-XRNOD/XHM)
! res is the result of integrating FUNCG from upsnod to upsprev
11 CALL QROMO (FUNCG, upsnod, upsprev,res,MIDPNT)
beta = ABS(res/(upsprev-upsnod))
delta = beta - alpha
if (abs(delta) .gt. tol) then
if(delta .gt. 0) then ! move outward towards 1.0
xrtemp1 = xrnod
xrnod =(xrnod+xrtemp2)/2.
UPSNOD = ACOS(XHP/XHM-XRNOD/XHM)
goto 11
else !move inward towards 0.704...
xrtemp2 = xrnod
xrnod =(xrnod+xrtemp1)/2.
UPSNOD = ACOS(XHP/XHM-XRNOD/XHM)
goto 11
end if
end if
XRPREV = XRNOD
10 CONTINUE
RETURN
END
Would appreciate your comments and any suggestions you might have.
Regards.
e p chandler
>
>
> ! res is the result of integrating FUNCG from upsnod to upsprev
> 11 CALL QROMO (FUNCG, upsnod, upsprev,res,MIDPNT)
> beta = ABS(res/(upsprev-upsnod))
> delta = beta - alpha
>
> if (abs(delta) .gt. tol) then
> if(delta .gt. 0) then ! move outward towards 1.0
> xrtemp1 = xrnod
> xrnod =(xrnod+xrtemp2)/2.
> UPSNOD = ACOS(XHP/XHM-XRNOD/XHM)
> goto 11
> else !move inward towards 0.704...
> xrtemp2 = xrnod
> xrnod =(xrnod+xrtemp1)/2.
> UPSNOD = ACOS(XHP/XHM-XRNOD/XHM)
> goto 11
> end if
> end if
>
> XRPREV = XRNOD
### to here ###
> 10 CONTINUE
>
> RETURN
> END
>
> Would appreciate your comments and any suggestions you might have.
>
> Regards.
OK. Without looking at your specification, does restructuring the code as
follows make any difference?
I've left the code mostly un-indented to make it easier for me to post.
The loop is restructured to eliminate the GO TO. It exits when abs(delta) <=
TOL.
### from here ###
do
CALL QROMO (FUNCG, upsnod, upsprev,res,MIDPNT)
beta = ABS(res/(upsprev-upsnod))
delta = beta - alpha
if ( abs(delta) .LE. tol) exit ! within tol, quit the DO loop
if ( delta .gt. 0) then ! move outward towards 1.0
xrtemp1 = xrnod
xrnod =(xrnod+xrtemp2)/2.
UPSNOD = ACOS(XHP/XHM-XRNOD/XHM)
else ! move inward towards 0.704...
xrtemp2 = xrnod
xrnod =(xrnod+xrtemp1)/2.
UPSNOD = ACOS(XHP/XHM-XRNOD/XHM)
end if
end do
XRPREV = XRNOD
### to here ###
Please let us know if this change also fixes your problem. If not, then it's
time to see if the code I just posted meets your specification posted above.
Sorry I don't have time to read your spec in detail now.
Also I suggest you print out the variables of interest from within this
routine for debugging purposes.
If the code does not work with these changes, then I suspect you have not
coded it properly OR you have a numerical problem.
Again, please add print statements to dump out the intermediate variables,
run that program and post the output here.
Sorry, but your spec makes sense to you but it makes my head spin.
--- e
monir
Thank you for the attempt.
Your suggested modifications are in error since they represent two
serious misinterpretations of what the nested IF supposed to do:
A) If the condition is TRUE in the outermost IF, then exit the nested
IF block, and continue for the next "i" of the DO loop.
NOT "quit the DO loop".
B) In either branch of the inner IF, a new bisection attempt MUST be
followed by calling the integrator QROMO() to obtain the corresponding
new value of "res" before re-testing the condition.
You need GO TO 11 once you enter either branch of the inner IF.
You would need to review the simple description of the logic behind
the intended nested IF (Item 4. of my previous post, Dec 03, 2008,
11:35 am).
Regards.
Monir
e p chandler
I am confident that the EXIT only leaves the inner most of the nested
DO loops.
I do think my code does the same thing as your old code.
Perhaps someone else with a few more functioning brain cells can
figure it out. I guess being in management and having programmed for
many years have taken their toll.
- e
- e
robin
> e p chandler;
>
> Thank you for the attempt.
>
> Your suggested modifications are in error since they represent two
> serious misinterpretations of what the nested IF supposed to do:
> A) If the condition is TRUE in the outermost IF, then exit the nested
> IF block, and continue for the next "i" of the DO loop.
> NOT "quit the DO loop".
It exits the inner DO loop, not the outer one, and the code
appears to be equivalent to yours.
> B) In either branch of the inner IF, a new bisection attempt MUST be
> followed by calling the integrator QROMO() to obtain the corresponding
> new value of "res" before re-testing the condition.
> You need GO TO 11 once you enter either branch of the inner IF.
His code does NOT require GO TO 11,
because that is covered by the DO loop that he introduced.
> You would need to review the simple description of the logic behind
> the intended nested IF (Item 4. of my previous post, Dec 03, 2008,
> 11:35 am).
You need to actually DO what has been suggested.
robin
> Hi;
>
> 1) Robin: "Again, you've not provided the evidence requested. And
> just because your program has been "extensively tested" doesn't mean
> that it has no bugs."
> I believe I've all the correct declaration in the FULL version of the
> subroutine.
"Believing" is no subsitiute for a definitive check,
such as is obtained by the use of IMPLICIT NONE,
enabling all checks (subscript, uninitialized variables, etc).
By you own admission:-
"I've ... not [used] the full debugging and checking options
" available in g95 compiler".
Why do you expect help when you have NOT done
this basic first step?
> The posted (tested) version of Sub SHTSBD2b is a
> simplified almost self-contained version (apart from calling the
> integrator QROMO). Posting the FULL version of the Sub
If you expect any assistance, you need to do what has
been suggested.
> (together
> with the other 50 or 60 subroutines of the program)
All the other subrountines are not necessary.
> with all the
> declarations, data blocks, equivalence statements, 1000s of lines of
> code, etc., would be confusing to the reader and impractical.
The complete subroutine wouldn't go astray.
> 2) Robin: "By your own admission, the program "works" EXCEPT for that
> particular nested IF!"
> The program works perfectly for all other scenarios when that
> particular recently-added nested IF block is not required or not
> invoked in Subroutine SHTSBD2b.
Just as I said.
But that doesn't prove that it's the IF block is at fault.
What's plain is that you have not adequately tested your
code, by printing out values bvefore entering the
IF block, and during the block, in order to investigate
what is going on.
> 3) I've successfully tested an alternative concept using Excel VBA.
> So it's the time now to return to the Fortran nested IF structure in
> Sub SHTSBD2b where your help would be greatly appreciated.
> Let us focus only on the IF block, variables names, initial values,
> etc.
You haven't provided the declarations, have you! Why is that?
Why is pi given to 10 digits, whereas only 6 of these digits
will be used?
> XH = 0.20
> XHM = 0.5*(1.-XH)
> XHP = 0.5*(1.+XH)
> tol = 1.E-7
>
> XRPREV = 1.0
> XRGMAX = 0.70432264
>
>
> UPSPREV = pi
> NFTIP = 15
> GMAX = -0.033839315
> Gstep = abs(0.995*GMAX/(NFTIP-1))
>
>
> DO 10 I=2, NFTIP + 1
> if (i .eq. nftip+1) then
> alpha = (nftip-1)*Gstep
> else
> alpha = (i-1.5)*Gstep
> end if
> xrtemp1=xrgmax
> xrtemp2=xrprev
> XRNOD = XRGMAX
> UPSNOD = ACOS(XHP/XHM-XRNOD/XHM)
>
>
>
> ! res is the result of integrating FUNCG from upsnod to upsprev
> 11 CALL QROMO (FUNCG, upsnod, upsprev,res,MIDPNT)
> beta = ABS(res/(upsprev-upsnod))
> delta = beta - alpha
What are the values of beta and delta here?
xrnod? xrtemp2? xrtemp1? tol? (repeat, TOL ! )