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 |
Sign up for free
to join this conversation on GitHub.
Already have an account?
Sign in to comment
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
axpycalls that could be both done with a call togemm.Edit: just noticed that mrocklin is the author of the blogpost 😄