Kalman filter in BLAS/LAPACK Fortran

kalman.f90

subroutine f(mu, Sigma, H, INFO, R, Sigmavar_2, data, muvar_2, k, n)
implicit none

integer, intent(in) :: k
integer, intent(in) :: n
real*8, intent(in) :: Sigma(n, n)        !  Sigma
real*8, intent(in) :: H(k, n)            !  H
real*8, intent(in) :: mu(n)              !  mu
real*8, intent(in) :: R(k, k)            !  R, H*Sigma*H' + R
real*8, intent(in) :: data(k)            !  (H*Sigma*H' + R)^-1*((-1)*data + H*mu), data, (-1)*   data + H*mu
integer, intent(out) :: INFO             !  INFO
real*8, intent(out) :: muvar_2(n)        !  mu, Sigma*H'*(H*Sigma*H' + R)^-1*((-1)*data + H*  mu) + mu
real*8, intent(out) :: Sigmavar_2(n, n)  !  Sigma, (-1)*Sigma*H'*(H*Sigma*H' + R)^-1*H* Sigma + Sigma
real*8 :: var_17(n, k)                   !  Sigma*H', 0
real*8 :: Hvar_2(k, n)                   !  (H*Sigma*H' + R)^-1*H, H
real*8 :: var_11(n)                      !  0, H'*(H*Sigma*H' + R)^-1*((-1)*data + H*mu)
real*8 :: var_19(n, n)                   !  0, H'*(H*Sigma*H' + R)^-1*H
real*8 :: var_5(n, n)                    !  0
real*8 :: var_20(n, n)                   !  H'*(H*Sigma*H' + R)^-1*H*Sigma, 0

call dcopy(n**2, var_5, 1, var_20, 1)
call dsymm('L', 'U', n, k, 1, Sigma, n, H, k, 0, var_17, n)
call dgemm('N', 'N', k, k, n, 1, H, k, var_17, n, 1, R, k)
call dcopy(n**2, var_5, 1, var_19, 1)
call dcopy(n, mu, 1, muvar_2, 1)
call dcopy(n**2, Sigma, 1, Sigmavar_2, 1)
call dcopy(k*n, H, 1, Hvar_2, 1)
call dgemm('N', 'N', k, 1, n, 1, H, k, mu, n, -1, data, k)
call dposv('U', k, n, R, k, Hvar_2, k, INFO)
call dposv('U', k, 1, R, k, data, k, INFO)
call dgemm('N', 'N', n, n, k, 1, H, k, Hvar_2, k, 0, var_19, n)
call dgemm('N', 'N', n, 1, k, 1, H, k, data, k, 0, var_11, n)
call dsymm('L', 'U', n, n, 1, var_19, n, Sigma, n, 0, var_20, n)
call dsymm('L', 'U', n, 1, 1, Sigma, n, var_11, n, 1, muvar_2, n)
call dsymm('L', 'U', n, n, -1, Sigmavar_2, n, var_20, n, 1, Sigmavar_2, n)

RETURN
END

Could I have more information about the main.f90?

Could you vertify this code please?I could not run this code properly.

I'm no expert by any means but (a) to use this code you'll need to have external dgemm etc. somewhere (maybe main.f90?) so that Fortran knows where to find the symbols in BLAS/LAPACK. Also, I'm pretty sure the dimensions for the dsymm call on line 22 are wrong - should be H^T (n,k) not H (k,n) and dsymm doesn't offer a way to set a transpose option for its second argument. I also don't know if the last line (35) is correct, pointer aliasing is not permitted in Fortran (https://en.wikipedia.org/wiki/Pointer_aliasing): Sigmavar_2 is used twice and code for an implementation of dsymm (http://www.netlib.org/lapack/explore-html/d8/db0/dsymm_8f_source.html) where the problem lines are 300 and 302 when i=j; the right answer is highly dependent on how those expressions are compiled.

In this blogpost the author uses sympy to symbolically derive a Kalman filter implementation using BLAS and LAPACK. It might be worth comparing the two (note that the symbolic route uses some unneccesary axpy calls that could be both done with a call to gemm.

Edit: just noticed that mrocklin is the author of the blogpost 😄