write(*,*) 'Hello, world!'
end

subroutine sub()
    implicit none
    write(*,*) "This is subroutine"
end subroutine
 
function trace(a, n)
    implicit none
    real(8), dimension(n, n) :: a
    integer :: n
    real(8) :: trace
    integer :: i
    real(8) :: s
    s = 0
    do i = 1, n
        s = s + a(i, i)
    end do
    trace = s
end

! И. В. Проскуряков.
! Сборник задач по линейной алгебре. 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

! Д. К. Фаддеев, И. С. Соминский.
! Сборник задач по высшей алгебре. 1977.
! № 279. Вычислить определитель.
! Ответ: 900.
 
    integer, parameter :: n=4
    real*8 :: a(n,n) = &
! Обратите внимание, что матрица в фортране записывается по столбцам.
    (/ 1, -2,  3,  4, &
       2,  1, -4,  3, &
       3, -4, -1, -2, &
       4,  3,  2, -1  /)
    integer :: ia(n) = (/ (i, i=1,n) /)
    real*8 det
    print *, det(a, n, ia, n)
    end
    
    recursive real*8 function det(a, n, ia, k)
!    Функция вычисляет минор матрицы a размера n x n, стоящий на пересечении
!    k последних строк и k столбцов, номера которых перечислены в массиве ia.
    integer k,n
    real*8, dimension (n,n) :: a
    integer, dimension(k) :: ia
    integer i, j, sign
    integer, dimension(k-1) :: iar
    real*8 :: s
    if(k==1) then
        det=a(n,ia(1))
        return
    end if
    if(k==n) ia = (/ (i, i=1,n) /)
    s=0.0
    sign = 1
    do i = 1,k
        iar=(/ (ia(j), j=1,i-1), (ia(j), j=i+1,k) /)
        s = s + sign * a(n-k+1,ia(i)) * det(a, n, iar, k-1)
        sign = -sign
    end do
    det = s
    end function det

integer*2 s1,s2
real*4 r
print *,s1,s2
call randu(s1,s2,r)
print *,s1,s2
print *,r
end

integer*2 s1,s2,tm(4)
real*4 r
integer i
CALL TIME(tm)
s1=tm(3)
s2=tm(4)
do i=1,20
    call randu(s1,s2,r)
    print *,r
end do
end

integer*2 s1,s2,t1,t2
real*4 r
integer n
s1=1234
s2=4321
do i=1,100000000
    t1=s1
    t2=s2
    call randu(s1,s2,r)
    n=r
    if(n==1) then
        print *,t1,t2
        exit
    end if
end do
end

integer*2 s1,s2
real*4 r
integer n
s1=14570
s2=21825
call randu(s1,s2,r)
n=r
if(r.eq.1.0) then
    print *,r,n
end if
end

integer*2 s1,s2,tm(4)
real*4 r
integer a(100) ! Массив, который нужно заполнить
integer lb,ub ! Нижняя и верхняя граница диапазона
integer i,p,imx,imn
call time(tm)
s1=tm(3)
s2=tm(4)
lb=-100
ub=100
! Генерация
do i=1,100
    call randu(s1,s2,r)
    p = (ub-lb+1)*r
    p = lb+p
    if(p>ub) p=ub ! проверка, нужная при r=1
    a(i)=p
end do
! Считаем размах
imx=a(1)
imn=a(1)
do i=2,100
    if(a(i)>imx) imx=a(i)
    if(a(i)<imn) imn=a(i)
end do
print *,a ! печать массива
print *,imx-imn ! печать размаха
end

integer a(100) ! Массив, который нужно заполнить
integer lb,ub ! Нижняя и верхняя граница диапазона
integer i,p
call random_seed;
lb=-100
ub=100
! Генерация
do i=1,100
    call random_number(r)
    p = (ub-lb+1)*r
    p = lb+p
    if(p>ub) p=ub ! проверка, нужная при r=1
    a(i)=p
end do
print *,a ! печать массива
end

integer*2 s1,s2,tm(4)
real*4 r
integer a(36) ! Массив, который нужно заполнить
integer lb,ub ! Нижняя и верхняя граница диапазона
integer i,p,t
call time(tm)
s1=tm(3)
s2=tm(4)
lb=1
ub=36
! Начальное заполнение (в стиле фортрана 90)
forall(i=1:36) a(i)=i
! или в стиле фортрана 77
!   do i=1,36
!       a(i)=i
!   end do
! Тасование
do i=36,2,-1
    call randu(s1,s2,r)
    p = (ub-lb+1)*r
    p = lb+p
    if(p>ub) p=ub ! проверка, нужная при r=1
    if(p.ne.i) then
        t=a(i)
        a(i)=a(p)
        a(p)=t
    end if
    ub=ub-1
end do
print *,a ! печать массива
end

integer*2 :: s1=1234, s2=4321
real*4 r
real*8 :: Mean=3.0, Sigma=0.5 ! Требуемые средняя и стандартное отклонение.
real*8 vm, vs ! Выборочная средняя и выборочное стандартное отклонение
real*8 s
integer i, k
integer, parameter :: n=10000 ! объем выборки
 
vm=0
vs=0
do i=1,n
    ! начало вычисления нормальной величины
    s=0.0
    do k=1,12
        call randu(s1,s2,r)
        s=s+r
    end do
    s=(s-6.0)*Sigma+Mean
    ! конец вычисления нормальной величины
    vm=vm+s   ! накопление суммы
    vs=vs+s*s ! накопление суммы квадратов
end do
vm=vm/n ! выборочная средняя
vs=sqrt((vs/n-vm*vm)*n/(n-1)) ! выборочное стандарное отклонение
print *,vm, vs
end