Skip to content

Instantly share code, notes, and snippets.

@pozdneev

pozdneev/chol-no-pivoting.f90

Created March 11, 2011 22:31
Show Gist options
  • Save pozdneev/866702 to your computer and use it in GitHub Desktop.
Save pozdneev/866702 to your computer and use it in GitHub Desktop.
Download ZIP
Cholesky decomposition without pivoting: straightforward formulas, Fortran 90/95 features and BLAS routines
Raw
chol-no-pivoting.f90
! Cholesky decomposition without pivoting: straightforward formulas,
! Fortran 90/95 features 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 matrix
integer, parameter :: n = 5
! 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 Fortran 90/95 features
call Generate_Data()
call Fortran9x_Decomposition()
call Fortran9x_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 :: i, j, k
elimination_start = omp_get_wtime()
do k = 1, n
a(k, k) = sqrt(a(k, k))
do i = k 1, n
a(i, k) = a(i, k) / a(k, k)
end do
do j = k 1, n
do i = j, n
a(i, j) = a(i, j) - a(i, k) * a(j, k)
end do
end do
end do
elimination_finish = omp_get_wtime()
end subroutine Straightforward_Decomposition
!!---------------------------------------------------------------------------
subroutine Fortran9x_Decomposition()
integer :: k, j, i
elimination_start = omp_get_wtime()
do k = 1, n
a(k, k) = sqrt(a(k, k))
a(k 1: n, k) = a(k 1: n, k) / a(k, k)
do j = k 1, n
do i = j, n
a(i, j) = a(i, j) - a(i, k) * a(j, k)
end do
end do
end do
elimination_finish = omp_get_wtime()
end subroutine Fortran9x_Decomposition
!!---------------------------------------------------------------------------
subroutine BLAS_Decomposition()
integer :: k
elimination_start = omp_get_wtime()
do k = 1, n
a(k, k) = sqrt(a(k, k))
! x = a*x
call sscal(n-k, 1.0 / a(k, k), a(k 1, k), 1)
! A := alpha*x*x' A
call ssyr('L', n-k, -1.0, a(k 1, k), 1, a(k 1, k 1), 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 Fortran9x_Backsubstition()
integer :: i
backsubstition_start = omp_get_wtime()
! L x = f => x = L \ f
do i = 1, n
y(i) = y(i) - dot_product(a(i, 1: i-1), y(1: i-1))
y(i) = y(i) / a(i, i)
end do
! L^T y = x => y = L^T \ x
do i = n, 1, -1
y(i) = y(i) - dot_product(a(i 1: n, i), y(i 1: n))
y(i) = y(i) / a(i, i)
end do
backsubstition_finish = omp_get_wtime()
end subroutine Fortran9x_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
!!---------------------------------------------------------------------------
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment