#:include "common.fypp"
#:for rk,rt,ri in REAL_KINDS_TYPES
#:if not rk in ["sp","dp"]
module stdlib_linalg_blas_${ri}$
use stdlib_linalg_constants
use stdlib_linalg_blas_aux
use stdlib_linalg_blas_s
use stdlib_linalg_blas_c
use stdlib_linalg_blas_d
use stdlib_linalg_blas_z
implicit none(type,external)
private
public :: sp,dp,${rk}$,lk,ilp
public :: stdlib_${ri}$asum
public :: stdlib_${ri}$axpy
public :: stdlib_${ri}$copy
public :: stdlib_${ri}$dot
public :: stdlib_${ri}$gbmv
public :: stdlib_${ri}$gemm
public :: stdlib_${ri}$gemv
public :: stdlib_${ri}$ger
public :: stdlib_${ri}$nrm2
public :: stdlib_${ri}$rot
public :: stdlib_${ri}$rotg
public :: stdlib_${ri}$rotm
public :: stdlib_${ri}$rotmg
public :: stdlib_${ri}$sbmv
public :: stdlib_${ri}$scal
public :: stdlib_${ri}$sdot
public :: stdlib_${ri}$spmv
public :: stdlib_${ri}$spr
public :: stdlib_${ri}$spr2
public :: stdlib_${ri}$swap
public :: stdlib_${ri}$symm
public :: stdlib_${ri}$symv
public :: stdlib_${ri}$syr
public :: stdlib_${ri}$syr2
public :: stdlib_${ri}$syr2k
public :: stdlib_${ri}$syrk
public :: stdlib_${ri}$tbmv
public :: stdlib_${ri}$tbsv
public :: stdlib_${ri}$tpmv
public :: stdlib_${ri}$tpsv
public :: stdlib_${ri}$trmm
public :: stdlib_${ri}$trmv
public :: stdlib_${ri}$trsm
public :: stdlib_${ri}$trsv
public :: stdlib_${ri}$zasum
public :: stdlib_${ri}$znrm2
! 128-bit real constants
real(${rk}$), parameter, private :: negone = -1.00_${rk}$
real(${rk}$), parameter, private :: zero = 0.00_${rk}$
real(${rk}$), parameter, private :: half = 0.50_${rk}$
real(${rk}$), parameter, private :: one = 1.00_${rk}$
real(${rk}$), parameter, private :: two = 2.00_${rk}$
real(${rk}$), parameter, private :: three = 3.00_${rk}$
real(${rk}$), parameter, private :: four = 4.00_${rk}$
real(${rk}$), parameter, private :: eight = 8.00_${rk}$
real(${rk}$), parameter, private :: ten = 10.00_${rk}$
! 128-bit complex constants
complex(${rk}$), parameter, private :: czero = ( 0.0_${rk}$,0.0_${rk}$)
complex(${rk}$), parameter, private :: chalf = ( 0.5_${rk}$,0.0_${rk}$)
complex(${rk}$), parameter, private :: cone = ( 1.0_${rk}$,0.0_${rk}$)
complex(${rk}$), parameter, private :: cnegone = (-1.0_${rk}$,0.0_${rk}$)
! 128-bit scaling constants
integer, parameter, private :: maxexp = maxexponent(zero)
integer, parameter, private :: minexp = minexponent(zero)
real(${rk}$), parameter, private :: rradix = real(radix(zero),${rk}$)
real(${rk}$), parameter, private :: ulp = epsilon(zero)
real(${rk}$), parameter, private :: eps = ulp*half
real(${rk}$), parameter, private :: safmin = rradix**max(minexp-1,1-maxexp)
real(${rk}$), parameter, private :: safmax = one/safmin
real(${rk}$), parameter, private :: smlnum = safmin/ulp
real(${rk}$), parameter, private :: bignum = safmax*ulp
real(${rk}$), parameter, private :: rtmin = sqrt(smlnum)
real(${rk}$), parameter, private :: rtmax = sqrt(bignum)
! 128-bit Blue's scaling constants
! ssml>=1/s and sbig==1/S with s,S as defined in https://doi.org/10.1145/355769.355771
real(${rk}$), parameter, private :: tsml = rradix**ceiling((minexp-1)*half)
real(${rk}$), parameter, private :: tbig = rradix**floor((maxexp-digits(zero)+1)*half)
real(${rk}$), parameter, private :: ssml = rradix**(-floor((minexp-digits(zero))*half))
real(${rk}$), parameter, private :: sbig = rradix**(-ceiling((maxexp+digits(zero)-1)*half))
contains
pure real(${rk}$) function stdlib_${ri}$asum(n,dx,incx)
!! DASUM: takes the sum of the absolute values.
! -- reference blas level1 routine --
! -- reference blas is a software package provided by univ. of tennessee, --
! -- univ. of california berkeley, univ. of colorado denver and nag ltd..--
! Scalar Arguments
integer(ilp), intent(in) :: incx, n
! Array Arguments
real(${rk}$), intent(in) :: dx(*)
! =====================================================================
! Local Scalars
real(${rk}$) :: dtemp
integer(ilp) :: i, m, mp1, nincx
! Intrinsic Functions
intrinsic :: abs,mod
stdlib_${ri}$asum = zero
dtemp = zero
if (n<=0 .or. incx<=0) return
if (incx==1) then
! code for increment equal to 1
! clean-up loop
m = mod(n,6)
if (m/=0) then
do i = 1,m
dtemp = dtemp + abs(dx(i))
end do
if (n<6) then
stdlib_${ri}$asum = dtemp
return
end if
end if
mp1 = m + 1
do i = mp1,n,6
dtemp = dtemp + abs(dx(i)) + abs(dx(i+1)) +abs(dx(i+2)) + abs(dx(i+3)) +abs(dx(i+&
4)) + abs(dx(i+5))
end do
else
! code for increment not equal to 1
nincx = n*incx
do i = 1,nincx,incx
dtemp = dtemp + abs(dx(i))
end do
end if
stdlib_${ri}$asum = dtemp
return
end function stdlib_${ri}$asum
pure subroutine stdlib_${ri}$axpy(n,da,dx,incx,dy,incy)
!! DAXPY: constant times a vector plus a vector.
!! uses unrolled loops for increments equal to one.
! -- reference blas level1 routine --
! -- reference blas is a software package provided by univ. of tennessee, --
! -- univ. of california berkeley, univ. of colorado denver and nag ltd..--
! Scalar Arguments
real(${rk}$), intent(in) :: da
integer(ilp), intent(in) :: incx, incy, n
! Array Arguments
real(${rk}$), intent(in) :: dx(*)
real(${rk}$), intent(inout) :: dy(*)
! =====================================================================
! Local Scalars
integer(ilp) :: i, ix, iy, m, mp1
! Intrinsic Functions
intrinsic :: mod
if (n<=0) return
if (da==0.0_${rk}$) return
if (incx==1 .and. incy==1) then
! code for both increments equal to 1
! clean-up loop
m = mod(n,4)
if (m/=0) then
do i = 1,m
dy(i) = dy(i) + da*dx(i)
end do
end if
if (n<4) return
mp1 = m + 1
do i = mp1,n,4
dy(i) = dy(i) + da*dx(i)
dy(i+1) = dy(i+1) + da*dx(i+1)
dy(i+2) = dy(i+2) + da*dx(i+2)
dy(i+3) = dy(i+3) + da*dx(i+3)
end do
else
! code for unequal increments or equal increments
! not equal to 1
ix = 1
iy = 1
if (incx<0) ix = (-n+1)*incx + 1
if (incy<0) iy = (-n+1)*incy + 1
do i = 1,n
dy(iy) = dy(iy) + da*dx(ix)
ix = ix + incx
iy = iy + incy
end do
end if
return
end subroutine stdlib_${ri}$axpy
pure subroutine stdlib_${ri}$copy(n,dx,incx,dy,incy)
!! DCOPY: copies a vector, x, to a vector, y.
!! uses unrolled loops for increments equal to 1.
! -- reference blas level1 routine --
! -- reference blas is a software package provided by univ. of tennessee, --
! -- univ. of california berkeley, univ. of colorado denver and nag ltd..--
! Scalar Arguments
integer(ilp), intent(in) :: incx, incy, n
! Array Arguments
real(${rk}$), intent(in) :: dx(*)
real(${rk}$), intent(out) :: dy(*)
! =====================================================================
! Local Scalars
integer(ilp) :: i, ix, iy, m, mp1
! Intrinsic Functions
intrinsic :: mod
if (n<=0) return
if (incx==1 .and. incy==1) then
! code for both increments equal to 1
! clean-up loop
m = mod(n,7)
if (m/=0) then
do i = 1,m
dy(i) = dx(i)
end do
if (n<7) return
end if
mp1 = m + 1
do i = mp1,n,7
dy(i) = dx(i)
dy(i+1) = dx(i+1)
dy(i+2) = dx(i+2)
dy(i+3) = dx(i+3)
dy(i+4) = dx(i+4)
dy(i+5) = dx(i+5)
dy(i+6) = dx(i+6)
end do
else
! code for unequal increments or equal increments
! not equal to 1
ix = 1
iy = 1
if (incx<0) ix = (-n+1)*incx + 1
if (incy<0) iy = (-n+1)*incy + 1
do i = 1,n
dy(iy) = dx(ix)
ix = ix + incx
iy = iy + incy
end do
end if
return
end subroutine stdlib_${ri}$copy
pure real(${rk}$) function stdlib_${ri}$dot(n,dx,incx,dy,incy)
!! DDOT: forms the dot product of two vectors.
!! uses unrolled loops for increments equal to one.
! -- reference blas level1 routine --
! -- reference blas is a software package provided by univ. of tennessee, --
! -- univ. of california berkeley, univ. of colorado denver and nag ltd..--
! Scalar Arguments
integer(ilp), intent(in) :: incx, incy, n
! Array Arguments
real(${rk}$), intent(in) :: dx(*), dy(*)
! =====================================================================
! Local Scalars
real(${rk}$) :: dtemp
integer(ilp) :: i, ix, iy, m, mp1
! Intrinsic Functions
intrinsic :: mod
stdlib_${ri}$dot = zero
dtemp = zero
if (n<=0) return
if (incx==1 .and. incy==1) then
! code for both increments equal to 1
! clean-up loop
m = mod(n,5)
if (m/=0) then
do i = 1,m
dtemp = dtemp + dx(i)*dy(i)
end do
if (n<5) then
stdlib_${ri}$dot=dtemp
return
end if
end if
mp1 = m + 1
do i = mp1,n,5
dtemp = dtemp + dx(i)*dy(i) + dx(i+1)*dy(i+1) +dx(i+2)*dy(i+2) + dx(i+3)*dy(i+3) + &
dx(i+4)*dy(i+4)
end do
else
! code for unequal increments or equal increments
! not equal to 1
ix = 1
iy = 1
if (incx<0) ix = (-n+1)*incx + 1
if (incy<0) iy = (-n+1)*incy + 1
do i = 1,n
dtemp = dtemp + dx(ix)*dy(iy)
ix = ix + incx
iy = iy + incy
end do
end if
stdlib_${ri}$dot = dtemp
return
end function stdlib_${ri}$dot
pure subroutine stdlib_${ri}$gbmv(trans,m,n,kl,ku,alpha,a,lda,x,incx,beta,y,incy)
!! DGBMV: 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 band matrix, with kl sub-diagonals and ku super-diagonals.
! -- reference blas level2 routine --
! -- reference blas is a software package provided by univ. of tennessee, --
! -- univ. of california berkeley, univ. of colorado denver and nag ltd..--
! Scalar Arguments
real(${rk}$), intent(in) :: alpha, beta
integer(ilp), intent(in) :: incx, incy, kl, ku, lda, m, n
character, intent(in) :: trans
! Array Arguments
real(${rk}$), intent(in) :: a(lda,*), x(*)
real(${rk}$), intent(inout) :: y(*)
! =====================================================================
! Local Scalars
real(${rk}$) :: temp
integer(ilp) :: i, info, ix, iy, j, jx, jy, k, kup1, kx, ky, lenx, leny
! Intrinsic Functions
intrinsic :: max,min
! test the input parameters.
info = 0
if (.not.stdlib_lsame(trans,'N') .and. .not.stdlib_lsame(trans,'T') &
.and..not.stdlib_lsame(trans,'C')) then
info = 1
else if (m<0) then
info = 2
else if (n<0) then
info = 3
else if (kl<0) then
info = 4
else if (ku<0) then
info = 5
else if (lda< (kl+ku+1)) then
info = 8
else if (incx==0) then
info = 10
else if (incy==0) then
info = 13
end if
if (info/=0) then
call stdlib_xerbla('DGBMV ',info)
return
end if
! quick return if possible.
if ((m==0) .or. (n==0) .or.((alpha==zero).and. (beta==one))) return
! set lenx and leny, the lengths of the vectors x and y, and set
! up the start points in x and y.
if (stdlib_lsame(trans,'N')) then
lenx = n
leny = m
else
lenx = m
leny = n
end if
if (incx>0) then
kx = 1
else
kx = 1 - (lenx-1)*incx
end if
if (incy>0) then
ky = 1
else
ky = 1 - (leny-1)*incy
end if
! start the operations. in this version the elements of a are
! accessed sequentially with one pass through the band part of a.
! first form y := beta*y.
if (beta/=one) then
if (incy==1) then
if (beta==zero) then
do i = 1,leny
y(i) = zero
end do
else
do i = 1,leny
y(i) = beta*y(i)
end do
end if
else
iy = ky
if (beta==zero) then
do i = 1,leny
y(iy) = zero
iy = iy + incy
end do
else
do i = 1,leny
y(iy) = beta*y(iy)
iy = iy + incy
end do
end if
end if
end if
if (alpha==zero) return
kup1 = ku + 1
if (stdlib_lsame(trans,'N')) then
! form y := alpha*a*x + y.
jx = kx
if (incy==1) then
do j = 1,n
temp = alpha*x(jx)
k = kup1 - j
do i = max(1,j-ku),min(m,j+kl)
y(i) = y(i) + temp*a(k+i,j)
end do
jx = jx + incx
end do
else
do j = 1,n
temp = alpha*x(jx)
iy = ky
k = kup1 - j
do i = max(1,j-ku),min(m,j+kl)
y(iy) = y(iy) + temp*a(k+i,j)
iy = iy + incy
end do
jx = jx + incx
if (j>ku) ky = ky + incy
end do
end if
else
! form y := alpha*a**t*x + y.
jy = ky
if (incx==1) then
do j = 1,n
temp = zero
k = kup1 - j
do i = max(1,j-ku),min(m,j+kl)
temp = temp + a(k+i,j)*x(i)
end do
y(jy) = y(jy) + alpha*temp
jy = jy + incy
end do
else
do j = 1,n
temp = zero
ix = kx
k = kup1 - j
do i = max(1,j-ku),min(m,j+kl)
temp = temp + a(k+i,j)*x(ix)
ix = ix + incx
end do
y(jy) = y(jy) + alpha*temp
jy = jy + incy
if (j>ku) kx = kx + incx
end do
end if
end if
return
end subroutine stdlib_${ri}$gbmv
pure subroutine stdlib_${ri}$gemm(transa,transb,m,n,k,alpha,a,lda,b,ldb,beta,c,ldc)
!! DGEMM: performs one of the matrix-matrix operations
!! C := alpha*op( A )*op( B ) + beta*C,
!! where op( X ) is one of
!! op( X ) = X or op( X ) = X**T,
!! alpha and beta are scalars, and A, B and C are matrices, with op( A )
!! an m by k matrix, op( B ) a k by n matrix and C an m by n matrix.
! -- reference blas level3 routine --
! -- reference blas is a software package provided by univ. of tennessee, --
! -- univ. of california berkeley, univ. of colorado denver and nag ltd..--
! Scalar Arguments
real(${rk}$), intent(in) :: alpha, beta
integer(ilp), intent(in) :: k, lda, ldb, ldc, m, n
character, intent(in) :: transa, transb
! Array Arguments
real(${rk}$), intent(in) :: a(lda,*), b(ldb,*)
real(${rk}$), intent(inout) :: c(ldc,*)
! =====================================================================
! Intrinsic Functions
intrinsic :: max
! Local Scalars
real(${rk}$) :: temp
integer(ilp) :: i, info, j, l, nrowa, nrowb
logical(lk) :: nota, notb
! set nota and notb as true if a and b respectively are not
! transposed and set nrowa and nrowb as the number of rows of a
! and b respectively.
nota = stdlib_lsame(transa,'N')
notb = stdlib_lsame(transb,'N')
if (nota) then
nrowa = m
else
nrowa = k
end if
if (notb) then
nrowb = k
else
nrowb = n
end if
! test the input parameters.
info = 0
if ((.not.nota) .and. (.not.stdlib_lsame(transa,'C')) .and.(.not.stdlib_lsame(transa,&
'T'))) then
info = 1
else if ((.not.notb) .and. (.not.stdlib_lsame(transb,'C')) .and.(.not.stdlib_lsame(&
transb,'T'))) then
info = 2
else if (m<0) then
info = 3
else if (n<0) then
info = 4
else if (k<0) then
info = 5
else if (lda<max(1,nrowa)) then
info = 8
else if (ldb<max(1,nrowb)) then
info = 10
else if (ldc<max(1,m)) then
info = 13
end if
if (info/=0) then
call stdlib_xerbla('DGEMM ',info)
return
end if
! quick return if possible.
if ((m==0) .or. (n==0) .or.(((alpha==zero).or. (k==0)).and. (beta==one))) &
return
! and if alpha.eq.zero.
if (alpha==zero) then
if (beta==zero) then
do j = 1,n
do i = 1,m
c(i,j) = zero
end do
end do
else
do j = 1,n
do i = 1,m
c(i,j) = beta*c(i,j)
end do
end do
end if
return
end if
! start the operations.
if (notb) then
if (nota) then
! form c := alpha*a*b + beta*c.
do j = 1,n
if (beta==zero) then
do i = 1,m
c(i,j) = zero
end do
else if (beta/=one) then
do i = 1,m
c(i,j) = beta*c(i,j)
end do
end if
do l = 1,k
temp = alpha*b(l,j)
do i = 1,m
c(i,j) = c(i,j) + temp*a(i,l)
end do
end do
end do
else
! form c := alpha*a**t*b + beta*c
do j = 1,n
do i = 1,m
temp = zero
do l = 1,k
temp = temp + a(l,i)*b(l,j)
end do
if (beta==zero) then
c(i,j) = alpha*temp
else
c(i,j) = alpha*temp + beta*c(i,j)
end if
end do
end do
end if
else
if (nota) then
! form c := alpha*a*b**t + beta*c
do j = 1,n
if (beta==zero) then
do i = 1,m
c(i,j) = zero
end do
else if (beta/=one) then
do i = 1,m
c(i,j) = beta*c(i,j)
end do
end if
do l = 1,k
temp = alpha*b(j,l)
do i = 1,m
c(i,j) = c(i,j) + temp*a(i,l)
end do
end do
end do
else
! form c := alpha*a**t*b**t + beta*c
do j = 1,n
do i = 1,m
temp = zero
do l = 1,k
temp = temp + a(l,i)*b(j,l)
end do
if (beta==zero) then
c(i,j) = alpha*temp
else
c(i,j) = alpha*temp + beta*c(i,j)
end if
end do
end do
end if
end if
return
end subroutine stdlib_${ri}$gemm
pure subroutine stdlib_${ri}$gemv(trans,m,n,alpha,a,lda,x,incx,beta,y,incy)
!! DGEMV: 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.
! -- reference blas level2 routine --
! -- reference blas is a software package provided by univ. of tennessee, --
! -- univ. of california berkeley, univ. of colorado denver and nag ltd..--
! Scalar Arguments
real(${rk}$), intent(in) :: alpha, beta
integer(ilp), intent(in) :: incx, incy, lda, m, n
character, intent(in) :: trans
! Array Arguments
real(${rk}$), intent(in) :: a(lda,*), x(*)
real(${rk}$), intent(inout) :: y(*)
! =====================================================================
! Local Scalars
real(${rk}$) :: temp
integer(ilp) :: i, info, ix, iy, j, jx, jy, kx, ky, lenx, leny
! Intrinsic Functions
intrinsic :: max
! test the input parameters.
info = 0
if (.not.stdlib_lsame(trans,'N') .and. .not.stdlib_lsame(trans,'T') &
.and..not.stdlib_lsame(trans,'C')) then
info = 1
else if (m<0) then
info = 2
else if (n<0) then
info = 3
else if (lda<max(1,m)) then
info = 6
else if (incx==0) then
info = 8
else if (incy==0) then
info = 11
end if
if (info/=0) then
call stdlib_xerbla('DGEMV ',info)
return
end if
! quick return if possible.
if ((m==0) .or. (n==0) .or.((alpha==zero).and. (beta==one))) return
! set lenx and leny, the lengths of the vectors x and y, and set
! up the start points in x and y.
if (stdlib_lsame(trans,'N')) then
lenx = n
leny = m
else
lenx = m
leny = n
end if
if (incx>0) then
kx = 1
else
kx = 1 - (lenx-1)*incx
end if
if (incy>0) then
ky = 1
else
ky = 1 - (leny-1)*incy
end if
! start the operations. in this version the elements of a are
! accessed sequentially with one pass through a.
! first form y := beta*y.
if (beta/=one) then
if (incy==1) then
if (beta==zero) then
do i = 1,leny
y(i) = zero
end do
else
do i = 1,leny
y(i) = beta*y(i)
end do
end if
else
iy = ky
if (beta==zero) then
do i = 1,leny
y(iy) = zero
iy = iy + incy
end do
else
do i = 1,leny
y(iy) = beta*y(iy)
iy = iy + incy
end do
end if
end if
end if
if (alpha==zero) return
if (stdlib_lsame(trans,'N')) then
! form y := alpha*a*x + y.
jx = kx
if (incy==1) then
do j = 1,n
temp = alpha*x(jx)
do i = 1,m
y(i) = y(i) + temp*a(i,j)
end do
jx = jx + incx
end do
else
do j = 1,n
temp = alpha*x(jx)
iy = ky
do i = 1,m
y(iy) = y(iy) + temp*a(i,j)
iy = iy + incy
end do
jx = jx + incx
end do
end if
else
! form y := alpha*a**t*x + y.
jy = ky
if (incx==1) then
do j = 1,n
temp = zero
do i = 1,m
temp = temp + a(i,j)*x(i)
end do
y(jy) = y(jy) + alpha*temp
jy = jy + incy
end do
else
do j = 1,n
temp = zero
ix = kx
do i = 1,m
temp = temp + a(i,j)*x(ix)
ix = ix + incx
end do
y(jy) = y(jy) + alpha*temp
jy = jy + incy
end do
end if
end if
return
end subroutine stdlib_${ri}$gemv
pure subroutine stdlib_${ri}$ger(m,n,alpha,x,incx,y,incy,a,lda)
!! DGER: performs the rank 1 operation
!! A := alpha*x*y**T + A,
!! where alpha is a scalar, x is an m element vector, y is an n element
!! vector and A is an m by n matrix.
! -- reference blas level2 routine --
! -- reference blas is a software package provided by univ. of tennessee, --
! -- univ. of california berkeley, univ. of colorado denver and nag ltd..--
! Scalar Arguments
real(${rk}$), intent(in) :: alpha
integer(ilp), intent(in) :: incx, incy, lda, m, n
! Array Arguments
real(${rk}$), intent(inout) :: a(lda,*)
real(${rk}$), intent(in) :: x(*), y(*)
! =====================================================================
! Local Scalars
real(${rk}$) :: temp
integer(ilp) :: i, info, ix, j, jy, kx
! Intrinsic Functions
intrinsic :: max
! test the input parameters.
info = 0
if (m<0) then
info = 1
else if (n<0) then
info = 2
else if (incx==0) then
info = 5
else if (incy==0) then
info = 7
else if (lda<max(1,m)) then
info = 9
end if
if (info/=0) then
call stdlib_xerbla('DGER ',info)
return
end if
! quick return if possible.
if ((m==0) .or. (n==0) .or. (alpha==zero)) return
! start the operations. in this version the elements of a are
! accessed sequentially with one pass through a.
if (incy>0) then
jy = 1
else
jy = 1 - (n-1)*incy
end if
if (incx==1) then
do j = 1,n
if (y(jy)/=zero) then
temp = alpha*y(jy)
do i = 1,m
a(i,j) = a(i,j) + x(i)*temp
end do
end if
jy = jy + incy
end do
else
if (incx>0) then
kx = 1
else
kx = 1 - (m-1)*incx
end if
do j = 1,n
if (y(jy)/=zero) then
temp = alpha*y(jy)
ix = kx
do i = 1,m
a(i,j) = a(i,j) + x(ix)*temp
ix = ix + incx
end do
end if
jy = jy + incy
end do
end if
return
end subroutine stdlib_${ri}$ger
pure function stdlib_${ri}$nrm2( n, x, incx )
!! DNRM2: returns the euclidean norm of a vector via the function
!! name, so that
!! DNRM2 := sqrt( x'*x )
real(${rk}$) :: stdlib_${ri}$nrm2
! -- reference blas level1 routine (version 3.9.1_${rk}$) --
! -- reference blas is a software package provided by univ. of tennessee, --
! -- univ. of california berkeley, univ. of colorado denver and nag ltd..--
! march 2021
! Constants
real(${rk}$), parameter :: maxn = huge(0.0_${rk}$)
! .. blue's scaling constants ..
! Scalar Arguments
integer(ilp), intent(in) :: incx, n
! Array Arguments
real(${rk}$), intent(in) :: x(*)
! Local Scalars
integer(ilp) :: i, ix
logical(lk) :: notbig
real(${rk}$) :: abig, amed, asml, ax, scl, sumsq, ymax, ymin
! quick return if possible
stdlib_${ri}$nrm2 = zero
if( n <= 0 ) return
scl = one
sumsq = zero
! compute the sum of squares in 3 accumulators:
! abig -- sums of squares scaled down to avoid overflow
! asml -- sums of squares scaled up to avoid underflow
! amed -- sums of squares that do not require scaling
! the thresholds and multipliers are
! tbig -- values bigger than this are scaled down by sbig
! tsml -- values smaller than this are scaled up by ssml
notbig = .true.
asml = zero
amed = zero
abig = zero
ix = 1
if( incx < 0 ) ix = 1 - (n-1)*incx
do i = 1, n
ax = abs(x(ix))
if (ax > tbig) then
abig = abig + (ax*sbig)**2
notbig = .false.
else if (ax < tsml) then
if (notbig) asml = asml + (ax*ssml)**2
else
amed = amed + ax**2
end if
ix = ix + incx
end do
! combine abig and amed or amed and asml if more than one
! accumulator was used.
if (abig > zero) then
! combine abig and amed if abig > 0.
if ( (amed > zero) .or. (amed > maxn) .or. (amed /= amed) ) then
abig = abig + (amed*sbig)*sbig
end if
scl = one / sbig
sumsq = abig
else if (asml > zero) then
! combine amed and asml if asml > 0.
if ( (amed > zero) .or. (amed > maxn) .or. (amed /= amed) ) then
amed = sqrt(amed)
asml = sqrt(asml) / ssml
if (asml > amed) then
ymin = amed
ymax = asml
else
ymin = asml
ymax = amed
end if
scl = one
sumsq = ymax**2*( one + (ymin/ymax)**2 )
else
scl = one / ssml
sumsq = asml
end if
else
! otherwise all values are mid-range
scl = one
sumsq = amed
end if
stdlib_${ri}$nrm2 = scl*sqrt( sumsq )
return
end function stdlib_${ri}$nrm2
pure subroutine stdlib_${ri}$rot(n,dx,incx,dy,incy,c,s)
!! DROT: applies a plane rotation.
! -- reference blas level1 routine --
! -- reference blas is a software package provided by univ. of tennessee, --
! -- univ. of california berkeley, univ. of colorado denver and nag ltd..--
! Scalar Arguments
real(${rk}$), intent(in) :: c, s
integer(ilp), intent(in) :: incx, incy, n
! Array Arguments
real(${rk}$), intent(inout) :: dx(*), dy(*)
! =====================================================================
! Local Scalars
real(${rk}$) :: dtemp
integer(ilp) :: i, ix, iy
if (n<=0) return
if (incx==1 .and. incy==1) then
! code for both increments equal to 1
do i = 1,n
dtemp = c*dx(i) + s*dy(i)
dy(i) = c*dy(i) - s*dx(i)
dx(i) = dtemp
end do
else
! code for unequal increments or equal increments not equal
! to 1
ix = 1
iy = 1
if (incx<0) ix = (-n+1)*incx + 1
if (incy<0) iy = (-n+1)*incy + 1
do i = 1,n
dtemp = c*dx(ix) + s*dy(iy)
dy(iy) = c*dy(iy) - s*dx(ix)
dx(ix) = dtemp
ix = ix + incx
iy = iy + incy
end do
end if
return
end subroutine stdlib_${ri}$rot
pure subroutine stdlib_${ri}$rotg( a, b, c, s )
!! The computation uses the formulas
!! sigma = sgn(a) if |a| > |b|
!! = sgn(b) if |b| >= |a|
!! r = sigma*sqrt( a**2 + b**2 )
!! c = 1; s = 0 if r = 0
!! c = a/r; s = b/r if r != 0
!! The subroutine also computes
!! z = s if |a| > |b|,
!! = 1/c if |b| >= |a| and c != 0
!! = 1 if c = 0
!! This allows c and s to be reconstructed from z as follows:
!! If z = 1, set c = 0, s = 1.
!! If |z| < 1, set c = sqrt(1 - z**2) and s = z.
!! If |z| > 1, set c = 1/z and s = sqrt( 1 - c**2).
! -- reference blas level1 routine --
! -- reference blas is a software package provided by univ. of tennessee, --
! -- univ. of california berkeley, univ. of colorado denver and nag ltd..--
! Scaling Constants
! Scalar Arguments
real(${rk}$), intent(inout) :: a, b
real(${rk}$), intent(out) :: c, s
! Local Scalars
real(${rk}$) :: anorm, bnorm, scl, sigma, r, z
anorm = abs(a)
bnorm = abs(b)
if( bnorm == zero ) then
c = one
s = zero
b = zero
else if( anorm == zero ) then
c = zero
s = one
a = b
b = one
else
scl = min( safmax, max( safmin, anorm, bnorm ) )
if( anorm > bnorm ) then
sigma = sign(one,a)
else
sigma = sign(one,b)
end if
r = sigma*( scl*sqrt((a/scl)**2 + (b/scl)**2) )
c = a/r
s = b/r
if( anorm > bnorm ) then
z = s
else if( c /= zero ) then
z = one/c
else
z = one
end if
a = r
b = z
end if
return
end subroutine stdlib_${ri}$rotg
pure subroutine stdlib_${ri}$rotm(n,dx,incx,dy,incy,dparam)
!! QROTM applies the modified Givens transformation, \(H\), to the 2-by-N matrix
!! $$ \left[ \begin{array}{c}DX^T\\DY^T\\ \end{array} \right], $$
!! where \(^T\) indicates transpose. The elements of \(DX\) are in
!! DX(LX+I*INCX), I = 0:N-1, where LX = 1 if INCX >= 0, else LX = (-INCX)*N,
!! and similarly for DY using LY and INCY.
!! With DPARAM(1)=DFLAG, \(H\) has one of the following forms:
!! $$ H=\underbrace{\begin{bmatrix}DH_{11} & DH_{12}\\DH_{21} & DH_{22}\end{bmatrix}}_{DFLAG=-1},
!! \underbrace{\begin{bmatrix}1 & DH_{12}\\DH_{21} & 1\end{bmatrix}}_{DFLAG=0},
!! \underbrace{\begin{bmatrix}DH_{11} & 1\\-1 & DH_{22}\end{bmatrix}}_{DFLAG=1},
!! \underbrace{\begin{bmatrix}1 & 0\\0 & 1\end{bmatrix}}_{DFLAG=-2}. $$
!! See QROTMG for a description of data storage in DPARAM.
! -- reference blas level1 routine --
! -- reference blas is a software package provided by univ. of tennessee, --
! -- univ. of california berkeley, univ. of colorado denver and nag ltd..--
! Scalar Arguments
integer(ilp), intent(in) :: incx, incy, n
! Array Arguments
real(${rk}$), intent(in) :: dparam(5)
real(${rk}$), intent(inout) :: dx(*), dy(*)
! =====================================================================
! Local Scalars
real(${rk}$) :: dflag, dh11, dh12, dh21, dh22, two, w, z, zero
integer(ilp) :: i, kx, ky, nsteps
! Data Statements
zero = 0.0_${rk}$
two = 2.0_${rk}$
dflag = dparam(1)
if (n<=0 .or. (dflag+two==zero)) return
if (incx==incy.and.incx>0) then
nsteps = n*incx
if (dflag<zero) then
dh11 = dparam(2)
dh12 = dparam(4)
dh21 = dparam(3)
dh22 = dparam(5)
do i = 1,nsteps,incx
w = dx(i)
z = dy(i)
dx(i) = w*dh11 + z*dh12
dy(i) = w*dh21 + z*dh22
end do
else if (dflag==zero) then
dh12 = dparam(4)
dh21 = dparam(3)
do i = 1,nsteps,incx
w = dx(i)
z = dy(i)
dx(i) = w + z*dh12
dy(i) = w*dh21 + z
end do
else
dh11 = dparam(2)
dh22 = dparam(5)
do i = 1,nsteps,incx
w = dx(i)
z = dy(i)
dx(i) = w*dh11 + z
dy(i) = -w + dh22*z
end do
end if
else
kx = 1
ky = 1
if (incx<0) kx = 1 + (1-n)*incx
if (incy<0) ky = 1 + (1-n)*incy
if (dflag<zero) then
dh11 = dparam(2)
dh12 = dparam(4)
dh21 = dparam(3)
dh22 = dparam(5)
do i = 1,n
w = dx(kx)
z = dy(ky)
dx(kx) = w*dh11 + z*dh12
dy(ky) = w*dh21 + z*dh22
kx = kx + incx
ky = ky + incy
end do
else if (dflag==zero) then
dh12 = dparam(4)
dh21 = dparam(3)
do i = 1,n
w = dx(kx)
z = dy(ky)
dx(kx) = w + z*dh12
dy(ky) = w*dh21 + z
kx = kx + incx
ky = ky + incy
end do
else
dh11 = dparam(2)
dh22 = dparam(5)
do i = 1,n
w = dx(kx)
z = dy(ky)
dx(kx) = w*dh11 + z
dy(ky) = -w + dh22*z
kx = kx + incx
ky = ky + incy
end do
end if
end if
return
end subroutine stdlib_${ri}$rotm
pure subroutine stdlib_${ri}$rotmg(dd1,dd2,dx1,dy1,dparam)
!! QROTMG Constructs the modified Givens transformation matrix \(H\) which zeros the
!! second component of the 2-vector
!! $$ \left[ {\sqrt{DD_1}\cdot DX_1,\sqrt{DD_2}\cdot DY_2} \right]^T. $$
!! With DPARAM(1)=DFLAG, \(H\) has one of the following forms:
!! $$ H=\underbrace{\begin{bmatrix}DH_{11} & DH_{12}\\DH_{21} & DH_{22}\end{bmatrix}}_{DFLAG=-1},
!! \underbrace{\begin{bmatrix}1 & DH_{12}\\DH_{21} & 1\end{bmatrix}}_{DFLAG=0},
!! \underbrace{\begin{bmatrix}DH_{11} & 1\\-1 & DH_{22}\end{bmatrix}}_{DFLAG=1},
!! \underbrace{\begin{bmatrix}1 & 0\\0 & 1\end{bmatrix}}_{DFLAG=-2}. $$
!! Locations 2-4 of DPARAM contain DH11, DH21, DH12 and DH22 respectively.
!! (Values of 1.0, -1.0, or 0.0 implied by the value of DPARAM(1) are not stored in DPARAM.)
!! The values of parameters GAMSQ and RGAMSQ may be inexact. This is OK as they are only
!! used for testing the size of DD1 and DD2. All actual scaling of data is done using GAM.
! -- reference blas level1 routine --
! -- reference blas is a software package provided by univ. of tennessee, --
! -- univ. of california berkeley, univ. of colorado denver and nag ltd..--
! Scalar Arguments
real(${rk}$), intent(inout) :: dd1, dd2, dx1
real(${rk}$), intent(in) :: dy1
! Array Arguments
real(${rk}$), intent(out) :: dparam(5)
! =====================================================================
! Local Scalars
real(${rk}$) :: dflag, dh11, dh12, dh21, dh22, dp1, dp2, dq1, dq2, dtemp, du, gam, gamsq, &
one, rgamsq, two, zero
! Intrinsic Functions
intrinsic :: abs
! Data Statements
zero = 0.0_${rk}$
one = 1.0_${rk}$
two = 2.0_${rk}$
gam = 4096.0_${rk}$
gamsq = 16777216.0_${rk}$
rgamsq = 5.9604645e-8_${rk}$
if (dd1<zero) then
! go zero-h-d-and-dx1..
dflag = -one
dh11 = zero
dh12 = zero
dh21 = zero
dh22 = zero
dd1 = zero
dd2 = zero
dx1 = zero
else
! case-dd1-nonnegative
dp2 = dd2*dy1
if (dp2==zero) then
dflag = -two
dparam(1) = dflag
return
end if
! regular-case..
dp1 = dd1*dx1
dq2 = dp2*dy1
dq1 = dp1*dx1
if (abs(dq1)>abs(dq2)) then
dh21 = -dy1/dx1
dh12 = dp2/dp1
du = one - dh12*dh21
if (du>zero) then
dflag = zero
dd1 = dd1/du
dd2 = dd2/du
dx1 = dx1*du
else
! this code path if here for safety. we do not expect this
! condition to ever hold except in edge cases with rounding
! errors. see doi: 10.1145/355841.355847
dflag = -one
dh11 = zero
dh12 = zero
dh21 = zero
dh22 = zero
dd1 = zero
dd2 = zero
dx1 = zero
end if
else
if (dq2<zero) then
! go zero-h-d-and-dx1..
dflag = -one
dh11 = zero
dh12 = zero
dh21 = zero
dh22 = zero
dd1 = zero
dd2 = zero
dx1 = zero
else
dflag = one
dh11 = dp1/dp2
dh22 = dx1/dy1
du = one + dh11*dh22
dtemp = dd2/du
dd2 = dd1/du
dd1 = dtemp
dx1 = dy1*du
end if
end if
! procedure..scale-check
if (dd1/=zero) then
do while ((dd1<=rgamsq) .or. (dd1>=gamsq))
if (dflag==zero) then
dh11 = one
dh22 = one
dflag = -one
else
dh21 = -one
dh12 = one
dflag = -one
end if
if (dd1<=rgamsq) then
dd1 = dd1*gam**2
dx1 = dx1/gam
dh11 = dh11/gam
dh12 = dh12/gam
else
dd1 = dd1/gam**2
dx1 = dx1*gam
dh11 = dh11*gam
dh12 = dh12*gam
end if
enddo
end if
if (dd2/=zero) then
do while ( (abs(dd2)<=rgamsq) .or. (abs(dd2)>=gamsq) )
if (dflag==zero) then
dh11 = one
dh22 = one
dflag = -one
else
dh21 = -one
dh12 = one
dflag = -one
end if
if (abs(dd2)<=rgamsq) then
dd2 = dd2*gam**2
dh21 = dh21/gam
dh22 = dh22/gam
else
dd2 = dd2/gam**2
dh21 = dh21*gam
dh22 = dh22*gam
end if
end do
end if
end if
if (dflag<zero) then
dparam(2) = dh11
dparam(3) = dh21
dparam(4) = dh12
dparam(5) = dh22
else if (dflag==zero) then
dparam(3) = dh21
dparam(4) = dh12
else
dparam(2) = dh11
dparam(5) = dh22
end if
dparam(1) = dflag
return
end subroutine stdlib_${ri}$rotmg
pure subroutine stdlib_${ri}$sbmv(uplo,n,k,alpha,a,lda,x,incx,beta,y,incy)
!! DSBMV: performs the matrix-vector operation
!! y := alpha*A*x + beta*y,
!! where alpha and beta are scalars, x and y are n element vectors and
!! A is an n by n symmetric band matrix, with k super-diagonals.
! -- reference blas level2 routine --
! -- reference blas is a software package provided by univ. of tennessee, --
! -- univ. of california berkeley, univ. of colorado denver and nag ltd..--
! Scalar Arguments
real(${rk}$), intent(in) :: alpha, beta
integer(ilp), intent(in) :: incx, incy, k, lda, n
character, intent(in) :: uplo
! Array Arguments
real(${rk}$), intent(in) :: a(lda,*), x(*)
real(${rk}$), intent(inout) :: y(*)
! =====================================================================
! Local Scalars
real(${rk}$) :: temp1, temp2
integer(ilp) :: i, info, ix, iy, j, jx, jy, kplus1, kx, ky, l
! Intrinsic Functions
intrinsic :: max,min
! test the input parameters.
info = 0
if (.not.stdlib_lsame(uplo,'U') .and. .not.stdlib_lsame(uplo,'L')) then
info = 1
else if (n<0) then
info = 2
else if (k<0) then
info = 3
else if (lda< (k+1)) then
info = 6
else if (incx==0) then
info = 8
else if (incy==0) then
info = 11
end if
if (info/=0) then
call stdlib_xerbla('DSBMV ',info)
return
end if
! quick return if possible.
if ((n==0) .or. ((alpha==zero).and. (beta==one))) return
! set up the start points in x and y.
if (incx>0) then
kx = 1
else
kx = 1 - (n-1)*incx
end if
if (incy>0) then
ky = 1
else
ky = 1 - (n-1)*incy
end if
! start the operations. in this version the elements of the array a
! are accessed sequentially with one pass through a.
! first form y := beta*y.
if (beta/=one) then
if (incy==1) then
if (beta==zero) then
do i = 1,n
y(i) = zero
end do
else
do i = 1,n
y(i) = beta*y(i)
end do
end if
else
iy = ky
if (beta==zero) then
do i = 1,n
y(iy) = zero
iy = iy + incy
end do
else
do i = 1,n
y(iy) = beta*y(iy)
iy = iy + incy
end do
end if
end if
end if
if (alpha==zero) return
if (stdlib_lsame(uplo,'U')) then
! form y when upper triangle of a is stored.
kplus1 = k + 1
if ((incx==1) .and. (incy==1)) then
do j = 1,n
temp1 = alpha*x(j)
temp2 = zero
l = kplus1 - j
do i = max(1,j-k),j - 1
y(i) = y(i) + temp1*a(l+i,j)
temp2 = temp2 + a(l+i,j)*x(i)
end do
y(j) = y(j) + temp1*a(kplus1,j) + alpha*temp2
end do
else
jx = kx
jy = ky
do j = 1,n
temp1 = alpha*x(jx)
temp2 = zero
ix = kx
iy = ky
l = kplus1 - j
do i = max(1,j-k),j - 1
y(iy) = y(iy) + temp1*a(l+i,j)
temp2 = temp2 + a(l+i,j)*x(ix)
ix = ix + incx
iy = iy + incy
end do
y(jy) = y(jy) + temp1*a(kplus1,j) + alpha*temp2
jx = jx + incx
jy = jy + incy
if (j>k) then
kx = kx + incx
ky = ky + incy
end if
end do
end if
else
! form y when lower triangle of a is stored.
if ((incx==1) .and. (incy==1)) then
do j = 1,n
temp1 = alpha*x(j)
temp2 = zero
y(j) = y(j) + temp1*a(1,j)
l = 1 - j
do i = j + 1,min(n,j+k)
y(i) = y(i) + temp1*a(l+i,j)
temp2 = temp2 + a(l+i,j)*x(i)
end do
y(j) = y(j) + alpha*temp2
end do
else
jx = kx
jy = ky
do j = 1,n
temp1 = alpha*x(jx)
temp2 = zero
y(jy) = y(jy) + temp1*a(1,j)
l = 1 - j
ix = jx
iy = jy
do i = j + 1,min(n,j+k)
ix = ix + incx
iy = iy + incy
y(iy) = y(iy) + temp1*a(l+i,j)
temp2 = temp2 + a(l+i,j)*x(ix)
end do
y(jy) = y(jy) + alpha*temp2
jx = jx + incx
jy = jy + incy
end do
end if
end if
return
end subroutine stdlib_${ri}$sbmv
pure subroutine stdlib_${ri}$scal(n,da,dx,incx)
!! DSCAL: scales a vector by a constant.
!! uses unrolled loops for increment equal to 1.
! -- reference blas level1 routine --
! -- reference blas is a software package provided by univ. of tennessee, --
! -- univ. of california berkeley, univ. of colorado denver and nag ltd..--
! Scalar Arguments
real(${rk}$), intent(in) :: da
integer(ilp), intent(in) :: incx, n
! Array Arguments
real(${rk}$), intent(inout) :: dx(*)
! =====================================================================
! Local Scalars
integer(ilp) :: i, m, mp1, nincx
! Intrinsic Functions
intrinsic :: mod
if (n<=0 .or. incx<=0) return
if (incx==1) then
! code for increment equal to 1
! clean-up loop
m = mod(n,5)
if (m/=0) then
do i = 1,m
dx(i) = da*dx(i)
end do
if (n<5) return
end if
mp1 = m + 1
do i = mp1,n,5
dx(i) = da*dx(i)
dx(i+1) = da*dx(i+1)
dx(i+2) = da*dx(i+2)
dx(i+3) = da*dx(i+3)
dx(i+4) = da*dx(i+4)
end do
else
! code for increment not equal to 1
nincx = n*incx
do i = 1,nincx,incx
dx(i) = da*dx(i)
end do
end if
return
end subroutine stdlib_${ri}$scal
pure real(${rk}$) function stdlib_${ri}$sdot(n,sx,incx,sy,incy)
!! Compute the inner product of two vectors with extended
!! precision accumulation and result.
!! Returns D.P. dot product accumulated in D.P., for S.P. SX and SY
!! DSDOT: = sum for I = 0 to N-1 of SX(LX+I*INCX) * SY(LY+I*INCY),
!! where LX = 1 if INCX >= 0, else LX = 1+(1-N)*INCX, and LY is
!! defined in a similar way using INCY.
! -- reference blas level1 routine --
! -- reference blas is a software package provided by univ. of tennessee, --
! -- univ. of california berkeley, univ. of colorado denver and nag ltd..--
! Scalar Arguments
integer(ilp), intent(in) :: incx, incy, n
! Array Arguments
real(dp), intent(in) :: sx(*), sy(*)
! authors:
! ========
! lawson, c. l., (jpl), hanson, r. j., (snla),
! kincaid, d. r., (u. of texas), krogh, f. t., (jpl)
! =====================================================================
! Local Scalars
integer(ilp) :: i, kx, ky, ns
! Intrinsic Functions
intrinsic :: real
stdlib_${ri}$sdot = zero
if (n<=0) return
if (incx==incy .and. incx>0) then
! code for equal, positive, non-unit increments.
ns = n*incx
do i = 1,ns,incx
stdlib_${ri}$sdot = stdlib_${ri}$sdot + real(sx(i),KIND=${rk}$)*real(sy(i),KIND=${rk}$)
end do
else
! code for unequal or nonpositive increments.
kx = 1
ky = 1
if (incx<0) kx = 1 + (1-n)*incx
if (incy<0) ky = 1 + (1-n)*incy
do i = 1,n
stdlib_${ri}$sdot = stdlib_${ri}$sdot + real(sx(kx),KIND=${rk}$)*real(sy(ky),KIND=${rk}$)
kx = kx + incx
ky = ky + incy
end do
end if
return
end function stdlib_${ri}$sdot
pure subroutine stdlib_${ri}$spmv(uplo,n,alpha,ap,x,incx,beta,y,incy)
!! DSPMV: performs the matrix-vector operation
!! y := alpha*A*x + beta*y,
!! where alpha and beta are scalars, x and y are n element vectors and
!! A is an n by n symmetric matrix, supplied in packed form.
! -- reference blas level2 routine --
! -- reference blas is a software package provided by univ. of tennessee, --
! -- univ. of california berkeley, univ. of colorado denver and nag ltd..--
! Scalar Arguments
real(${rk}$), intent(in) :: alpha, beta
integer(ilp), intent(in) :: incx, incy, n
character, intent(in) :: uplo
! Array Arguments
real(${rk}$), intent(in) :: ap(*), x(*)
real(${rk}$), intent(inout) :: y(*)
! =====================================================================
! Local Scalars
real(${rk}$) :: temp1, temp2
integer(ilp) :: i, info, ix, iy, j, jx, jy, k, kk, kx, ky
! test the input parameters.
info = 0
if (.not.stdlib_lsame(uplo,'U') .and. .not.stdlib_lsame(uplo,'L')) then
info = 1
else if (n<0) then
info = 2
else if (incx==0) then
info = 6
else if (incy==0) then
info = 9
end if
if (info/=0) then
call stdlib_xerbla('DSPMV ',info)
return
end if
! quick return if possible.
if ((n==0) .or. ((alpha==zero).and. (beta==one))) return
! set up the start points in x and y.
if (incx>0) then
kx = 1
else
kx = 1 - (n-1)*incx
end if
if (incy>0) then
ky = 1
else
ky = 1 - (n-1)*incy
end if
! start the operations. in this version the elements of the array ap
! are accessed sequentially with one pass through ap.
! first form y := beta*y.
if (beta/=one) then
if (incy==1) then
if (beta==zero) then
do i = 1,n
y(i) = zero
end do
else
do i = 1,n
y(i) = beta*y(i)
end do
end if
else
iy = ky
if (beta==zero) then
do i = 1,n
y(iy) = zero
iy = iy + incy
end do
else
do i = 1,n
y(iy) = beta*y(iy)
iy = iy + incy
end do
end if
end if
end if
if (alpha==zero) return
kk = 1
if (stdlib_lsame(uplo,'U')) then
! form y when ap contains the upper triangle.
if ((incx==1) .and. (incy==1)) then
do j = 1,n
temp1 = alpha*x(j)
temp2 = zero
k = kk
do i = 1,j - 1
y(i) = y(i) + temp1*ap(k)
temp2 = temp2 + ap(k)*x(i)
k = k + 1
end do
y(j) = y(j) + temp1*ap(kk+j-1) + alpha*temp2
kk = kk + j
end do
else
jx = kx
jy = ky
do j = 1,n
temp1 = alpha*x(jx)
temp2 = zero
ix = kx
iy = ky
do k = kk,kk + j - 2
y(iy) = y(iy) + temp1*ap(k)
temp2 = temp2 + ap(k)*x(ix)
ix = ix + incx
iy = iy + incy
end do
y(jy) = y(jy) + temp1*ap(kk+j-1) + alpha*temp2
jx = jx + incx
jy = jy + incy
kk = kk + j
end do
end if
else
! form y when ap contains the lower triangle.
if ((incx==1) .and. (incy==1)) then
do j = 1,n
temp1 = alpha*x(j)
temp2 = zero
y(j) = y(j) + temp1*ap(kk)
k = kk + 1
do i = j + 1,n
y(i) = y(i) + temp1*ap(k)
temp2 = temp2 + ap(k)*x(i)
k = k + 1
end do
y(j) = y(j) + alpha*temp2
kk = kk + (n-j+1)
end do
else
jx = kx
jy = ky
do j = 1,n
temp1 = alpha*x(jx)
temp2 = zero
y(jy) = y(jy) + temp1*ap(kk)
ix = jx
iy = jy
do k = kk + 1,kk + n - j
ix = ix + incx
iy = iy + incy
y(iy) = y(iy) + temp1*ap(k)
temp2 = temp2 + ap(k)*x(ix)
end do
y(jy) = y(jy) + alpha*temp2
jx = jx + incx
jy = jy + incy
kk = kk + (n-j+1)
end do
end if
end if
return
end subroutine stdlib_${ri}$spmv
pure subroutine stdlib_${ri}$spr(uplo,n,alpha,x,incx,ap)
!! DSPR: performs the symmetric rank 1 operation
!! A := alpha*x*x**T + A,
!! where alpha is a real scalar, x is an n element vector and A is an
!! n by n symmetric matrix, supplied in packed form.
! -- reference blas level2 routine --
! -- reference blas is a software package provided by univ. of tennessee, --
! -- univ. of california berkeley, univ. of colorado denver and nag ltd..--
! Scalar Arguments
real(${rk}$), intent(in) :: alpha
integer(ilp), intent(in) :: incx, n
character, intent(in) :: uplo
! Array Arguments
real(${rk}$), intent(inout) :: ap(*)
real(${rk}$), intent(in) :: x(*)
! =====================================================================
! Local Scalars
real(${rk}$) :: temp
integer(ilp) :: i, info, ix, j, jx, k, kk, kx
! test the input parameters.
info = 0
if (.not.stdlib_lsame(uplo,'U') .and. .not.stdlib_lsame(uplo,'L')) then
info = 1
else if (n<0) then
info = 2
else if (incx==0) then
info = 5
end if
if (info/=0) then
call stdlib_xerbla('DSPR ',info)
return
end if
! quick return if possible.
if ((n==0) .or. (alpha==zero)) return
! set the start point in x if the increment is not unity.
if (incx<=0) then
kx = 1 - (n-1)*incx
else if (incx/=1) then
kx = 1
end if
! start the operations. in this version the elements of the array ap
! are accessed sequentially with one pass through ap.
kk = 1
if (stdlib_lsame(uplo,'U')) then
! form a when upper triangle is stored in ap.
if (incx==1) then
do j = 1,n
if (x(j)/=zero) then
temp = alpha*x(j)
k = kk
do i = 1,j
ap(k) = ap(k) + x(i)*temp
k = k + 1
end do
end if
kk = kk + j
end do
else
jx = kx
do j = 1,n
if (x(jx)/=zero) then
temp = alpha*x(jx)
ix = kx
do k = kk,kk + j - 1
ap(k) = ap(k) + x(ix)*temp
ix = ix + incx
end do
end if
jx = jx + incx
kk = kk + j
end do
end if
else
! form a when lower triangle is stored in ap.
if (incx==1) then
do j = 1,n
if (x(j)/=zero) then
temp = alpha*x(j)
k = kk
do i = j,n
ap(k) = ap(k) + x(i)*temp
k = k + 1
end do
end if
kk = kk + n - j + 1
end do
else
jx = kx
do j = 1,n
if (x(jx)/=zero) then
temp = alpha*x(jx)
ix = jx
do k = kk,kk + n - j
ap(k) = ap(k) + x(ix)*temp
ix = ix + incx
end do
end if
jx = jx + incx
kk = kk + n - j + 1
end do
end if
end if
return
end subroutine stdlib_${ri}$spr
pure subroutine stdlib_${ri}$spr2(uplo,n,alpha,x,incx,y,incy,ap)
!! DSPR2: performs the symmetric rank 2 operation
!! A := alpha*x*y**T + alpha*y*x**T + A,
!! where alpha is a scalar, x and y are n element vectors and A is an
!! n by n symmetric matrix, supplied in packed form.
! -- reference blas level2 routine --
! -- reference blas is a software package provided by univ. of tennessee, --
! -- univ. of california berkeley, univ. of colorado denver and nag ltd..--
! Scalar Arguments
real(${rk}$), intent(in) :: alpha
integer(ilp), intent(in) :: incx, incy, n
character, intent(in) :: uplo
! Array Arguments
real(${rk}$), intent(inout) :: ap(*)
real(${rk}$), intent(in) :: x(*), y(*)
! =====================================================================
! Local Scalars
real(${rk}$) :: temp1, temp2
integer(ilp) :: i, info, ix, iy, j, jx, jy, k, kk, kx, ky
! test the input parameters.
info = 0
if (.not.stdlib_lsame(uplo,'U') .and. .not.stdlib_lsame(uplo,'L')) then
info = 1
else if (n<0) then
info = 2
else if (incx==0) then
info = 5
else if (incy==0) then
info = 7
end if
if (info/=0) then
call stdlib_xerbla('DSPR2 ',info)
return
end if
! quick return if possible.
if ((n==0) .or. (alpha==zero)) return
! set up the start points in x and y if the increments are not both
! unity.
if ((incx/=1) .or. (incy/=1)) then
if (incx>0) then
kx = 1
else
kx = 1 - (n-1)*incx
end if
if (incy>0) then
ky = 1
else
ky = 1 - (n-1)*incy
end if
jx = kx
jy = ky
end if
! start the operations. in this version the elements of the array ap
! are accessed sequentially with one pass through ap.
kk = 1
if (stdlib_lsame(uplo,'U')) then
! form a when upper triangle is stored in ap.
if ((incx==1) .and. (incy==1)) then
do j = 1,n
if ((x(j)/=zero) .or. (y(j)/=zero)) then
temp1 = alpha*y(j)
temp2 = alpha*x(j)
k = kk
do i = 1,j
ap(k) = ap(k) + x(i)*temp1 + y(i)*temp2
k = k + 1
end do
end if
kk = kk + j
end do
else
do j = 1,n
if ((x(jx)/=zero) .or. (y(jy)/=zero)) then
temp1 = alpha*y(jy)
temp2 = alpha*x(jx)
ix = kx
iy = ky
do k = kk,kk + j - 1
ap(k) = ap(k) + x(ix)*temp1 + y(iy)*temp2
ix = ix + incx
iy = iy + incy
end do
end if
jx = jx + incx
jy = jy + incy
kk = kk + j
end do
end if
else
! form a when lower triangle is stored in ap.
if ((incx==1) .and. (incy==1)) then
do j = 1,n
if ((x(j)/=zero) .or. (y(j)/=zero)) then
temp1 = alpha*y(j)
temp2 = alpha*x(j)
k = kk
do i = j,n
ap(k) = ap(k) + x(i)*temp1 + y(i)*temp2
k = k + 1
end do
end if
kk = kk + n - j + 1
end do
else
do j = 1,n
if ((x(jx)/=zero) .or. (y(jy)/=zero)) then
temp1 = alpha*y(jy)
temp2 = alpha*x(jx)
ix = jx
iy = jy
do k = kk,kk + n - j
ap(k) = ap(k) + x(ix)*temp1 + y(iy)*temp2
ix = ix + incx
iy = iy + incy
end do
end if
jx = jx + incx
jy = jy + incy
kk = kk + n - j + 1
end do
end if
end if
return
end subroutine stdlib_${ri}$spr2
pure subroutine stdlib_${ri}$swap(n,dx,incx,dy,incy)
!! DSWAP: interchanges two vectors.
!! uses unrolled loops for increments equal to 1.
! -- reference blas level1 routine --
! -- reference blas is a software package provided by univ. of tennessee, --
! -- univ. of california berkeley, univ. of colorado denver and nag ltd..--
! Scalar Arguments
integer(ilp), intent(in) :: incx, incy, n
! Array Arguments
real(${rk}$), intent(inout) :: dx(*), dy(*)
! =====================================================================
! Local Scalars
real(${rk}$) :: dtemp
integer(ilp) :: i, ix, iy, m, mp1
! Intrinsic Functions
intrinsic :: mod
if (n<=0) return
if (incx==1 .and. incy==1) then
! code for both increments equal to 1
! clean-up loop
m = mod(n,3)
if (m/=0) then
do i = 1,m
dtemp = dx(i)
dx(i) = dy(i)
dy(i) = dtemp
end do
if (n<3) return
end if
mp1 = m + 1
do i = mp1,n,3
dtemp = dx(i)
dx(i) = dy(i)
dy(i) = dtemp
dtemp = dx(i+1)
dx(i+1) = dy(i+1)
dy(i+1) = dtemp
dtemp = dx(i+2)
dx(i+2) = dy(i+2)
dy(i+2) = dtemp
end do
else
! code for unequal increments or equal increments not equal
! to 1
ix = 1
iy = 1
if (incx<0) ix = (-n+1)*incx + 1
if (incy<0) iy = (-n+1)*incy + 1
do i = 1,n
dtemp = dx(ix)
dx(ix) = dy(iy)
dy(iy) = dtemp
ix = ix + incx
iy = iy + incy
end do
end if
return
end subroutine stdlib_${ri}$swap
pure subroutine stdlib_${ri}$symm(side,uplo,m,n,alpha,a,lda,b,ldb,beta,c,ldc)
!! DSYMM: performs one of the matrix-matrix operations
!! C := alpha*A*B + beta*C,
!! or
!! C := alpha*B*A + beta*C,
!! where alpha and beta are scalars, A is a symmetric matrix and B and
!! C are m by n matrices.
! -- reference blas level3 routine --
! -- reference blas is a software package provided by univ. of tennessee, --
! -- univ. of california berkeley, univ. of colorado denver and nag ltd..--
! Scalar Arguments
real(${rk}$), intent(in) :: alpha, beta
integer(ilp), intent(in) :: lda, ldb, ldc, m, n
character, intent(in) :: side, uplo
! Array Arguments
real(${rk}$), intent(in) :: a(lda,*), b(ldb,*)
real(${rk}$), intent(inout) :: c(ldc,*)
! =====================================================================
! Intrinsic Functions
intrinsic :: max
! Local Scalars
real(${rk}$) :: temp1, temp2
integer(ilp) :: i, info, j, k, nrowa
logical(lk) :: upper
! set nrowa as the number of rows of a.
if (stdlib_lsame(side,'L')) then
nrowa = m
else
nrowa = n
end if
upper = stdlib_lsame(uplo,'U')
! test the input parameters.
info = 0
if ((.not.stdlib_lsame(side,'L')) .and. (.not.stdlib_lsame(side,'R'))) then
info = 1
else if ((.not.upper) .and. (.not.stdlib_lsame(uplo,'L'))) then
info = 2
else if (m<0) then
info = 3
else if (n<0) then
info = 4
else if (lda<max(1,nrowa)) then
info = 7
else if (ldb<max(1,m)) then
info = 9
else if (ldc<max(1,m)) then
info = 12
end if
if (info/=0) then
call stdlib_xerbla('DSYMM ',info)
return
end if
! quick return if possible.
if ((m==0) .or. (n==0) .or.((alpha==zero).and. (beta==one))) return
! and when alpha.eq.zero.
if (alpha==zero) then
if (beta==zero) then
do j = 1,n
do i = 1,m
c(i,j) = zero
end do
end do
else
do j = 1,n
do i = 1,m
c(i,j) = beta*c(i,j)
end do
end do
end if
return
end if
! start the operations.
if (stdlib_lsame(side,'L')) then
! form c := alpha*a*b + beta*c.
if (upper) then
do j = 1,n
do i = 1,m
temp1 = alpha*b(i,j)
temp2 = zero
do k = 1,i - 1
c(k,j) = c(k,j) + temp1*a(k,i)
temp2 = temp2 + b(k,j)*a(k,i)
end do
if (beta==zero) then
c(i,j) = temp1*a(i,i) + alpha*temp2
else
c(i,j) = beta*c(i,j) + temp1*a(i,i) +alpha*temp2
end if
end do
end do
else
do j = 1,n
do i = m,1,-1
temp1 = alpha*b(i,j)
temp2 = zero
do k = i + 1,m
c(k,j) = c(k,j) + temp1*a(k,i)
temp2 = temp2 + b(k,j)*a(k,i)
end do
if (beta==zero) then
c(i,j) = temp1*a(i,i) + alpha*temp2
else
c(i,j) = beta*c(i,j) + temp1*a(i,i) +alpha*temp2
end if
end do
end do
end if
else
! form c := alpha*b*a + beta*c.
loop_170: do j = 1,n
temp1 = alpha*a(j,j)
if (beta==zero) then
do i = 1,m
c(i,j) = temp1*b(i,j)
end do
else
do i = 1,m
c(i,j) = beta*c(i,j) + temp1*b(i,j)
end do
end if
do k = 1,j - 1
if (upper) then
temp1 = alpha*a(k,j)
else
temp1 = alpha*a(j,k)
end if
do i = 1,m
c(i,j) = c(i,j) + temp1*b(i,k)
end do
end do
do k = j + 1,n
if (upper) then
temp1 = alpha*a(j,k)
else
temp1 = alpha*a(k,j)
end if
do i = 1,m
c(i,j) = c(i,j) + temp1*b(i,k)
end do
end do
end do loop_170
end if
return
end subroutine stdlib_${ri}$symm
pure subroutine stdlib_${ri}$symv(uplo,n,alpha,a,lda,x,incx,beta,y,incy)
!! DSYMV: performs the matrix-vector operation
!! y := alpha*A*x + beta*y,
!! where alpha and beta are scalars, x and y are n element vectors and
!! A is an n by n symmetric matrix.
! -- reference blas level2 routine --
! -- reference blas is a software package provided by univ. of tennessee, --
! -- univ. of california berkeley, univ. of colorado denver and nag ltd..--
! Scalar Arguments
real(${rk}$), intent(in) :: alpha, beta
integer(ilp), intent(in) :: incx, incy, lda, n
character, intent(in) :: uplo
! Array Arguments
real(${rk}$), intent(in) :: a(lda,*), x(*)
real(${rk}$), intent(inout) :: y(*)
! =====================================================================
! Local Scalars
real(${rk}$) :: temp1, temp2
integer(ilp) :: i, info, ix, iy, j, jx, jy, kx, ky
! Intrinsic Functions
intrinsic :: max
! test the input parameters.
info = 0
if (.not.stdlib_lsame(uplo,'U') .and. .not.stdlib_lsame(uplo,'L')) then
info = 1
else if (n<0) then
info = 2
else if (lda<max(1,n)) then
info = 5
else if (incx==0) then
info = 7
else if (incy==0) then
info = 10
end if
if (info/=0) then
call stdlib_xerbla('DSYMV ',info)
return
end if
! quick return if possible.
if ((n==0) .or. ((alpha==zero).and. (beta==one))) return
! set up the start points in x and y.
if (incx>0) then
kx = 1
else
kx = 1 - (n-1)*incx
end if
if (incy>0) then
ky = 1
else
ky = 1 - (n-1)*incy
end if
! start the operations. in this version the elements of a are
! accessed sequentially with one pass through the triangular part
! of a.
! first form y := beta*y.
if (beta/=one) then
if (incy==1) then
if (beta==zero) then
do i = 1,n
y(i) = zero
end do
else
do i = 1,n
y(i) = beta*y(i)
end do
end if
else
iy = ky
if (beta==zero) then
do i = 1,n
y(iy) = zero
iy = iy + incy
end do
else
do i = 1,n
y(iy) = beta*y(iy)
iy = iy + incy
end do
end if
end if
end if
if (alpha==zero) return
if (stdlib_lsame(uplo,'U')) then
! form y when a is stored in upper triangle.
if ((incx==1) .and. (incy==1)) then
do j = 1,n
temp1 = alpha*x(j)
temp2 = zero
do i = 1,j - 1
y(i) = y(i) + temp1*a(i,j)
temp2 = temp2 + a(i,j)*x(i)
end do
y(j) = y(j) + temp1*a(j,j) + alpha*temp2
end do
else
jx = kx
jy = ky
do j = 1,n
temp1 = alpha*x(jx)
temp2 = zero
ix = kx
iy = ky
do i = 1,j - 1
y(iy) = y(iy) + temp1*a(i,j)
temp2 = temp2 + a(i,j)*x(ix)
ix = ix + incx
iy = iy + incy
end do
y(jy) = y(jy) + temp1*a(j,j) + alpha*temp2
jx = jx + incx
jy = jy + incy
end do
end if
else
! form y when a is stored in lower triangle.
if ((incx==1) .and. (incy==1)) then
do j = 1,n
temp1 = alpha*x(j)
temp2 = zero
y(j) = y(j) + temp1*a(j,j)
do i = j + 1,n
y(i) = y(i) + temp1*a(i,j)
temp2 = temp2 + a(i,j)*x(i)
end do
y(j) = y(j) + alpha*temp2
end do
else
jx = kx
jy = ky
do j = 1,n
temp1 = alpha*x(jx)
temp2 = zero
y(jy) = y(jy) + temp1*a(j,j)
ix = jx
iy = jy
do i = j + 1,n
ix = ix + incx
iy = iy + incy
y(iy) = y(iy) + temp1*a(i,j)
temp2 = temp2 + a(i,j)*x(ix)
end do
y(jy) = y(jy) + alpha*temp2
jx = jx + incx
jy = jy + incy
end do
end if
end if
return
end subroutine stdlib_${ri}$symv
pure subroutine stdlib_${ri}$syr(uplo,n,alpha,x,incx,a,lda)
!! DSYR: performs the symmetric rank 1 operation
!! A := alpha*x*x**T + A,
!! where alpha is a real scalar, x is an n element vector and A is an
!! n by n symmetric matrix.
! -- reference blas level2 routine --
! -- reference blas is a software package provided by univ. of tennessee, --
! -- univ. of california berkeley, univ. of colorado denver and nag ltd..--
! Scalar Arguments
real(${rk}$), intent(in) :: alpha
integer(ilp), intent(in) :: incx, lda, n
character, intent(in) :: uplo
! Array Arguments
real(${rk}$), intent(inout) :: a(lda,*)
real(${rk}$), intent(in) :: x(*)
! =====================================================================
! Local Scalars
real(${rk}$) :: temp
integer(ilp) :: i, info, ix, j, jx, kx
! Intrinsic Functions
intrinsic :: max
! test the input parameters.
info = 0
if (.not.stdlib_lsame(uplo,'U') .and. .not.stdlib_lsame(uplo,'L')) then
info = 1
else if (n<0) then
info = 2
else if (incx==0) then
info = 5
else if (lda<max(1,n)) then
info = 7
end if
if (info/=0) then
call stdlib_xerbla('DSYR ',info)
return
end if
! quick return if possible.
if ((n==0) .or. (alpha==zero)) return
! set the start point in x if the increment is not unity.
if (incx<=0) then
kx = 1 - (n-1)*incx
else if (incx/=1) then
kx = 1
end if
! start the operations. in this version the elements of a are
! accessed sequentially with one pass through the triangular part
! of a.
if (stdlib_lsame(uplo,'U')) then
! form a when a is stored in upper triangle.
if (incx==1) then
do j = 1,n
if (x(j)/=zero) then
temp = alpha*x(j)
do i = 1,j
a(i,j) = a(i,j) + x(i)*temp
end do
end if
end do
else
jx = kx
do j = 1,n
if (x(jx)/=zero) then
temp = alpha*x(jx)
ix = kx
do i = 1,j
a(i,j) = a(i,j) + x(ix)*temp
ix = ix + incx
end do
end if
jx = jx + incx
end do
end if
else
! form a when a is stored in lower triangle.
if (incx==1) then
do j = 1,n
if (x(j)/=zero) then
temp = alpha*x(j)
do i = j,n
a(i,j) = a(i,j) + x(i)*temp
end do
end if
end do
else
jx = kx
do j = 1,n
if (x(jx)/=zero) then
temp = alpha*x(jx)
ix = jx
do i = j,n
a(i,j) = a(i,j) + x(ix)*temp
ix = ix + incx
end do
end if
jx = jx + incx
end do
end if
end if
return
end subroutine stdlib_${ri}$syr
pure subroutine stdlib_${ri}$syr2(uplo,n,alpha,x,incx,y,incy,a,lda)
!! DSYR2: performs the symmetric rank 2 operation
!! A := alpha*x*y**T + alpha*y*x**T + A,
!! where alpha is a scalar, x and y are n element vectors and A is an n
!! by n symmetric matrix.
! -- reference blas level2 routine --
! -- reference blas is a software package provided by univ. of tennessee, --
! -- univ. of california berkeley, univ. of colorado denver and nag ltd..--
! Scalar Arguments
real(${rk}$), intent(in) :: alpha
integer(ilp), intent(in) :: incx, incy, lda, n
character, intent(in) :: uplo
! Array Arguments
real(${rk}$), intent(inout) :: a(lda,*)
real(${rk}$), intent(in) :: x(*), y(*)
! =====================================================================
! Local Scalars
real(${rk}$) :: temp1, temp2
integer(ilp) :: i, info, ix, iy, j, jx, jy, kx, ky
! Intrinsic Functions
intrinsic :: max
! test the input parameters.
info = 0
if (.not.stdlib_lsame(uplo,'U') .and. .not.stdlib_lsame(uplo,'L')) then
info = 1
else if (n<0) then
info = 2
else if (incx==0) then
info = 5
else if (incy==0) then
info = 7
else if (lda<max(1,n)) then
info = 9
end if
if (info/=0) then
call stdlib_xerbla('DSYR2 ',info)
return
end if
! quick return if possible.
if ((n==0) .or. (alpha==zero)) return
! set up the start points in x and y if the increments are not both
! unity.
if ((incx/=1) .or. (incy/=1)) then
if (incx>0) then
kx = 1
else
kx = 1 - (n-1)*incx
end if
if (incy>0) then
ky = 1
else
ky = 1 - (n-1)*incy
end if
jx = kx
jy = ky
end if
! start the operations. in this version the elements of a are
! accessed sequentially with one pass through the triangular part
! of a.
if (stdlib_lsame(uplo,'U')) then
! form a when a is stored in the upper triangle.
if ((incx==1) .and. (incy==1)) then
do j = 1,n
if ((x(j)/=zero) .or. (y(j)/=zero)) then
temp1 = alpha*y(j)
temp2 = alpha*x(j)
do i = 1,j
a(i,j) = a(i,j) + x(i)*temp1 + y(i)*temp2
end do
end if
end do
else
do j = 1,n
if ((x(jx)/=zero) .or. (y(jy)/=zero)) then
temp1 = alpha*y(jy)
temp2 = alpha*x(jx)
ix = kx
iy = ky
do i = 1,j
a(i,j) = a(i,j) + x(ix)*temp1 + y(iy)*temp2
ix = ix + incx
iy = iy + incy
end do
end if
jx = jx + incx
jy = jy + incy
end do
end if
else
! form a when a is stored in the lower triangle.
if ((incx==1) .and. (incy==1)) then
do j = 1,n
if ((x(j)/=zero) .or. (y(j)/=zero)) then
temp1 = alpha*y(j)
temp2 = alpha*x(j)
do i = j,n
a(i,j) = a(i,j) + x(i)*temp1 + y(i)*temp2
end do
end if
end do
else
do j = 1,n
if ((x(jx)/=zero) .or. (y(jy)/=zero)) then
temp1 = alpha*y(jy)
temp2 = alpha*x(jx)
ix = jx
iy = jy
do i = j,n
a(i,j) = a(i,j) + x(ix)*temp1 + y(iy)*temp2
ix = ix + incx
iy = iy + incy
end do
end if
jx = jx + incx
jy = jy + incy
end do
end if
end if
return
end subroutine stdlib_${ri}$syr2
pure subroutine stdlib_${ri}$syr2k(uplo,trans,n,k,alpha,a,lda,b,ldb,beta,c,ldc)
!! DSYR2K: performs one of the symmetric rank 2k operations
!! C := alpha*A*B**T + alpha*B*A**T + beta*C,
!! or
!! C := alpha*A**T*B + alpha*B**T*A + beta*C,
!! where alpha and beta are scalars, C is an n by n symmetric matrix
!! and A and B are n by k matrices in the first case and k by n
!! matrices in the second case.
! -- reference blas level3 routine --
! -- reference blas is a software package provided by univ. of tennessee, --
! -- univ. of california berkeley, univ. of colorado denver and nag ltd..--
! Scalar Arguments
real(${rk}$), intent(in) :: alpha, beta
integer(ilp), intent(in) :: k, lda, ldb, ldc, n
character, intent(in) :: trans, uplo
! Array Arguments
real(${rk}$), intent(in) :: a(lda,*), b(ldb,*)
real(${rk}$), intent(inout) :: c(ldc,*)
! =====================================================================
! Intrinsic Functions
intrinsic :: max
! Local Scalars
real(${rk}$) :: temp1, temp2
integer(ilp) :: i, info, j, l, nrowa
logical(lk) :: upper
! test the input parameters.
if (stdlib_lsame(trans,'N')) then
nrowa = n
else
nrowa = k
end if
upper = stdlib_lsame(uplo,'U')
info = 0
if ((.not.upper) .and. (.not.stdlib_lsame(uplo,'L'))) then
info = 1
else if ((.not.stdlib_lsame(trans,'N')) .and.(.not.stdlib_lsame(trans,'T')) .and.(&
.not.stdlib_lsame(trans,'C'))) then
info = 2
else if (n<0) then
info = 3
else if (k<0) then
info = 4
else if (lda<max(1,nrowa)) then
info = 7
else if (ldb<max(1,nrowa)) then
info = 9
else if (ldc<max(1,n)) then
info = 12
end if
if (info/=0) then
call stdlib_xerbla('DSYR2K',info)
return
end if
! quick return if possible.
if ((n==0) .or. (((alpha==zero).or.(k==0)).and. (beta==one))) return
! and when alpha.eq.zero.
if (alpha==zero) then
if (upper) then
if (beta==zero) then
do j = 1,n
do i = 1,j
c(i,j) = zero
end do
end do
else
do j = 1,n
do i = 1,j
c(i,j) = beta*c(i,j)
end do
end do
end if
else
if (beta==zero) then
do j = 1,n
do i = j,n
c(i,j) = zero
end do
end do
else
do j = 1,n
do i = j,n
c(i,j) = beta*c(i,j)
end do
end do
end if
end if
return
end if
! start the operations.
if (stdlib_lsame(trans,'N')) then
! form c := alpha*a*b**t + alpha*b*a**t + c.
if (upper) then
do j = 1,n
if (beta==zero) then
do i = 1,j
c(i,j) = zero
end do
else if (beta/=one) then
do i = 1,j
c(i,j) = beta*c(i,j)
end do
end if
do l = 1,k
if ((a(j,l)/=zero) .or. (b(j,l)/=zero)) then
temp1 = alpha*b(j,l)
temp2 = alpha*a(j,l)
do i = 1,j
c(i,j) = c(i,j) + a(i,l)*temp1 +b(i,l)*temp2
end do
end if
end do
end do
else
do j = 1,n
if (beta==zero) then
do i = j,n
c(i,j) = zero
end do
else if (beta/=one) then
do i = j,n
c(i,j) = beta*c(i,j)
end do
end if
do l = 1,k
if ((a(j,l)/=zero) .or. (b(j,l)/=zero)) then
temp1 = alpha*b(j,l)
temp2 = alpha*a(j,l)
do i = j,n
c(i,j) = c(i,j) + a(i,l)*temp1 +b(i,l)*temp2
end do
end if
end do
end do
end if
else
! form c := alpha*a**t*b + alpha*b**t*a + c.
if (upper) then
do j = 1,n
do i = 1,j
temp1 = zero
temp2 = zero
do l = 1,k
temp1 = temp1 + a(l,i)*b(l,j)
temp2 = temp2 + b(l,i)*a(l,j)
end do
if (beta==zero) then
c(i,j) = alpha*temp1 + alpha*temp2
else
c(i,j) = beta*c(i,j) + alpha*temp1 +alpha*temp2
end if
end do
end do
else
do j = 1,n
do i = j,n
temp1 = zero
temp2 = zero
do l = 1,k
temp1 = temp1 + a(l,i)*b(l,j)
temp2 = temp2 + b(l,i)*a(l,j)
end do
if (beta==zero) then
c(i,j) = alpha*temp1 + alpha*temp2
else
c(i,j) = beta*c(i,j) + alpha*temp1 +alpha*temp2
end if
end do
end do
end if
end if
return
end subroutine stdlib_${ri}$syr2k
pure subroutine stdlib_${ri}$syrk(uplo,trans,n,k,alpha,a,lda,beta,c,ldc)
!! DSYRK: performs one of the symmetric rank k operations
!! C := alpha*A*A**T + beta*C,
!! or
!! C := alpha*A**T*A + beta*C,
!! where alpha and beta are scalars, C is an n by n symmetric matrix
!! and A is an n by k matrix in the first case and a k by n matrix
!! in the second case.
! -- reference blas level3 routine --
! -- reference blas is a software package provided by univ. of tennessee, --
! -- univ. of california berkeley, univ. of colorado denver and nag ltd..--
! Scalar Arguments
real(${rk}$), intent(in) :: alpha, beta
integer(ilp), intent(in) :: k, lda, ldc, n
character, intent(in) :: trans, uplo
! Array Arguments
real(${rk}$), intent(in) :: a(lda,*)
real(${rk}$), intent(inout) :: c(ldc,*)
! =====================================================================
! Intrinsic Functions
intrinsic :: max
! Local Scalars
real(${rk}$) :: temp
integer(ilp) :: i, info, j, l, nrowa
logical(lk) :: upper
! test the input parameters.
if (stdlib_lsame(trans,'N')) then
nrowa = n
else
nrowa = k
end if
upper = stdlib_lsame(uplo,'U')
info = 0
if ((.not.upper) .and. (.not.stdlib_lsame(uplo,'L'))) then
info = 1
else if ((.not.stdlib_lsame(trans,'N')) .and.(.not.stdlib_lsame(trans,'T')) .and.(&
.not.stdlib_lsame(trans,'C'))) then
info = 2
else if (n<0) then
info = 3
else if (k<0) then
info = 4
else if (lda<max(1,nrowa)) then
info = 7
else if (ldc<max(1,n)) then
info = 10
end if
if (info/=0) then
call stdlib_xerbla('DSYRK ',info)
return
end if
! quick return if possible.
if ((n==0) .or. (((alpha==zero).or.(k==0)).and. (beta==one))) return
! and when alpha.eq.zero.
if (alpha==zero) then
if (upper) then
if (beta==zero) then
do j = 1,n
do i = 1,j
c(i,j) = zero
end do
end do
else
do j = 1,n
do i = 1,j
c(i,j) = beta*c(i,j)
end do
end do
end if
else
if (beta==zero) then
do j = 1,n
do i = j,n
c(i,j) = zero
end do
end do
else
do j = 1,n
do i = j,n
c(i,j) = beta*c(i,j)
end do
end do
end if
end if
return
end if
! start the operations.
if (stdlib_lsame(trans,'N')) then
! form c := alpha*a*a**t + beta*c.
if (upper) then
do j = 1,n
if (beta==zero) then
do i = 1,j
c(i,j) = zero
end do
else if (beta/=one) then
do i = 1,j
c(i,j) = beta*c(i,j)
end do
end if
do l = 1,k
if (a(j,l)/=zero) then
temp = alpha*a(j,l)
do i = 1,j
c(i,j) = c(i,j) + temp*a(i,l)
end do
end if
end do
end do
else
do j = 1,n
if (beta==zero) then
do i = j,n
c(i,j) = zero
end do
else if (beta/=one) then
do i = j,n
c(i,j) = beta*c(i,j)
end do
end if
do l = 1,k
if (a(j,l)/=zero) then
temp = alpha*a(j,l)
do i = j,n
c(i,j) = c(i,j) + temp*a(i,l)
end do
end if
end do
end do
end if
else
! form c := alpha*a**t*a + beta*c.
if (upper) then
do j = 1,n
do i = 1,j
temp = zero
do l = 1,k
temp = temp + a(l,i)*a(l,j)
end do
if (beta==zero) then
c(i,j) = alpha*temp
else
c(i,j) = alpha*temp + beta*c(i,j)
end if
end do
end do
else
do j = 1,n
do i = j,n
temp = zero
do l = 1,k
temp = temp + a(l,i)*a(l,j)
end do
if (beta==zero) then
c(i,j) = alpha*temp
else
c(i,j) = alpha*temp + beta*c(i,j)
end if
end do
end do
end if
end if
return
end subroutine stdlib_${ri}$syrk
pure subroutine stdlib_${ri}$tbmv(uplo,trans,diag,n,k,a,lda,x,incx)
!! DTBMV: performs one of the matrix-vector operations
!! x := A*x, or x := A**T*x,
!! where x is an n element vector and A is an n by n unit, or non-unit,
!! upper or lower triangular band matrix, with ( k + 1 ) diagonals.
! -- reference blas level2 routine --
! -- reference blas is a software package provided by univ. of tennessee, --
! -- univ. of california berkeley, univ. of colorado denver and nag ltd..--
! Scalar Arguments
integer(ilp), intent(in) :: incx, k, lda, n
character, intent(in) :: diag, trans, uplo
! Array Arguments
real(${rk}$), intent(in) :: a(lda,*)
real(${rk}$), intent(inout) :: x(*)
! =====================================================================
! Local Scalars
real(${rk}$) :: temp
integer(ilp) :: i, info, ix, j, jx, kplus1, kx, l
logical(lk) :: nounit
! Intrinsic Functions
intrinsic :: max,min
! test the input parameters.
info = 0
if (.not.stdlib_lsame(uplo,'U') .and. .not.stdlib_lsame(uplo,'L')) then
info = 1
else if (.not.stdlib_lsame(trans,'N') .and. .not.stdlib_lsame(trans,'T') &
.and..not.stdlib_lsame(trans,'C')) then
info = 2
else if (.not.stdlib_lsame(diag,'U') .and. .not.stdlib_lsame(diag,'N')) then
info = 3
else if (n<0) then
info = 4
else if (k<0) then
info = 5
else if (lda< (k+1)) then
info = 7
else if (incx==0) then
info = 9
end if
if (info/=0) then
call stdlib_xerbla('DTBMV ',info)
return
end if
! quick return if possible.
if (n==0) return
nounit = stdlib_lsame(diag,'N')
! set up the start point in x if the increment is not unity. this
! will be ( n - 1 )*incx too small for descending loops.
if (incx<=0) then
kx = 1 - (n-1)*incx
else if (incx/=1) then
kx = 1
end if
! start the operations. in this version the elements of a are
! accessed sequentially with one pass through a.
if (stdlib_lsame(trans,'N')) then
! form x := a*x.
if (stdlib_lsame(uplo,'U')) then
kplus1 = k + 1
if (incx==1) then
do j = 1,n
if (x(j)/=zero) then
temp = x(j)
l = kplus1 - j
do i = max(1,j-k),j - 1
x(i) = x(i) + temp*a(l+i,j)
end do
if (nounit) x(j) = x(j)*a(kplus1,j)
end if
end do
else
jx = kx
do j = 1,n
if (x(jx)/=zero) then
temp = x(jx)
ix = kx
l = kplus1 - j
do i = max(1,j-k),j - 1
x(ix) = x(ix) + temp*a(l+i,j)
ix = ix + incx
end do
if (nounit) x(jx) = x(jx)*a(kplus1,j)
end if
jx = jx + incx
if (j>k) kx = kx + incx
end do
end if
else
if (incx==1) then
do j = n,1,-1
if (x(j)/=zero) then
temp = x(j)
l = 1 - j
do i = min(n,j+k),j + 1,-1
x(i) = x(i) + temp*a(l+i,j)
end do
if (nounit) x(j) = x(j)*a(1,j)
end if
end do
else
kx = kx + (n-1)*incx
jx = kx
do j = n,1,-1
if (x(jx)/=zero) then
temp = x(jx)
ix = kx
l = 1 - j
do i = min(n,j+k),j + 1,-1
x(ix) = x(ix) + temp*a(l+i,j)
ix = ix - incx
end do
if (nounit) x(jx) = x(jx)*a(1,j)
end if
jx = jx - incx
if ((n-j)>=k) kx = kx - incx
end do
end if
end if
else
! form x := a**t*x.
if (stdlib_lsame(uplo,'U')) then
kplus1 = k + 1
if (incx==1) then
do j = n,1,-1
temp = x(j)
l = kplus1 - j
if (nounit) temp = temp*a(kplus1,j)
do i = j - 1,max(1,j-k),-1
temp = temp + a(l+i,j)*x(i)
end do
x(j) = temp
end do
else
kx = kx + (n-1)*incx
jx = kx
do j = n,1,-1
temp = x(jx)
kx = kx - incx
ix = kx
l = kplus1 - j
if (nounit) temp = temp*a(kplus1,j)
do i = j - 1,max(1,j-k),-1
temp = temp + a(l+i,j)*x(ix)
ix = ix - incx
end do
x(jx) = temp
jx = jx - incx
end do
end if
else
if (incx==1) then
do j = 1,n
temp = x(j)
l = 1 - j
if (nounit) temp = temp*a(1,j)
do i = j + 1,min(n,j+k)
temp = temp + a(l+i,j)*x(i)
end do
x(j) = temp
end do
else
jx = kx
do j = 1,n
temp = x(jx)
kx = kx + incx
ix = kx
l = 1 - j
if (nounit) temp = temp*a(1,j)
do i = j + 1,min(n,j+k)
temp = temp + a(l+i,j)*x(ix)
ix = ix + incx
end do
x(jx) = temp
jx = jx + incx
end do
end if
end if
end if
return
end subroutine stdlib_${ri}$tbmv
pure subroutine stdlib_${ri}$tbsv(uplo,trans,diag,n,k,a,lda,x,incx)
!! DTBSV: solves one of the systems of equations
!! A*x = b, or A**T*x = b,
!! where b and x are n element vectors and A is an n by n unit, or
!! non-unit, upper or lower triangular band matrix, with ( k + 1 )
!! diagonals.
!! No test for singularity or near-singularity is included in this
!! routine. Such tests must be performed before calling this routine.
! -- reference blas level2 routine --
! -- reference blas is a software package provided by univ. of tennessee, --
! -- univ. of california berkeley, univ. of colorado denver and nag ltd..--
! Scalar Arguments
integer(ilp), intent(in) :: incx, k, lda, n
character, intent(in) :: diag, trans, uplo
! Array Arguments
real(${rk}$), intent(in) :: a(lda,*)
real(${rk}$), intent(inout) :: x(*)
! =====================================================================
! Local Scalars
real(${rk}$) :: temp
integer(ilp) :: i, info, ix, j, jx, kplus1, kx, l
logical(lk) :: nounit
! Intrinsic Functions
intrinsic :: max,min
! test the input parameters.
info = 0
if (.not.stdlib_lsame(uplo,'U') .and. .not.stdlib_lsame(uplo,'L')) then
info = 1
else if (.not.stdlib_lsame(trans,'N') .and. .not.stdlib_lsame(trans,'T') &
.and..not.stdlib_lsame(trans,'C')) then
info = 2
else if (.not.stdlib_lsame(diag,'U') .and. .not.stdlib_lsame(diag,'N')) then
info = 3
else if (n<0) then
info = 4
else if (k<0) then
info = 5
else if (lda< (k+1)) then
info = 7
else if (incx==0) then
info = 9
end if
if (info/=0) then
call stdlib_xerbla('DTBSV ',info)
return
end if
! quick return if possible.
if (n==0) return
nounit = stdlib_lsame(diag,'N')
! set up the start point in x if the increment is not unity. this
! will be ( n - 1 )*incx too small for descending loops.
if (incx<=0) then
kx = 1 - (n-1)*incx
else if (incx/=1) then
kx = 1
end if
! start the operations. in this version the elements of a are
! accessed by sequentially with one pass through a.
if (stdlib_lsame(trans,'N')) then
! form x := inv( a )*x.
if (stdlib_lsame(uplo,'U')) then
kplus1 = k + 1
if (incx==1) then
do j = n,1,-1
if (x(j)/=zero) then
l = kplus1 - j
if (nounit) x(j) = x(j)/a(kplus1,j)
temp = x(j)
do i = j - 1,max(1,j-k),-1
x(i) = x(i) - temp*a(l+i,j)
end do
end if
end do
else
kx = kx + (n-1)*incx
jx = kx
do j = n,1,-1
kx = kx - incx
if (x(jx)/=zero) then
ix = kx
l = kplus1 - j
if (nounit) x(jx) = x(jx)/a(kplus1,j)
temp = x(jx)
do i = j - 1,max(1,j-k),-1
x(ix) = x(ix) - temp*a(l+i,j)
ix = ix - incx
end do
end if
jx = jx - incx
end do
end if
else
if (incx==1) then
do j = 1,n
if (x(j)/=zero) then
l = 1 - j
if (nounit) x(j) = x(j)/a(1,j)
temp = x(j)
do i = j + 1,min(n,j+k)
x(i) = x(i) - temp*a(l+i,j)
end do
end if
end do
else
jx = kx
do j = 1,n
kx = kx + incx
if (x(jx)/=zero) then
ix = kx
l = 1 - j
if (nounit) x(jx) = x(jx)/a(1,j)
temp = x(jx)
do i = j + 1,min(n,j+k)
x(ix) = x(ix) - temp*a(l+i,j)
ix = ix + incx
end do
end if
jx = jx + incx
end do
end if
end if
else
! form x := inv( a**t)*x.
if (stdlib_lsame(uplo,'U')) then
kplus1 = k + 1
if (incx==1) then
do j = 1,n
temp = x(j)
l = kplus1 - j
do i = max(1,j-k),j - 1
temp = temp - a(l+i,j)*x(i)
end do
if (nounit) temp = temp/a(kplus1,j)
x(j) = temp
end do
else
jx = kx
do j = 1,n
temp = x(jx)
ix = kx
l = kplus1 - j
do i = max(1,j-k),j - 1
temp = temp - a(l+i,j)*x(ix)
ix = ix + incx
end do
if (nounit) temp = temp/a(kplus1,j)
x(jx) = temp
jx = jx + incx
if (j>k) kx = kx + incx
end do
end if
else
if (incx==1) then
do j = n,1,-1
temp = x(j)
l = 1 - j
do i = min(n,j+k),j + 1,-1
temp = temp - a(l+i,j)*x(i)
end do
if (nounit) temp = temp/a(1,j)
x(j) = temp
end do
else
kx = kx + (n-1)*incx
jx = kx
do j = n,1,-1
temp = x(jx)
ix = kx
l = 1 - j
do i = min(n,j+k),j + 1,-1
temp = temp - a(l+i,j)*x(ix)
ix = ix - incx
end do
if (nounit) temp = temp/a(1,j)
x(jx) = temp
jx = jx - incx
if ((n-j)>=k) kx = kx - incx
end do
end if
end if
end if
return
end subroutine stdlib_${ri}$tbsv
pure subroutine stdlib_${ri}$tpmv(uplo,trans,diag,n,ap,x,incx)
!! DTPMV: performs one of the matrix-vector operations
!! x := A*x, or x := A**T*x,
!! where x is an n element vector and A is an n by n unit, or non-unit,
!! upper or lower triangular matrix, supplied in packed form.
! -- reference blas level2 routine --
! -- reference blas is a software package provided by univ. of tennessee, --
! -- univ. of california berkeley, univ. of colorado denver and nag ltd..--
! Scalar Arguments
integer(ilp), intent(in) :: incx, n
character, intent(in) :: diag, trans, uplo
! Array Arguments
real(${rk}$), intent(in) :: ap(*)
real(${rk}$), intent(inout) :: x(*)
! =====================================================================
! Local Scalars
real(${rk}$) :: temp
integer(ilp) :: i, info, ix, j, jx, k, kk, kx
logical(lk) :: nounit
! test the input parameters.
info = 0
if (.not.stdlib_lsame(uplo,'U') .and. .not.stdlib_lsame(uplo,'L')) then
info = 1
else if (.not.stdlib_lsame(trans,'N') .and. .not.stdlib_lsame(trans,'T') &
.and..not.stdlib_lsame(trans,'C')) then
info = 2
else if (.not.stdlib_lsame(diag,'U') .and. .not.stdlib_lsame(diag,'N')) then
info = 3
else if (n<0) then
info = 4
else if (incx==0) then
info = 7
end if
if (info/=0) then
call stdlib_xerbla('DTPMV ',info)
return
end if
! quick return if possible.
if (n==0) return
nounit = stdlib_lsame(diag,'N')
! set up the start point in x if the increment is not unity. this
! will be ( n - 1 )*incx too small for descending loops.
if (incx<=0) then
kx = 1 - (n-1)*incx
else if (incx/=1) then
kx = 1
end if
! start the operations. in this version the elements of ap are
! accessed sequentially with one pass through ap.
if (stdlib_lsame(trans,'N')) then
! form x:= a*x.
if (stdlib_lsame(uplo,'U')) then
kk = 1
if (incx==1) then
do j = 1,n
if (x(j)/=zero) then
temp = x(j)
k = kk
do i = 1,j - 1
x(i) = x(i) + temp*ap(k)
k = k + 1
end do
if (nounit) x(j) = x(j)*ap(kk+j-1)
end if
kk = kk + j
end do
else
jx = kx
do j = 1,n
if (x(jx)/=zero) then
temp = x(jx)
ix = kx
do k = kk,kk + j - 2
x(ix) = x(ix) + temp*ap(k)
ix = ix + incx
end do
if (nounit) x(jx) = x(jx)*ap(kk+j-1)
end if
jx = jx + incx
kk = kk + j
end do
end if
else
kk = (n* (n+1))/2
if (incx==1) then
do j = n,1,-1
if (x(j)/=zero) then
temp = x(j)
k = kk
do i = n,j + 1,-1
x(i) = x(i) + temp*ap(k)
k = k - 1
end do
if (nounit) x(j) = x(j)*ap(kk-n+j)
end if
kk = kk - (n-j+1)
end do
else
kx = kx + (n-1)*incx
jx = kx
do j = n,1,-1
if (x(jx)/=zero) then
temp = x(jx)
ix = kx
do k = kk,kk - (n- (j+1)),-1
x(ix) = x(ix) + temp*ap(k)
ix = ix - incx
end do
if (nounit) x(jx) = x(jx)*ap(kk-n+j)
end if
jx = jx - incx
kk = kk - (n-j+1)
end do
end if
end if
else
! form x := a**t*x.
if (stdlib_lsame(uplo,'U')) then
kk = (n* (n+1))/2
if (incx==1) then
do j = n,1,-1
temp = x(j)
if (nounit) temp = temp*ap(kk)
k = kk - 1
do i = j - 1,1,-1
temp = temp + ap(k)*x(i)
k = k - 1
end do
x(j) = temp
kk = kk - j
end do
else
jx = kx + (n-1)*incx
do j = n,1,-1
temp = x(jx)
ix = jx
if (nounit) temp = temp*ap(kk)
do k = kk - 1,kk - j + 1,-1
ix = ix - incx
temp = temp + ap(k)*x(ix)
end do
x(jx) = temp
jx = jx - incx
kk = kk - j
end do
end if
else
kk = 1
if (incx==1) then
do j = 1,n
temp = x(j)
if (nounit) temp = temp*ap(kk)
k = kk + 1
do i = j + 1,n
temp = temp + ap(k)*x(i)
k = k + 1
end do
x(j) = temp
kk = kk + (n-j+1)
end do
else
jx = kx
do j = 1,n
temp = x(jx)
ix = jx
if (nounit) temp = temp*ap(kk)
do k = kk + 1,kk + n - j
ix = ix + incx
temp = temp + ap(k)*x(ix)
end do
x(jx) = temp
jx = jx + incx
kk = kk + (n-j+1)
end do
end if
end if
end if
return
end subroutine stdlib_${ri}$tpmv
pure subroutine stdlib_${ri}$tpsv(uplo,trans,diag,n,ap,x,incx)
!! DTPSV: solves one of the systems of equations
!! A*x = b, or A**T*x = b,
!! where b and x are n element vectors and A is an n by n unit, or
!! non-unit, upper or lower triangular matrix, supplied in packed form.
!! No test for singularity or near-singularity is included in this
!! routine. Such tests must be performed before calling this routine.
! -- reference blas level2 routine --
! -- reference blas is a software package provided by univ. of tennessee, --
! -- univ. of california berkeley, univ. of colorado denver and nag ltd..--
! Scalar Arguments
integer(ilp), intent(in) :: incx, n
character, intent(in) :: diag, trans, uplo
! Array Arguments
real(${rk}$), intent(in) :: ap(*)
real(${rk}$), intent(inout) :: x(*)
! =====================================================================
! Local Scalars
real(${rk}$) :: temp
integer(ilp) :: i, info, ix, j, jx, k, kk, kx
logical(lk) :: nounit
! test the input parameters.
info = 0
if (.not.stdlib_lsame(uplo,'U') .and. .not.stdlib_lsame(uplo,'L')) then
info = 1
else if (.not.stdlib_lsame(trans,'N') .and. .not.stdlib_lsame(trans,'T') &
.and..not.stdlib_lsame(trans,'C')) then
info = 2
else if (.not.stdlib_lsame(diag,'U') .and. .not.stdlib_lsame(diag,'N')) then
info = 3
else if (n<0) then
info = 4
else if (incx==0) then
info = 7
end if
if (info/=0) then
call stdlib_xerbla('DTPSV ',info)
return
end if
! quick return if possible.
if (n==0) return
nounit = stdlib_lsame(diag,'N')
! set up the start point in x if the increment is not unity. this
! will be ( n - 1 )*incx too small for descending loops.
if (incx<=0) then
kx = 1 - (n-1)*incx
else if (incx/=1) then
kx = 1
end if
! start the operations. in this version the elements of ap are
! accessed sequentially with one pass through ap.
if (stdlib_lsame(trans,'N')) then
! form x := inv( a )*x.
if (stdlib_lsame(uplo,'U')) then
kk = (n* (n+1))/2
if (incx==1) then
do j = n,1,-1
if (x(j)/=zero) then
if (nounit) x(j) = x(j)/ap(kk)
temp = x(j)
k = kk - 1
do i = j - 1,1,-1
x(i) = x(i) - temp*ap(k)
k = k - 1
end do
end if
kk = kk - j
end do
else
jx = kx + (n-1)*incx
do j = n,1,-1
if (x(jx)/=zero) then
if (nounit) x(jx) = x(jx)/ap(kk)
temp = x(jx)
ix = jx
do k = kk - 1,kk - j + 1,-1
ix = ix - incx
x(ix) = x(ix) - temp*ap(k)
end do
end if
jx = jx - incx
kk = kk - j
end do
end if
else
kk = 1
if (incx==1) then
do j = 1,n
if (x(j)/=zero) then
if (nounit) x(j) = x(j)/ap(kk)
temp = x(j)
k = kk + 1
do i = j + 1,n
x(i) = x(i) - temp*ap(k)
k = k + 1
end do
end if
kk = kk + (n-j+1)
end do
else
jx = kx
do j = 1,n
if (x(jx)/=zero) then
if (nounit) x(jx) = x(jx)/ap(kk)
temp = x(jx)
ix = jx
do k = kk + 1,kk + n - j
ix = ix + incx
x(ix) = x(ix) - temp*ap(k)
end do
end if
jx = jx + incx
kk = kk + (n-j+1)
end do
end if
end if
else
! form x := inv( a**t )*x.
if (stdlib_lsame(uplo,'U')) then
kk = 1
if (incx==1) then
do j = 1,n
temp = x(j)
k = kk
do i = 1,j - 1
temp = temp - ap(k)*x(i)
k = k + 1
end do
if (nounit) temp = temp/ap(kk+j-1)
x(j) = temp
kk = kk + j
end do
else
jx = kx
do j = 1,n
temp = x(jx)
ix = kx
do k = kk,kk + j - 2
temp = temp - ap(k)*x(ix)
ix = ix + incx
end do
if (nounit) temp = temp/ap(kk+j-1)
x(jx) = temp
jx = jx + incx
kk = kk + j
end do
end if
else
kk = (n* (n+1))/2
if (incx==1) then
do j = n,1,-1
temp = x(j)
k = kk
do i = n,j + 1,-1
temp = temp - ap(k)*x(i)
k = k - 1
end do
if (nounit) temp = temp/ap(kk-n+j)
x(j) = temp
kk = kk - (n-j+1)
end do
else
kx = kx + (n-1)*incx
jx = kx
do j = n,1,-1
temp = x(jx)
ix = kx
do k = kk,kk - (n- (j+1)),-1
temp = temp - ap(k)*x(ix)
ix = ix - incx
end do
if (nounit) temp = temp/ap(kk-n+j)
x(jx) = temp
jx = jx - incx
kk = kk - (n-j+1)
end do
end if
end if
end if
return
end subroutine stdlib_${ri}$tpsv
pure subroutine stdlib_${ri}$trmm(side,uplo,transa,diag,m,n,alpha,a,lda,b,ldb)
!! DTRMM: performs one of the matrix-matrix operations
!! B := alpha*op( A )*B, or B := alpha*B*op( A ),
!! where alpha is a scalar, B is an m by n matrix, A is a unit, or
!! non-unit, upper or lower triangular matrix and op( A ) is one of
!! op( A ) = A or op( A ) = A**T.
! -- reference blas level3 routine --
! -- reference blas is a software package provided by univ. of tennessee, --
! -- univ. of california berkeley, univ. of colorado denver and nag ltd..--
! Scalar Arguments
real(${rk}$), intent(in) :: alpha
integer(ilp), intent(in) :: lda, ldb, m, n
character, intent(in) :: diag, side, transa, uplo
! Array Arguments
real(${rk}$), intent(in) :: a(lda,*)
real(${rk}$), intent(inout) :: b(ldb,*)
! =====================================================================
! Intrinsic Functions
intrinsic :: max
! Local Scalars
real(${rk}$) :: temp
integer(ilp) :: i, info, j, k, nrowa
logical(lk) :: lside, nounit, upper
! test the input parameters.
lside = stdlib_lsame(side,'L')
if (lside) then
nrowa = m
else
nrowa = n
end if
nounit = stdlib_lsame(diag,'N')
upper = stdlib_lsame(uplo,'U')
info = 0
if ((.not.lside) .and. (.not.stdlib_lsame(side,'R'))) then
info = 1
else if ((.not.upper) .and. (.not.stdlib_lsame(uplo,'L'))) then
info = 2
else if ((.not.stdlib_lsame(transa,'N')) .and.(.not.stdlib_lsame(transa,'T')) .and.(&
.not.stdlib_lsame(transa,'C'))) then
info = 3
else if ((.not.stdlib_lsame(diag,'U')) .and. (.not.stdlib_lsame(diag,'N'))) &
then
info = 4
else if (m<0) then
info = 5
else if (n<0) then
info = 6
else if (lda<max(1,nrowa)) then
info = 9
else if (ldb<max(1,m)) then
info = 11
end if
if (info/=0) then
call stdlib_xerbla('DTRMM ',info)
return
end if
! quick return if possible.
if (m==0 .or. n==0) return
! and when alpha.eq.zero.
if (alpha==zero) then
do j = 1,n
do i = 1,m
b(i,j) = zero
end do
end do
return
end if
! start the operations.
if (lside) then
if (stdlib_lsame(transa,'N')) then
! form b := alpha*a*b.
if (upper) then
do j = 1,n
do k = 1,m
if (b(k,j)/=zero) then
temp = alpha*b(k,j)
do i = 1,k - 1
b(i,j) = b(i,j) + temp*a(i,k)
end do
if (nounit) temp = temp*a(k,k)
b(k,j) = temp
end if
end do
end do
else
do j = 1,n
do k = m,1,-1
if (b(k,j)/=zero) then
temp = alpha*b(k,j)
b(k,j) = temp
if (nounit) b(k,j) = b(k,j)*a(k,k)
do i = k + 1,m
b(i,j) = b(i,j) + temp*a(i,k)
end do
end if
end do
end do
end if
else
! form b := alpha*a**t*b.
if (upper) then
do j = 1,n
do i = m,1,-1
temp = b(i,j)
if (nounit) temp = temp*a(i,i)
do k = 1,i - 1
temp = temp + a(k,i)*b(k,j)
end do
b(i,j) = alpha*temp
end do
end do
else
do j = 1,n
do i = 1,m
temp = b(i,j)
if (nounit) temp = temp*a(i,i)
do k = i + 1,m
temp = temp + a(k,i)*b(k,j)
end do
b(i,j) = alpha*temp
end do
end do
end if
end if
else
if (stdlib_lsame(transa,'N')) then
! form b := alpha*b*a.
if (upper) then
do j = n,1,-1
temp = alpha
if (nounit) temp = temp*a(j,j)
do i = 1,m
b(i,j) = temp*b(i,j)
end do
do k = 1,j - 1
if (a(k,j)/=zero) then
temp = alpha*a(k,j)
do i = 1,m
b(i,j) = b(i,j) + temp*b(i,k)
end do
end if
end do
end do
else
do j = 1,n
temp = alpha
if (nounit) temp = temp*a(j,j)
do i = 1,m
b(i,j) = temp*b(i,j)
end do
do k = j + 1,n
if (a(k,j)/=zero) then
temp = alpha*a(k,j)
do i = 1,m
b(i,j) = b(i,j) + temp*b(i,k)
end do
end if
end do
end do
end if
else
! form b := alpha*b*a**t.
if (upper) then
do k = 1,n
do j = 1,k - 1
if (a(j,k)/=zero) then
temp = alpha*a(j,k)
do i = 1,m
b(i,j) = b(i,j) + temp*b(i,k)
end do
end if
end do
temp = alpha
if (nounit) temp = temp*a(k,k)
if (temp/=one) then
do i = 1,m
b(i,k) = temp*b(i,k)
end do
end if
end do
else
do k = n,1,-1
do j = k + 1,n
if (a(j,k)/=zero) then
temp = alpha*a(j,k)
do i = 1,m
b(i,j) = b(i,j) + temp*b(i,k)
end do
end if
end do
temp = alpha
if (nounit) temp = temp*a(k,k)
if (temp/=one) then
do i = 1,m
b(i,k) = temp*b(i,k)
end do
end if
end do
end if
end if
end if
return
end subroutine stdlib_${ri}$trmm
pure subroutine stdlib_${ri}$trmv(uplo,trans,diag,n,a,lda,x,incx)
!! DTRMV: performs one of the matrix-vector operations
!! x := A*x, or x := A**T*x,
!! where x is an n element vector and A is an n by n unit, or non-unit,
!! upper or lower triangular matrix.
! -- reference blas level2 routine --
! -- reference blas is a software package provided by univ. of tennessee, --
! -- univ. of california berkeley, univ. of colorado denver and nag ltd..--
! Scalar Arguments
integer(ilp), intent(in) :: incx, lda, n
character, intent(in) :: diag, trans, uplo
! Array Arguments
real(${rk}$), intent(in) :: a(lda,*)
real(${rk}$), intent(inout) :: x(*)
! =====================================================================
! Local Scalars
real(${rk}$) :: temp
integer(ilp) :: i, info, ix, j, jx, kx
logical(lk) :: nounit
! Intrinsic Functions
intrinsic :: max
! test the input parameters.
info = 0
if (.not.stdlib_lsame(uplo,'U') .and. .not.stdlib_lsame(uplo,'L')) then
info = 1
else if (.not.stdlib_lsame(trans,'N') .and. .not.stdlib_lsame(trans,'T') &
.and..not.stdlib_lsame(trans,'C')) then
info = 2
else if (.not.stdlib_lsame(diag,'U') .and. .not.stdlib_lsame(diag,'N')) then
info = 3
else if (n<0) then
info = 4
else if (lda<max(1,n)) then
info = 6
else if (incx==0) then
info = 8
end if
if (info/=0) then
call stdlib_xerbla('DTRMV ',info)
return
end if
! quick return if possible.
if (n==0) return
nounit = stdlib_lsame(diag,'N')
! set up the start point in x if the increment is not unity. this
! will be ( n - 1 )*incx too small for descending loops.
if (incx<=0) then
kx = 1 - (n-1)*incx
else if (incx/=1) then
kx = 1
end if
! start the operations. in this version the elements of a are
! accessed sequentially with one pass through a.
if (stdlib_lsame(trans,'N')) then
! form x := a*x.
if (stdlib_lsame(uplo,'U')) then
if (incx==1) then
do j = 1,n
if (x(j)/=zero) then
temp = x(j)
do i = 1,j - 1
x(i) = x(i) + temp*a(i,j)
end do
if (nounit) x(j) = x(j)*a(j,j)
end if
end do
else
jx = kx
do j = 1,n
if (x(jx)/=zero) then
temp = x(jx)
ix = kx
do i = 1,j - 1
x(ix) = x(ix) + temp*a(i,j)
ix = ix + incx
end do
if (nounit) x(jx) = x(jx)*a(j,j)
end if
jx = jx + incx
end do
end if
else
if (incx==1) then
do j = n,1,-1
if (x(j)/=zero) then
temp = x(j)
do i = n,j + 1,-1
x(i) = x(i) + temp*a(i,j)
end do
if (nounit) x(j) = x(j)*a(j,j)
end if
end do
else
kx = kx + (n-1)*incx
jx = kx
do j = n,1,-1
if (x(jx)/=zero) then
temp = x(jx)
ix = kx
do i = n,j + 1,-1
x(ix) = x(ix) + temp*a(i,j)
ix = ix - incx
end do
if (nounit) x(jx) = x(jx)*a(j,j)
end if
jx = jx - incx
end do
end if
end if
else
! form x := a**t*x.
if (stdlib_lsame(uplo,'U')) then
if (incx==1) then
do j = n,1,-1
temp = x(j)
if (nounit) temp = temp*a(j,j)
do i = j - 1,1,-1
temp = temp + a(i,j)*x(i)
end do
x(j) = temp
end do
else
jx = kx + (n-1)*incx
do j = n,1,-1
temp = x(jx)
ix = jx
if (nounit) temp = temp*a(j,j)
do i = j - 1,1,-1
ix = ix - incx
temp = temp + a(i,j)*x(ix)
end do
x(jx) = temp
jx = jx - incx
end do
end if
else
if (incx==1) then
do j = 1,n
temp = x(j)
if (nounit) temp = temp*a(j,j)
do i = j + 1,n
temp = temp + a(i,j)*x(i)
end do
x(j) = temp
end do
else
jx = kx
do j = 1,n
temp = x(jx)
ix = jx
if (nounit) temp = temp*a(j,j)
do i = j + 1,n
ix = ix + incx
temp = temp + a(i,j)*x(ix)
end do
x(jx) = temp
jx = jx + incx
end do
end if
end if
end if
return
end subroutine stdlib_${ri}$trmv
pure subroutine stdlib_${ri}$trsm(side,uplo,transa,diag,m,n,alpha,a,lda,b,ldb)
!! DTRSM: solves one of the matrix equations
!! op( A )*X = alpha*B, or X*op( A ) = alpha*B,
!! where alpha is a scalar, X and B are m by n matrices, A is a unit, or
!! non-unit, upper or lower triangular matrix and op( A ) is one of
!! op( A ) = A or op( A ) = A**T.
!! The matrix X is overwritten on B.
! -- reference blas level3 routine --
! -- reference blas is a software package provided by univ. of tennessee, --
! -- univ. of california berkeley, univ. of colorado denver and nag ltd..--
! Scalar Arguments
real(${rk}$), intent(in) :: alpha
integer(ilp), intent(in) :: lda, ldb, m, n
character, intent(in) :: diag, side, transa, uplo
! Array Arguments
real(${rk}$), intent(in) :: a(lda,*)
real(${rk}$), intent(inout) :: b(ldb,*)
! =====================================================================
! Intrinsic Functions
intrinsic :: max
! Local Scalars
real(${rk}$) :: temp
integer(ilp) :: i, info, j, k, nrowa
logical(lk) :: lside, nounit, upper
! test the input parameters.
lside = stdlib_lsame(side,'L')
if (lside) then
nrowa = m
else
nrowa = n
end if
nounit = stdlib_lsame(diag,'N')
upper = stdlib_lsame(uplo,'U')
info = 0
if ((.not.lside) .and. (.not.stdlib_lsame(side,'R'))) then
info = 1
else if ((.not.upper) .and. (.not.stdlib_lsame(uplo,'L'))) then
info = 2
else if ((.not.stdlib_lsame(transa,'N')) .and.(.not.stdlib_lsame(transa,'T')) .and.(&
.not.stdlib_lsame(transa,'C'))) then
info = 3
else if ((.not.stdlib_lsame(diag,'U')) .and. (.not.stdlib_lsame(diag,'N'))) &
then
info = 4
else if (m<0) then
info = 5
else if (n<0) then
info = 6
else if (lda<max(1,nrowa)) then
info = 9
else if (ldb<max(1,m)) then
info = 11
end if
if (info/=0) then
call stdlib_xerbla('DTRSM ',info)
return
end if
! quick return if possible.
if (m==0 .or. n==0) return
! and when alpha.eq.zero.
if (alpha==zero) then
do j = 1,n
do i = 1,m
b(i,j) = zero
end do
end do
return
end if
! start the operations.
if (lside) then
if (stdlib_lsame(transa,'N')) then
! form b := alpha*inv( a )*b.
if (upper) then
do j = 1,n
if (alpha/=one) then
do i = 1,m
b(i,j) = alpha*b(i,j)
end do
end if
do k = m,1,-1
if (b(k,j)/=zero) then
if (nounit) b(k,j) = b(k,j)/a(k,k)
do i = 1,k - 1
b(i,j) = b(i,j) - b(k,j)*a(i,k)
end do
end if
end do
end do
else
do j = 1,n
if (alpha/=one) then
do i = 1,m
b(i,j) = alpha*b(i,j)
end do
end if
do k = 1,m
if (b(k,j)/=zero) then
if (nounit) b(k,j) = b(k,j)/a(k,k)
do i = k + 1,m
b(i,j) = b(i,j) - b(k,j)*a(i,k)
end do
end if
end do
end do
end if
else
! form b := alpha*inv( a**t )*b.
if (upper) then
do j = 1,n
do i = 1,m
temp = alpha*b(i,j)
do k = 1,i - 1
temp = temp - a(k,i)*b(k,j)
end do
if (nounit) temp = temp/a(i,i)
b(i,j) = temp
end do
end do
else
do j = 1,n
do i = m,1,-1
temp = alpha*b(i,j)
do k = i + 1,m
temp = temp - a(k,i)*b(k,j)
end do
if (nounit) temp = temp/a(i,i)
b(i,j) = temp
end do
end do
end if
end if
else
if (stdlib_lsame(transa,'N')) then
! form b := alpha*b*inv( a ).
if (upper) then
do j = 1,n
if (alpha/=one) then
do i = 1,m
b(i,j) = alpha*b(i,j)
end do
end if
do k = 1,j - 1
if (a(k,j)/=zero) then
do i = 1,m
b(i,j) = b(i,j) - a(k,j)*b(i,k)
end do
end if
end do
if (nounit) then
temp = one/a(j,j)
do i = 1,m
b(i,j) = temp*b(i,j)
end do
end if
end do
else
do j = n,1,-1
if (alpha/=one) then
do i = 1,m
b(i,j) = alpha*b(i,j)
end do
end if
do k = j + 1,n
if (a(k,j)/=zero) then
do i = 1,m
b(i,j) = b(i,j) - a(k,j)*b(i,k)
end do
end if
end do
if (nounit) then
temp = one/a(j,j)
do i = 1,m
b(i,j) = temp*b(i,j)
end do
end if
end do
end if
else
! form b := alpha*b*inv( a**t ).
if (upper) then
do k = n,1,-1
if (nounit) then
temp = one/a(k,k)
do i = 1,m
b(i,k) = temp*b(i,k)
end do
end if
do j = 1,k - 1
if (a(j,k)/=zero) then
temp = a(j,k)
do i = 1,m
b(i,j) = b(i,j) - temp*b(i,k)
end do
end if
end do
if (alpha/=one) then
do i = 1,m
b(i,k) = alpha*b(i,k)
end do
end if
end do
else
do k = 1,n
if (nounit) then
temp = one/a(k,k)
do i = 1,m
b(i,k) = temp*b(i,k)
end do
end if
do j = k + 1,n
if (a(j,k)/=zero) then
temp = a(j,k)
do i = 1,m
b(i,j) = b(i,j) - temp*b(i,k)
end do
end if
end do
if (alpha/=one) then
do i = 1,m
b(i,k) = alpha*b(i,k)
end do
end if
end do
end if
end if
end if
return
end subroutine stdlib_${ri}$trsm
pure subroutine stdlib_${ri}$trsv(uplo,trans,diag,n,a,lda,x,incx)
!! DTRSV: solves one of the systems of equations
!! A*x = b, or A**T*x = b,
!! where b and x are n element vectors and A is an n by n unit, or
!! non-unit, upper or lower triangular matrix.
!! No test for singularity or near-singularity is included in this
!! routine. Such tests must be performed before calling this routine.
! -- reference blas level1 routine --
! -- reference blas is a software package provided by univ. of tennessee, --
! -- univ. of california berkeley, univ. of colorado denver and nag ltd..--
! Scalar Arguments
integer(ilp), intent(in) :: incx, lda, n
character, intent(in) :: diag, trans, uplo
! Array Arguments
real(${rk}$), intent(in) :: a(lda,*)
real(${rk}$), intent(inout) :: x(*)
! =====================================================================
! Local Scalars
real(${rk}$) :: temp
integer(ilp) :: i, info, ix, j, jx, kx
logical(lk) :: nounit
! Intrinsic Functions
intrinsic :: max
! test the input parameters.
info = 0
if (.not.stdlib_lsame(uplo,'U') .and. .not.stdlib_lsame(uplo,'L')) then
info = 1
else if (.not.stdlib_lsame(trans,'N') .and. .not.stdlib_lsame(trans,'T') &
.and..not.stdlib_lsame(trans,'C')) then
info = 2
else if (.not.stdlib_lsame(diag,'U') .and. .not.stdlib_lsame(diag,'N')) then
info = 3
else if (n<0) then
info = 4
else if (lda<max(1,n)) then
info = 6
else if (incx==0) then
info = 8
end if
if (info/=0) then
call stdlib_xerbla('DTRSV ',info)
return
end if
! quick return if possible.
if (n==0) return
nounit = stdlib_lsame(diag,'N')
! set up the start point in x if the increment is not unity. this
! will be ( n - 1 )*incx too small for descending loops.
if (incx<=0) then
kx = 1 - (n-1)*incx
else if (incx/=1) then
kx = 1
end if
! start the operations. in this version the elements of a are
! accessed sequentially with one pass through a.
if (stdlib_lsame(trans,'N')) then
! form x := inv( a )*x.
if (stdlib_lsame(uplo,'U')) then
if (incx==1) then
do j = n,1,-1
if (x(j)/=zero) then
if (nounit) x(j) = x(j)/a(j,j)
temp = x(j)
do i = j - 1,1,-1
x(i) = x(i) - temp*a(i,j)
end do
end if
end do
else
jx = kx + (n-1)*incx
do j = n,1,-1
if (x(jx)/=zero) then
if (nounit) x(jx) = x(jx)/a(j,j)
temp = x(jx)
ix = jx
do i = j - 1,1,-1
ix = ix - incx
x(ix) = x(ix) - temp*a(i,j)
end do
end if
jx = jx - incx
end do
end if
else
if (incx==1) then
do j = 1,n
if (x(j)/=zero) then
if (nounit) x(j) = x(j)/a(j,j)
temp = x(j)
do i = j + 1,n
x(i) = x(i) - temp*a(i,j)
end do
end if
end do
else
jx = kx
do j = 1,n
if (x(jx)/=zero) then
if (nounit) x(jx) = x(jx)/a(j,j)
temp = x(jx)
ix = jx
do i = j + 1,n
ix = ix + incx
x(ix) = x(ix) - temp*a(i,j)
end do
end if
jx = jx + incx
end do
end if
end if
else
! form x := inv( a**t )*x.
if (stdlib_lsame(uplo,'U')) then
if (incx==1) then
do j = 1,n
temp = x(j)
do i = 1,j - 1
temp = temp - a(i,j)*x(i)
end do
if (nounit) temp = temp/a(j,j)
x(j) = temp
end do
else
jx = kx
do j = 1,n
temp = x(jx)
ix = kx
do i = 1,j - 1
temp = temp - a(i,j)*x(ix)
ix = ix + incx
end do
if (nounit) temp = temp/a(j,j)
x(jx) = temp
jx = jx + incx
end do
end if
else
if (incx==1) then
do j = n,1,-1
temp = x(j)
do i = n,j + 1,-1
temp = temp - a(i,j)*x(i)
end do
if (nounit) temp = temp/a(j,j)
x(j) = temp
end do
else
kx = kx + (n-1)*incx
jx = kx
do j = n,1,-1
temp = x(jx)
ix = kx
do i = n,j + 1,-1
temp = temp - a(i,j)*x(ix)
ix = ix - incx
end do
if (nounit) temp = temp/a(j,j)
x(jx) = temp
jx = jx - incx
end do
end if
end if
end if
return
end subroutine stdlib_${ri}$trsv
pure real(${rk}$) function stdlib_${ri}$zasum(n,zx,incx)
!! DZASUM: takes the sum of the (|Re(.)| + |Im(.)|)'s of a complex vector and
!! returns a quad precision result.
! -- reference blas level1 routine --
! -- reference blas is a software package provided by univ. of tennessee, --
! -- univ. of california berkeley, univ. of colorado denver and nag ltd..--
! Scalar Arguments
integer(ilp), intent(in) :: incx, n
! Array Arguments
complex(${rk}$), intent(in) :: zx(*)
! =====================================================================
! Local Scalars
real(${rk}$) :: stemp
integer(ilp) :: i, nincx
stdlib_${ri}$zasum = zero
stemp = zero
if (n<=0 .or. incx<=0) return
if (incx==1) then
! code for increment equal to 1
do i = 1,n
stemp = stemp + stdlib_cabs1(zx(i))
end do
else
! code for increment not equal to 1
nincx = n*incx
do i = 1,nincx,incx
stemp = stemp + stdlib_cabs1(zx(i))
end do
end if
stdlib_${ri}$zasum = stemp
return
end function stdlib_${ri}$zasum
pure function stdlib_${ri}$znrm2( n, x, incx )
!! DZNRM2: returns the euclidean norm of a vector via the function
!! name, so that
!! DZNRM2 := sqrt( x**H*x )
real(${rk}$) :: stdlib_${ri}$znrm2
! -- reference blas level1 routine (version 3.9.1_${rk}$) --
! -- reference blas is a software package provided by univ. of tennessee, --
! -- univ. of california berkeley, univ. of colorado denver and nag ltd..--
! march 2021
! Constants
real(${rk}$), parameter :: maxn = huge(0.0_${rk}$)
! .. blue's scaling constants ..
! Scalar Arguments
integer(ilp), intent(in) :: incx, n
! Array Arguments
complex(${rk}$), intent(in) :: x(*)
! Local Scalars
integer(ilp) :: i, ix
logical(lk) :: notbig
real(${rk}$) :: abig, amed, asml, ax, scl, sumsq, ymax, ymin
! quick return if possible
stdlib_${ri}$znrm2 = zero
if( n <= 0 ) return
scl = one
sumsq = zero
! compute the sum of squares in 3 accumulators:
! abig -- sums of squares scaled down to avoid overflow
! asml -- sums of squares scaled up to avoid underflow
! amed -- sums of squares that do not require scaling
! the thresholds and multipliers are
! tbig -- values bigger than this are scaled down by sbig
! tsml -- values smaller than this are scaled up by ssml
notbig = .true.
asml = zero
amed = zero
abig = zero
ix = 1
if( incx < 0 ) ix = 1 - (n-1)*incx
do i = 1, n
ax = abs(real(x(ix),KIND=${rk}$))
if (ax > tbig) then
abig = abig + (ax*sbig)**2
notbig = .false.
else if (ax < tsml) then
if (notbig) asml = asml + (ax*ssml)**2
else
amed = amed + ax**2
end if
ax = abs(aimag(x(ix)))
if (ax > tbig) then
abig = abig + (ax*sbig)**2
notbig = .false.
else if (ax < tsml) then
if (notbig) asml = asml + (ax*ssml)**2
else
amed = amed + ax**2
end if
ix = ix + incx
end do
! combine abig and amed or amed and asml if more than one
! accumulator was used.
if (abig > zero) then
! combine abig and amed if abig > 0.
if ( (amed > zero) .or. (amed > maxn) .or. (amed /= amed) ) then
abig = abig + (amed*sbig)*sbig
end if
scl = one / sbig
sumsq = abig
else if (asml > zero) then
! combine amed and asml if asml > 0.
if ( (amed > zero) .or. (amed > maxn) .or. (amed /= amed) ) then
amed = sqrt(amed)
asml = sqrt(asml) / ssml
if (asml > amed) then
ymin = amed
ymax = asml
else
ymin = asml
ymax = amed
end if
scl = one
sumsq = ymax**2*( one + (ymin/ymax)**2 )
else
scl = one / ssml
sumsq = asml
end if
else
! otherwise all values are mid-range
scl = one
sumsq = amed
end if
stdlib_${ri}$znrm2 = scl*sqrt( sumsq )
return
end function stdlib_${ri}$znrm2
end module stdlib_linalg_blas_${ri}$
#:endif
#:endfor
