program lincombtest
    use iso_fortran_env, only: int64
    implicit none

    integer, parameter :: n = 100

    real :: f
    real :: var(10*n,10*n,n,2), res(10*n,10*n,n)

    integer(int64) :: t0, t1, rate
    integer :: ntot = size(res)

    call random_number(f)
    call random_number(var)

    ! Naive approach
    call system_clock(t0, rate)
    res = f*var(:,:,:,1) + (1-f)*var(:,:,:,2)
    call system_clock(t1, rate)

    write(*,"(a,f0.4)") "Simple expression: t1 - t0 = ", (t1 - t0)/(1.0d0*rate)

    ! (Open)BLAS
    call system_clock(t0, rate)
    call sscal(ntot, f, var(:,:,:,1), 1)
    call saxpy(ntot, 1-f, var(:,:,:,2), 1, var(:,:,:,1), 1)
    call system_clock(t1, rate)

    write(*,"(a,f0.4)") "BLAS: t1 - t0 = ", (t1 - t0)/(1.0d0*rate)
    write(*,"(a,x,f0.6,x,f0.6)") "Comparing results:", var(1,1,1,1), res(1,1,1)

end program

! gfortran -Ofast -march=native blah.f90 -lblas && ./a.out
! Simple expression: t1 - t0 = .1500
! BLAS: t1 - t0 = .0398
! Comparing results: .477540 .477540