- Mark as New
- Bookmark
- Subscribe
- Mute
- Subscribe to RSS Feed
- Permalink
- Report Inappropriate Content
NESTED DO loops
Suppose we have the following simple program:
program test
parameter (m=4)
dimension ib(m),ie(m),is(m),i(m)
data ib/1,1,1,1/
data ie/3,4,-3,2/
data is/1,2,-2,1/
write(*,*) 'Enter a number >=1,and <=',m
read(*,*)n
if (n.gt.m) stop
if (n.lt.1) stop
call donest(ib,ie,is,n,i)
end
We want to construct: subroutine donest(ib,ie,is,n,i) such as to implement a variable number of NESTEDDO-loops.
For example for n=2 we should get the output that the following code yields:
EQUIVELANT CODE:
DO 1 I1=IB(1),IE(1),IS(1)
DO 2 I2=IB(2),IE(2),IS(2)
WRITE(*,*)I1,I2
2 CONTINUE
1 CONTINUE
OUTPUT:
1 1
1 3
2 1
2 3
3 1
3 3
Similarly for n=3 we have:
EQUIVELANT CODE:
DO 1 I1=IB(1),IE(1),IS(1)
DO 2 I2=IB(2),IE(2),IS(2)
DO 3 I3=IB(3),IE(3),IS(3)
WRITE(*,*)I1,I2,I3
3 CONTINUE
2 CONTINUE
1 CONTINUE
OUTPUT:
1 1 1
1 1 -1
1 1 -3
1 3 1
1 3 -1
1 3 -3
2 1 1
2 1 -1
2 1 -3
...
3 1 -3
Correspondigly for n=4, etc
Sample Solution
subroutine donest(ib,ie,is,n,i)
dimension ib(n),ie(n),is(n),i(n)
* Determine if the loops are doable and the exact right bounds
do 2 j=1,n
if (is(j).eq.0) return
k = (ie(j)-ib(j))/is(j)
if ( ib(j).le.ie(j) .and. is(j).gt.0 ) then
ie(j)=ib(j)+ k*is(j)
elseif ( ib(j).ge.ie(j) .and. is(j).lt.0 ) then
ie(j)=ib(j)+ k*is(j)
else
return
endif
2 continue
* Set initial values
do 3 j = 1,n
i(j) = ib(j)
3 continue
* Print all values now (running the innermost loop
4 continue
do 1 in = ib(n),ie(n),is(n)
i(n) = in
write(*,*)(i(j),j=1,n)
1 continue
*
do 10 j = n-1,1,-1
if ( i(j).ne.ie(j) ) then
i(j) = i(j)+is(j)
go to 4
else
i(j) = ib(j)
endif
10 continue
end
*
program test
parameter (m=4)
dimension ib(m),ie(m),is(m),i(m)
data ib/1,1,1,1/
data ie/3,4,-3,2/
data is/1,2,-2,1/
write(*,*) 'Enter a number >=1,and <=',m
read(*,*)n
if (n.gt.m) stop
if (n.lt.1) stop
call donest(ib,ie,is,n,i)
endWinston and Horn, Lisp 3rd Edition, page 72 has an excellent algorithm for problems like yours and a sample using the Fibonacci function that is 5 lines for the function, add 3 for the program and you have an 8 line program.
Fortran was not designed from the start for recursion, Lisp was, a second text is
https://www.amazon.com/Structure-Interpretation-Computer-Programs-Engineering/dp/0262510871 which is where any computer language degree should start, instead at my place they rely on Python.
Of course if you are going to do it in Fortran use the Winston method.
Link Copied
- Mark as New
- Bookmark
- Subscribe
- Mute
- Subscribe to RSS Feed
- Permalink
- Report Inappropriate Content
Lisp is also not strongly typed, so you have to be careful, but you do not run into the Fortran problem that
if(st .eq. 0.1) then
fails if st is double precision and 0.1 is taken as a real and yet st is set to 0.1.
Ran into that little problem yesterday on updating an old FE Model generator.
- Mark as New
- Bookmark
- Subscribe
- Mute
- Subscribe to RSS Feed
- Permalink
- Report Inappropriate Content
There is an old woodwork saying that if all you have is a hammer, then every screw is a nail. We should not see Fortran in these head lights.
- Mark as New
- Bookmark
- Subscribe
- Mute
- Subscribe to RSS Feed
- Permalink
- Report Inappropriate Content
John, you can do better than that:
! NestedDoLoops.f90
module mod_donest
! Note, while donest is recursive, as written it is not reentrant (thread-safe)
! if you wish to make this thread-safe, declare iLevel as threadprivate
integer, allocatable :: iLevel(:)
contains
recursive subroutine donest(ib,ie,is,n,i)
implicit none
! dummy arguments
integer :: ib(:), ie(:), is(:), n
integer, value :: i
! local variables
integer :: idx
! sanity checks (remove in production code)
if(n > size(ib)) stop "(n > size(ib))"
if(n > size(ie)) stop "(n > size(ie))"
if(n > size(is)) stop "(n > size(is))"
if(i < 1 .or. i > n) stop "(i < 1 .or. i > n)"
! code
if(allocated(iLevel)) then
if(n > size(iLevel)) then
deallocate(iLevel)
allocate(iLevel(n))
endif
else
allocate(iLevel(n))
endif
do idx = ib(i),ie(i),is(i)
iLevel(i) = idx
if(i == n) then
write(*,*) iLevel(1:n)
else
call donest(ib,ie,is,n,i+1)
endif
end do
end subroutine donest
end module mod_donest
program NestedDoLoops
use mod_donest
parameter (m=4)
dimension ib(m),ie(m),is(m),i(m)
data ib/1,1,1,1/
data ie/3,4,-3,2/
data is/1,2,-2,1/
do
write(*,*) 'Enter a number >=1,and <=',m
read(*,*)n
if (n.gt.m) stop
if (n.lt.1) stop
call donest(ib,ie,is,n,1)
end do
end program NestedDoLoops Enter a number >=1,and <= 4
1
1
2
3
Enter a number >=1,and <= 4
2
1 1
1 3
2 1
2 3
3 1
3 3
Enter a number >=1,and <= 4
3
1 1 1
1 1 -1
1 1 -3
1 3 1
1 3 -1
1 3 -3
2 1 1
2 1 -1
2 1 -3
2 3 1
2 3 -1
2 3 -3
3 1 1
3 1 -1
3 1 -3
3 3 1
3 3 -1
3 3 -3
Enter a number >=1,and <= 4
4
1 1 1 1
1 1 1 2
1 1 -1 1
1 1 -1 2
1 1 -3 1
1 1 -3 2
1 3 1 1
1 3 1 2
1 3 -1 1
1 3 -1 2
1 3 -3 1
1 3 -3 2
2 1 1 1
2 1 1 2
2 1 -1 1
2 1 -1 2
2 1 -3 1
2 1 -3 2
2 3 1 1
2 3 1 2
2 3 -1 1
2 3 -1 2
2 3 -3 1
2 3 -3 2
3 1 1 1
3 1 1 2
3 1 -1 1
3 1 -1 2
3 1 -3 1
3 1 -3 2
3 3 1 1
3 3 1 2
3 3 -1 1
3 3 -1 2
3 3 -3 1
3 3 -3 2
Enter a number >=1,and <= 4The reentrant code can do any levels.
As stated in the code, you can remove the sanity checks (and rely on runtime checks for array bounds issues in debug build).
Jim Dempsey
- Mark as New
- Bookmark
- Subscribe
- Mute
- Subscribe to RSS Feed
- Permalink
- Report Inappropriate Content
I did not write the code, a couple of people complained that someone had added this code to the end of a 3 year old very long thread. I just moved it to the forefront to make it easy for all.
It has been a few years since I wrote code of that quality, using the term quality loosely.
PS In one other post I used one instead of won, apologies.
- Mark as New
- Bookmark
- Subscribe
- Mute
- Subscribe to RSS Feed
- Permalink
- Report Inappropriate Content
- Mark as New
- Bookmark
- Subscribe
- Mute
- Subscribe to RSS Feed
- Permalink
- Report Inappropriate Content
I am amazed that you do these changes at 4 am in the morning. I was up at 4am doing FE modelling, but looking here, no. Only @jimdempseyatthecove is that dedicated and lovable, like Animal in the Muppets.
If the distance between two real numbers labelled X1 and X2 is d, and you split d in two and add it to X1 and subtract from X2 you should get X3 which is the same point so to speak, except in blasted Fortran in a FE model. I am going insane.
My daughter thanked me yesterday, she said Dad you can be really annoying sometimes, and then she looked at me through 16 yr. old eyes and said no all the time.
Got to love the honesty of a teenager.
- Mark as New
- Bookmark
- Subscribe
- Mute
- Subscribe to RSS Feed
- Permalink
- Report Inappropriate Content
@JohnNichols, temporarily working unusual hours/days from EDT, instead of my usual days/hours PDT. 4am is a bit early for FEM in ANY timezone!
Glad your daughter keeps you on your toes!
- Mark as New
- Bookmark
- Subscribe
- Mute
- Subscribe to RSS Feed
- Permalink
- Report Inappropriate Content
>>If the distance between two real numbers labelled X1 and X2 is d, and you split d in two and add it to X1 and subtract from X2 you should get X3 which is the same point so to speak, except in blasted Fortran in a FE model. I am going insane.
With infinite precision that assumption might be true. But alas, we do not have infinite precision, we must accept (and code for) limited precision and that all of are calculations are to be presumed approximate. (might.... I cannot say if infinite precision has round-off errors or not)
I think you are old enough to have used a slide-rule (I am), about 3 digits of precision.
Calculate X3 once (X1 + d/2 .OR. X2 - d/2, and do not assume the two are equal)
I am sure you know this, but other readers might not be aware of this. One should be aware of code that can accumulate errors
real, parameter :: DeltaT = 0.1 ! repeading fraction in floating point
real :: T
...
T = 0.0
do
(code)
if(mod(T, 60.) == 0.0) call report()
T = T + DeltaT
end do
You may find that the above program does not report as it should.
Worse yet, the accumulated values of T will drift either above or below the expected values of T
The better way to code this is:
real, parameter :: DeltaT = 0.1 ! repeating fraction in floating point
real :: T
integer :: interval
...
interval = 0
do
T = DeltaT * interval
(code)
if(mod(interval, 600) == 0) call report()
interval = interval + 1
end do
Or some variation on this:
real, parameter :: DeltaT = 0.1 ! repeating fraction in floating point
real :: T
integer :: interval
...
interval = 0
do
T = real(interval / 10) + DeltaT * mod(interval,10)
(code)
if(mod(interval, 600) == 0) call report()
interval = interval + 1
end doThe actual solution will depend on the interval values required. Often it is not as simple as above.
Jim Dempsey
- Mark as New
- Bookmark
- Subscribe
- Mute
- Subscribe to RSS Feed
- Permalink
- Report Inappropriate Content
@jimdempseyatthecove , I have no defence for such a lousy if statement, I have known better since 1979, but it was old code and I am in a desperate hurry with daily phone calls, you know the drill. I did not expect a change from sp to dp to break it.
I will return and fix it, just not till I get the real stuff out. The real people pay money, they do not care except for the pretty pictures.
@Barbara_P_Intel , if you are running FE models and they take 3 hours to run, you get up at midnight and 3am to restart them. If you are downloading updates to a BR1 Cell Phone modem in Washington, you get up hourly to check it at 7KB/sec
Life for a Fortran programmer is not easy.
Someone rewrote the Abelson and Sussman book in Javascript. Why?
- Subscribe to RSS Feed
- Mark Topic as New
- Mark Topic as Read
- Float this Topic for Current User
- Bookmark
- Subscribe
- Printer Friendly Page