Accelerating seismic computations using customized number representations on FPGAs

Number Representation Exploration

Based on all the above modules, we can now perform exploration of different number representations for the FPGA implementation of a specific piece of Fortran code.

The current tools support two different number representations: fixed-point and floating-point numbers (the value simulator for LNS is still in progress). For all the different number formats, the users can also specify arbitrary bit-widths for each different variable.

There are usually a large number of different variables involved in one circuit design. In our previous work, we usually apply heuristic algorithms, such as ASA (Ingber, 2004), to find out a close-to-optimal set of multiple values for the bit-widths of different variables. The heuristic algorithms may require millions of test runs to check whether a specific set of values meet the constraints or not. This is acceptable when the test run is only a simple error function and can be processed in nanoseconds. In our seismic processing application, depending on the problem size, it takes half an hour to several days to run one test set. Thus, heuristic algorithms become impractical.

A simple and straightforward method to solve the problem is to use uniform bit-width over all the different variables to either iterate over a set of possible values or use a binary search algorithm to jump to an appropriate bit-width value.

Based on the range information and the internal behavior of the program, we can also try to divide the variables in the target Fortran code into several different groups, and assign a different uniform bit-width for each different group. For instance, in the `complex exponential' function, there is a clear boundary that the first half performs square, square root and division operations to calculate an integer result, and the second half uses the integer result as a table index, and performs sine, cosine and complex multiplications to get the final result. Thus, in the hardware circuit design, we divide the variables into two groups based on which half they belong to. Furthermore, in the second half of the function, some of the variables are trigonometric values in the range of $[-1, 1]$, while the other variables represent the seismic image data and scale up to $10^6$. Thus they can be further divided into two groups and assigned bit-widths separately.

! generation of table step%ctable

do i=1,size(step%ctable)
    k=ko*step%dstep*dsr%phase(i)
    step%ctable(i)=dsr%amp(i)*cmplx(cos(k),sin(k))
end do

! the core part of function wei_wem

do i4=1,size(wfld,4)
    do i3=1,size(wfld,3)
        do i2=1,size(wfld,2)
            do i1=1,size(wfld,1)

                k = sqrt(step%kx(i1,i3)**2 + step%ky(i2,i4)**2)
                itable =max(1, min(int(1 + k/ko / dsr%d) , dsr%n))
                wfld(i1,i2,i3,i4,i5)=wfld(i1,i2,i3,i4,i5)*step%ctable(itable)

            end do
        end do
    end do
end do

Accelerating seismic computations using customized number representations on FPGAs