0

I have an array with multiple dimensions (the goal is to allow for about 100) and each dimension has a size of about 2^10 and I only need to store in it about 1000 double precision coefficients. I don't need to do any operation with this array aside from reading and writing into it. The code is written in Fortran 90.

I assume that if I a library like one of the ones mentioned in this answer I would be able to store the do this, but would this be optimized for the simple reading and writing operations? Is there a library that would be most efficient for that purpose?

Edit: By "simple reading and writing operations" I mean the following. Suppose

REAL(8), DIMENSION(1000)   :: coeff1
INTEGER, DIMENSION(1000,5) :: index

I want to define coeff2 to store the values in coeff1 and then read itat the indices in index, that is

DO i = 1,1000
     index(i,:) = [something]
     coeff1(i)  = [another something]
     coeff2(index(i,1),index(i,2),index(i,3),index(i,4),index(i,5)) = coeff1(i)
ENDDO

Then, for any i I would like to access the value of 

coeff2(index(i,1),index(i,2),index(i,3),index(i,4),index(i,5))

as quickly as possible. Being able to do this fast is what I mean by "efficient".

Since the indices in [something] are at most 2^10 I am currently defining coeff2 as follows:

REAL(8), DIMENSION(2**10,2**10,2**10,2**10,2**10) :: coeff2

but this is too wasteful of memory specially since I need to increase the number of dimensions, now 5, to the order of 100 and most elements of this array are equal to 0. So, another measure of efficiency that is relevant to me is that the memory necessary to store coeff2 should not explode as I increase the number of dimensions.

mzp
  • 181
  • 9
  • Questions asking us to **recommend or find a book, tool, software library, tutorial or other off-site resource** are off-topic for Stack Overflow as they tend to attract opinionated answers and spam. Instead, describe the problem and what has been done so far to solve it. – Vladimir F Героям слава May 24 '18 at 16:14
  • Even if there is a non-opinionated answer, you still have to say what kind of operations you have in mind and in which sense they should be optimized. What are *"the"* simple operations you want? Should they be optimized for speed, for storage size or for what? – Vladimir F Героям слава May 24 '18 at 16:15
  • @VladimirF I attempted to clarify the points you raised. Thank you for your help. – mzp May 24 '18 at 17:09
  • 1
    Note that using a magic constant such as 8 in `real(8)` is ugly and not portable. – Vladimir F Героям слава May 24 '18 at 20:50

1 Answers1

2

Well, It's still not totally clear to me the nature of your data and the way you want to use it.

If what you need is indexed data, whose indices are not consecutive, Sparse matrix can be an answer, and there are many solutions already implemented over the internet (as shown in the link you provided). But maybe it would be overkill for what I think you are trying to do. Maybe a simple datatype could serve your purpose, like this:

program indexed_values
  implicit none

  type :: indexed
    integer :: index
    real(8) :: value
  end type

  integer, parameter :: n_coeffs = 1000
  integer, parameter :: n_indices = 5
  integer :: i

  real(8), dimension(n_coeffs) :: coeff1
  integer, dimension(n_coeffs, n_indices) :: index
  type(indexed), dimension(n_coeffs, n_indices) :: coeff2
  type(indexed) :: var

  do i = 1, n_coeffs
    index(i, :) = [1, 2, 4, 16, 32] * i ! your calc here
    coeff1(i) = real(i * 3, 8) ! more calc here
    coeff2(i, :)%index = index(i, :)
    coeff2(i, :)%value = coeff1(i)
  end do

  ! that's how you fetch the indices and values by stored position
  var = coeff2(500, 2)
  print*, var%index, var%value ! outputs: 1000   1500.0

  ! that's how you fetch a value by its index
  print*, fetch_by_index(coeff2(500, :), 1000) ! outputs: 1500.0

contains

  real(8) function fetch_by_index(indexed_pairs, index)
    type(indexed), dimension(:) :: indexed_pairs
    integer, intent(in) :: index
    integer :: i

    do i=1, size(indexed_pairs)
      if(index == indexed_pairs(i)%index) then
        fetch_by_index = indexed_pairs(i)%value
        return
      end if
    end do
    stop "No value stored for this index"
  end
end

The provided function for fetching values by its indices could be improved if your indices will be alwyas stored in ascending order (no need to traverse the whole list to fail). Moreover, if you will assing a constant result of coeff1 to all the indices at each row, you could do even better and just not having a coeff2 array at all, just have coeff1 for values and index for the indices, and correlate them by position.

Rodrigo Rodrigues
  • 7,545
  • 1
  • 24
  • 36
  • Thank you very much for your answer, it is very very helpful! This is almost exactly what I need, but I realize now I did not explain well enough one key part. I need a function similar to the one you wrote down, but that would take an index and output the associated coefficient. For example `function([1, 2, 4, 16, 32])` should output `3.0` in your example. Is there a way to adjust the function above to be able to do this? – mzp May 24 '18 at 22:14
  • The key is that I need to get the value in `coeff2` without knowing `i` but only the 5 dimensional index (and without having to go through all the `i`'s to find the associated 5 dimensional index). – mzp May 24 '18 at 22:38
  • Again, I got lost in what you need. The only way you would't need to go through all the list to find the location of a specific index list if knowing its location (like, 1 for [1, 2, 4, 16, 32] or 500 for [500, 1000, 2000, 8000, 16000]). And, if you already know the location, you can use the location itself to fetch the value. – Rodrigo Rodrigues May 25 '18 at 04:02
  • I see. Notice that if I define the `coeff2` as in the question, it is possible to obtain the coefficient form the 5 dimensional array. Though, even if I define `coeff2` only in the relevant indexes this solution would not work if the dimension of the index list is more than 7. – mzp May 25 '18 at 08:40
  • I accepted the answer because I think it answers the question as I wrote it very well actually. But I attempted to ask it in a different way [here](https://stackoverflow.com/q/50534287/3553449), in case you are interested. – mzp May 25 '18 at 19:02