array operation in fortran

user_4571350 0
Viewed 214

I am writing a code with a lot of 2D arrays and manipulation of them. I would like the code to be as concise as possible, for that I would like to use as many 'implicit' operation on array as possible but I don't really know how to write them for 2D arrays. For axample:

DO  J=1,N
   DO  I=1,M
    A(I,J)=B(J)*A(I,J)
   ENDDO
ENDDO

become easily:

DO  J=1,N
 A(:,J)=B(J)*A(:,J)
ENDDO

Is there a way to reduce also the loop J?

Thanks

fortran
3 Answers

For brevity and clarity, you could wrap these operations in a derived type. I wrote a minimal example which is not so concise because I need to initialise the objects, but once this initialisation is done, manipulating your arrays becomes very concise and elegant.

I stored in arrays_module.f90 a derived type arrays2d_T which can hold the array coefficients, plus useful information (number of rows and columns). This type contains procedures for initialisation, and the operation you are trying to perform.

module arrays_module

implicit none

integer, parameter :: dp = kind(0.d0) !double precision definition

 type :: arrays2d_T
   real(kind=dp), allocatable :: dat(:,:)
   integer :: nRow, nCol     
  contains
   procedure :: kindOfMultiply => array_kindOfMuliply_vec
   procedure :: init           => initialize_with_an_allocatable
 end type

 contains

 subroutine initialize_with_an_allocatable(self, source_dat, nRow, nCol)
   class(arrays2d_t), intent(inOut) :: self
   real(kind=dp), allocatable, intent(in) :: source_dat(:,:)
   integer, intent(in) :: nRow, nCol
   allocate (self%dat(nRow, nCol), source=source_dat)
   self%nRow = nRow
   self%nCol = nCol
 end subroutine

 subroutine array_kindOfMuliply_vec(self, vec)
   class(arrays2d_t), intent(inOut) :: self
   real(kind=dp), allocatable, intent(in) :: vec(:)
   integer :: iRow, jCol
   do jCol = 1, self%nCol
   do iRow = 1, self%nRow
      self%dat(iRow, jCol) = vec(jCol)*self%dat(iRow, jCol)
   end do
   end do
 end subroutine

end module arrays_module

Then, in main.f90, I check the behaviour of this multiplication on a simple example:

program main
use arrays_module
implicit none

type(arrays2d_T) :: A
real(kind=dp), allocatable :: B(:)

! auxilliary variables that are only useful for initialization
real(kind=dp), allocatable :: Aux_array(:,:)
integer :: M = 3
integer :: N = 2

! initialise the 2d array
allocate(Aux_array(M,N))
Aux_array(:,1) = [2._dp, -1.4_dp, 0.3_dp]
Aux_array(:,2) = [4._dp, -3.4_dp, 2.3_dp]
call A%init(aux_array, M, N)

! initialise vector
allocate (B(N))
B = [0.3_dp, -2._dp]

! compute the product
call A%kindOfMultiply(B)

print *, A%dat(:,1)
print *, A%dat(:,2)
end program main

Compilation can be as simple as gfortran -c arrays_module.f90 && gfortran -c main.f90 && gfortran -o main.out main.o arrays_module.o

Once this object-oriented machinery exists, call A%kindOfMultiply(B) is much clearer than a FORALL approach (and much less error prone).

user_45861550

A one-line solution can be achieved by using FORALL:

FORALL(J=1:N) A(:,J) = B(J)*A(:,J)

Note that FORALL is deprecated in the most recent versions of the standard, but as far as I know, that is the only way you can perform that operation as a single line of code.

user_47028700

No one has mentioned do concurrent construct here, which has the potential to automatically parallelize and speed up your code,

do concurrent(j=1:n); A(:,j)=B(j)*A(:,j); end do
user_20886940
Related