Intel® Fortran Compiler
Build applications that can scale for the future with optimized code designed for Intel® Xeon® and compatible processors.

Nested do loops

JohnNichols
Valued Contributor III
875 Views
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)
      end

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

 

9 Replies
JohnNichols
Valued Contributor III
872 Views

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.  

0 Kudos
JohnNichols
Valued Contributor III
862 Views

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. 

 

0 Kudos
jimdempseyatthecove
Honored Contributor III
854 Views

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

The 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

 

0 Kudos
JohnNichols
Valued Contributor III
814 Views

@jimdempseyatthecove 

 

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.  

 

0 Kudos
Barbara_P_Intel
Employee
757 Views

@ndd21@JohnNichols started a new thread for you with your issue. Thanks, John!

0 Kudos
JohnNichols
Valued Contributor III
742 Views

@Barbara_P_Intel 

 

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.  

 

maxresdefault.jpg

0 Kudos
Barbara_P_Intel
Employee
731 Views

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

0 Kudos
jimdempseyatthecove
Honored Contributor III
690 Views

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

The actual solution will depend on the interval values required. Often it is not as simple as above.

Jim Dempsey

 

0 Kudos
JohnNichols
Valued Contributor III
681 Views

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

 

 

0 Kudos
Reply