Hauptprogramm

module mysubs

contains

  subroutine init(a, b)
    use kinds
    implicit none

    real(kind=REAL8), dimension(:,:), intent(out) :: a, b
    integer :: n, i, j
    
    n = size(a,1)
    do i=1, n
       do j=1, n
          a(i,j) = i + j
          b(i,j) = i - j
       enddo
       b(i,i) = 1.0
    enddo
    
  end subroutine init


  subroutine test(c)
    use kinds
    implicit none
    
    real(kind=REAL8), dimension(:,:), intent(in) :: c
    integer          :: n, i
    real(kind=REAL8) :: spur, spur1, cth
  
    n = size(c,1)
    spur1 = n * (n+1)
    print *, "Spur theoretisch: ", spur1
    spur = 0.0
    do i = 1, n
       spur = spur + c(i,i)
    enddo
    print *, "Spur            : ", spur
  
    cth = 333331003.0D0
    print *, "c(1,2) theor.   : ", cth
    print *, "c(1,2)          : ", c(1,2)
  
  end subroutine test

end module mysubs


program matmulblas

  ! Multiplikation von Matrizen

  use kinds
  use blas
  use mysubs
  implicit none
  
  integer, parameter                :: n = 1000 
  real(kind=REAL8), dimension(n, n) :: a, b, c

  ! Anfangswerte setzen
  call init(a, b)

  ! Matrixmultiplikation
  call gemm(a, b, c)

  ! Testdaten ausgeben
  call test(c)

end program matmulblas