Last active
August 27, 2023 20:21
-
-
Save mrocklin/5144149 to your computer and use it in GitHub Desktop.
Kalman filter in BLAS/LAPACK Fortran
This file contains hidden or bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
| 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 |
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 😄
Sign up for free
to join this conversation on GitHub.
Already have an account?
Sign in to comment
I'm no expert by any means but (a) to use this code you'll need to have
external dgemmetc. 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.