stdlib_chpr Subroutine

public pure subroutine stdlib_chpr(uplo, n, alpha, x, incx, ap)

CHPR performs the hermitian rank 1 operation A := alphaxx**H + A, where alpha is a real scalar, x is an n element vector and A is an n by n hermitian matrix, supplied in packed form.

Arguments

Type IntentOptional Attributes Name
character(len=1), intent(in) :: uplo
integer(kind=ilp), intent(in) :: n
real(kind=sp), intent(in) :: alpha
complex(kind=sp), intent(in) :: x(*)
integer(kind=ilp), intent(in) :: incx
complex(kind=sp), intent(inout) :: ap(*)

Source Code

     pure subroutine stdlib_chpr(uplo,n,alpha,x,incx,ap)
     !! CHPR performs the hermitian rank 1 operation
     !! A := alpha*x*x**H + A,
     !! where alpha is a real scalar, x is an n element vector and A is an
     !! n by n hermitian 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(sp), intent(in) :: alpha
           integer(ilp), intent(in) :: incx, n
           character, intent(in) :: uplo
           ! Array Arguments 
           complex(sp), intent(inout) :: ap(*)
           complex(sp), intent(in) :: x(*)
        ! =====================================================================
           
           ! Local Scalars 
           complex(sp) :: temp
           integer(ilp) :: i, info, ix, j, jx, k, kk, kx
           ! Intrinsic Functions 
           intrinsic :: conjg,real
           ! 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('CHPR  ',info)
               return
           end if
           ! quick return if possible.
           if ((n==0) .or. (alpha==real(czero,KIND=sp))) 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 cone 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)/=czero) then
                           temp = alpha*conjg(x(j))
                           k = kk
                           do i = 1,j - 1
                               ap(k) = ap(k) + x(i)*temp
                               k = k + 1
                           end do
                           ap(kk+j-1) = real(ap(kk+j-1),KIND=sp) + real(x(j)*temp,KIND=sp)
                       else
                           ap(kk+j-1) = real(ap(kk+j-1),KIND=sp)
                       end if
                       kk = kk + j
                   end do
               else
                   jx = kx
                   do j = 1,n
                       if (x(jx)/=czero) then
                           temp = alpha*conjg(x(jx))
                           ix = kx
                           do k = kk,kk + j - 2
                               ap(k) = ap(k) + x(ix)*temp
                               ix = ix + incx
                           end do
                           ap(kk+j-1) = real(ap(kk+j-1),KIND=sp) + real(x(jx)*temp,KIND=sp)
                                     
                       else
                           ap(kk+j-1) = real(ap(kk+j-1),KIND=sp)
                       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)/=czero) then
                           temp = alpha*conjg(x(j))
                           ap(kk) = real(ap(kk),KIND=sp) + real(temp*x(j),KIND=sp)
                           k = kk + 1
                           do i = j + 1,n
                               ap(k) = ap(k) + x(i)*temp
                               k = k + 1
                           end do
                       else
                           ap(kk) = real(ap(kk),KIND=sp)
                       end if
                       kk = kk + n - j + 1
                   end do
               else
                   jx = kx
                   do j = 1,n
                       if (x(jx)/=czero) then
                           temp = alpha*conjg(x(jx))
                           ap(kk) = real(ap(kk),KIND=sp) + real(temp*x(jx),KIND=sp)
                           ix = jx
                           do k = kk + 1,kk + n - j
                               ix = ix + incx
                               ap(k) = ap(k) + x(ix)*temp
                           end do
                       else
                           ap(kk) = real(ap(kk),KIND=sp)
                       end if
                       jx = jx + incx
                       kk = kk + n - j + 1
                   end do
               end if
           end if
           return
     end subroutine stdlib_chpr