Умножение матрицы Лапака на тождество Fortran DGEMM - PullRequest
Новая
       139

Умножение матрицы Лапака на тождество Fortran DGEMM

0 голосов
/ 30 марта 2020

Я столкнулся с проблемой, когда умножение ненулевой матрицы на единичной матрице через Лапака дает мне нулевую матрицу. Все матрицы с положительными элементами 

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

матрицы идентичности DM

Я пытаюсь получить умножение D * W '* M, где W' - транспонирование W Вот код Фортрана, использующий операцию DGEMM 

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

Ответы последовательной печати: True, True, False. Итак, первое умножение работает, и я получаю не нулевую матрицу, но во втором все элементы равны 0.

Я знаю, что мне не нужно умножение матриц на единичных матрицах, но я хочу выяснить, в чем проблема если я сделаю это.

Еще один вопрос, могу ли я сделать это эффективным с точки зрения памяти без временных матриц Main и C?

EDIT

Разобрался с загрузкой полученной матрицы после первого умножения что все элементы нулевые. Не могу понять, почему ЛЮБОЙ (C> 0) является Истиной на втором этапе.

1 Ответ

2 голосов
/ 31 марта 2020

Во-первых, обратите внимание, что DGEMM является частью библиотеки BLAS, а не LAPACK - последняя является библиотекой более высокого уровня.

У вас есть несколько проблем с вашими вызовами в DGEMM

  1. Настоящие константы имеют вид, и вы должны предоставить правильный вид. В частности, DGEMM ожидает, что в старомодном Фортране называлось Double Precision. Заданные вами константы являются реальными по умолчанию. Это приведет к ошибкам, и я настоятельно рекомендую предоставить вид для ВСЕХ реальных констант, если точность, требуемая вашей программой, не является точностью по умолчанию
  2. Из первых трех целых чисел первые два всегда являются формой результата матрица. C равно n1xm1, следовательно, изменение первого вызова DGEMM
  3. Начальное измерение матрицы составляет 99,999% от времени, когда первое измерение было выделено / объявлено - оно не имеет ничего общего с математикой в все, это чисто так, что DGEMM может определить, как матрица размещена в памяти. Это было неправильно для W в первой DGEMM

. Я бы также предположил, что в качестве правильной проверки результатов достаточно просто использовать MATUL intrinsi c, а не половину. тот, который вы делаете в оригинале, и избегайте использования единичных матриц для тестов, если вы действительно не хотите использовать единичные матрицы, так как высокая симметрия и большое количество нулей могут легко замаскировать ошибки.

Объединение всего этого вместе нерелевантные части и изменение вашей программы в чуть более современной форме, я получаю 

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
...