Cholesky partitioned decomposition without pivoting: straightforward formulas and BLAS routines

chol-partitioned.f90

! Cholesky partitioned decomposition without pivoting: straightforward
! formulas and BLAS routines
!!---------------------------------------------------------------------------
! gfortran -ffree-form -Wall -Wextra -lblas -fopenmp
!!---------------------------------------------------------------------------
program chol
!!---------------------------------------------------------------------------
    use omp_lib, only: omp_get_wtime

    implicit none

    ! Matrix of coefficients; the one is "random" symmetric positive-definite
    ! Only the lower triangular part is significant
    real, dimension(:, :), allocatable :: A

    ! "Analytical" solution; the one is filled by random_number()
    real, dimension(:), allocatable :: u

    ! 1. Right-hand side (RHS); the one is calculated as f = A * u
    ! 2. Numerical solution (NS) of the equation A y = f
    ! #. RHS is overwritten by NS
    real, dimension(:), allocatable :: y

    ! Size of block
    integer, parameter :: b = 4

    ! Size of matrix; must be a multiple of b
    integer, parameter :: n = b * 3

    ! Time marks
    real(kind(0.d0)) :: elimination_start, elimination_finish
    real(kind(0.d0)) :: backsubstition_start, backsubstition_finish

    ! Allocate memory
    allocate(A(1: n, 1: n))
    allocate(u(1: n))
    allocate(y(1: n))

    ! Algorithm uses straightforward formulas
    call Generate_Data()
    call Straightforward_Decomposition()
    call Straightforward_Backsubstition()
    call Print_Norms()
    call Print_Times()

    ! Algorithm uses BLAS
    call Generate_Data()
    call BLAS_Decomposition()
    call BLAS_Backsubstition()
    call Print_Norms()
    call Print_Times()

    ! Free memory
    deallocate(A)
    deallocate(u)
    deallocate(y)
!!---------------------------------------------------------------------------
contains
!!---------------------------------------------------------------------------
subroutine Print_Norms()
    write (*, *) "Norms:", maxval(abs(u)), maxval(abs(y - u))
end subroutine Print_Norms
!!---------------------------------------------------------------------------
subroutine Print_Times()
    write (*, *) "Times:", &
            elimination_finish - elimination_start, &
            backsubstition_finish - backsubstition_start
end subroutine Print_Times
!!---------------------------------------------------------------------------
! This version is a simplified modification of
! http://gcc.gnu.org/onlinedocs/gfortran/RANDOM_005fSEED.html
subroutine Init_Random_Seed()
    integer :: i, n
    integer, dimension(:), allocatable :: seed

    call random_seed(size = n)
    allocate(seed(n))

    seed = 37 * (/ (i - 1, i = 1, n) /)
    call random_seed(put = seed)

    deallocate(seed)
end subroutine Init_Random_Seed
!!---------------------------------------------------------------------------
subroutine Generate_Data()
    call Init_Random_Seed()

    call random_number(A)
    A = matmul(A, transpose(A))

    call random_number(u)

    y = matmul(A, u)
end subroutine Generate_Data
!!---------------------------------------------------------------------------
subroutine Straightforward_Decomposition()
    integer :: ib
    integer :: i, j, k
    integer :: p1b, p1e, p2b, p2e

    elimination_start = omp_get_wtime()

    do ib = 0, n/b - 1
        p1b = ib * b + 1
        p1e = (ib + 1) * b
        p2b = p1e + 1
        p2e = n

        ! L_{11} L_{11}^T = A_{11}
        do k = p1b, p1e
            a(k, k) = sqrt(a(k, k))

            do i = k+1, p1e
                a(i, k) = a(i, k) / a(k, k)
            end do

            do j = k+1, p1e
                do i = j, p1e
                    a(i, j) = a(i, j) - a(i, k) * a(j, k)
                end do
            end do
        end do

        ! L_{21} = A_{21} L_{11}^T^{-1}
        do k = p1b, p1e
            do i = p2b, p2e
                do j = p1b, k-1
                    a(i, k) = a(i, k) - a(i, j) * a(k, j)
                end do

                a(i, k) = a(i, k) / a(k, k)
            end do
        end do

        ! A_{22} -= L_{21} L_{21}^T
        do j = p2b, p2e
            do k = p1b, p1e
                do i = j, p2e
                    a(i, j) = a(i, j) - a(i, k) * a(j, k)
                end do
            end do
        end do
    end do

    elimination_finish = omp_get_wtime()
end subroutine Straightforward_Decomposition
!!---------------------------------------------------------------------------
subroutine BLAS_Decomposition()
    integer :: ib
    integer :: k
    integer :: p1b, p1e, p2b, p2e

    elimination_start = omp_get_wtime()

    do ib = 0, n/b - 1
        p1b = ib * b + 1
        p1e = (ib + 1) * b
        p2b = p1e + 1
        p2e = n

        ! L_{11} L_{11}^T = A_{11}
        do k = p1b, p1e
            a(k, k) = sqrt(a(k, k))

            ! x = a*x
            call sscal(p1e-k, 1.0 / a(k, k), a(k+1, k), 1)

            ! A := alpha*x*x' + A
            call ssyr('L', p1e-k, -1.0, a(k+1, k), 1, a(k+1, k+1), n)
        end do

        ! L_{21} L_{11}^T = A_{21}  =>  L_{21} = A_{21} L_{11}^T^{-1}
        ! X*op(A) = alpha*B
        call strsm('R', 'L', 'T', 'N', p2e-p2b+1, b, 1.0, &
                a(p1b, p1b), n, a(p2b, p1b), n)

        ! A_{22} -= L_{21} L_{21}^T
        ! C := alpha*A*A' + beta*C,
        call ssyrk('L', 'N', p2e-p2b+1, b, -1.0, a(p2b, p1b), n, 1.0, &
                a(p2b, p2b), n)
    end do

    elimination_finish = omp_get_wtime()
end subroutine BLAS_Decomposition
!!---------------------------------------------------------------------------
subroutine Straightforward_Backsubstition()
    integer :: i, j

    backsubstition_start = omp_get_wtime()

    ! L x = f  =>  x = L \ f
    do i = 1, n
        do j = 1, i-1
            y(i) = y(i) - a(i, j) * y(j)
        end do

        y(i) = y(i) / a(i, i)
    end do

    ! L^T y = x  =>  y = L^T \ x
    do i = n, 1, -1
        do j = i+1, n
            y(i) = y(i) - a(j, i) * y(j)
        end do

        y(i) = y(i) / a(i, i)
    end do

    backsubstition_finish = omp_get_wtime()
end subroutine Straightforward_Backsubstition
!!---------------------------------------------------------------------------
subroutine BLAS_Backsubstition()
    backsubstition_start = omp_get_wtime()

    ! L x = f  =>  x = L \ f
    ! ... A * x = b
    call strsv('L', 'N', 'N', n, a, n, y, 1)

    ! L^T y = x  =>  y = L^T \ x
    ! ... A * x = b
    call strsv('L', 'T', 'N', n, a, n, y, 1)

    backsubstition_finish = omp_get_wtime()
end subroutine BLAS_Backsubstition
!!---------------------------------------------------------------------------
end program chol
!!---------------------------------------------------------------------------