Работа с массивами в современном фортране

! И. В. Проскуряков.
! Сборник задач по линейной алгебре. 2005.
! № 857. Найти обратную матрицу к матрице
! |3  3 -4 -3|
! |0  6  1  1|
! |5  4  2  1|
! |2  3  3  2|
! В задачнике имеется следующий ответ
! | -7   5  12  -19|
! |  3  -2  -5    8|
! | 41 -30 -69  111|
! |-59  43  99 -159|
real*8 :: a(4,4) = &
! (Обратите внимание, что матрица в фортране записывается по столбцам.)
    (/ 3, 0, 5, 2, &
       3, 6, 4, 3, &
      -4, 1, 2, 3, &
      -3, 1, 1, 2 /)
interface
    function inv2(a)
    real(8), dimension(2,2) :: a, inv2
    end function inv2
end interface
integer i,j
real*8 :: a11_inv(2,2), d_inv(2,2)
 
a11_inv=inv2(a(1:2,1:2))
d_inv=inv2(a(3:4,3:4)-matmul(a(3:4,1:2),matmul(a11_inv,a(1:2,3:4))))
a(1:2,1:2)=a11_inv+matmul(a11_inv,matmul(a(1:2,3:4), &
                            matmul(d_inv,matmul(a(3:4,1:2),a11_inv))))
a(1:2,3:4)=-matmul(a11_inv,matmul(a(1:2,3:4),d_inv))
a(3:4,1:2)=-matmul(d_inv,matmul(a(3:4,1:2),a11_inv))
a(3:4,3:4)=d_inv
! Результат
print '(4F7.1)', ((a(i,j), j=1,4), i=1,4)
end
    
function inv2(a)
! Обращение матрицы порядка 2
real(8) :: a(2,2), inv2(2,2), det
inv2(1,1)=a(2,2)
inv2(1,2)=-a(1,2)
inv2(2,1)=-a(2,1)
inv2(2,2)=a(1,1)
det=1.0/(a(1,1)*a(2,2)-a(1,2)*a(2,1))
inv2=inv2*det
end function inv2