p_d_matmul_f Function

function p_d_matmul_f(A, B) result(C)

Arguments

Type IntentOptional Attributes Name
real(kind=8), intent(in), dimension(:,:) :: A
real(kind=8), intent(in), dimension(:,:) :: B

Return Value real(kind=8), dimension(size(A,1),size(B,2))


Calls

proc~~p_d_matmul_f~2~~CallsGraph proc~p_d_matmul_f~2 p_d_matmul_f blacs_gridinfo blacs_gridinfo proc~p_d_matmul_f~2->blacs_gridinfo descinit descinit proc~p_d_matmul_f~2->descinit distribute_blacs distribute_blacs proc~p_d_matmul_f~2->distribute_blacs gather_blacs gather_blacs proc~p_d_matmul_f~2->gather_blacs numroc numroc proc~p_d_matmul_f~2->numroc pdgemm pdgemm proc~p_d_matmul_f~2->pdgemm shift_dw shift_dw proc~p_d_matmul_f~2->shift_dw

Source Code

function p_d_matmul_f(A,B) result(C)
  real(8),dimension(:,:),intent(in)      :: A ![N,K]
  real(8),dimension(:,:),intent(in)      :: B ![K,M]
  real(8),dimension(size(A,1),size(B,2)) :: C ![N,M]
  real(8)                                :: alfa_,beta_
  integer                                :: Nb
  integer                                :: N,K,M
  integer                                :: Nrows,Kcols,Krows,Mcols
  !
  integer                                :: i,j,lda,info
  integer,external                       :: numroc,indxG2L,indxG2P,indxL2G
  !
  integer,dimension(9)                   :: descAloc,descBloc,descCloc,descC
  real(8),dimension(:,:),allocatable     :: Aloc,Bloc,Cloc
  integer                                :: rankX,rankY
  integer                                :: sendR,sendC,Nipiv
  integer                                :: Nrow,Ncol
  integer                                :: myi,myj,unit
  real(8)                                :: t_stop,t_start
  logical                                :: master
  !
  ! C = alfa*A*B + beta*C
  !
  alfa_     = 1d0 
  beta_     = 0d0 
  !
  N = size(A,1)
  K = size(A,2) !==size(B,1)
  M = size(B,2)
  if(any(shape(B)/=[K,M]))stop "p_d_matmul error: B has illegal shape"
  if(any(shape(C)/=[N,M]))stop "p_d_matmul error: C has illegal shape"
  !
  !INIT SCALAPACK TREATMENT:
  !< Get coordinate of the processes
  call blacs_gridinfo( p_context, p_Nx, p_Ny, rankX, rankY)
  master = (rankX==0).AND.(rankY==0)
  if(rankX<0.AND.rankY<0)goto 201
  !
  I  = min(N,K,M)
  Nb = min(shift_dw(I)/2,64)
  !
  Nrows = numroc(N, Nb, rankX, 0, p_Nx)
  Kcols = numroc(K, Nb, rankY, 0, p_Ny)
  Krows = numroc(K, Nb, rankX, 0, p_Nx)
  Mcols = numroc(M, Nb, rankY, 0, p_Ny)
  !
  !< allocate local distributed A
  !
  call descinit( descAloc, N, K, Nb, Nb, 0, 0, p_context, Nrows, info )
  call descinit( descBloc, K, M, Nb, Nb, 0, 0, p_context, Krows, info )
  call descinit( descCloc, N, M, Nb, Nb, 0, 0, p_context, Nrows, info )
  !
  !< Distribute A,B; allocate C
  allocate(Aloc(Nrows,Kcols));Aloc=0d0
  allocate(Bloc(Krows,Mcols));Bloc=0d0
  allocate(Cloc(Nrows,Mcols));Cloc=0d0
  call Distribute_BLACS(A,Aloc,descAloc)
  call Distribute_BLACS(B,Bloc,descBloc)
  !
  call PDGEMM( 'No transpose', 'No transpose', N, M, K, &
       alfa_, &
       Aloc, 1, 1, descAloc, &
       Bloc, 1, 1, descBloc, &
       beta_, &
       Cloc, 1, 1, descCloc )
  !
  C=0d0
  call descinit( descC, N, M, Nb, Nb, 0, 0, p_context, Nrows, info )
  call Gather_BLACS(Cloc,C,descC)
201 continue
  return
  !
end function p_d_matmul_f