Как правильно позвонить SGEMV в Фортран? - PullRequest
Новая
       117

Как правильно позвонить SGEMV в Фортран?

0 голосов
/ 14 января 2019

Я хочу выполнить продукт Matrix-Vector на фортране, используя подпрограмму SGEMV от BLAS. У меня есть код, который похож на это: 

program test
integer, parameter :: DP = selected_real_kind(15)
real(kind=DP), dimension (3,3) :: A
real(kind=DP), dimension (3) :: XP,YP
call sgemv(A,XP,YP)

A - матрица 3x3, XP и YP - векторы. Во включенном модуле можно увидеть следующий код: 

PURE SUBROUTINE SGEMV_F95(A,X,Y,ALPHA,BETA,TRANS)
    ! Fortran77 call:
    ! SGEMV(TRANS,M,N,ALPHA,A,LDA,X,INCX,BETA,Y,INCY)
    USE F95_PRECISION, ONLY: WP => SP
    REAL(WP), INTENT(IN), OPTIONAL :: ALPHA
    REAL(WP), INTENT(IN), OPTIONAL :: BETA
    CHARACTER(LEN=1), INTENT(IN), OPTIONAL :: TRANS
    REAL(WP), INTENT(IN) :: A(:,:)
    REAL(WP), INTENT(IN) :: X(:)
    REAL(WP), INTENT(INOUT) :: Y(:)
END SUBROUTINE SGEMV_F95

Я понимаю, что некоторые параметры являются необязательными, поэтому где я ошибаюсь при вызове метода?

Ответы [ 3 ]

0 голосов
/ 14 января 2019

Точности несовместимы. Вы вызываете sgemv, который принимает аргументы одинарной точности, но вы передаете массивы и векторы двойной точности.

0 голосов
/ 14 января 2019

Когда вы смотрите на процедуры BLAS или LAPACK, вы всегда должны смотреть на первую букву:

  • S: одинарная точность
  • D: двойная точность
  • C: комплекс одинарной точности
  • Z: комплекс двойной точности

Вы определили свою матрицу A, а также векторы XP и YP как число двойной точности, используя выражение: 

integer, parameter :: DP = selected_real_kind(15)

Так что для этого вам нужно использовать dgemv или определить свою точность как одинарную точность.

Существует также разница между звонками dgemv и dgemv_f95. dgemv_f95 является частью Intel MKL и не является общепринятым наименованием. Из соображений переносимости я бы не стал использовать это обозначение, а придерживался классического dgemv, который также является частью Intel MKL.

DGEMV выполняет одну из матрично-векторных операций 

y := alpha*A*x + beta*y,   or   y := alpha*A**T*x + beta*y,

, где alpha и beta - скаляры, x и y - векторы, а A - это m на n матрица.

Если вы хотите знать, как вызывать функцию, я предлагаю посмотреть здесь , но в итоге она должна выглядеть примерно так: 

call DGEMV('N',3,3,ALPHA,A,3,XP,1,BETA,YP,1)
0 голосов
/ 14 января 2019

Возможно, требуется параметр trans? 

trans: Must be 'N', 'C', or 'T'.

(Согласно примечанию внизу Справочник разработчика для библиотеки Intel® Math Kernel - Fortran .)

...