I am working on a problem where I have a symmetric function afnc(x,x'), whose values I need to compute for a rectangular grid and store in a suitable matrix a(i,j).

Currently, I am using a normal square matrix. I evaluate half of the matrix and then I populate the other half using the property a(i,j)=a(j,i), as shown below:

```
program main
implicit none
integer, parameter ::nc = 1000
real :: x(nc), a(nc, nc)
integer :: i, j
! Make a grid
do i = 1, nc
x(i) = 1.*i
end do
! Method 1
do j = 1, nc
do i = j, nc
a(i, j) = afnc(x(i), x(j))
end do
end do
do j = 2, nc
do i = 1, j - 1
a(i, j) = a(j, i)
end do
end do
contains
pure real function afnc(x1, x2) result(res)
!! "Very complex" kernel. Here made simple on purpose.
real, intent(in) :: x1, x2
res = x1*x2
end function
end program main
```

From a computational perspective, IMO, this seems already not-too-bad, because there are no redundant function calls and retrieving the values from the matrix is a â€śdirectâ€ť operation for all index combinations. On the flip side, I am wasting half of the matrix (although we are only talking about a small ~(1000^2)/2 memory waste).

So, I was curious how I could improve this implementation without actually making things worse. For example, I could imagine storing the info in a vector v(k) (by flattening the â€śhalf-arrayâ€ť), but I would then have to wrap it into a function to compute the index k based on the array indexes (i,j), which would involve an if-else construct on the indexes, making it slower to access.

`lapack`

probably has a nice solution to this (?), but as far as I know `fpm`

cannot yet install it automatically, which somehow makes me a bit reluctant to go down that road â€¦

Any thoughts are very welcome!