The product terms.

By operations count and direct measurement about $95$% of the cpu time required to execute this wavelet-Galerkin algorithm is spent in evaluating the product terms.

Consider the product terms for the three term connection coefficients. These would appear in the acoustic wave equation with $c^2=1$.

$$\Theta(p,q)=\Omega_{j,l}^{101}\Psi(l+p,m+q) K(j+p,n+q)\Omega_{n,m}^{000}+$$

$$\Omega_{j,l}^{000}\Psi(l+p,m+q) K(j+p,n+q\Omega_{n,m}^{101}$$

Let $\Omega^{000}$ and $\Omega^{101}$ be the real $(m,m)$ matricies of coefficients that depend on the wavelet basis. Let $\hat K(p,q)$ and $\hat \Psi(p,q)$ be the $(m,m)$ matrices of coefficients centered on the $(p,q)$ element. Then

$$\Theta(p,q) = \Omega^{101}\hat \Psi(p,q).^*\hat K(p,q) \Omega^{000}$$

$$+\Omega^{000}\hat \Psi(p,q).^*\hat K(p,q) \Omega^{101}$$

where $.^*$ is the term by term product and sum for matrices.

The unoptimized evaluation of these expressions requires

$$(8m^3 + 4m^2)N^2$$

floating point operations for a complete evaluation of the product terms.

If we know $\Omega\hat \Psi(p,q)$ then

$$\Omega\hat \Psi(p,q+1)=\Omega\{\hat \Psi(p,q) \hat D + \hat S(p,q+1)\}.$$

$\hat D$ is the matrix that shifts the columns to the left by one column and assigns the null column to last one and $\hat S(p,q+1)$ is the matrix with null columns except the last one which is the last column of $\hat \Psi(p,q+1)$. This reduces the operation count to

$$8m^2N^2 + 8m^3N.$$

A similar optimization exists for the four term connection coefficients.

This type of procedure is suitable for parallel implementations since the data structures are, more or less, local. However, even with this optimization the four term evaluations are expensive. Further optimizations take into account the structure of the connection coefficients applied to piecewise constant $k$ and $c$.