Search code examples
matrixfortranhpclapackblas

Lapack Matrix multiplication on identity Fortran DGEMM

I faced a problem, where multiplication not null matrix on identity matrix via Lapack gives me a null matrix. All matrices are with positive elements

Dimensions of the matrices:
W = M1*N1
D = N1*N1
M = M1*M1

D M -identity matrices

What I try to do is to get multiplication of D*W'*M , where W' is a transpose of W. Here is a Fortran code using DGEMM operation

PROGRAM SIRT
DOUBLE PRECISION,ALLOCATABLE,DIMENSION(:,:) :: W, D, M, C, MAIN
DOUBLE PRECISION,ALLOCATABLE,DIMENSION(:) :: b, x
DOUBLE PRECISION :: pho
INTEGER, PARAMETER :: M1 = 27000
INTEGER, PARAMETER :: N1 = 1000
INTEGER, PARAMETER :: num_iterations = 200
INTEGER :: i, j, k

allocate(W(1:M1,1:N1))
allocate(D(1:N1,1:N1))
allocate(M(1:M1,1:M1))
allocate(C(1:N1,1:M1))
allocate(MAIN(1:N1,1:M1))
allocate(b(1:M1))
allocate(x(1:N1))

D = 0
M = 0

DO i=1,N1
    D(i,i) = 1
END DO

DO i=1,M1
    M(i,i) = 1
END DO

OPEN(UNIT=11, FILE="Wmpi.txt")
DO i = 1,M1
    READ(11,*) (W(i,j),j=1,N1)
END DO
print *,ANY(W>0)
CLOSE (11, STATUS='KEEP') 

OPEN(UNIT=11, FILE="bmpi.txt")
DO i = 1,M1


READ(11,*) b(i)
END DO
CLOSE (11, STATUS='KEEP') 

CALL DGEMM('N', 'T', N1, N1, N1, 1.0, D, N1, W, N1, 0.0, C, N1)
print *,ANY(C>0)
CALL DGEMM('N', 'N', N1, M1, M1, 1.0, C, N1, M, M1, 0.0, MAIN, N1)
print *,ANY(MAIN>0)
pho = DLANGE('F', N1, M1, C, N1, x)
END PROGRAM SIRT

The answers of sequential print are True,True,False. So the first multiplication works and I get not null matrix, but in the second all elements are 0.

I know that I don't need matrix multiplication on identity matrices, but I want to figure out what's the problem if I do that.

Another question can I do it memory efficient without temp matrices Main and C?

EDIT

Figured out downloading a resulting matrix after first multiplication that all elements are null. Can't understand why ANY(C>0) is True at second stage.

Solution

  • Firstly note that DGEMM is part of the BLAS library, not LAPACK - the latter is a higher level library.

    You have a few problems with your calls to DGEMM

    1. Real constant have kind, and you must provide the correct kind. In particular DGEMM expects what in old fashioned Fortran was called Double Precision. The constants you provided are default kind real. This will cause errors, and I would strongly recommend providing a kind for ALL real constants if the precision required by your program is not the default precision
    2. Of the first three integers the first two are always the shape of the result matrix. C is n1xm1, hence the change to the first call to DGEMM
    3. The leading dimension of a matrix is 99.999% of the time the first dimension as allocated/declared - it has nothing to do with the maths at all, it is purely so DGEMM can work out how the matrix is laid out in memory. It was wrong for W in the first DGEMM

    I would also suggest that as a proper check of the results is fairly easy to do by use of the MATUL intrinsic do that rather than the half hearted one you do in the original, and avoid using unit matrices for tests unless you really want unit matrices, as the high symmetry and large numbers of zeros can easily mask errors.

    Pulling this all together, cutting our the irrelevant parts and modifying your program to a slightly more modern form I get

    Program sirt

  Integer, Parameter :: wp = Kind( 1.0d0 )

  Real( wp ),Allocatable,Dimension(:,:) :: w, d, m, c, main, main_compare
  Real( wp ),Allocatable,Dimension(:) :: b, x
  ! Change problem size to something manageable on my laptop
!!$  Integer, Parameter :: m1 = 27000
!!$  Integer, Parameter :: n1 = 1000
  Integer, Parameter :: m1 = 2700
  Integer, Parameter :: n1 = 100
!!$  Integer :: i

  Allocate(w(1:m1,1:n1))
  Allocate(d(1:n1,1:n1))
  Allocate(m(1:m1,1:m1))
  Allocate(c(1:n1,1:m1))
  Allocate(main(1:n1,1:m1))
  Allocate(main_compare(1:n1,1:m1))
  Allocate(b(1:m1))
  Allocate(x(1:n1))

  d = 0.0_wp
  m = 0.0_wp

  ! Don't trust unit matrices for tests, too much symmetry, too many zeros - use Random numbers
!!$  Do i=1,n1
!!$     d(i,i) = 1.0_wp
!!$  End Do
  Call Random_number( d )
!!$
!!$  Do i=1,m1
!!$     m(i,i) = 1.0_wp
!!$  End Do
  Call Random_number( m )

  ! Don't have the file - use random numbers
!!$  open(unit=11, file="wmpi.txt")
!!$  do i = 1,m1
!!$     read(11,*) (w(i,j),j=1,n1)
!!$  end do
!!$  print *,any(w>0)
!!$  close (11, status='keep') 
  Call random_Number( w )

  ! 1) Real constant have kind, and you must provide the correct kind
  ! 2) Of the first three constant the first two are always the shape
  !    of the result matrix. C is n1xm1, hece the change
  ! 3) The leading dimension of a matrix is 99.999% of the time
  !    the first dimension as allocated/declared  -it has nothing
  !    to do with the maths at all. It was wrong for W in the
  !    first matmul
  ! 4) DGEMM is part of the BLAS library - not LAPACK
!!$  CALL DGEMM('N', 'T', N1, N1, N1, 1.0, D, N1, W, N1, 0.0, C, N1)
!!$  CALL DGEMM('N', 'N', N1, M1, M1, 1.0, C, N1, M, M1, 0.0, MAIN, N1)
  Call dgemm('n', 't', n1, m1, n1, 1.0_wp, d, n1, w, m1, 0.0_wp, c, n1)
  Call dgemm('n', 'n', n1, m1, m1, 1.0_wp, c, n1, m, m1, 0.0_wp, main, n1)

  ! 5) Don't do half hearted checks on the results when proper checks are easy
  main_compare = Matmul( Matmul( d, Transpose( w ) ), m )
  Write( *, * ) 'Max error ', Maxval( Abs( main - main_compare ) )

  ! Check we haven't somehow managed all zeros in both matrices ...
  Write( *, * ) main( 1:3, 1 )
  Write( *, * ) main_compare( 1:3, 1 )

End Program sirt
ian@eris:~/work/stack$ gfortran-8  -std=f2008 -fcheck=all -Wall -Wextra -pedantic -O matmul.f90 -lblas
/usr/bin/ld: warning: libgfortran.so.4, needed by //usr/lib/x86_64-linux-gnu/libopenblas.so.0, may conflict with libgfortran.so.5
ian@eris:~/work/stack$ ./a.out
 Max error    3.6379788070917130E-011
   38576.055405987529        33186.640731082334        33818.909332709263     
   38576.055405987536        33186.640731082334        33818.909332709263     
ian@eris:~/work/stack$ ./a.out
 Max error    2.9103830456733704E-011
   34303.739077708480        34227.623080598998        34987.143088270866     
   34303.739077708473        34227.623080598998        34987.143088270859     
ian@eris:~/work/stack$ ./a.out
 Max error    3.2741809263825417E-011
   35968.603030053979        34778.110740682620        32732.657800858586     
   35968.603030053971        34778.110740682612        32732.657800858586     
ian@eris:~/work/stack$ ./a.out
 Max error    2.9103830456733704E-011
   31575.076511213174        35879.913361891951        35278.030249048912     
   31575.076511213178        35879.913361891951        35278.030249048912