6

When we deal with large arrays, it may be important to consider the cost of change of rank and shape of arrays specially when it happens a couple of times in multiple subroutines/functions.

The main purpose of my question is to change rank of arrays from 2nd to 1st and vice versa...

To do this one can use :

  1. RESHAPE statement

  2. POINTER variables. Below is a code which shows how to utilize pointer variables:

    program test
    real, pointer :: a(:)
    real, pointer :: b(:,:)

    allocate(a(6))
    a = (/1,2,3,4,5,6/)
    b (1:2, 1:3) => a

    WRITE(*,*) 'b=',b(3,1)
end program test

Questions : 1. Which method is faster? 2. Is there any other faster method? 3. Any other suggestion to do this work?

Thanks...

Vahid
  • 153
  • 1
  • 2
  • 15
  • 1
    You could try this yourself... In general, this is implementation specific and not defined by the standard. I would suspect that pointers are usually faster than transpositions, as no copying is required. I know that `gfortran` tries to avoid/delay copies when transposing arrays, but I guess there are scenarios where this cannot be avoided. – Alexander Vogt Jul 01 '15 at 19:20
  • 1
    Since the *main purpose of [the] question is to change rank of arrays from 2nd to 1st*, why don't you simply work on array slices? As long as the data is contiguous, this would be fastest (I guess). – Alexander Vogt Jul 01 '15 at 19:23
  • @AlexanderVogt I assume by slicing you mean taking the rows or columns. I know matrix to vector multiplications are based on slicing the matrix and using a loop. I think the compiler should deal with tens of row/columns slicing in that case. by using pointer, there is no need to copy the elements of array in memory, which might be faster. – Vahid Jul 01 '15 at 19:55
  • Why don't you do some benchmarks and report the results? – Alexander Vogt Jul 01 '15 at 20:03
  • Faster to what, to create it? To access it later? – Vladimir F Героям слава Jul 01 '15 at 20:18
  • @HighPerformanceMark I agree that the best is to not move the data at all. I tried **EQUIVALENCE** statement. I think this is the fastest way. However, when I tried to implement it to my code I encounter the DUMMY attribute conflict error. It is inside a subroutine... – Vahid Jul 01 '15 at 23:13
  • There are other ways of remapping ranks - such as via sequence and storage association. All methods have upsides and downsides. The answer to this question completely depends on what you want to do with the arrays. – IanH Jul 02 '15 at 00:14

1 Answers1

4

Well, Fortran is designed to be the language of MATH. I dug a little bit and I came across the below points in Fortran.

Some explanation before explaining the points :

My subroutine must work with 1st rank arrays. I am calling 2nd rank arrays as input at the beginning of the subroutine. Then, I need to change the rank from 2nd to 1st. Later in the subroutine, I need to change the rank back to 2. This changing of ranks is happening 3-4 times in the code.

  1. USING EQUIVALENCE statement:

This is the fastest method. Nothing changes in the memory and I thought it is the best. However, it ends with attribute conflict errors for my problem, because I am working inside a subroutine.

  1. Using Pointer attribute :

I have tried pointer then. But, it seems that it is not possible to remap a 2nd rank array into a 1st rank array. Remapping 1st rank to 2nd rank arrays works fine.

The simple code, I wrote to remap a 1st rank array into a 2nd rank array :

program ptrtest
real, pointer :: a(:)
real, pointer :: b(:,:)

allocate(b(1:2,1:3))
b = transpose(reshape((/ 1, 2, 3, 4, 5, 6 /), shape(b)))
a(1:6) => b(:,:) 

WRITE(*,*) a(4), b(2,2) ! to see if the remapped elements are same?
end program ptrtest

The error that I have received :

gfortran -Wall -o "POINTER" "POINTER.f90" (in directory: /home/vahid/Desktop)
POINTER.f90:12.14:
a(1:6) => b(:,:)
          1
Error: Rank remapping target must be rank 1 or simply contiguous at (1)
Compilation failed.
  1. RESHAPE statement:

The slowest method which is able to do any type of transition. Basically, it assigns another memory locations for the transitioned element which is expensive, considering both memory efficiency and processing costs.

As a consequence, Fortran 2003 manual states that: (section 2.2.4 and 5.2)

Data objects may be dynamic in size, shape, type, or length type parameters, but not rank or kind type parameters.

I don't know the consequences, but I think the rank of arrays should be dynamic also. Please correct if any part is wrong.

Vahid
  • 153
  • 1
  • 2
  • 15
  • 2
    Regarding point (2), if you do as the error suggest and give `b` the attribute `contiguous` then it works. – casey Jul 02 '15 at 01:40