Create a gist now

Instantly share code, notes, and snippets.

What would you like to do?
Program to test solution of cyclic tridiagonal matrix using the Sherman-Morrison formula
!
! Program to test solution of cyclic tridiagonal matrix
!
program testCyclicTridiag
integer, parameter :: n=5
integer :: i, j
double precision :: a(5), b(5), c(5), d(5), x(5)
!
! make up a particular tridiagonal matrix; x is the solution;
! a, b, and c are the "tridiagonal" matrix elements; d is the RHS,
! and x is the solution vector
!
a = 0.
b = 0.
c = 0.
d = 0.
x = 0.
do i = 1, n
a(i) = 1
b(i) = -2
c(i) = 1
d(i) = dfloat(i)
enddo
a(1) = 2.0
c(n) = -1.0
call cyclicTriDiag(a, b, c, d, x)
! call TriDiag(a, b, c, d, x) ! can test the tridiag routine separately
write(*,*) "Solution x ="
write(*,*)
do i = 1, n
write(*,'(f8.4)') x(i)
enddo
write(*,*)
contains
!--------------------------------------------------------------------------------------------------
subroutine cyclicTriDiag(a, b, c, r, x)
! Solves for a vector u of length n the modified tridiagonal linear set
! m u = r, where a, b and c are the three main diagonals of matrix
! m(n,n), the other terms are 0. r is the right side vector.
!
! m(1,n) = a(1) (or beta from NR); and m(n,1) = c(n) (or alpha from NR)
!--------------------------------------------------------------------------------------------------
implicit none
double precision, intent(in) :: a(:), b(:), c(:), r(:)
double precision, intent(out) :: x(:)
!
! Local Variables
!
double precision, dimension(size(x)) :: bb, u, z
integer :: n, i
double precision :: fact, gamma
n = size(x)
gamma = -b(1) ! gamma can be arbit; -b1 is recommended
bb(1) = b(1) - gamma
bb(n) = b(n) - a(1) * c(n)/gamma
bb(2:n-1) = b(2:n-1)
call tridiag(a, bb, c, r, x) ! Solve Ax = r
u(1) = gamma ! Set up vector u
u(n) = c(n) !a(1)
u(2:n-1) = 0.d0
call tridiag(a, bb, c, u, z) ! Solve Az = u
!
! Form v.x/(1 + v.z)
!
fact = (x(1) + a(1)*x(n)/gamma)/(1.0 + z(1) + a(1) * z(n)/gamma)
x = x - fact * z ! get solution vector z
end subroutine cyclicTriDiag
!--------------------------------------------------------------------------------------------------
subroutine tridiag(a, b, c, r, u)
! Solves for a vector u of length n the tridiagonal linear set ! from numerical recipes
! m u = r, where a, b and c are the three main diagonals of matrix !
! m(n,n), the other terms are 0. r is the right side vector. !
!
! a(1) and c(n) are not used in the calculation
!--------------------------------------------------------------------------------------------------
implicit none
double precision, intent(in) :: a(:), b(:), c(:), r(:)
double precision, intent(out) :: u(:)
!
! Local Variables
!
double precision, dimension(size(b)) :: gam
integer :: n, j
double precision :: bet
n = size(a)
bet = b(1)
if (bet == 0.0) then
write(*,*) "TriDiag Error"
stop
endif
u(1) = r(1)/bet
do j = 2, n
gam(j) = c(j-1) / bet
bet = b(j) - a(j) * gam(j)
if (bet == 0.0) then
write(*,*) "TriDiag Error"
stop
endif
u(j) = ( r(j) - a(j) * u(j-1) )/bet
end do
do j = n-1, 1, -1
u(j) = u(j) - gam(j+1) * u(j+1)
end do
end subroutine tridiag
end program testCyclicTridiag
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment