! This file is part of dftd4. ! SPDX-Identifier: LGPL-3.0-or-later ! ! dftd4 is free software: you can redistribute it and/or modify it under ! the terms of the Lesser GNU General Public License as published by ! the Free Software Foundation, either version 3 of the License, or ! (at your option) any later version. ! ! dftd4 is distributed in the hope that it will be useful, ! but WITHOUT ANY WARRANTY; without even the implied warranty of ! MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the ! Lesser GNU General Public License for more details. ! ! You should have received a copy of the Lesser GNU General Public License ! along with dftd4. If not, see <https://www.gnu.org/licenses/>. !> Interface to BLAS library #ifndef IK #define IK i4 #endif module dftd4_blas use mctc_env, only : sp, dp, ik => IK implicit none private public :: d4_gemv, blas_gemv !> Performs one of the matrix-vector operations !> !> y := alpha*A*x + beta*y, or y := alpha*A**T*x + beta*y, !> !> where alpha and beta are scalars, x and y are vectors and A is an !> m by n matrix. interface d4_gemv module procedure :: d4_sgemv module procedure :: d4_dgemv module procedure :: d4_sgemv312 module procedure :: d4_sgemv321 module procedure :: d4_dgemv312 module procedure :: d4_dgemv321 end interface d4_gemv !> Performs one of the matrix-vector operations !> !> y := alpha*A*x + beta*y, or y := alpha*A**T*x + beta*y, !> !> where alpha and beta are scalars, x and y are vectors and A is an !> m by n matrix. interface blas_gemv pure subroutine sgemv(trans, m, n, alpha, a, lda, x, incx, beta, y, incy) import :: sp, ik integer(ik), intent(in) :: lda real(sp), intent(in) :: a(lda, *) real(sp), intent(in) :: x(*) real(sp), intent(inout) :: y(*) real(sp), intent(in) :: alpha real(sp), intent(in) :: beta character(len=1), intent(in) :: trans integer(ik), intent(in) :: incx integer(ik), intent(in) :: incy integer(ik), intent(in) :: m integer(ik), intent(in) :: n end subroutine sgemv pure subroutine dgemv(trans, m, n, alpha, a, lda, x, incx, beta, y, incy) import :: dp, ik integer(ik), intent(in) :: lda real(dp), intent(in) :: a(lda, *) real(dp), intent(in) :: x(*) real(dp), intent(inout) :: y(*) real(dp), intent(in) :: alpha real(dp), intent(in) :: beta character(len=1), intent(in) :: trans integer(ik), intent(in) :: incx integer(ik), intent(in) :: incy integer(ik), intent(in) :: m integer(ik), intent(in) :: n end subroutine dgemv end interface blas_gemv contains subroutine d4_sgemv312(amat, xvec, yvec, alpha, beta, trans) real(sp), intent(in), contiguous, target :: amat(:, :, :) real(sp), intent(in) :: xvec(:) real(sp), intent(inout), contiguous, target :: yvec(:, :) real(sp), intent(in), optional :: alpha real(sp), intent(in), optional :: beta character(len=1), intent(in), optional :: trans real(sp), pointer :: aptr(:, :), yptr(:) character(len=1) :: tra if (present(trans)) then tra = trans else tra = 'n' end if if (any(tra == ['n', 'N'])) then aptr(1:size(amat, 1)*size(amat, 2), 1:size(amat, 3)) => amat yptr(1:size(yvec, 1)*size(yvec, 2)) => yvec else aptr(1:size(amat, 1), 1:size(amat, 2)*size(amat, 3)) => amat yptr(1:size(yvec, 1) * size(yvec, 2)) => yvec end if call d4_gemv(aptr, xvec, yptr, alpha, beta, tra) end subroutine d4_sgemv312 subroutine d4_sgemv321(amat, xvec, yvec, alpha, beta, trans) real(sp), intent(in), contiguous, target :: amat(:, :, :) real(sp), intent(in), contiguous, target :: xvec(:, :) real(sp), intent(inout) :: yvec(:) real(sp), intent(in), optional :: alpha real(sp), intent(in), optional :: beta character(len=1), intent(in), optional :: trans real(sp), pointer :: aptr(:, :), xptr(:) character(len=1) :: tra if (present(trans)) then tra = trans else tra = 'n' end if if (any(tra == ['n', 'N'])) then aptr(1:size(amat, 1), 1:size(amat, 2)*size(amat, 3)) => amat xptr(1:size(xvec, 1)*size(xvec, 2)) => xvec else aptr(1:size(amat, 1)*size(amat, 2), 1:size(amat, 3)) => amat xptr(1:size(xvec, 1) * size(xvec, 2)) => xvec end if call d4_gemv(aptr, xptr, yvec, alpha, beta, tra) end subroutine d4_sgemv321 subroutine d4_dgemv312(amat, xvec, yvec, alpha, beta, trans) real(dp), intent(in), contiguous, target :: amat(:, :, :) real(dp), intent(in) :: xvec(:) real(dp), intent(inout), contiguous, target :: yvec(:, :) real(dp), intent(in), optional :: alpha real(dp), intent(in), optional :: beta character(len=1), intent(in), optional :: trans real(dp), pointer :: aptr(:, :), yptr(:) character(len=1) :: tra if (present(trans)) then tra = trans else tra = 'n' end if if (any(tra == ['n', 'N'])) then aptr(1:size(amat, 1)*size(amat, 2), 1:size(amat, 3)) => amat yptr(1:size(yvec, 1)*size(yvec, 2)) => yvec else aptr(1:size(amat, 1), 1:size(amat, 2)*size(amat, 3)) => amat yptr(1:size(yvec, 1) * size(yvec, 2)) => yvec end if call d4_gemv(aptr, xvec, yptr, alpha, beta, tra) end subroutine d4_dgemv312 subroutine d4_dgemv321(amat, xvec, yvec, alpha, beta, trans) real(dp), intent(in), contiguous, target :: amat(:, :, :) real(dp), intent(in), contiguous, target :: xvec(:, :) real(dp), intent(inout) :: yvec(:) real(dp), intent(in), optional :: alpha real(dp), intent(in), optional :: beta character(len=1), intent(in), optional :: trans real(dp), pointer :: aptr(:, :), xptr(:) character(len=1) :: tra if (present(trans)) then tra = trans else tra = 'n' end if if (any(tra == ['n', 'N'])) then aptr(1:size(amat, 1), 1:size(amat, 2)*size(amat, 3)) => amat xptr(1:size(xvec, 1)*size(xvec, 2)) => xvec else aptr(1:size(amat, 1)*size(amat, 2), 1:size(amat, 3)) => amat xptr(1:size(xvec, 1) * size(xvec, 2)) => xvec end if call d4_gemv(aptr, xptr, yvec, alpha, beta, tra) end subroutine d4_dgemv321 pure subroutine d4_sgemv(amat, xvec, yvec, alpha, beta, trans) real(sp), intent(in) :: amat(:, :) real(sp), intent(in) :: xvec(:) real(sp), intent(inout) :: yvec(:) real(sp), intent(in), optional :: alpha real(sp), intent(in), optional :: beta character(len=1), intent(in), optional :: trans real(sp) :: a, b character(len=1) :: tra integer(ik) :: incx, incy, m, n, lda if (present(alpha)) then a = alpha else a = 1.0_sp end if if (present(beta)) then b = beta else b = 0 end if if (present(trans)) then tra = trans else tra = 'n' end if incx = 1 incy = 1 lda = max(1, size(amat, 1)) m = size(amat, 1) n = size(amat, 2) call blas_gemv(tra, m, n, a, amat, lda, xvec, incx, b, yvec, incy) end subroutine d4_sgemv pure subroutine d4_dgemv(amat, xvec, yvec, alpha, beta, trans) real(dp), intent(in) :: amat(:, :) real(dp), intent(in) :: xvec(:) real(dp), intent(inout) :: yvec(:) real(dp), intent(in), optional :: alpha real(dp), intent(in), optional :: beta character(len=1), intent(in), optional :: trans real(dp) :: a, b character(len=1) :: tra integer(ik) :: incx, incy, m, n, lda if (present(alpha)) then a = alpha else a = 1.0_dp end if if (present(beta)) then b = beta else b = 0 end if if (present(trans)) then tra = trans else tra = 'n' end if incx = 1 incy = 1 lda = max(1, size(amat, 1)) m = size(amat, 1) n = size(amat, 2) call blas_gemv(tra, m, n, a, amat, lda, xvec, incx, b, yvec, incy) end subroutine d4_dgemv end module dftd4_blas