I am trying to optimize the execution time for this particular piece of fortran code.

Details: i_gc is a (ngpts, 3) array of containing (i,j,k) indices for each grid point. This is a subset of the full NX x NY x NZ grid. The subset is connected in the sense that density of sampling varies, but the points aren't scattered randomly throughout the domain.

Essentially I have a non-local operator on the grid which is a tensor product of 3 one dimensional operators. If I had the data on the full 3d grid it would be pretty easy to write this code without all the if statements in the inner loop. Also, I would like to avoid storing a set {i},{j},{k} indices for each grid point, but maybe this is the way to go?

I am wondering if there is a way to do much better than this? I do have some freedom to reorder the points if needed. Maybe I can exploit that somehow? I suppose I could try something like a Hilbert curve, but I'm not sure how well that could apply to my case where I have an oddly shaped subset of a regular grid.

edit: mathematically the operation is like this

$Y_{ijk} = \sum^{N^x}_{q=1} M^x_{qi} X_{qjk} + \sum^{N^y}_{q=1} M^y_{qj} X_{iqk} + \sum^{N^z}_{q=1} M^z_{qk} X_{ijq} $

where the sums are only over points in the subset (i_gc)

subroutine apply_lf4(ngpts, i_gc, mx,my,mz,x,y)
  implicit none
  integer, intent(in) :: ngpts
  integer, intent(in), dimension(:,:) :: i_gc
  real(8), intent(in), dimension(:,:) :: mx,my,mz
  real(8), intent(in), dimension(ngpts) :: x
  real(8), intent(out), dimension(ngpts) :: y

  integer :: ig,jg,i1,i2,i3,j1,j2,j3


  do ig=1, ngpts
     i1=i_gc(ig,1); i2=i_gc(ig,2); i3=i_gc(ig,3)
    inner: do jg=1, ngpts
        if (j1==i1) then
           if (j2==i2) then
              j3 = i_gc(jg,3)
              y(ig) = y(ig) + mz(j3,i3)*x(jg)
              cycle inner
           end if
           if (j3==i3) then
              y(ig) = y(ig) + my(j2,i2)*x(jg)
              cycle inner
           end if
           j2 = i_gc(jg,2)
           if (j2 == i2) then
              j3 = i_gc(jg,3)
              if (j3 == i3) then
                 y(ig) = y(ig) + mx(j1,i1)*x(jg)
                 cycle inner
              end if
           end if
        end if                        
     end do inner
  end do             

end subroutine apply_lf4
  • $\begingroup$ How large is ngpts relative to NX, NY, NZ? $\endgroup$ – Jesse Chan May 6 '15 at 0:21
  • $\begingroup$ What is this computing in mathematical terms? I'm a bit lost. $\endgroup$ – Bill Barth May 6 '15 at 0:22
  • $\begingroup$ @JesseChan a typical NX=NY=NZ=400, Ngpts=10^6 $\endgroup$ – user1984528 May 6 '15 at 2:15
  • $\begingroup$ @BillBarth I have updated the question to include some mathematical terms. Hopefully that is clear enough? $\endgroup$ – user1984528 May 6 '15 at 2:22

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Browse other questions tagged or ask your own question.