| 1 | ! If using gfortran, force a call to a linked BLAS library from the matmaul() |
|---|
| 2 | ! intrinsic by using the following flags: |
|---|
| 3 | ! -fexternal-blas -fblas-matmul-limit=1 |
|---|
| 4 | ! (See gfortran man page for details.) |
|---|
| 5 | ! |
|---|
| 6 | ! Otherwise replace matmul with a BLAS equivalent (e.g., dgemm). |
|---|
| 7 | |
|---|
| 8 | implicit none |
|---|
| 9 | |
|---|
| 10 | character(len=50) :: arg |
|---|
| 11 | integer :: n, iostat, i, j |
|---|
| 12 | real(8), dimension(:,:), allocatable :: a, b |
|---|
| 13 | |
|---|
| 14 | n = 1000 |
|---|
| 15 | if (command_argument_count() > 0) then |
|---|
| 16 | call get_command_argument(1, arg) |
|---|
| 17 | read(arg,*,iostat=iostat) n |
|---|
| 18 | if (iostat /= 0) error stop 'First argument must be an integer n' |
|---|
| 19 | if (n <= 0) error stop 'n must be > 0' |
|---|
| 20 | endif |
|---|
| 21 | |
|---|
| 22 | allocate(a(n,n), b(n,n)) |
|---|
| 23 | |
|---|
| 24 | a = reshape([((i+j, j=1,n), i=1,n)], [n,n]) |
|---|
| 25 | b = a |
|---|
| 26 | |
|---|
| 27 | a = matmul(a, b) |
|---|
| 28 | |
|---|
| 29 | end |
|---|
