Program to test solution of cyclic tridiagonal matrix using the Sherman-Morrison formula

testCyclicTridiag.f90

!
! 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