program main

    character :: uplo, trans, diag
    integer :: n
    double precision, allocatable :: ah(:,:), al(:,:)
    integer :: lda
    double precision, allocatable :: xh(:), xl(:)
    integer :: incx

    integer :: i, j

    uplo = 'u'
    trans = 'n'
    diag = 'u'
    n = 3
    incx = 1

    lda = max(1, n)
    allocate(ah(lda,n))
    allocate(al(lda,n))
    allocate(xh(1+(n-1)*abs(incx)))
    allocate(xl(1+(n-1)*abs(incx)))

    call random_number(ah)
    ah(1,2) = 0.1d0
    ah(1,3) = 0.2d0
    ah(2,3) = 0.3d0
    al = 0d0

    call random_number(xh)
    xl = 0d0

    print *, "INPUT"
    call print_matrix("a", ah, al)
    call print_vector("x", xh, xl)

    call ddtrsv(uplo, trans, diag, n, ah, al, lda, xh, xl, incx)

    print *, "OUTPUT"
    call print_vector("x", xh, xl)

    contains

    subroutine print_matrix(vname, ah, al)

        character(*), intent(in) :: vname
        double precision, intent(in) :: ah(:,:), al(:,:)
        integer :: i, j

        do i = 1, size(ah, 1)
            do j = 1, size(ah, 2)
                print '(a,"(",i0,",",i0,") =",2e23.15)', &
                    & vname, i, j, ah(i,j), al(i,j)
            end do
        end do

    end subroutine

    subroutine print_vector(vname, xh, xl)

        character(*), intent(in) :: vname
        double precision, intent(in) :: xh(:), xl(:)
        integer :: i

        do i = 1, size(xh, 1)
            print '(a,"(",i0,") =",2e23.15)', vname, i, xh(i), xl(i)
        end do

    end subroutine

end program