how to initialize a large array in Fortran?

Question

I have a Fortran function in which I would like to initialize a large array at compile time. A simplified working example is below, where the parameter coeff in fill_coefficients has been reduced in size greatly.

How do I write similar code when coeff is large, without exceeding the maximum of 255 continuation lines, or the maximum of 132 characters per line? Here fill_coefficients should really be PURE, which probably makes it impossible to read coeff from a file once during runtime, and then store the result.

The file "main.f03":

    PROGRAM main
        USE coefficients
        IMPLICIT NONE

        REAL(dp), ALLOCATABLE, DIMENSION(:,:) :: matrix

        CALL fill_coefficients(matrix,2)

        PRINT *, "The first row of 'matrix':"
        PRINT *, matrix(1,:)
    END PROGRAM main

The file "coefficients.f03":

    MODULE coefficients
        USE iso_fortran_env
        IMPLICIT NONE

        INTEGER, PARAMETER :: dp = REAL64

    CONTAINS
        PURE SUBROUTINE fill_coefficients(my_coefficients, n)
            IMPLICIT NONE
            REAL(dp), ALLOCATABLE, DIMENSION(:,:), INTENT(OUT) :: my_coefficients
            INTEGER, INTENT(IN) :: n

            ! The size of the following array would be roughly 200 x 200 = 40.000.
            REAL(dp), DIMENSION(3,3), PARAMETER :: coeff = &
                RESHAPE ( &
                [ + 10.6770782520313112108115239655957106_dp, &
                - 854.166260162504896864921917247656850_dp, &
                - 85.4166260162504896864921917247656850_dp, &
                + 16250.5130995916556628551394756366716_dp, &
                + 6747.91345528378868523288314625648912_dp, &
                + 106.770782520313112108115239655957106_dp, &
                - 123256.191341449456617608232658836883_dp, &
                - 8328.12103658442274443298869316465429_dp, &
                + 500381.272281447399894682070647642979_dp ], &
                [3,3] )

            IF (ALLOCATED(my_coefficients)) DEALLOCATE(my_coefficients)
            ALLOCATE(my_coefficients(n,n))

            my_coefficients = coeff(1:n,1:n)
        END SUBROUTINE fill_coefficients
    END MODULE coefficients

The output:

The first row of 'matrix': 10.677078252031311 16250.513099591655

If it is that large, can't you just read the thing from a file? Do you need it to be `parameter`? — Vladimir F Героям слава, Feb 15 '16 at 15:27
I can't find the answer/question I thought existed, but you can also see http://stackoverflow.com/q/25811884. You can effectively build up a named constant from other smaller named constants in a similar way to that linked answer. — francescalus, Feb 15 '16 at 15:27
Further to @VladimirF's comment: you say "should really be pure". Do you really need it to be pure (which invalidates my first suggested duplicate) or are you asking about that requirement also? That is, can we make suggestions based on not being pure? — francescalus, Feb 15 '16 at 15:30
I would like `fill_coefficients` to be PURE, because in about a year from now I will probably need to do run my FEM multi-core, and this routine is in the innermost loop. I believe this makes it impossible to read `coeff` from a file once during runtime, and then store the result. I could perhaps use the 'trick' of building up the array from several other smaller array, or do this in many steps. Besides the ugliness of that solution, wouldn't I end up storing all the elements several times? — knia, Feb 15 '16 at 15:37
I thought of the [DATA statement](http://www.idris.fr/data/cours/lang/fortran/F77.html#p7.3.1), but I don't see how that splits my statement into many lines. — knia, Feb 15 '16 at 15:40
You can't use `data` statements for a named constant, and can't use them for a local variable of the procedure (that conflicts with the pureness). So, for something with the `parameter` attribute the "building up" is the way to go. That probably won't "store" things twice. — francescalus, Feb 15 '16 at 15:41
As for "ugliness", I think storing 40000 values in the source code is going to be pretty horrifying with any approach. — francescalus, Feb 15 '16 at 15:43
Regarding your intended use and the requirement for `pure`: I suspect you'd probably want to think about approaches which don't involve that huge array assignment in the inner loop. Perhaps it may be possible to read it from a file and then just do things with array sections (like how we did things back with static memory layout in the F77 days). [But that's off-topic: if you really want to require `pure` I'm not going to argue - there are times when that is what is wanted.] — francescalus, Feb 15 '16 at 15:46
I should perhaps add that this is not really source code, as I am using python to generate it. It produces fairly elegant Fortran code. It doesn't really copy a huge array, but it needs segments of `coeff` (like `r = matmul(coeff(1:20,1:20),vector)`). — knia, Feb 15 '16 at 15:54

score 3 · Accepted Answer · answered Feb 15 '16 at 21:05

From a maintenance perspective (and as perhaps suggested in the comments), I would read the data into a module variable in a separate non-pure subroutine that is called once at program start-up. fill_coefficients then becomes a simple assignment from that module variable and can still be PURE.

MODULE coefficients
  IMPLICIT NONE
  ...
  ! Could be PUBLIC, PROTECTED, then you could directly 
  ! assign from it and dispense with fill_coefficients 
  ! altogether.
  REAL(dp), PRIVATE :: coeff(200,200)
CONTAINS
  SUBROUTINE init
    INTEGER :: unit 
    OPEN( NEWUNIT=unit,  &
          FILE='lots-of-numbers.bin',  &
          FORM='UNFORMATTED',  &
!         ACCESS='STREAM',  &    ! Maybe - depending on how you write it.
          STATUS='OLD' )
    READ (unit) coeff
    CLOSE(unit)
  END SUBROUTINE init

  PURE SUBROUTINE fill_coefficients(my_coefficients, n)
    ! implicit none already in force due to the statement in 
    ! the specification part of the host module.
    ! IMPLICIT NONE    
    REAL(dp), ALLOCATABLE, DIMENSION(:,:), INTENT(OUT) :: my_coefficients
    INTEGER, INTENT(IN) :: n

    ! This test is redundant - my_coefficients is INTENT(OUT) so 
    ! it must be not allocated at this point.
    ! IF (ALLOCATED(my_coefficients)) DEALLOCATE(my_coefficients)

    ! This allocate statement is redundant - allocation will 
    ! happen automatically under F2003 with the assignment.
    ! ALLOCATE(my_coefficients(n,n))

    my_coefficients = coeff(1:n,1:n)
  END SUBROUTINE fill_coefficients
END MODULE coefficients

If you must have coeff as a compile time parameter, then assemble it in source manageable chunks - perhaps column by column. Your limits per declaration are line length (132) and number of continuation lines (255).

REAL(dp), PARAMETER :: column_1(200) = [  &
     + 10.6770782520313112108115239655957106_dp, &
     - 854.166260162504896864921917247656850_dp, &
     - 85.4166260162504896864921917247656850_dp, &
     ... ]
REAL(dp), PARAMETER :: column_2(200) = [ ... ]
...
REAL(dp), PARAMETER :: column_200(200) = [ ... ]

REAL(dp), PARAMETER :: coeff(200,200) = RESHAPE( [  &
       column_1, column_2, ..., column_200 ],  &
     SHAPE=[200,200] )

Things declared with PARAMETER are named constants. Conceptually these only exist at compile time - depending on what you do with a named constant the compiler may or may not set aside storage in the executable image for the constants.

Large named constants may result in the compiler having issues compiling the file.

I'm auto-generating to Fortran code to be used on any computer, which makes me even more skeptical of `FORM='UNFORMATTED'` than I usually already am. But this is certainly the right idea, though it's a little trickier to auto-generate the whole Fortran code from Python this way. — knia, Feb 16 '16 at 08:35
You can write and read the file in as human readable text if you want, it will just take a little bit longer for the IO operations, but on the Fortran side they only need to happen once per image. — IanH, Feb 16 '16 at 09:55
For the record, I'll mention my final solution, which is a little more straightforward to implement when auto-generating the code. Instead of reading from a file, I just do `SUBROUTINE init()` `coeff(1,1)=1.0_dp ` `coeff(2,1)=2.0_dp ` `.... ` `coeff(199,200)=0.0_dp ` `coeff(200,200)=-43.7_dp ` `END SUBROUTINE init ` The slightly bloated size of the program is no issue for me, and the module as a whole is very simple, even if there is this long list at the beginning. Thanks for all the help! — knia, Feb 16 '16 at 12:10

how to initialize a large array in Fortran?

1 Answers1

Linked

Related