invert_block Subroutine

public recursive subroutine invert_block(a, n, m_in)

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.

Arguments

Type IntentOptional Attributes Name
complex, intent(inout), dimension(n, n), target :: a
integer, intent(in) :: n
integer, intent(in), optional :: m_in

Contents

Source Code


Source Code

  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