0

Consider that you have a square matrix M(N,N) and you want to sum all pairs such that:

M(i,j)*M(i+1,j)+M(i,j)*M(i,j+1)

to do this, the easiest way is to compute:

INTEGER i,j,N, SUMT
INTEGER M(100,100), c(101)

N=100
SUMT = 0
do j=1,N   
  c(j) = j
end do
c(N+1)=1

do j=1,N   
   do i=1,N
      SUMT = SUMT + M(i,j)*M(c(i+1),j)+M(i,j)*M(i,c(j*1))
   end do 
end do

NOTE: c is a fast way to apply periodic boundary condition.

In my problem, for a 3D system { M(N,N,N) } I should do the following:

M(i,j,k)*M(i+1,j,k)+M(i,j,k)*M(i,j+1,k)+M(i,j,k)*M(i,j,k+1)

So the code is:

INTEGER i,j,k,N, SUMT
INTEGER M(100,100), c(101)
SUMT = 0
do j=1,N   
  c(j) = j
end do

c(N+1)=1 
N=100

do j=1,N   
   do i=1,N
      SUMT = SUMT +M(i,j,k)*M(c(i+1),j,k)+M(i,j,k)*M(i,c(j+1),k)+M(i,j,k)*M(i,j,c(k+1))
   end do 
end do

At this point, my question is:

Is there any way to compute this problem with nested loops such dimension of M matrix is a parameter? I mean, I could do:

INTEGER i,j,k,l,m,n,....
INTEGER N, SUMT, D
PARAMETER (N=100)
PARAMETER (D=3) !DIMENSION

INTEGER M(N**D), c(N+1)

if (dim=1) then
  do i=1,N
else if (dim=2) then
  do j=1,N
    do i=1,N
else if (dim=3) then
  do k=1,N
    do j=1,N
       do i=1,N
 ....

but do you think about a more elegant solution in fortran 77?

I was thinking in accessing the M matrix with dim D as if it would have only one dimension with N**D spaces but I think that if apply if-instructions inside the to control limit of N it will work very slow. Any good idea or I should consider the nasty if-do loops?

Learning from masters
  • 2,032
  • 3
  • 29
  • 42
  • 4
    If you care about elegance or brevity, forget Fortran **77**, it is 2015. But I fail to see the problem. Why do you need all three branches when it is parameter? The compiler will discard them during optimizations anyway. – Vladimir F Героям слава Nov 16 '15 at 18:54

1 Answers1

0

I don't know fortran and not sure if this is you want, but here's some pseudocode that might help

INTEGER N, D, m
INTEGER I(D)     ! indices in each dimension
do j=1,N**D      ! loop over whole matrix
  m = j          
  do k=1,D           ! for each dimension
    I(k) = MOD(m,N)  ! get index in dimension k
    m = m / N        ! and move to next most-significant "register"
  end do
  ! now you have an array of indices for this item
  ! do stuff here
end do

e.g. D=3, N=100 means j goes from 1 to 1000000
I successively becomes the array [1,1,1], [2,1,1], [3,1,1], ... [100,100,100]

Sanjay Manohar
  • 6,920
  • 3
  • 35
  • 58