f90ファイル実行時のエラー

実現したいこと

ある.f90ファイルを実行したい

発生している問題・分からないこと

.f90ファイルをターミナルでコンパイルした後、./a.outで実行しようとするとエラーが表示された。

エラーメッセージ

error
1Program received signal SIGSEGV: Segmentation fault - invalid memory reference.
2
3Backtrace for this error:
4#0  0x7f0bad623960 in ???
5#1  0x7f0bad622ac5 in ???
6#2  0x7f0bad24251f in ???
7	at ./signal/../sysdeps/unix/sysv/linux/x86_64/libc_sigaction.c:0
8#3  0x55598b12e736 in ???
9#4  0x55598b12f0b8 in ???
10#5  0x55598b12f6b6 in ???
11#6  0x7f0bad229d8f in __libc_start_call_main
12	at ../sysdeps/nptl/libc_start_call_main.h:58
13#7  0x7f0bad229e3f in __libc_start_main_impl
14	at ../csu/libc-start.c:392
15#8  0x55598b12e154 in ???
16#9  0xffffffffffffffff in ???
17Segmentation fault
18

該当のソースコード

fortran
1module subprogs
2  implicit none
3contains
4  subroutine inverse_mat_gjp(a0, ai, n)
5    integer, intent(in) :: n
6    real(8), intent(in) :: a0(n, n)
7    real(8), intent(out) :: ai(n, n)
8    integer i, k, is, m, n2
9    real(8) ar, am
10    real(8), allocatable :: a(:, :), w(:)
11    n2 = 2 * n
12    allocate(a(n, n2), w(n2), stat = is)
13    if (is /= 0) stop 'cannot allocate (n is too large)'
14    a(1:n, 1:n) = a0(1:n, 1:n)
15    a(1:n, n+1:2*n) = 0.0d0
16    do i = 1, n
17       a(i, i+n) = 1.0d0
18    end do
19    do k = 1, n
20       m = k
21       am = abs(a(k, k))
22       do i = k+1, n
23          if (abs(a(i, k)) > am) then
24             am = abs(a(i, k))
25             m = i
26          endif
27       end do
28       if (am == 0.0d0) stop 'A is singular'
29       if (k /= m) then
30          w(k:n2) = a(k, k:n2)
31          a(k, k:n2) = a(m, k:n2)
32          a(m, k:n2) = w(k:n2)
33       end if
34       ar = 1.0d0 / a(k, k)
35       a(k, k) = 1.0d0
36       a(k, k+1:n2) = ar * a(k, k+1:n2)
37       do i = 1, n
38          if (i /= k) then
39             a(i, k+1:n2) = a(i, k+1:n2) - a(i, k) * a(k, k+1:n2)
40             a(i, k) = 0.0d0
41          end if
42       end do
43    end do
44    ai(1:n, 1:n) = a(1:n, n+1:n2)
45  end subroutine inverse_mat_gjp
46end module subprogs
47
48program main
49  use subprogs
50  implicit none
51  real(8), allocatable :: a(:, :), ai(:, :)
52  integer n, m, i
53  integer :: fi = 10, fo = 11
54  open(fi, file = 'input.d')
55  open(fo, file = 'output.d')
56  read(fi, *) a(:, :)
57  n = size(a, 1)
58  m = size(a, 2)
59  if (n /= m) stop 'not square'
60  call inverse_mat_gjp(a, ai, n)
61  do i = 1, n
62     write(fo, '(10e12.4)') ai(i, 1:n)
63  end do
64  a = matmul(a, ai)
65  close(fi)
66  close(fo)
67end program main
68