Automatic Implementation of an 8 x 8 DCT

The DCT (discrete cosine transform) is defined as the matrix

where
We want to show, how a fast implementation of this transform can be generated using AREP and SPL (formerly TPL) according to the following diagram.
The first two steps are performed with the GAP share package AREP. GAP is a computer algebra system for computational group theory, AREP. is a package for symbolic computation with group representations. The third step uses SPL (signal processing language) and the SPL compiler. SPL is a domain specific language for linear computations and allows to express fast algorithms given as sparse matrix products. The SPL compiler translates an SPL program into an efficient FORTRAN program for a matrix-vector multplication with the represented matrix.

Step 1: Computing a sparse matrix factorization

First we factorize the matrix into a product of sparse matrices using the function MatrixDecompositionByPermIrredSymmetry from AREP. The matrix is analyzed for symmetry and decomposed (explanation of the procedure). We start with the transpose of the DCT, usually called DCT, type III, and transpose the result.

gap> M := DCT_III(8);
gap> A := MatrixDecompositionByPermIrredSymmetry(M);
gap> A := TransposedAMat(A);
The function MatrixDecompositionByPermIrredSymmetry returns an AMat, a recursive data structure for representing sparse structured matrix products. Click here to see, what A looks like.

Step 2: Translating into an SPL program

Next we want to create a SPL program from the AMat A. First, A is translated into an SPL object, an interim data structure to represent SPL programs within GAP. Second, the SPL object is exported to a file as an SPL program.

gap> T := SPLAMat(A);
gap> T := PrepareForExportingSPL(T);
gap> ExportSPL("tmp.spl", T);
The file tmp.spl now contains an SPL program representing the sparse matrix factorization computed above.

Step 3: Compiling to a FORTRAN program

The SPL compiler translates the SPL program into a FORTRAN progam performing different optimizations. The option "-R" creates a program for real input vectors.

spl -R tmp.spl > tmp.f
See, how tmp.f looks like.

Step 4: Compiling the FORTRAN program

Any standard FORTRAN compiler can be used to compile tmp.f:

f77 -c tmp.f