CODEWITHSUNDEEP
fortranfortran90
Maarten Dhondt
Maarten Dhondt

Reputation: 657

Is it possible to eliminate do loop

As we know, more recent versions of Fortran support array operations, which can eliminate many loops. So I was wondering if it would be possible to eliminate even the last remaining loop in following code snippet (as to make it a one-liner):

subroutine test(n,x,lambda)
   integer, intent(in) :: n
   real, dimension(:), intent(in) :: x
   real, dimension(:), intent(out) :: lambda
   real :: eps
   integer :: i

   do i=1,n
      lambda(i) = product(x(i)-x, mask=(abs(x(i)-x) > epsilon(eps)))
   enddo
end subroutine

Its intention is to calculate n lambda(i) values in which

lambda(i) = (x(i)-x(1))*(x(i)-x(2))*...*(x(i)-x(i-1)*(x(i)-x(i+1))*...*(x(i)-x(n))

Upvotes: 0

Views: 101

Answers (3)

High Performance Mark
High Performance Mark

Reputation: 78324

OK, try this

lambda = product(max(spread(x, dim=1, ncopies=size(x)) - &
            spread(x, dim=2, ncopies=size(x)), eps), dim=2)

That's a one-liner. It's also rather wasteful of memory and much less comprehensible than the original.

Upvotes: 2

Alexander Vogt
Alexander Vogt

Reputation: 18098

Yes, you can shorten this, product can use 2D arrays:

  1. You would first need to set up a matrix of the differences: 
   do i=1,n
     mat(:,i) = x(i) - x
   enddo

or, as a one-liner: 

   forall ( i=1:n ) mat(:,i) = x(i) - x
  1. Now you can do the product along the second dimension: 
   lambda = product(mat, dim=2, mask=(abs(mat) > epsilon(eps)))

Whole program:

program test
   integer, parameter  :: n = 3
   real, dimension(n) :: x
   real, dimension(n) :: lambda
   real, dimension(n,n) :: mat
   real :: eps = 1.
   integer :: i

   call random_number( x )

   do i=1,n
      lambda(i) = product(x(i)-x, mask=(abs(x(i)-x) > epsilon(eps)))
   enddo
   print *,lambda

   forall ( i=1:n ) mat(:,i) = x(i) - x
   lambda = product(mat, dim=2, mask=(abs(mat) > epsilon(eps)))
   print *,lambda

end program

Upvotes: 1

WWhisperer
WWhisperer

Reputation: 97

Have you tried it with an implied do-loop in the array creation? something like real, dimension(:), intent(out):: lambda =(/product(x(i)-x, mask=(abs(x(i)-x)>epsilon(eps))), i=1, n/) ... I am not sure about the syntax here, but something like that might work. You might even be able to create the array without calling the subroutine and do it in your main program, if your x-array is available.

Hope it helps.

Upvotes: 1

Related Questions