program main character :: transa, transb integer :: m, n, k double precision :: alphah, alphal double precision, allocatable :: ah(:,:), al(:,:) integer :: lda double precision, allocatable :: bh(:,:), bl(:,:) integer :: ldb double precision :: betah, betal double precision, allocatable :: ch(:,:), cl(:,:) integer :: ldc integer :: i, j transa = 'n' transb = 'n' m = 2 n = 3 k = 4 lda = max(1, m) ldb = max(1, k) ldc = max(1, m) allocate(ah(lda,k)) allocate(al(lda,k)) allocate(bh(ldb,n)) allocate(bl(ldb,n)) allocate(ch(ldc,n)) allocate(cl(ldc,n)) call random_number(ah) al = 0d0 call random_number(bh) bl = 0d0 call random_number(ch) cl = 0d0 call random_number(alphah) alphal = 0d0 call random_number(betah) betal = 0d0 print *, "INPUT" call print_scalar("alpha", alphah, alphal) call print_matrix("a", ah, al) call print_matrix("b", bh, bl) call print_scalar("beta", betah, betal) call print_matrix("c", ch, cl) call ddgemm(transa, transb, m, n, k, alphah, alphal, ah, al, lda, & & bh, bl, ldb, betah, betal, ch, cl, ldc) print *, "OUTPUT" call print_matrix("c", ch, cl) contains subroutine print_matrix(vname, ah, al) character(*), intent(in) :: vname double precision, intent(in) :: ah(:,:), al(:,:) integer :: i, j do i = 1, size(ah, 1) do j = 1, size(ah, 2) print '(a,"(",i0,",",i0,") =",2e23.15)', & & vname, i, j, ah(i,j), al(i,j) end do end do end subroutine subroutine print_scalar(vname, xh, xl) character(*), intent(in) :: vname double precision, intent(in) :: xh, xl print '(a," =",2e23.15)', vname, xh, xl end subroutine end program