Умножение матриц с LAPACK, BLAS, dgemm, целочисленным типом параметра

Я тестирую разные варианты умножения матриц с разными типами параметров для матриц. Одним из них является процедура dgemm в BLAS. Когда я хотел создать матрицу, определяемую как целое число (вид = 1) с размером 1000x1000 (nxp), это привело к сбою с dgemm, но matmul делает это хорошо. Когда я уменьшаю размер матриц до 500x500, обе работают хорошо. Кроме того, я определил все матрицы как действительные (8), обе вычисляли матричное произведение, но результаты были разными. Код, который я использую:

    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

Я компилирую код с оператором make, который запускает следующий код (я назвал его 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  

Кто-нибудь знает, в чем проблема с подпрограммой dgemm для больших матриц, определенных как целые числа, и что может быть причиной разных результатов с matmul/dgemm?