Serial Block inverse
We split the matrix into four sub-blocks:
| A B |
| C D |
where A is m x m and D is (n-m) x (n-m). The inverse can then be written as
| A B |^-1 == | M N |
| C D | | O P |
With
M = [A-B(D^-1)C]^-1
N = -MB(D^-1)
O = -(D^-1)CM
P = (D^-1) + (D^-1)CMB(D^-1)
In practice if we define
T = -B(D^-1)
S = -(D^-1)C
we can then write
M = [A + TC]^-1
N = MT
O = SM
P = D^-1 + SN
We can therefore find the inverse of our matrix using two inversions of matrices of size ~ m x m ~ (n/2)x(n/2) two matrix additions and six matrix multiplies.
The performance of this approach is dependent on an efficient matrix multiplication implementation. Most compilers will offer decent implementations and the option to use an external library such as blas.
Type | Intent | Optional | Attributes | Name | ||
---|---|---|---|---|---|---|
complex, | intent(inout), | dimension(n, n), target | :: | a | ||
integer, | intent(in) | :: | n | |||
integer, | intent(in), | optional | :: | m_in |
recursive subroutine invert_block(a,n,m_in)
implicit none
integer, intent(in):: n
complex, dimension(n, n), intent(in out), target :: a
integer, intent(in), optional :: m_in
integer :: m
complex, dimension(:, :), allocatable :: tmp
complex, dimension(:, :), pointer :: aa, bb, cc, dd
! If the matrix is small then use direct inversion
if (n <= 4) then
call invert_gauss_jordan(a, n)
return
end if
! Ensure the sub-matrix size is set. We split the length n into
! two regions, 1:m and m+1:n.
if(present(m_in)) then
m = m_in
else
m = max(n/2, 1)
end if
! Here we use pointers to give more useful names to the sub-blocks of a
! An associate block could be another option, but currently this leads to
! errors in the result.
aa => a(1:m, 1:m)
bb => a(1:m, m+1:n)
cc => a(m+1:n, 1:m)
dd => a(m+1:n, m+1:n)
! First calculate D^-1 and store in place
call invert_block(dd, n-m)
!Calculate T = -B(D^-1)
tmp = -matmul(bb, dd)
! Calculate [A+TC]
aa = aa + matmul(tmp, cc)
! Invert [A+TC] and store in place. M is now done
call invert_block(aa, m)
! Calculate N and store.
bb = matmul(aa, tmp)
! Calculate S = -(D^-1)C
tmp = -matmul(dd, cc)
! Calculate O
cc = matmul(tmp, aa)
! Calculate P
dd = dd + matmul(tmp, bb)
end subroutine invert_block