What ChatGPT says about Fortran

AI
38

In Computing the mean: accuracy and speed @urbanjost asked,

I am wondering how far away AI is from being able to refactor netlib into modern code structure

If you want a modern Fortran procedure that solves a linear algebra problem, ChatGPT-4 will write a LAPACK wrapper for you. I asked it to

Write a self-contained Fortran to solve the set of linear equations A*x = b, where matrix A and vector b are given. Fill in the following code:

module solve_lapack_mod
implicit none
integer, parameter :: dp = kind(1.0d0)
contains
subroutine solve_eq(a, b, x, ierr)
real(kind=dp), intent(in)  :: a(:,:) ! (nobs, ncol) -- LHS matrix
real(kind=dp), intent(in)  :: b(:)   ! (nobs)       -- RHS vector
real(kind=dp), intent(out) :: x(:)   ! (ncol)       -- solution vector
integer      , intent(out) :: ierr   ! returned nonzero if error encountered
end subroutine solve_eq
end module solve_lapack_mod

It produced

module solve_lapack_mod
    use, intrinsic :: iso_fortran_env, only : dp => real64
    implicit none
    contains
    subroutine solve_eq(a, b, x, ierr)
        real(dp), intent(in)    :: a(:,:)  ! LHS matrix
        real(dp), intent(in)    :: b(:)    ! RHS vector
        real(dp), intent(out)   :: x(:)    ! solution vector
        integer   , intent(out) :: ierr    ! returned nonzero if error encountered
        integer :: n, nrhs, lda, ldb, info
        integer, allocatable :: ipiv(:)
        
        n = size(a, dim=1)   ! dimension of the matrix A
        nrhs = 1             ! number of right hand side vectors
        lda = size(a, dim=1) ! leading dimension of the matrix A
        ldb = max(1, n)      ! leading dimension of the matrix B
        allocate(ipiv(n))    ! allocate space for pivot indices

        x = b                ! copy b into x as dgesv overwrites b with the solution
        a = a                ! copy a into a as dgesv overwrites a with the LU factors

        ! Call LAPACK routine to solve the system
        call dgesv(n, nrhs, a, lda, ipiv, x, ldb, info)

        if (info /= 0) then
            ierr = info
        else
            ierr = 0
        end if

        deallocate(ipiv)      ! deallocate pivot indices array
    end subroutine solve_eq
end module solve_lapack_mod

which is invalid because of the line a = a since a is intent(in). When fed the error message of gfortran, it produces

module solve_lapack_mod
    use kind_mod, only: dp
    implicit none
    contains
subroutine solve_eq(a, b, x, ierr)
    use, intrinsic :: iso_fortran_env, only: dp => real64
    real(kind=dp), intent(in) :: a(:,:), b(:)
    real(kind=dp), intent(out) :: x(:)
    integer, intent(out) :: ierr
    integer :: n, lda, nrhs, ldb
    integer, allocatable :: ipiv(:)
    real(kind=dp), allocatable :: a_copy(:,:), b_copy(:)
    n = size(a, 1)
    lda = size(a, 1)
    nrhs = 1
    ldb = size(b)
    allocate(ipiv(n))
    allocate(a_copy(n, n))
    allocate(b_copy(n))

    a_copy = a
    b_copy = b
    call dgesv(n, nrhs, a_copy, lda, ipiv, b_copy, ldb, ierr)
    x = b_copy

    deallocate(ipiv, a_copy, b_copy)
end subroutine solve_eq
end module solve_lapack_mod

which compiles and runs with the main program

program main
use solve_lapack_mod, only: solve_eq
integer, parameter :: n = 1000, dp = kind(1.0d0)
real(kind=dp) :: a(n, n), x(n), b(n), bchk(n)
integer :: ierr
logical, parameter :: print_b = .false.
call random_seed()
call random_number(a)
call random_number(b)
call solve_eq(a, b, x, ierr)
if (print_b) then
   bchk = matmul(a, x)
   print*,b
   print*,bchk
   print*,b - bchk
end if
print*,maxval(abs(matmul(a,x) - b))
end program main