Difficulty with Relatively Simple IF Structure
Author Message
Difficulty with Relatively Simple IF Structure

Hello;

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

Sun, 08 May 2011 02:49:48 GMT
Difficulty with Relatively Simple IF Structure

Quote:
> Hello;

> 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

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

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

Sun, 08 May 2011 10:46:55 GMT
Difficulty with Relatively Simple IF Structure

Quote:

> > Hello;

> > 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

> 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

> 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- Hide quoted text -

> - Show quoted text -

Hi e p chandler;

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

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
variable names in the latest abbreviated version below.

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

Mon, 09 May 2011 10:11:14 GMT
Difficulty with Relatively Simple IF Structure

Quote:

> > > Hello;

> > > 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

> > 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

> > 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- Hide quoted text -

> > - Show quoted text -

> Hi e p chandler;

> 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
> variable names in the latest abbreviated version below.

> 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- Hide quoted text -

> - 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
...

Mon, 09 May 2011 11:41:58 GMT
Difficulty with Relatively Simple IF Structure

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
--------------- 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

Mon, 09 May 2011 13:43:22 GMT
Difficulty with Relatively Simple IF Structure

[big snip]

Quote:
> 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.

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.

Your help would be greatly appreciated.

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

Mon, 09 May 2011 13:50:17 GMT
Difficulty with Relatively Simple IF Structure

Quote:

> [big snip]
>> 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.

> 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.

> Your help would be greatly appreciated.

> 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

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/

Mon, 09 May 2011 17:15:24 GMT
Difficulty with Relatively Simple IF Structure

Quote:

> > [big snip]
> >> 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.

> > 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

> > ------------- 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

> 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 Dayhttp://apod.nasa.gov/apod/- Hide quoted text -

> - 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

Mon, 09 May 2011 21:12:51 GMT
Difficulty with Relatively Simple IF Structure

Quote:
>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?

Quote:
>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.

Thu, 12 May 2011 08:11:38 GMT
Difficulty with Relatively Simple IF Structure

Quote:

> >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.

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

Tue, 17 May 2011 23:05:15 GMT
Difficulty with Relatively Simple IF Structure

~>

~> >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.

for that particular nested IF!

Sun, 22 May 2011 20:07:17 GMT
Difficulty with Relatively Simple IF Structure
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.  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

Regards.

Mon, 23 May 2011 00:35:28 GMT
Difficulty with Relatively Simple IF Structure

Quote:
> 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.  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)

### from here ###

> ! 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 ###

Quote:
> 10 CONTINUE

> RETURN
> END

> 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.

run that program and post the output here.

Sorry, but your spec makes sense to you but it makes my head spin.

--- e

Mon, 23 May 2011 01:45:34 GMT
Difficulty with Relatively Simple IF Structure
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".

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

Mon, 23 May 2011 03:32:21 GMT
Difficulty with Relatively Simple IF Structure

Quote:
> 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".

> 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

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

Mon, 23 May 2011 10:11:47 GMT

 Page 1 of 2 [ 17 post ] Go to page: [1] [2]

Relevant Pages