PROGRAM Main
    IMPLICIT NONE
    INTEGER, PARAMETER :: dp = SELECTED_REAL_KIND(12)
    REAL(KIND=dp), DIMENSION(:, :), ALLOCATABLE :: a, a1, b, b1
    INTEGER :: size, info, i
    CHARACTER(LEN=100) :: filename

! Read the file name from the argument and set up the data.
    IF (COMMAND_ARGUMENT_COUNT() /= 1) THEN
        WRITE (*,'("Wrong number of arguments")')
        STOP 1
    END IF
    CALL GET_COMMAND_ARGUMENT(1,filename)
    OPEN (42,FILE=filename,ACTION='READ',FORM='UNFORMATTED')
    READ (42) size
    IF (size < 1 .OR. size > 1000000) THEN
        WRITE (*,'("Invalid value of size: ",I0)') size
        STOP 1
    END IF
    ALLOCATE(a(size,size),b(size,size),a1(size,size),b1(size,size))
    READ (42) a
    READ (42) b
    CLOSE(42)
    a1 = a
    b1 = b

! Time and check the use of LAPACK.
    CALL Times('')
    CALL DPOTRF('l',size,a,size,info)
    IF (info /= 0) CALL Fail
    DO i = 2, size
        a(:i-1,i) = 0.0_dp
    END DO
    CALL Times('LAPACK Cholesky')
    WRITE (*,'("Error in result = ",ES10.2)') &
        MAXVAL(ABS(a1 - MATMUL(a, TRANSPOSE(a))))
    CALL Times('')
    CALL DPOTRS('l',size,size,a,size,b,size,info)
    IF (info /= 0) CALL Fail
    CALL Times('LAPACK solver')
    WRITE (*,'("Error in result = ",ES10.2)') &
        MAXVAL(ABS(b1 - MATMUL(a1, b)))
 
! Time and check the use of the LAPACK code converted to Fortran 90.
    a = a1
    CALL Times('')
    CALL Cholesky(a)
    CALL Times('Coded Cholesky')
    WRITE (*,'("Error in result = ",ES10.2)') &
        MAXVAL(ABS(a1 - MATMUL(a, TRANSPOSE(a))))
    b = b1
    CALL Times('')
    CALL Solve(a, b)
    CALL Times('Coded solver')
    WRITE (*,'("Error in result = ",ES10.2)') &
        MAXVAL(ABS(b1 - MATMUL(a1, b)))
    
CONTAINS

    SUBROUTINE Fail
        WRITE (*,*) 'Error in calling LAPACK'
        STOP 1
    END SUBROUTINE Fail

    SUBROUTINE Times (which)
! If which is not empty, print the times since the previous call.
        CHARACTER(LEN=*), INTENT(IN) :: which
        REAL(KIND=dp), SAVE :: last_wall = 0.0_dp, last_cpu = 0.0_dp
        REAL(KIND=dp) :: wall, cpu
        INTEGER :: m, n

        wall = last_wall
        cpu = last_cpu
        CALL SYSTEM_CLOCK(m,n)
        last_wall = m/REAL(n,KIND=dp)
        CALL CPU_TIME(last_cpu)
        IF (LEN(which) > 0) THEN
            wall = last_wall-wall
            cpu = last_cpu-cpu
            WRITE (*,'(A," time = ",F0.2," seconds, CPU = ",F0.2," seconds")') &
                which,wall,cpu
        END IF
    END SUBROUTINE Times

    SUBROUTINE Cholesky(a)
        REAL(KIND=dp), INTENT(INOUT) :: a(:, :)
        INTEGER :: i, j, n
        REAL(KIND=dp) :: x

        n = UBOUND(a,1)
        DO j = 1, n
! Expanding the following:
!     a(:j-1,j) = 0.0_dp
!     x = SQRT(a(j,j)-DOT_PRODUCT(a(j, :j-1), a(j, :j-1)))
!     a(j,j) = x
!     a(j+1:, j) = (a(j+1:, j) - &
!     MATMUL(a(j+1:, :j-1), a(j, :j-1))) / x
            x = 0.0_dp
            DO i = 1, j-1
                a(i,j) = 0.0_dp
                x = x+a(j,i)*a(j,i)
            END DO
            x = SQRT(a(j,j)-x)
            a(j,j) = x
            DO i = j+1, n
                a(i, j) = (a(i, j) - &
                    DOT_PRODUCT(a(i, :j-1), a(j, :j-1))) / x
            END DO
! End of expansion
        END DO
    END SUBROUTINE Cholesky

    SUBROUTINE Solve(a,b)
        REAL(KIND=dp), INTENT(IN) :: a(:, :)
        REAL(KIND=dp), INTENT(INOUT) :: b(:, :)
        INTEGER :: i, j, k, n

        n = UBOUND(a,1)
        DO j = 1,n
            DO i = 1,n
                b(i,j) = b(i,j)/a(i,i)
! Expanding the following:
!     b(i+1:,j) = b(i+1:,j)-b(i,j)*a(i+1:,i)
                DO k = i+1, n
                    b(k,j) = b(k,j)-b(i,j)*a(k,i)
                END DO
! End of expansion
            END DO
        END DO
        DO j = 1,n
            DO i = n,1,-1
! Expanding the following:
!     b(i,j) = (b(i,j)-DOT_PRODUCT(a(i+1:,i),b(i+1:,j)))/a(i,i)
                DO k = i+1, n
                    b(i,j) = b(i,j)-a(k,i)*b(k,j)
                END DO
                b(i,j) = b(i,j)/a(i,i)
! End of expansion
            END DO
        END DO
    END SUBROUTINE Solve

END PROGRAM Main
