Instantly share code, notes, and snippets.

# shane5ul/testCyclicTridiag.f90 Created Feb 17, 2014

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