Difficulty with Relatively Simple IF Structure
Author 
Message 
moni #1 / 17

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} / (upsprevupsnod)  (i1) *Gstep or..abs(res/(upsprevupsnod))(i1)*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 > (i1)*Gstep} or  move outward {if the current step < (i1)*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.E7 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/XHMXRNOD/XHM) 11 CALL QROMO (FUNCG, upsnod, upsprev,res,MIDPNT) if (abs(abs(res/(upsprevupsnod))(i1)*Gstep) .gt. tol) then if (abs(res/(upsprevupsnod)) .gt. (i1)*Gstep) then ! secant inward towards 1.0 xrtemp = xrnod xrnod =(xrnod+xrprev)/2. UPSNOD = ACOS(XHP/XHMXRNOD/XHM) goto 11 else !secant outward xrnod =(xrnod+xrtemp)/2. UPSNOD = ACOS(XHP/XHMXRNOD/XHM) goto 11 end if end if GNOD = SIGN(real(res/(upsprevupsnod)),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 


e p chandle #2 / 17

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} / (upsprevupsnod)  (i1) > *Gstep > or..abs(res/(upsprevupsnod))(i1)*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 > (i1)*Gstep} or >  move outward {if the current step < (i1)*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.E7 > 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/XHMXRNOD/XHM) > 11 CALL QROMO (FUNCG, upsnod, upsprev,res,MIDPNT) > ? ? ? ? if (abs(abs(res/(upsprevupsnod))(i1)*Gstep) .gt. tol) then > ? ? ? ? ? ? ? if (abs(res/(upsprevupsnod)) .gt. (i1)*Gstep) then ?! > secant inward towards 1.0 > ? ? ? ? ? ? ? ? ?xrtemp = xrnod > ? ? ? ? ? ? ? ? ?xrnod =(xrnod+xrprev)/2. > ? ? ? ? ? ? ? ? ?UPSNOD = ACOS(XHP/XHMXRNOD/XHM) > ? ? ? ? ? ? ? ? ?goto 11 > ? ? ? ? ? ? ?else ? !secant > outward > ? ? ? ? ? ? ? ? ?xrnod =(xrnod+xrtemp)/2. > ? ? ? ? ? ? ? ? ?UPSNOD = ACOS(XHP/XHMXRNOD/XHM) > ? ? ? ? ? ? ? ? ?goto 11 > ? ? ? ? ? ? ?end if > ? ? ? ? end if > ? ? ? ? ? ?GNOD = SIGN(real(res/(upsprevupsnod)),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 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/(upsprevupsnod)),1.) This statement DOES NOT do what you think it does. As you use it, it returns ABS(real(res/(upsprevupsnod)). 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 (AH,OZ) as REAL, hence the REAL() is not needed. see 1a too. 3. I don't know exactly how your twoway branch works, but it should be symmetric in some sense. 4. (i1)*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/(upsprevupsnod)) 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/XHMXRNOD/XHM) ALPHA = (I1)*Gstep !## 11 CALL QROMO (FUNCG, upsnod, upsprev,res,MIDPNT) BETA = ABS(RES/(UPSPREVUPSNOD)) !## DELTA = ALPHA  BETA !## ! if (abs(abs(res/(upsprevupsnod))(i1)*Gstep) .gt. tol) then ! if (abs(res/(upsprevupsnod)) .gt. (i1)*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/XHMXRNOD/XHM) goto 11 else IF (DELTA .LT. 0) THEN !secant outward xrnod =(xrnod+xrtemp)/2. UPSNOD = ACOS(XHP/XHMXRNOD/XHM) goto 11 else WHAT DO YOU WANT HERE IF DELTA .EQ. 0? end if end if !#### GNOD = SIGN(real(res/(upsprevupsnod)),1.) WHAT YOU REALLY WANT HERE GPREV = GNOD XRPREV = XRNOD 10 CONTINUE See this diagram for conditions: BETATOL BETA BETA+TOL ALPHA: < > ALPHA: <> HTH  e

Sun, 08 May 2011 10:46:55 GMT 


moni #3 / 17

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} / (upsprevupsnod)  (i1) > > *Gstep > > or..abs(res/(upsprevupsnod))(i1)*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 > (i1)*Gstep} or > >  move outward {if the current step < (i1)*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.E7 > > 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/XHMXRNOD/XHM) > > 11 CALL QROMO (FUNCG, upsnod, upsprev,res,MIDPNT) > > ? ? ? ? if (abs(abs(res/(upsprevupsnod))(i1)*Gstep) .gt. tol) then > > ? ? ? ? ? ? ? if (abs(res/(upsprevupsnod)) .gt. (i1)*Gstep) then ?! > > secant inward towards 1.0 > > ? ? ? ? ? ? ? ? ?xrtemp = xrnod > > ? ? ? ? ? ? ? ? ?xrnod =(xrnod+xrprev)/2. > > ? ? ? ? ? ? ? ? ?UPSNOD = ACOS(XHP/XHMXRNOD/XHM) > > ? ? ? ? ? ? ? ? ?goto 11 > > ? ? ? ? ? ? ?else ? !secant > > outward > > ? ? ? ? ? ? ? ? ?xrnod =(xrnod+xrtemp)/2. > > ? ? ? ? ? ? ? ? ?UPSNOD = ACOS(XHP/XHMXRNOD/XHM) > > ? ? ? ? ? ? ? ? ?goto 11 > > ? ? ? ? ? ? ?end if > > ? ? ? ? end if > > ? ? ? ? ? ?GNOD = SIGN(real(res/(upsprevupsnod)),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 > 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/(upsprevupsnod)),1.) > This statement DOES NOT do what you think it does. As you use it, it > returns ABS(real(res/(upsprevupsnod)). 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 (AH,OZ) as REAL, > hence the REAL() is not needed. see 1a too. > 3. I don't know exactly how your twoway branch works, but it should > be symmetric in some sense. > 4. (i1)*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/(upsprevupsnod)) 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/XHMXRNOD/XHM) > ? ? ?ALPHA = (I1)*Gstep !## > 11 CALL QROMO (FUNCG, upsnod, upsprev,res,MIDPNT) > ? ? ?BETA = ABS(RES/(UPSPREVUPSNOD)) !## > ? ? ?DELTA = ALPHA  BETA !## > ! ? ? ? if (abs(abs(res/(upsprevupsnod))(i1)*Gstep) .gt. tol) then > ! ? ? ? ? ? ? if (abs(res/(upsprevupsnod)) .gt. (i1)*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/XHMXRNOD/XHM) > ? ? ? ? ? ? ? ? ?goto 11 > ? ? ? ? ? ? ?else IF (DELTA .LT. 0) THEN !secant outward > ? ? ? ? ? ? ? ? ?xrnod =(xrnod+xrtemp)/2. > ? ? ? ? ? ? ? ? ?UPSNOD = ACOS(XHP/XHMXRNOD/XHM) > ? ? ? ? ? ? ? ? ?goto 11 > ? ? ? ? ? ? ?else > ? ? ? ? ? ? ? ? ?WHAT DO YOU WANT HERE IF DELTA .EQ. 0? > ? ? ? ? ? ? ?end if > ? ? ? ? end if > !#### ? ? ? GNOD = SIGN(real(res/(upsprevupsnod)),1.) > ? ? ? ? WHAT YOU REALLY WANT HERE > ? ? GPREV = GNOD > ? ? XRPREV = XRNOD > 10 CONTINUE > See this diagram for conditions: > ? ? ? ? ? ? ? ?BETATOL ? ? BETA ? ?BETA+TOL > ALPHA: ? ? ? < ? ? ? ? ? ? ? ? ? ?> > ALPHA: ? ? ? ? <> > HTH >  e Hide quoted text  >  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 selfcontained (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 twoway 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.E7 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 = (i1)*Gstep xrtemp1=xrgmax xrtemp2=xrprev XRNOD = XRGMAX XRNOD = XRGMAX UPSNOD = ACOS(XHP/XHMXRNOD/XHM) ! res is the result of integrating FUNCG from upsnod to upsprev 11 CALL QROMO (FUNCG, upsnod, upsprev,res,MIDPNT) beta = ABS(res/(upsprevupsnod)) 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/XHMXRNOD/XHM) goto 11 else !secant outward xrtemp2 = xrnod xrnod =(xrnod+xrtemp1)/2. UPSNOD = ACOS(XHP/XHMXRNOD/XHM) goto 11 end if end if GNOD = SIGN(real(res/(upsprevupsnod)),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 


moni #4 / 17

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} / (upsprevupsnod)  (i1) > > > *Gstep > > > or..abs(res/(upsprevupsnod))(i1)*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 > (i1)*Gstep} or > > >  move outward {if the current step < (i1)*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.E7 > > > 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/XHMXRNOD/XHM) > > > 11 CALL QROMO (FUNCG, upsnod, upsprev,res,MIDPNT) > > > ? ? ? ? if (abs(abs(res/(upsprevupsnod))(i1)*Gstep) .gt. tol) then > > > ? ? ? ? ? ? ? if (abs(res/(upsprevupsnod)) .gt. (i1)*Gstep) then ?! > > > secant inward towards 1.0 > > > ? ? ? ? ? ? ? ? ?xrtemp = xrnod > > > ? ? ? ? ? ? ? ? ?xrnod =(xrnod+xrprev)/2. > > > ? ? ? ? ? ? ? ? ?UPSNOD = ACOS(XHP/XHMXRNOD/XHM) > > > ? ? ? ? ? ? ? ? ?goto 11 > > > ? ? ? ? ? ? ?else ? !secant > > > outward > > > ? ? ? ? ? ? ? ? ?xrnod =(xrnod+xrtemp)/2. > > > ? ? ? ? ? ? ? ? ?UPSNOD = ACOS(XHP/XHMXRNOD/XHM) > > > ? ? ? ? ? ? ? ? ?goto 11 > > > ? ? ? ? ? ? ?end if > > > ? ? ? ? end if > > > ? ? ? ? ? ?GNOD = SIGN(real(res/(upsprevupsnod)),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 > > 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/(upsprevupsnod)),1.) > > This statement DOES NOT do what you think it does. As you use it, it > > returns ABS(real(res/(upsprevupsnod)). 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 (AH,OZ) as REAL, > > hence the REAL() is not needed. see 1a too. > > 3. I don't know exactly how your twoway branch works, but it should > > be symmetric in some sense. > > 4. (i1)*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/(upsprevupsnod)) 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/XHMXRNOD/XHM) > > ? ? ?ALPHA = (I1)*Gstep !## > > 11 CALL QROMO (FUNCG, upsnod, upsprev,res,MIDPNT) > > ? ? ?BETA = ABS(RES/(UPSPREVUPSNOD)) !## > > ? ? ?DELTA = ALPHA  BETA !## > > ! ? ? ? if (abs(abs(res/(upsprevupsnod))(i1)*Gstep) .gt. tol) then > > ! ? ? ? ? ? ? if (abs(res/(upsprevupsnod)) .gt. (i1)*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/XHMXRNOD/XHM) > > ? ? ? ? ? ? ? ? ?goto 11 > > ? ? ? ? ? ? ?else IF (DELTA .LT. 0) THEN !secant outward > > ? ? ? ? ? ? ? ? ?xrnod =(xrnod+xrtemp)/2. > > ? ? ? ? ? ? ? ? ?UPSNOD = ACOS(XHP/XHMXRNOD/XHM) > > ? ? ? ? ? ? ? ? ?goto 11 > > ? ? ? ? ? ? ?else > > ? ? ? ? ? ? ? ? ?WHAT DO YOU WANT HERE IF DELTA .EQ. 0? > > ? ? ? ? ? ? ?end if > > ? ? ? ? end if > > !#### ? ? ? GNOD = SIGN(real(res/(upsprevupsnod)),1.) > > ? ? ? ? WHAT YOU REALLY WANT HERE > > ? ? GPREV = GNOD > > ? ? XRPREV = XRNOD > > 10 CONTINUE > > See this diagram for conditions: > > ? ? ? ? ? ? ? ?BETATOL ? ? BETA ? ?BETA+TOL > > ALPHA: ? ? ? < ? ? ? ? ? ? ? ? ? ?> > > ALPHA: ? ? ? ? <> > > HTH > >  e Hide quoted text  > >  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 selfcontained (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 twoway 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.E7 > 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 = (i1)*Gstep > ? ? ?xrtemp1=xrgmax > ? ? ?xrtemp2=xrprev > ? ? ?XRNOD = XRGMAX > ? ? ?XRNOD = XRGMAX > ? ? ?UPSNOD = ACOS(XHP/XHMXRNOD/XHM) > ! res is the result of integrating FUNCG from upsnod to upsprev > 11 CALL QROMO (FUNCG, upsnod, upsprev,res,MIDPNT) > ? ? ?beta = ABS(res/(upsprevupsnod)) > ? ? ?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/XHMXRNOD/XHM) > ? ? ? ? ? ? ? ? ? ?goto 11 > ? ? ? ? ? ? ?else ? !secant > outward > ? ? ? ? ? ? ? ? ? ?xrtemp2 = xrnod > ? ? ? ? ? ? ? ? ? ?xrnod =(xrnod+xrtemp1)/2. > ? ? ? ? ? ? ? ? ? ?UPSNOD = ACOS(XHP/XHMXRNOD/XHM) > ? ? ? ? ? ? ? ? ? ?goto 11 > ? ? ? ? ? ? ?end if > ? ? ? ? end if > ? ? ? ? ? ?GNOD = SIGN(real(res/(upsprevupsnod)),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) = .329564E01 GM( 2) = 0.457676E02 GM( 3) = .265454E03 GM( 4) = 0.101691E03
... read more »

Mon, 09 May 2011 11:41:58 GMT 


e p chandle #5 / 17

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.E7 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 = (i1)*Gstep xrtemp1=xrgmax xrtemp2=xrprev XRNOD = XRGMAX XRNOD = XRGMAX UPSNOD = ACOS(XHP/XHMXRNOD/XHM) ! res is the result of integrating FUNCG from upsnod to upsprev 11 CALL QROMO (FUNCG, upsnod, upsprev,res,MIDPNT) beta = ABS(res/(upsprevupsnod)) 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/XHMXRNOD/XHM) goto 11 else !secant outward xrtemp2 = xrnod xrnod =(xrnod+xrtemp1)/2. UPSNOD = ACOS(XHP/XHMXRNOD/XHM) goto 11 end if end if GNOD = SIGN(real(res/(upsprevupsnod)),1.) GPREV = GNOD XRPREV = XRNOD 10 CONTINUE RETURN END  end of quoted text  Well you've corrected the problem with xrtemp being uninitialized. 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 


e p chandle #6 / 17

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.E7 > 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 = (i1)*Gstep > xrtemp1=xrgmax > xrtemp2=xrprev > XRNOD = XRGMAX > XRNOD = XRGMAX > UPSNOD = ACOS(XHP/XHMXRNOD/XHM) > ! res is the result of integrating FUNCG from upsnod to upsprev > 11 CALL QROMO (FUNCG, upsnod, upsprev,res,MIDPNT) > beta = ABS(res/(upsprevupsnod)) > 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/XHMXRNOD/XHM) > goto 11 > else !secant > outward > xrtemp2 = xrnod > xrnod =(xrnod+xrtemp1)/2. > UPSNOD = ACOS(XHP/XHMXRNOD/XHM) > goto 11 > end if > end if > GNOD = SIGN(real(res/(upsprevupsnod)),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) = .329564E01 GM( 2) = 0.457676E02 GM( 3) = .265454E03 GM( 4) = 0.101691E03 GM( 5) = 0.523956E04 GM( 6) = 0.112262E04 GM( 7) = .596746E06 GM( 8) = 0.863187E05 GM( 9) = 0.398432E05 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/XHMXR/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.51E05 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.51E05 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 rangeofinterest. 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 nonhomogeneous 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

Mon, 09 May 2011 13:50:17 GMT 


Georg #7 / 17

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.E7 >> 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 = (i1)*Gstep >> xrtemp1=xrgmax >> xrtemp2=xrprev >> XRNOD = XRGMAX >> XRNOD = XRGMAX >> UPSNOD = ACOS(XHP/XHMXRNOD/XHM) >> ! res is the result of integrating FUNCG from upsnod to upsprev >> 11 CALL QROMO (FUNCG, upsnod, upsprev,res,MIDPNT) >> beta = ABS(res/(upsprevupsnod)) >> 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/XHMXRNOD/XHM) >> goto 11 >> else !secant >> outward >> xrtemp2 = xrnod >> xrnod =(xrnod+xrtemp1)/2. >> UPSNOD = ACOS(XHP/XHMXRNOD/XHM) >> goto 11 >> end if >> end if >> GNOD = SIGN(real(res/(upsprevupsnod)),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) = .329564E01 > GM( 2) = 0.457676E02 > GM( 3) = .265454E03 > GM( 4) = 0.101691E03 > GM( 5) = 0.523956E04 > GM( 6) = 0.112262E04 > GM( 7) = .596746E06 > GM( 8) = 0.863187E05 > GM( 9) = 0.398432E05 > 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/XHMXR/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.51E05 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.51E05 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 rangeofinterest. 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 nonhomogeneous 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 Day http://apod.nasa.gov/apod/

Mon, 09 May 2011 17:15:24 GMT 


moni #8 / 17

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.E7 > >> 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 = (i1)*Gstep > >> xrtemp1=xrgmax > >> xrtemp2=xrprev > >> XRNOD = XRGMAX > >> XRNOD = XRGMAX > >> UPSNOD = ACOS(XHP/XHMXRNOD/XHM) > >> ! res is the result of integrating FUNCG from upsnod to upsprev > >> 11 CALL QROMO (FUNCG, upsnod, upsprev,res,MIDPNT) > >> beta = ABS(res/(upsprevupsnod)) > >> 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/XHMXRNOD/XHM) > >> goto 11 > >> else !secant > >> outward > >> xrtemp2 = xrnod > >> xrnod =(xrnod+xrtemp1)/2. > >> UPSNOD = ACOS(XHP/XHMXRNOD/XHM) > >> goto 11 > >> end if > >> end if > >> GNOD = SIGN(real(res/(upsprevupsnod)),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) = .329564E01 > > ? GM( 2) = 0.457676E02 > > ? GM( 3) = .265454E03 > > ? GM( 4) = 0.101691E03 > > ? GM( 5) = 0.523956E04 > > ? GM( 6) = 0.112262E04 > > ? GM( 7) = .596746E06 > > ? GM( 8) = 0.863187E05 > > ? GM( 9) = 0.398432E05 > > 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/XHMXR/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.51E05 ? ? ? 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.51E05 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 rangeofinterest. ?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 nonhomogeneous 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 rangeofinterest. 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.numanalysis DG. Regards. Monir

Mon, 09 May 2011 21:12:51 GMT 


robi #9 / 17

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 


moni #10 / 17

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 welldefined 1D functions over the XRrangeofinterest. 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 


robi #11 / 17

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. By your own admission, the program "works" EXCEPT for that particular nested IF!

Sun, 22 May 2011 20:07:17 GMT 


moni #12 / 17

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 selfcontained 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 recentlyadded 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.E7 XRPREV = 1.0 XRGMAX = 0.70432264 UPSPREV = pi NFTIP = 15 GMAX = 0.033839315 Gstep = abs(0.995*GMAX/(NFTIP1)) DO 10 I=2, NFTIP + 1 if (i .eq. nftip+1) then alpha = (nftip1)*Gstep else alpha = (i1.5)*Gstep end if xrtemp1=xrgmax xrtemp2=xrprev XRNOD = XRGMAX UPSNOD = ACOS(XHP/XHMXRNOD/XHM) ! res is the result of integrating FUNCG from upsnod to upsprev 11 CALL QROMO (FUNCG, upsnod, upsprev,res,MIDPNT) beta = ABS(res/(upsprevupsnod)) 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/XHMXRNOD/XHM) goto 11 else !move inward towards 0.704... xrtemp2 = xrnod xrnod =(xrnod+xrtemp1)/2. UPSNOD = ACOS(XHP/XHMXRNOD/XHM) goto 11 end if end if XRPREV = XRNOD 10 CONTINUE RETURN END Would appreciate your comments and any suggestions you might have. Regards.

Mon, 23 May 2011 00:35:28 GMT 


e p chandle #13 / 17

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 selfcontained 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 recentlyadded 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.E7 > XRPREV = 1.0 > XRGMAX = 0.70432264 > UPSPREV = pi > NFTIP = 15 > GMAX = 0.033839315 > Gstep = abs(0.995*GMAX/(NFTIP1)) > DO 10 I=2, NFTIP + 1 > if (i .eq. nftip+1) then > alpha = (nftip1)*Gstep > else > alpha = (i1.5)*Gstep > end if > xrtemp1=xrgmax > xrtemp2=xrprev > XRNOD = XRGMAX > UPSNOD = ACOS(XHP/XHMXRNOD/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/(upsprevupsnod)) > 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/XHMXRNOD/XHM) > goto 11 > else !move inward towards 0.704... > xrtemp2 = xrnod > xrnod =(xrnod+xrtemp1)/2. > UPSNOD = ACOS(XHP/XHMXRNOD/XHM) > goto 11 > end if > end if > XRPREV = XRNOD
### to here ### Quote: > 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 unindented 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/(upsprevupsnod)) 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/XHMXRNOD/XHM) else ! move inward towards 0.704... xrtemp2 = xrnod xrnod =(xrnod+xrtemp1)/2. UPSNOD = ACOS(XHP/XHMXRNOD/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

Mon, 23 May 2011 01:45:34 GMT 


moni #14 / 17

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


e p chandle #15 / 17

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