Matrix multiplication with LAPACK, BLAS, dgemm, intiger parameter type

Question

I am testing different options for matrix multiplication with different parameter types for matrices. One of them is dgemm routine within BLAS. When I wanted to make a matrix defined as integer(kind=1) with a size of 1000x1000 (nxp) it crashed with dgemm but matmul does it well. When I decrease the size of the matrices to 500x500, both work well. Furthermore, I defined all the matrices as real(8), both computed the matrix product but the results were different. The code that I am using is:

    program test
    implicit none

    real(8), allocatable :: x(:,:),xi(:,:),xt(:,:)
    integer(kind=1), allocatable :: z(:,:)
    integer :: i,j,n,p
    real(8):: u,myone= 1.d0

    n=1000
    p=1000
    allocate(x(n,n),z(n,p),xi(n,n),xt(n,n))
    do i=1,n
      do j=1,p
        call random_number(u)
        z(i,j)=real(floor(u*3),8)
      enddo
    enddo

    print*,"matmul"
    x=matmul(z,transpose(z))

    do i=1,min(10,n)
       write(*,'(10(g10.3,x))') x(i,1:min(10,n))
    enddo 

    print*,"dgemm"
    call dgemm('n'   ,'t'   ,n,n,p,myone  ,Z,n  ,Z,n  ,myone ,X,n)

    do i=1,min(10,n)
       write(*,'(10(g10.3,x))') x(i,1:min(10,n))
    enddo 

    end program test

I compile the code with make statement which runs the following code (I have named it: Makefile):

    f90=ifort
    optf90=-O0 -heap-arrays
    optdir=-I
    mkl=-L/opt/intel/mkl/lib -lmkl_intel_lp64 -lmkl_intel_thread -lmkl_core -openmp -lpthread
    prog=test
    dir= .
    a.out:  $(prog).o
            $(f90)  $(optf90) \
                    $(prog).o  \
                    $(mkl)  $(libf77)

    $(prog).o:      $(prog).f90
            $(f90)  $(optdir)$(dir)  -c $(optf90) $(prog).f90  

Does anyone know what is the problem with dgemm routine for big matrices defined as intigers and what might be the reason for different outcomes with matmul/dgemm?

Answer 1

DGEMM works on double precision real numbers, not on integers (of any kind).

I would be surprised if you got any (correct) results at all when using integer numbers with DGEMM.

MATMUL, on the other hand, accepts integers as input.