Talk:Algorithm Implementation/Linear Algebra/Tridiagonal matrix algorithm

Are we sure that the fortran implementation gives the right answers? It seems to be wrong, there is a minimal example showing this below:

code:

program tdmatester implicit none ! the extent of the arrays integer :: n=6 call runner(n) end program tdmatester

subroutine runner(n) integer :: n, i   real, dimension(n,n) :: matrix real, dimension(n) :: lower_digonal, upper_diagonal real, dimension(n) :: main_diagonal real, dimension(n) :: rhs, rhs_2 real, dimension(n) :: solution

do i = 1, n       lower_digonal(i) = 1.0 + rand ! these rands can be commented out main_diagonal(i) = 2.0 + rand ! the same problem remains upper_diagonal(i) = 3.0 + rand rhs(i) = 4.0 + rand end do

write(*,*) "Original right hand side:" write(*,*) rhs call make_matrix(lower_digonal, main_diagonal, upper_diagonal, n, matrix)

write(*,*) "Matrix:" do i = 1, n       print *, matrix(i,:) end do

call tdma(n, lower_digonal, main_diagonal, upper_diagonal, rhs, solution)

rhs_2 = matmul(matrix, solution) write(*,*) "Reconstructed right hand side matmul(matrix, solution):" write(*,*) rhs_2

end subroutine

subroutine tdma(K,a,b,c,d,x) implicit none integer, intent(in) :: K   real, intent(in) :: a(K), c(K) real, intent(inout), dimension(K) :: b, d   real, intent(out) :: x(K) ! --- Local variables --- integer :: i   real :: q        !  --- Elimination --- do i = 2,K q = a(i)/b(i - 1) b(i) = b(i) - c(i - 1)*q d(i) = d(i) - d(i - 1)*q end do       ! --- Backsubstitution --- q = d(K)/b(K) x(K) = q   do i = K - 1,1,-1 q = (d(i) - c(i)*q)/b(i) x(i) = q   end do    return end SUBROUTINE tdma

subroutine make_matrix(lower_digonal,main_diagonal,upper_diagonal, n, matrix) !Constructs a tri diagonl matrix writes the output to matrix implicit none integer, intent(in) :: n   real, intent(out), dimension(n, n) :: matrix real, intent(in) :: lower_digonal(n), upper_diagonal(n) real, intent(in) :: main_diagonal(n)

integer :: i

matrix = 0.0 do i = 1, n       matrix(i,i) = main_diagonal(i) matrix(i, i+1) = upper_diagonal(i) if (i>1) then matrix(i, i-1) = lower_digonal(i) end if   end do

end subroutine

Console output (I've added * to show the results of interest):

Original right hand side: 4.45865011      4.67886448       4.51941633       *4.52970028*       4.06684208       4.93043613 Matrix: 2.13153768      3.75560522       0.00000000       0.00000000       0.00000000       0.00000000   1.53276706       2.21895909       3.04704452       0.00000000       0.00000000       0.00000000   0.00000000       1.67929626       2.93469286       3.38350201       0.00000000       0.00000000   0.00000000       0.00000000       1.83096528       2.03457189       3.05346155       0.00000000   0.00000000       0.00000000       0.00000000       1.67114925       2.00769806       3.38341546   0.00000000       0.00000000       0.00000000       0.00000000       1.41748595       2.68677258 Reconstructed right hand side matmul(matrix, solution): 4.45864964      4.67886400       4.51941586       *3.58897662*       4.06684208       4.93043613

You can see that the right hand side of the equation is not reconstructed correctly, so the found solution is not correct.

Demonstrating the Tridiagonal matrix algorithm in MAST?
It is helpful to see the examples of the Tridiagonal matrix algorithm in Python and C++, but what about modeling languages such as MAST or other HDL? 198.208.47.88 (discuss) 12:59, 6 October 2022 (UTC)