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!