#:include "common.fypp" submodule(stdlib_lapack_base) stdlib_lapack_blas_like_base implicit none contains #:for ik,it,ii in LINALG_INT_KINDS_TYPES pure module subroutine stdlib${ii}$_slaset( uplo, m, n, alpha, beta, a, lda ) !! SLASET initializes an m-by-n matrix A to BETA on the diagonal and !! ALPHA on the offdiagonals. ! -- lapack auxiliary routine -- ! -- lapack is a software package provided by univ. of tennessee, -- ! -- univ. of california berkeley, univ. of colorado denver and nag ltd..-- ! Scalar Arguments character, intent(in) :: uplo integer(${ik}$), intent(in) :: lda, m, n real(sp), intent(in) :: alpha, beta ! Array Arguments real(sp), intent(out) :: a(lda,*) ! ===================================================================== ! Local Scalars integer(${ik}$) :: i, j ! Intrinsic Functions ! Executable Statements if( stdlib_lsame( uplo, 'U' ) ) then ! set the strictly upper triangular or trapezoidal part of the ! array to alpha. do j = 2, n do i = 1, min( j-1, m ) a( i, j ) = alpha end do end do else if( stdlib_lsame( uplo, 'L' ) ) then ! set the strictly lower triangular or trapezoidal part of the ! array to alpha. do j = 1, min( m, n ) do i = j + 1, m a( i, j ) = alpha end do end do else ! set the leading m-by-n submatrix to alpha. do j = 1, n do i = 1, m a( i, j ) = alpha end do end do end if ! set the first min(m,n) diagonal elements to beta. do i = 1, min( m, n ) a( i, i ) = beta end do return end subroutine stdlib${ii}$_slaset pure module subroutine stdlib${ii}$_dlaset( uplo, m, n, alpha, beta, a, lda ) !! DLASET initializes an m-by-n matrix A to BETA on the diagonal and !! ALPHA on the offdiagonals. ! -- lapack auxiliary routine -- ! -- lapack is a software package provided by univ. of tennessee, -- ! -- univ. of california berkeley, univ. of colorado denver and nag ltd..-- ! Scalar Arguments character, intent(in) :: uplo integer(${ik}$), intent(in) :: lda, m, n real(dp), intent(in) :: alpha, beta ! Array Arguments real(dp), intent(out) :: a(lda,*) ! ===================================================================== ! Local Scalars integer(${ik}$) :: i, j ! Intrinsic Functions ! Executable Statements if( stdlib_lsame( uplo, 'U' ) ) then ! set the strictly upper triangular or trapezoidal part of the ! array to alpha. do j = 2, n do i = 1, min( j-1, m ) a( i, j ) = alpha end do end do else if( stdlib_lsame( uplo, 'L' ) ) then ! set the strictly lower triangular or trapezoidal part of the ! array to alpha. do j = 1, min( m, n ) do i = j + 1, m a( i, j ) = alpha end do end do else ! set the leading m-by-n submatrix to alpha. do j = 1, n do i = 1, m a( i, j ) = alpha end do end do end if ! set the first min(m,n) diagonal elements to beta. do i = 1, min( m, n ) a( i, i ) = beta end do return end subroutine stdlib${ii}$_dlaset #:for rk,rt,ri in REAL_KINDS_TYPES #:if not rk in ["sp","dp"] pure module subroutine stdlib${ii}$_${ri}$laset( uplo, m, n, alpha, beta, a, lda ) !! DLASET: initializes an m-by-n matrix A to BETA on the diagonal and !! ALPHA on the offdiagonals. ! -- lapack auxiliary routine -- ! -- lapack is a software package provided by univ. of tennessee, -- ! -- univ. of california berkeley, univ. of colorado denver and nag ltd..-- ! Scalar Arguments character, intent(in) :: uplo integer(${ik}$), intent(in) :: lda, m, n real(${rk}$), intent(in) :: alpha, beta ! Array Arguments real(${rk}$), intent(out) :: a(lda,*) ! ===================================================================== ! Local Scalars integer(${ik}$) :: i, j ! Intrinsic Functions ! Executable Statements if( stdlib_lsame( uplo, 'U' ) ) then ! set the strictly upper triangular or trapezoidal part of the ! array to alpha. do j = 2, n do i = 1, min( j-1, m ) a( i, j ) = alpha end do end do else if( stdlib_lsame( uplo, 'L' ) ) then ! set the strictly lower triangular or trapezoidal part of the ! array to alpha. do j = 1, min( m, n ) do i = j + 1, m a( i, j ) = alpha end do end do else ! set the leading m-by-n submatrix to alpha. do j = 1, n do i = 1, m a( i, j ) = alpha end do end do end if ! set the first min(m,n) diagonal elements to beta. do i = 1, min( m, n ) a( i, i ) = beta end do return end subroutine stdlib${ii}$_${ri}$laset #:endif #:endfor pure module subroutine stdlib${ii}$_claset( uplo, m, n, alpha, beta, a, lda ) !! CLASET initializes a 2-D array A to BETA on the diagonal and !! ALPHA on the offdiagonals. ! -- lapack auxiliary routine -- ! -- lapack is a software package provided by univ. of tennessee, -- ! -- univ. of california berkeley, univ. of colorado denver and nag ltd..-- ! Scalar Arguments character, intent(in) :: uplo integer(${ik}$), intent(in) :: lda, m, n complex(sp), intent(in) :: alpha, beta ! Array Arguments complex(sp), intent(out) :: a(lda,*) ! ===================================================================== ! Local Scalars integer(${ik}$) :: i, j ! Intrinsic Functions ! Executable Statements if( stdlib_lsame( uplo, 'U' ) ) then ! set the diagonal to beta and the strictly upper triangular ! part of the array to alpha. do j = 2, n do i = 1, min( j-1, m ) a( i, j ) = alpha end do end do do i = 1, min( n, m ) a( i, i ) = beta end do else if( stdlib_lsame( uplo, 'L' ) ) then ! set the diagonal to beta and the strictly lower triangular ! part of the array to alpha. do j = 1, min( m, n ) do i = j + 1, m a( i, j ) = alpha end do end do do i = 1, min( n, m ) a( i, i ) = beta end do else ! set the array to beta on the diagonal and alpha on the ! offdiagonal. do j = 1, n do i = 1, m a( i, j ) = alpha end do end do do i = 1, min( m, n ) a( i, i ) = beta end do end if return end subroutine stdlib${ii}$_claset pure module subroutine stdlib${ii}$_zlaset( uplo, m, n, alpha, beta, a, lda ) !! ZLASET initializes a 2-D array A to BETA on the diagonal and !! ALPHA on the offdiagonals. ! -- lapack auxiliary routine -- ! -- lapack is a software package provided by univ. of tennessee, -- ! -- univ. of california berkeley, univ. of colorado denver and nag ltd..-- ! Scalar Arguments character, intent(in) :: uplo integer(${ik}$), intent(in) :: lda, m, n complex(dp), intent(in) :: alpha, beta ! Array Arguments complex(dp), intent(out) :: a(lda,*) ! ===================================================================== ! Local Scalars integer(${ik}$) :: i, j ! Intrinsic Functions ! Executable Statements if( stdlib_lsame( uplo, 'U' ) ) then ! set the diagonal to beta and the strictly upper triangular ! part of the array to alpha. do j = 2, n do i = 1, min( j-1, m ) a( i, j ) = alpha end do end do do i = 1, min( n, m ) a( i, i ) = beta end do else if( stdlib_lsame( uplo, 'L' ) ) then ! set the diagonal to beta and the strictly lower triangular ! part of the array to alpha. do j = 1, min( m, n ) do i = j + 1, m a( i, j ) = alpha end do end do do i = 1, min( n, m ) a( i, i ) = beta end do else ! set the array to beta on the diagonal and alpha on the ! offdiagonal. do j = 1, n do i = 1, m a( i, j ) = alpha end do end do do i = 1, min( m, n ) a( i, i ) = beta end do end if return end subroutine stdlib${ii}$_zlaset #:for ck,ct,ci in CMPLX_KINDS_TYPES #:if not ck in ["sp","dp"] pure module subroutine stdlib${ii}$_${ci}$laset( uplo, m, n, alpha, beta, a, lda ) !! ZLASET: initializes a 2-D array A to BETA on the diagonal and !! ALPHA on the offdiagonals. ! -- lapack auxiliary routine -- ! -- lapack is a software package provided by univ. of tennessee, -- ! -- univ. of california berkeley, univ. of colorado denver and nag ltd..-- ! Scalar Arguments character, intent(in) :: uplo integer(${ik}$), intent(in) :: lda, m, n complex(${ck}$), intent(in) :: alpha, beta ! Array Arguments complex(${ck}$), intent(out) :: a(lda,*) ! ===================================================================== ! Local Scalars integer(${ik}$) :: i, j ! Intrinsic Functions ! Executable Statements if( stdlib_lsame( uplo, 'U' ) ) then ! set the diagonal to beta and the strictly upper triangular ! part of the array to alpha. do j = 2, n do i = 1, min( j-1, m ) a( i, j ) = alpha end do end do do i = 1, min( n, m ) a( i, i ) = beta end do else if( stdlib_lsame( uplo, 'L' ) ) then ! set the diagonal to beta and the strictly lower triangular ! part of the array to alpha. do j = 1, min( m, n ) do i = j + 1, m a( i, j ) = alpha end do end do do i = 1, min( n, m ) a( i, i ) = beta end do else ! set the array to beta on the diagonal and alpha on the ! offdiagonal. do j = 1, n do i = 1, m a( i, j ) = alpha end do end do do i = 1, min( m, n ) a( i, i ) = beta end do end if return end subroutine stdlib${ii}$_${ci}$laset #:endif #:endfor pure module subroutine stdlib${ii}$_slarnv( idist, iseed, n, x ) !! SLARNV returns a vector of n random real numbers from a uniform or !! normal distribution. ! -- lapack auxiliary routine -- ! -- lapack is a software package provided by univ. of tennessee, -- ! -- univ. of california berkeley, univ. of colorado denver and nag ltd..-- use stdlib_blas_constants_sp, only: one, two ! Scalar Arguments integer(${ik}$), intent(in) :: idist, n ! Array Arguments integer(${ik}$), intent(inout) :: iseed(4_${ik}$) real(sp), intent(out) :: x(*) ! ===================================================================== ! Parameters integer(${ik}$), parameter :: lv = 128_${ik}$ real(sp), parameter :: twopi = 6.28318530717958647692528676655900576839e+0_sp ! Local Scalars integer(${ik}$) :: i, il, il2, iv ! Local Arrays real(sp) :: u(lv) ! Intrinsic Functions ! Executable Statements do 40 iv = 1, n, lv / 2 il = min( lv / 2_${ik}$, n-iv+1 ) if( idist==3_${ik}$ ) then il2 = 2_${ik}$*il else il2 = il end if ! call stdlib${ii}$_slaruv to generate il2 numbers from a uniform (0,1) ! distribution (il2 <= lv) call stdlib${ii}$_slaruv( iseed, il2, u ) if( idist==1_${ik}$ ) then ! copy generated numbers do i = 1, il x( iv+i-1 ) = u( i ) end do else if( idist==2_${ik}$ ) then ! convert generated numbers to uniform (-1,1) distribution do i = 1, il x( iv+i-1 ) = two*u( i ) - one end do else if( idist==3_${ik}$ ) then ! convert generated numbers to normal (0,1) distribution do i = 1, il x( iv+i-1 ) = sqrt( -two*log( u( 2_${ik}$*i-1 ) ) )*cos( twopi*u( 2_${ik}$*i ) ) end do end if 40 continue return end subroutine stdlib${ii}$_slarnv pure module subroutine stdlib${ii}$_dlarnv( idist, iseed, n, x ) !! DLARNV returns a vector of n random real numbers from a uniform or !! normal distribution. ! -- lapack auxiliary routine -- ! -- lapack is a software package provided by univ. of tennessee, -- ! -- univ. of california berkeley, univ. of colorado denver and nag ltd..-- use stdlib_blas_constants_dp, only: one, two ! Scalar Arguments integer(${ik}$), intent(in) :: idist, n ! Array Arguments integer(${ik}$), intent(inout) :: iseed(4_${ik}$) real(dp), intent(out) :: x(*) ! ===================================================================== ! Parameters integer(${ik}$), parameter :: lv = 128_${ik}$ real(dp), parameter :: twopi = 6.28318530717958647692528676655900576839e+0_dp ! Local Scalars integer(${ik}$) :: i, il, il2, iv ! Local Arrays real(dp) :: u(lv) ! Intrinsic Functions ! Executable Statements do 40 iv = 1, n, lv / 2 il = min( lv / 2_${ik}$, n-iv+1 ) if( idist==3_${ik}$ ) then il2 = 2_${ik}$*il else il2 = il end if ! call stdlib${ii}$_dlaruv to generate il2 numbers from a uniform (0,1) ! distribution (il2 <= lv) call stdlib${ii}$_dlaruv( iseed, il2, u ) if( idist==1_${ik}$ ) then ! copy generated numbers do i = 1, il x( iv+i-1 ) = u( i ) end do else if( idist==2_${ik}$ ) then ! convert generated numbers to uniform (-1,1) distribution do i = 1, il x( iv+i-1 ) = two*u( i ) - one end do else if( idist==3_${ik}$ ) then ! convert generated numbers to normal (0,1) distribution do i = 1, il x( iv+i-1 ) = sqrt( -two*log( u( 2_${ik}$*i-1 ) ) )*cos( twopi*u( 2_${ik}$*i ) ) end do end if 40 continue return end subroutine stdlib${ii}$_dlarnv #:for rk,rt,ri in REAL_KINDS_TYPES #:if not rk in ["sp","dp"] pure module subroutine stdlib${ii}$_${ri}$larnv( idist, iseed, n, x ) !! DLARNV: returns a vector of n random real numbers from a uniform or !! normal distribution. ! -- lapack auxiliary routine -- ! -- lapack is a software package provided by univ. of tennessee, -- ! -- univ. of california berkeley, univ. of colorado denver and nag ltd..-- use stdlib_blas_constants_${rk}$, only: one, two ! Scalar Arguments integer(${ik}$), intent(in) :: idist, n ! Array Arguments integer(${ik}$), intent(inout) :: iseed(4_${ik}$) real(${rk}$), intent(out) :: x(*) ! ===================================================================== ! Parameters integer(${ik}$), parameter :: lv = 128_${ik}$ real(${rk}$), parameter :: twopi = 6.28318530717958647692528676655900576839e+0_${rk}$ ! Local Scalars integer(${ik}$) :: i, il, il2, iv ! Local Arrays real(${rk}$) :: u(lv) ! Intrinsic Functions ! Executable Statements do 40 iv = 1, n, lv / 2 il = min( lv / 2_${ik}$, n-iv+1 ) if( idist==3_${ik}$ ) then il2 = 2_${ik}$*il else il2 = il end if ! call stdlib${ii}$_${ri}$laruv to generate il2 numbers from a uniform (0,1) ! distribution (il2 <= lv) call stdlib${ii}$_${ri}$laruv( iseed, il2, u ) if( idist==1_${ik}$ ) then ! copy generated numbers do i = 1, il x( iv+i-1 ) = u( i ) end do else if( idist==2_${ik}$ ) then ! convert generated numbers to uniform (-1,1) distribution do i = 1, il x( iv+i-1 ) = two*u( i ) - one end do else if( idist==3_${ik}$ ) then ! convert generated numbers to normal (0,1) distribution do i = 1, il x( iv+i-1 ) = sqrt( -two*log( u( 2_${ik}$*i-1 ) ) )*cos( twopi*u( 2_${ik}$*i ) ) end do end if 40 continue return end subroutine stdlib${ii}$_${ri}$larnv #:endif #:endfor pure module subroutine stdlib${ii}$_clarnv( idist, iseed, n, x ) !! CLARNV returns a vector of n random complex numbers from a uniform or !! normal distribution. ! -- lapack auxiliary routine -- ! -- lapack is a software package provided by univ. of tennessee, -- ! -- univ. of california berkeley, univ. of colorado denver and nag ltd..-- use stdlib_blas_constants_sp, only: zero, one, two ! Scalar Arguments integer(${ik}$), intent(in) :: idist, n ! Array Arguments integer(${ik}$), intent(inout) :: iseed(4_${ik}$) complex(sp), intent(out) :: x(*) ! ===================================================================== ! Parameters integer(${ik}$), parameter :: lv = 128_${ik}$ real(sp), parameter :: twopi = 6.28318530717958647692528676655900576839e+0_sp ! Local Scalars integer(${ik}$) :: i, il, iv ! Local Arrays real(sp) :: u(lv) ! Intrinsic Functions ! Executable Statements do 60 iv = 1, n, lv / 2 il = min( lv / 2_${ik}$, n-iv+1 ) ! call stdlib${ii}$_slaruv to generate 2*il realnumbers from a uniform (0,1,KIND=sp) ! distribution (2*il <= lv) call stdlib${ii}$_slaruv( iseed, 2_${ik}$*il, u ) if( idist==1_${ik}$ ) then ! copy generated numbers do i = 1, il x( iv+i-1 ) = cmplx( u( 2_${ik}$*i-1 ), u( 2_${ik}$*i ),KIND=sp) end do else if( idist==2_${ik}$ ) then ! convert generated numbers to uniform (-1,1) distribution do i = 1, il x( iv+i-1 ) = cmplx( two*u( 2_${ik}$*i-1 )-one,two*u( 2_${ik}$*i )-one,KIND=sp) end do else if( idist==3_${ik}$ ) then ! convert generated numbers to normal (0,1) distribution do i = 1, il x( iv+i-1 ) = sqrt( -two*log( u( 2_${ik}$*i-1 ) ) )*exp( cmplx( zero, twopi*u( 2_${ik}$*i ),& KIND=sp) ) end do else if( idist==4_${ik}$ ) then ! convert generated numbers to complex numbers uniformly ! distributed on the unit disk do i = 1, il x( iv+i-1 ) = sqrt( u( 2_${ik}$*i-1 ) )*exp( cmplx( zero, twopi*u( 2_${ik}$*i ),KIND=sp) ) end do else if( idist==5_${ik}$ ) then ! convert generated numbers to complex numbers uniformly ! distributed on the unit circle do i = 1, il x( iv+i-1 ) = exp( cmplx( zero, twopi*u( 2_${ik}$*i ),KIND=sp) ) end do end if 60 continue return end subroutine stdlib${ii}$_clarnv pure module subroutine stdlib${ii}$_zlarnv( idist, iseed, n, x ) !! ZLARNV returns a vector of n random complex numbers from a uniform or !! normal distribution. ! -- lapack auxiliary routine -- ! -- lapack is a software package provided by univ. of tennessee, -- ! -- univ. of california berkeley, univ. of colorado denver and nag ltd..-- use stdlib_blas_constants_dp, only: zero, one, two ! Scalar Arguments integer(${ik}$), intent(in) :: idist, n ! Array Arguments integer(${ik}$), intent(inout) :: iseed(4_${ik}$) complex(dp), intent(out) :: x(*) ! ===================================================================== ! Parameters integer(${ik}$), parameter :: lv = 128_${ik}$ real(dp), parameter :: twopi = 6.28318530717958647692528676655900576839e+0_dp ! Local Scalars integer(${ik}$) :: i, il, iv ! Local Arrays real(dp) :: u(lv) ! Intrinsic Functions ! Executable Statements do 60 iv = 1, n, lv / 2 il = min( lv / 2_${ik}$, n-iv+1 ) ! call stdlib${ii}$_dlaruv to generate 2*il realnumbers from a uniform (0,1,KIND=dp) ! distribution (2*il <= lv) call stdlib${ii}$_dlaruv( iseed, 2_${ik}$*il, u ) if( idist==1_${ik}$ ) then ! copy generated numbers do i = 1, il x( iv+i-1 ) = cmplx( u( 2_${ik}$*i-1 ), u( 2_${ik}$*i ),KIND=dp) end do else if( idist==2_${ik}$ ) then ! convert generated numbers to uniform (-1,1) distribution do i = 1, il x( iv+i-1 ) = cmplx( two*u( 2_${ik}$*i-1 )-one,two*u( 2_${ik}$*i )-one,KIND=dp) end do else if( idist==3_${ik}$ ) then ! convert generated numbers to normal (0,1) distribution do i = 1, il x( iv+i-1 ) = sqrt( -two*log( u( 2_${ik}$*i-1 ) ) )*exp( cmplx( zero, twopi*u( 2_${ik}$*i ),& KIND=dp) ) end do else if( idist==4_${ik}$ ) then ! convert generated numbers to complex numbers uniformly ! distributed on the unit disk do i = 1, il x( iv+i-1 ) = sqrt( u( 2_${ik}$*i-1 ) )*exp( cmplx( zero, twopi*u( 2_${ik}$*i ),KIND=dp) ) end do else if( idist==5_${ik}$ ) then ! convert generated numbers to complex numbers uniformly ! distributed on the unit circle do i = 1, il x( iv+i-1 ) = exp( cmplx( zero, twopi*u( 2_${ik}$*i ),KIND=dp) ) end do end if 60 continue return end subroutine stdlib${ii}$_zlarnv #:for ck,ct,ci in CMPLX_KINDS_TYPES #:if not ck in ["sp","dp"] pure module subroutine stdlib${ii}$_${ci}$larnv( idist, iseed, n, x ) !! ZLARNV: returns a vector of n random complex numbers from a uniform or !! normal distribution. ! -- lapack auxiliary routine -- ! -- lapack is a software package provided by univ. of tennessee, -- ! -- univ. of california berkeley, univ. of colorado denver and nag ltd..-- use stdlib_blas_constants_${ck}$, only: zero, one, two ! Scalar Arguments integer(${ik}$), intent(in) :: idist, n ! Array Arguments integer(${ik}$), intent(inout) :: iseed(4_${ik}$) complex(${ck}$), intent(out) :: x(*) ! ===================================================================== ! Parameters integer(${ik}$), parameter :: lv = 128_${ik}$ real(${ck}$), parameter :: twopi = 6.28318530717958647692528676655900576839e+0_${ck}$ ! Local Scalars integer(${ik}$) :: i, il, iv ! Local Arrays real(${ck}$) :: u(lv) ! Intrinsic Functions ! Executable Statements do 60 iv = 1, n, lv / 2 il = min( lv / 2_${ik}$, n-iv+1 ) ! call stdlib${ii}$_${c2ri(ci)}$laruv to generate 2*il realnumbers from a uniform (0,1,KIND=${ck}$) ! distribution (2*il <= lv) call stdlib${ii}$_${c2ri(ci)}$laruv( iseed, 2_${ik}$*il, u ) if( idist==1_${ik}$ ) then ! copy generated numbers do i = 1, il x( iv+i-1 ) = cmplx( u( 2_${ik}$*i-1 ), u( 2_${ik}$*i ),KIND=${ck}$) end do else if( idist==2_${ik}$ ) then ! convert generated numbers to uniform (-1,1) distribution do i = 1, il x( iv+i-1 ) = cmplx( two*u( 2_${ik}$*i-1 )-one,two*u( 2_${ik}$*i )-one,KIND=${ck}$) end do else if( idist==3_${ik}$ ) then ! convert generated numbers to normal (0,1) distribution do i = 1, il x( iv+i-1 ) = sqrt( -two*log( u( 2_${ik}$*i-1 ) ) )*exp( cmplx( zero, twopi*u( 2_${ik}$*i ),& KIND=${ck}$) ) end do else if( idist==4_${ik}$ ) then ! convert generated numbers to complex numbers uniformly ! distributed on the unit disk do i = 1, il x( iv+i-1 ) = sqrt( u( 2_${ik}$*i-1 ) )*exp( cmplx( zero, twopi*u( 2_${ik}$*i ),KIND=${ck}$) ) end do else if( idist==5_${ik}$ ) then ! convert generated numbers to complex numbers uniformly ! distributed on the unit circle do i = 1, il x( iv+i-1 ) = exp( cmplx( zero, twopi*u( 2_${ik}$*i ),KIND=${ck}$) ) end do end if 60 continue return end subroutine stdlib${ii}$_${ci}$larnv #:endif #:endfor pure module subroutine stdlib${ii}$_slaruv( iseed, n, x ) !! SLARUV returns a vector of n random real numbers from a uniform (0,1) !! distribution (n <= 128). !! This is an auxiliary routine called by SLARNV and CLARNV. ! -- lapack auxiliary routine -- ! -- lapack is a software package provided by univ. of tennessee, -- ! -- univ. of california berkeley, univ. of colorado denver and nag ltd..-- use stdlib_blas_constants_sp, only: one ! Scalar Arguments integer(${ik}$), intent(in) :: n ! Array Arguments integer(${ik}$), intent(inout) :: iseed(4_${ik}$) real(sp), intent(out) :: x(n) ! ===================================================================== ! Parameters integer(${ik}$), parameter :: lv = 128_${ik}$ integer(${ik}$), parameter :: ipw2 = 4096_${ik}$ real(sp), parameter :: r = one/ipw2 ! Local Scalars integer(${ik}$) :: i, i1, i2, i3, i4, it1, it2, it3, it4 ! Local Arrays integer(${ik}$) :: mm(lv,4_${ik}$) ! Intrinsic Functions ! Data Statements mm(1_${ik}$,1_${ik}$:4_${ik}$)=[494_${ik}$,322_${ik}$,2508_${ik}$,2549_${ik}$] mm(2_${ik}$,1_${ik}$:4_${ik}$)=[2637_${ik}$,789_${ik}$,3754_${ik}$,1145_${ik}$] mm(3_${ik}$,1_${ik}$:4_${ik}$)=[255_${ik}$,1440_${ik}$,1766_${ik}$,2253_${ik}$] mm(4_${ik}$,1_${ik}$:4_${ik}$)=[2008_${ik}$,752_${ik}$,3572_${ik}$,305_${ik}$] mm(5_${ik}$,1_${ik}$:4_${ik}$)=[1253_${ik}$,2859_${ik}$,2893_${ik}$,3301_${ik}$] mm(6_${ik}$,1_${ik}$:4_${ik}$)=[3344_${ik}$,123_${ik}$,307_${ik}$,1065_${ik}$] mm(7_${ik}$,1_${ik}$:4_${ik}$)=[4084_${ik}$,1848_${ik}$,1297_${ik}$,3133_${ik}$] mm(8_${ik}$,1_${ik}$:4_${ik}$)=[1739_${ik}$,643_${ik}$,3966_${ik}$,2913_${ik}$] mm(9_${ik}$,1_${ik}$:4_${ik}$)=[3143_${ik}$,2405_${ik}$,758_${ik}$,3285_${ik}$] mm(10_${ik}$,1_${ik}$:4_${ik}$)=[3468_${ik}$,2638_${ik}$,2598_${ik}$,1241_${ik}$] mm(11_${ik}$,1_${ik}$:4_${ik}$)=[688_${ik}$,2344_${ik}$,3406_${ik}$,1197_${ik}$] mm(12_${ik}$,1_${ik}$:4_${ik}$)=[1657_${ik}$,46_${ik}$,2922_${ik}$,3729_${ik}$] mm(13_${ik}$,1_${ik}$:4_${ik}$)=[1238_${ik}$,3814_${ik}$,1038_${ik}$,2501_${ik}$] mm(14_${ik}$,1_${ik}$:4_${ik}$)=[3166_${ik}$,913_${ik}$,2934_${ik}$,1673_${ik}$] mm(15_${ik}$,1_${ik}$:4_${ik}$)=[1292_${ik}$,3649_${ik}$,2091_${ik}$,541_${ik}$] mm(16_${ik}$,1_${ik}$:4_${ik}$)=[3422_${ik}$,339_${ik}$,2451_${ik}$,2753_${ik}$] mm(17_${ik}$,1_${ik}$:4_${ik}$)=[1270_${ik}$,3808_${ik}$,1580_${ik}$,949_${ik}$] mm(18_${ik}$,1_${ik}$:4_${ik}$)=[2016_${ik}$,822_${ik}$,1958_${ik}$,2361_${ik}$] mm(19_${ik}$,1_${ik}$:4_${ik}$)=[154_${ik}$,2832_${ik}$,2055_${ik}$,1165_${ik}$] mm(20_${ik}$,1_${ik}$:4_${ik}$)=[2862_${ik}$,3078_${ik}$,1507_${ik}$,4081_${ik}$] mm(21_${ik}$,1_${ik}$:4_${ik}$)=[697_${ik}$,3633_${ik}$,1078_${ik}$,2725_${ik}$] mm(22_${ik}$,1_${ik}$:4_${ik}$)=[1706_${ik}$,2970_${ik}$,3273_${ik}$,3305_${ik}$] mm(23_${ik}$,1_${ik}$:4_${ik}$)=[491_${ik}$,637_${ik}$,17_${ik}$,3069_${ik}$] mm(24_${ik}$,1_${ik}$:4_${ik}$)=[931_${ik}$,2249_${ik}$,854_${ik}$,3617_${ik}$] mm(25_${ik}$,1_${ik}$:4_${ik}$)=[1444_${ik}$,2081_${ik}$,2916_${ik}$,3733_${ik}$] mm(26_${ik}$,1_${ik}$:4_${ik}$)=[444_${ik}$,4019_${ik}$,3971_${ik}$,409_${ik}$] mm(27_${ik}$,1_${ik}$:4_${ik}$)=[3577_${ik}$,1478_${ik}$,2889_${ik}$,2157_${ik}$] mm(28_${ik}$,1_${ik}$:4_${ik}$)=[3944_${ik}$,242_${ik}$,3831_${ik}$,1361_${ik}$] mm(29_${ik}$,1_${ik}$:4_${ik}$)=[2184_${ik}$,481_${ik}$,2621_${ik}$,3973_${ik}$] mm(30_${ik}$,1_${ik}$:4_${ik}$)=[1661_${ik}$,2075_${ik}$,1541_${ik}$,1865_${ik}$] mm(31_${ik}$,1_${ik}$:4_${ik}$)=[3482_${ik}$,4058_${ik}$,893_${ik}$,2525_${ik}$] mm(32_${ik}$,1_${ik}$:4_${ik}$)=[657_${ik}$,622_${ik}$,736_${ik}$,1409_${ik}$] mm(33_${ik}$,1_${ik}$:4_${ik}$)=[3023_${ik}$,3376_${ik}$,3992_${ik}$,3445_${ik}$] mm(34_${ik}$,1_${ik}$:4_${ik}$)=[3618_${ik}$,812_${ik}$,787_${ik}$,3577_${ik}$] mm(35_${ik}$,1_${ik}$:4_${ik}$)=[1267_${ik}$,234_${ik}$,2125_${ik}$,77_${ik}$] mm(36_${ik}$,1_${ik}$:4_${ik}$)=[1828_${ik}$,641_${ik}$,2364_${ik}$,3761_${ik}$] mm(37_${ik}$,1_${ik}$:4_${ik}$)=[164_${ik}$,4005_${ik}$,2460_${ik}$,2149_${ik}$] mm(38_${ik}$,1_${ik}$:4_${ik}$)=[3798_${ik}$,1122_${ik}$,257_${ik}$,1449_${ik}$] mm(39_${ik}$,1_${ik}$:4_${ik}$)=[3087_${ik}$,3135_${ik}$,1574_${ik}$,3005_${ik}$] mm(40_${ik}$,1_${ik}$:4_${ik}$)=[2400_${ik}$,2640_${ik}$,3912_${ik}$,225_${ik}$] mm(41_${ik}$,1_${ik}$:4_${ik}$)=[2870_${ik}$,2302_${ik}$,1216_${ik}$,85_${ik}$] mm(42_${ik}$,1_${ik}$:4_${ik}$)=[3876_${ik}$,40_${ik}$,3248_${ik}$,3673_${ik}$] mm(43_${ik}$,1_${ik}$:4_${ik}$)=[1905_${ik}$,1832_${ik}$,3401_${ik}$,3117_${ik}$] mm(44_${ik}$,1_${ik}$:4_${ik}$)=[1593_${ik}$,2247_${ik}$,2124_${ik}$,3089_${ik}$] mm(45_${ik}$,1_${ik}$:4_${ik}$)=[1797_${ik}$,2034_${ik}$,2762_${ik}$,1349_${ik}$] mm(46_${ik}$,1_${ik}$:4_${ik}$)=[1234_${ik}$,2637_${ik}$,149_${ik}$,2057_${ik}$] mm(47_${ik}$,1_${ik}$:4_${ik}$)=[3460_${ik}$,1287_${ik}$,2245_${ik}$,413_${ik}$] mm(48_${ik}$,1_${ik}$:4_${ik}$)=[328_${ik}$,1691_${ik}$,166_${ik}$,65_${ik}$] mm(49_${ik}$,1_${ik}$:4_${ik}$)=[2861_${ik}$,496_${ik}$,466_${ik}$,1845_${ik}$] mm(50_${ik}$,1_${ik}$:4_${ik}$)=[1950_${ik}$,1597_${ik}$,4018_${ik}$,697_${ik}$] mm(51_${ik}$,1_${ik}$:4_${ik}$)=[617_${ik}$,2394_${ik}$,1399_${ik}$,3085_${ik}$] mm(52_${ik}$,1_${ik}$:4_${ik}$)=[2070_${ik}$,2584_${ik}$,190_${ik}$,3441_${ik}$] mm(53_${ik}$,1_${ik}$:4_${ik}$)=[3331_${ik}$,1843_${ik}$,2879_${ik}$,1573_${ik}$] mm(54_${ik}$,1_${ik}$:4_${ik}$)=[769_${ik}$,336_${ik}$,153_${ik}$,3689_${ik}$] mm(55_${ik}$,1_${ik}$:4_${ik}$)=[1558_${ik}$,1472_${ik}$,2320_${ik}$,2941_${ik}$] mm(56_${ik}$,1_${ik}$:4_${ik}$)=[2412_${ik}$,2407_${ik}$,18_${ik}$,929_${ik}$] mm(57_${ik}$,1_${ik}$:4_${ik}$)=[2800_${ik}$,433_${ik}$,712_${ik}$,533_${ik}$] mm(58_${ik}$,1_${ik}$:4_${ik}$)=[189_${ik}$,2096_${ik}$,2159_${ik}$,2841_${ik}$] mm(59_${ik}$,1_${ik}$:4_${ik}$)=[287_${ik}$,1761_${ik}$,2318_${ik}$,4077_${ik}$] mm(60_${ik}$,1_${ik}$:4_${ik}$)=[2045_${ik}$,2810_${ik}$,2091_${ik}$,721_${ik}$] mm(61_${ik}$,1_${ik}$:4_${ik}$)=[1227_${ik}$,566_${ik}$,3443_${ik}$,2821_${ik}$] mm(62_${ik}$,1_${ik}$:4_${ik}$)=[2838_${ik}$,442_${ik}$,1510_${ik}$,2249_${ik}$] mm(63_${ik}$,1_${ik}$:4_${ik}$)=[209_${ik}$,41_${ik}$,449_${ik}$,2397_${ik}$] mm(64_${ik}$,1_${ik}$:4_${ik}$)=[2770_${ik}$,1238_${ik}$,1956_${ik}$,2817_${ik}$] mm(65_${ik}$,1_${ik}$:4_${ik}$)=[3654_${ik}$,1086_${ik}$,2201_${ik}$,245_${ik}$] mm(66_${ik}$,1_${ik}$:4_${ik}$)=[3993_${ik}$,603_${ik}$,3137_${ik}$,1913_${ik}$] mm(67_${ik}$,1_${ik}$:4_${ik}$)=[192_${ik}$,840_${ik}$,3399_${ik}$,1997_${ik}$] mm(68_${ik}$,1_${ik}$:4_${ik}$)=[2253_${ik}$,3168_${ik}$,1321_${ik}$,3121_${ik}$] mm(69_${ik}$,1_${ik}$:4_${ik}$)=[3491_${ik}$,1499_${ik}$,2271_${ik}$,997_${ik}$] mm(70_${ik}$,1_${ik}$:4_${ik}$)=[2889_${ik}$,1084_${ik}$,3667_${ik}$,1833_${ik}$] mm(71_${ik}$,1_${ik}$:4_${ik}$)=[2857_${ik}$,3438_${ik}$,2703_${ik}$,2877_${ik}$] mm(72_${ik}$,1_${ik}$:4_${ik}$)=[2094_${ik}$,2408_${ik}$,629_${ik}$,1633_${ik}$] mm(73_${ik}$,1_${ik}$:4_${ik}$)=[1818_${ik}$,1589_${ik}$,2365_${ik}$,981_${ik}$] mm(74_${ik}$,1_${ik}$:4_${ik}$)=[688_${ik}$,2391_${ik}$,2431_${ik}$,2009_${ik}$] mm(75_${ik}$,1_${ik}$:4_${ik}$)=[1407_${ik}$,288_${ik}$,1113_${ik}$,941_${ik}$] mm(76_${ik}$,1_${ik}$:4_${ik}$)=[634_${ik}$,26_${ik}$,3922_${ik}$,2449_${ik}$] mm(77_${ik}$,1_${ik}$:4_${ik}$)=[3231_${ik}$,512_${ik}$,2554_${ik}$,197_${ik}$] mm(78_${ik}$,1_${ik}$:4_${ik}$)=[815_${ik}$,1456_${ik}$,184_${ik}$,2441_${ik}$] mm(79_${ik}$,1_${ik}$:4_${ik}$)=[3524_${ik}$,171_${ik}$,2099_${ik}$,285_${ik}$] mm(80_${ik}$,1_${ik}$:4_${ik}$)=[1914_${ik}$,1677_${ik}$,3228_${ik}$,1473_${ik}$] mm(81_${ik}$,1_${ik}$:4_${ik}$)=[516_${ik}$,2657_${ik}$,4012_${ik}$,2741_${ik}$] mm(82_${ik}$,1_${ik}$:4_${ik}$)=[164_${ik}$,2270_${ik}$,1921_${ik}$,3129_${ik}$] mm(83_${ik}$,1_${ik}$:4_${ik}$)=[303_${ik}$,2587_${ik}$,3452_${ik}$,909_${ik}$] mm(84_${ik}$,1_${ik}$:4_${ik}$)=[2144_${ik}$,2961_${ik}$,3901_${ik}$,2801_${ik}$] mm(85_${ik}$,1_${ik}$:4_${ik}$)=[3480_${ik}$,1970_${ik}$,572_${ik}$,421_${ik}$] mm(86_${ik}$,1_${ik}$:4_${ik}$)=[119_${ik}$,1817_${ik}$,3309_${ik}$,4073_${ik}$] mm(87_${ik}$,1_${ik}$:4_${ik}$)=[3357_${ik}$,676_${ik}$,3171_${ik}$,2813_${ik}$] mm(88_${ik}$,1_${ik}$:4_${ik}$)=[837_${ik}$,1410_${ik}$,817_${ik}$,2337_${ik}$] mm(89_${ik}$,1_${ik}$:4_${ik}$)=[2826_${ik}$,3723_${ik}$,3039_${ik}$,1429_${ik}$] mm(90_${ik}$,1_${ik}$:4_${ik}$)=[2332_${ik}$,2803_${ik}$,1696_${ik}$,1177_${ik}$] mm(91_${ik}$,1_${ik}$:4_${ik}$)=[2089_${ik}$,3185_${ik}$,1256_${ik}$,1901_${ik}$] mm(92_${ik}$,1_${ik}$:4_${ik}$)=[3780_${ik}$,184_${ik}$,3715_${ik}$,81_${ik}$] mm(93_${ik}$,1_${ik}$:4_${ik}$)=[1700_${ik}$,663_${ik}$,2077_${ik}$,1669_${ik}$] mm(94_${ik}$,1_${ik}$:4_${ik}$)=[3712_${ik}$,499_${ik}$,3019_${ik}$,2633_${ik}$] mm(95_${ik}$,1_${ik}$:4_${ik}$)=[150_${ik}$,3784_${ik}$,1497_${ik}$,2269_${ik}$] mm(96_${ik}$,1_${ik}$:4_${ik}$)=[2000_${ik}$,1631_${ik}$,1101_${ik}$,129_${ik}$] mm(97_${ik}$,1_${ik}$:4_${ik}$)=[3375_${ik}$,1925_${ik}$,717_${ik}$,1141_${ik}$] mm(98_${ik}$,1_${ik}$:4_${ik}$)=[1621_${ik}$,3912_${ik}$,51_${ik}$,249_${ik}$] mm(99_${ik}$,1_${ik}$:4_${ik}$)=[3090_${ik}$,1398_${ik}$,981_${ik}$,3917_${ik}$] mm(100_${ik}$,1_${ik}$:4_${ik}$)=[3765_${ik}$,1349_${ik}$,1978_${ik}$,2481_${ik}$] mm(101_${ik}$,1_${ik}$:4_${ik}$)=[1149_${ik}$,1441_${ik}$,1813_${ik}$,3941_${ik}$] mm(102_${ik}$,1_${ik}$:4_${ik}$)=[3146_${ik}$,2224_${ik}$,3881_${ik}$,2217_${ik}$] mm(103_${ik}$,1_${ik}$:4_${ik}$)=[33_${ik}$,2411_${ik}$,76_${ik}$,2749_${ik}$] mm(104_${ik}$,1_${ik}$:4_${ik}$)=[3082_${ik}$,1907_${ik}$,3846_${ik}$,3041_${ik}$] mm(105_${ik}$,1_${ik}$:4_${ik}$)=[2741_${ik}$,3192_${ik}$,3694_${ik}$,1877_${ik}$] mm(106_${ik}$,1_${ik}$:4_${ik}$)=[359_${ik}$,2786_${ik}$,1682_${ik}$,345_${ik}$] mm(107_${ik}$,1_${ik}$:4_${ik}$)=[3316_${ik}$,382_${ik}$,124_${ik}$,2861_${ik}$] mm(108_${ik}$,1_${ik}$:4_${ik}$)=[1749_${ik}$,37_${ik}$,1660_${ik}$,1809_${ik}$] mm(109_${ik}$,1_${ik}$:4_${ik}$)=[185_${ik}$,759_${ik}$,3997_${ik}$,3141_${ik}$] mm(110_${ik}$,1_${ik}$:4_${ik}$)=[2784_${ik}$,2948_${ik}$,479_${ik}$,2825_${ik}$] mm(111_${ik}$,1_${ik}$:4_${ik}$)=[2202_${ik}$,1862_${ik}$,1141_${ik}$,157_${ik}$] mm(112_${ik}$,1_${ik}$:4_${ik}$)=[2199_${ik}$,3802_${ik}$,886_${ik}$,2881_${ik}$] mm(113_${ik}$,1_${ik}$:4_${ik}$)=[1364_${ik}$,2423_${ik}$,3514_${ik}$,3637_${ik}$] mm(114_${ik}$,1_${ik}$:4_${ik}$)=[1244_${ik}$,2051_${ik}$,1301_${ik}$,1465_${ik}$] mm(115_${ik}$,1_${ik}$:4_${ik}$)=[2020_${ik}$,2295_${ik}$,3604_${ik}$,2829_${ik}$] mm(116_${ik}$,1_${ik}$:4_${ik}$)=[3160_${ik}$,1332_${ik}$,1888_${ik}$,2161_${ik}$] mm(117_${ik}$,1_${ik}$:4_${ik}$)=[2785_${ik}$,1832_${ik}$,1836_${ik}$,3365_${ik}$] mm(118_${ik}$,1_${ik}$:4_${ik}$)=[2772_${ik}$,2405_${ik}$,1990_${ik}$,361_${ik}$] mm(119_${ik}$,1_${ik}$:4_${ik}$)=[1217_${ik}$,3638_${ik}$,2058_${ik}$,2685_${ik}$] mm(120_${ik}$,1_${ik}$:4_${ik}$)=[1822_${ik}$,3661_${ik}$,692_${ik}$,3745_${ik}$] mm(121_${ik}$,1_${ik}$:4_${ik}$)=[1245_${ik}$,327_${ik}$,1194_${ik}$,2325_${ik}$] mm(122_${ik}$,1_${ik}$:4_${ik}$)=[2252_${ik}$,3660_${ik}$,20_${ik}$,3609_${ik}$] mm(123_${ik}$,1_${ik}$:4_${ik}$)=[3904_${ik}$,716_${ik}$,3285_${ik}$,3821_${ik}$] mm(124_${ik}$,1_${ik}$:4_${ik}$)=[2774_${ik}$,1842_${ik}$,2046_${ik}$,3537_${ik}$] mm(125_${ik}$,1_${ik}$:4_${ik}$)=[997_${ik}$,3987_${ik}$,2107_${ik}$,517_${ik}$] mm(126_${ik}$,1_${ik}$:4_${ik}$)=[2573_${ik}$,1368_${ik}$,3508_${ik}$,3017_${ik}$] mm(127_${ik}$,1_${ik}$:4_${ik}$)=[1148_${ik}$,1848_${ik}$,3525_${ik}$,2141_${ik}$] mm(128_${ik}$,1_${ik}$:4_${ik}$)=[545_${ik}$,2366_${ik}$,3801_${ik}$,1537_${ik}$] ! Executable Statements i1 = iseed( 1_${ik}$ ) i2 = iseed( 2_${ik}$ ) i3 = iseed( 3_${ik}$ ) i4 = iseed( 4_${ik}$ ) loop_10: do i = 1, min( n, lv ) 20 continue ! multiply the seed by i-th power of the multiplier modulo 2**48 it4 = i4*mm( i, 4_${ik}$ ) it3 = it4 / ipw2 it4 = it4 - ipw2*it3 it3 = it3 + i3*mm( i, 4_${ik}$ ) + i4*mm( i, 3_${ik}$ ) it2 = it3 / ipw2 it3 = it3 - ipw2*it2 it2 = it2 + i2*mm( i, 4_${ik}$ ) + i3*mm( i, 3_${ik}$ ) + i4*mm( i, 2_${ik}$ ) it1 = it2 / ipw2 it2 = it2 - ipw2*it1 it1 = it1 + i1*mm( i, 4_${ik}$ ) + i2*mm( i, 3_${ik}$ ) + i3*mm( i, 2_${ik}$ ) +i4*mm( i, 1_${ik}$ ) it1 = mod( it1, ipw2 ) ! convert 48-bit integer to a realnumber in the interval (0,1,KIND=sp) x( i ) = r*( real( it1,KIND=sp)+r*( real( it2,KIND=sp)+r*( real( it3,KIND=sp)+& r*real( it4,KIND=sp) ) ) ) if (x( i )==one) then ! if a real number has n bits of precision, and the first ! n bits of the 48-bit integer above happen to be all 1 (which ! will occur about once every 2**n calls), then x( i ) will ! be rounded to exactly one. in ieee single precision arithmetic, ! this will happen relatively often since n = 24. ! since x( i ) is not supposed to return exactly 0.0_sp or one, ! the statistically correct thing to do in this situation is ! simply to iterate again. ! n.b. the case x( i ) = 0.0_sp should not be possible. i1 = i1 + 2_${ik}$ i2 = i2 + 2_${ik}$ i3 = i3 + 2_${ik}$ i4 = i4 + 2_${ik}$ goto 20 end if end do loop_10 ! return final value of seed iseed( 1_${ik}$ ) = it1 iseed( 2_${ik}$ ) = it2 iseed( 3_${ik}$ ) = it3 iseed( 4_${ik}$ ) = it4 return end subroutine stdlib${ii}$_slaruv pure module subroutine stdlib${ii}$_dlaruv( iseed, n, x ) !! DLARUV returns a vector of n random real numbers from a uniform (0,1) !! distribution (n <= 128). !! This is an auxiliary routine called by DLARNV and ZLARNV. ! -- lapack auxiliary routine -- ! -- lapack is a software package provided by univ. of tennessee, -- ! -- univ. of california berkeley, univ. of colorado denver and nag ltd..-- use stdlib_blas_constants_dp, only: one ! Scalar Arguments integer(${ik}$), intent(in) :: n ! Array Arguments integer(${ik}$), intent(inout) :: iseed(4_${ik}$) real(dp), intent(out) :: x(n) ! ===================================================================== ! Parameters integer(${ik}$), parameter :: lv = 128_${ik}$ integer(${ik}$), parameter :: ipw2 = 4096_${ik}$ real(dp), parameter :: r = one/ipw2 ! Local Scalars integer(${ik}$) :: i, i1, i2, i3, i4, it1, it2, it3, it4 ! Local Arrays integer(${ik}$) :: mm(lv,4_${ik}$) ! Intrinsic Functions ! Data Statements mm(1_${ik}$,1_${ik}$:4_${ik}$)=[494_${ik}$,322_${ik}$,2508_${ik}$,2549_${ik}$] mm(2_${ik}$,1_${ik}$:4_${ik}$)=[2637_${ik}$,789_${ik}$,3754_${ik}$,1145_${ik}$] mm(3_${ik}$,1_${ik}$:4_${ik}$)=[255_${ik}$,1440_${ik}$,1766_${ik}$,2253_${ik}$] mm(4_${ik}$,1_${ik}$:4_${ik}$)=[2008_${ik}$,752_${ik}$,3572_${ik}$,305_${ik}$] mm(5_${ik}$,1_${ik}$:4_${ik}$)=[1253_${ik}$,2859_${ik}$,2893_${ik}$,3301_${ik}$] mm(6_${ik}$,1_${ik}$:4_${ik}$)=[3344_${ik}$,123_${ik}$,307_${ik}$,1065_${ik}$] mm(7_${ik}$,1_${ik}$:4_${ik}$)=[4084_${ik}$,1848_${ik}$,1297_${ik}$,3133_${ik}$] mm(8_${ik}$,1_${ik}$:4_${ik}$)=[1739_${ik}$,643_${ik}$,3966_${ik}$,2913_${ik}$] mm(9_${ik}$,1_${ik}$:4_${ik}$)=[3143_${ik}$,2405_${ik}$,758_${ik}$,3285_${ik}$] mm(10_${ik}$,1_${ik}$:4_${ik}$)=[3468_${ik}$,2638_${ik}$,2598_${ik}$,1241_${ik}$] mm(11_${ik}$,1_${ik}$:4_${ik}$)=[688_${ik}$,2344_${ik}$,3406_${ik}$,1197_${ik}$] mm(12_${ik}$,1_${ik}$:4_${ik}$)=[1657_${ik}$,46_${ik}$,2922_${ik}$,3729_${ik}$] mm(13_${ik}$,1_${ik}$:4_${ik}$)=[1238_${ik}$,3814_${ik}$,1038_${ik}$,2501_${ik}$] mm(14_${ik}$,1_${ik}$:4_${ik}$)=[3166_${ik}$,913_${ik}$,2934_${ik}$,1673_${ik}$] mm(15_${ik}$,1_${ik}$:4_${ik}$)=[1292_${ik}$,3649_${ik}$,2091_${ik}$,541_${ik}$] mm(16_${ik}$,1_${ik}$:4_${ik}$)=[3422_${ik}$,339_${ik}$,2451_${ik}$,2753_${ik}$] mm(17_${ik}$,1_${ik}$:4_${ik}$)=[1270_${ik}$,3808_${ik}$,1580_${ik}$,949_${ik}$] mm(18_${ik}$,1_${ik}$:4_${ik}$)=[2016_${ik}$,822_${ik}$,1958_${ik}$,2361_${ik}$] mm(19_${ik}$,1_${ik}$:4_${ik}$)=[154_${ik}$,2832_${ik}$,2055_${ik}$,1165_${ik}$] mm(20_${ik}$,1_${ik}$:4_${ik}$)=[2862_${ik}$,3078_${ik}$,1507_${ik}$,4081_${ik}$] mm(21_${ik}$,1_${ik}$:4_${ik}$)=[697_${ik}$,3633_${ik}$,1078_${ik}$,2725_${ik}$] mm(22_${ik}$,1_${ik}$:4_${ik}$)=[1706_${ik}$,2970_${ik}$,3273_${ik}$,3305_${ik}$] mm(23_${ik}$,1_${ik}$:4_${ik}$)=[491_${ik}$,637_${ik}$,17_${ik}$,3069_${ik}$] mm(24_${ik}$,1_${ik}$:4_${ik}$)=[931_${ik}$,2249_${ik}$,854_${ik}$,3617_${ik}$] mm(25_${ik}$,1_${ik}$:4_${ik}$)=[1444_${ik}$,2081_${ik}$,2916_${ik}$,3733_${ik}$] mm(26_${ik}$,1_${ik}$:4_${ik}$)=[444_${ik}$,4019_${ik}$,3971_${ik}$,409_${ik}$] mm(27_${ik}$,1_${ik}$:4_${ik}$)=[3577_${ik}$,1478_${ik}$,2889_${ik}$,2157_${ik}$] mm(28_${ik}$,1_${ik}$:4_${ik}$)=[3944_${ik}$,242_${ik}$,3831_${ik}$,1361_${ik}$] mm(29_${ik}$,1_${ik}$:4_${ik}$)=[2184_${ik}$,481_${ik}$,2621_${ik}$,3973_${ik}$] mm(30_${ik}$,1_${ik}$:4_${ik}$)=[1661_${ik}$,2075_${ik}$,1541_${ik}$,1865_${ik}$] mm(31_${ik}$,1_${ik}$:4_${ik}$)=[3482_${ik}$,4058_${ik}$,893_${ik}$,2525_${ik}$] mm(32_${ik}$,1_${ik}$:4_${ik}$)=[657_${ik}$,622_${ik}$,736_${ik}$,1409_${ik}$] mm(33_${ik}$,1_${ik}$:4_${ik}$)=[3023_${ik}$,3376_${ik}$,3992_${ik}$,3445_${ik}$] mm(34_${ik}$,1_${ik}$:4_${ik}$)=[3618_${ik}$,812_${ik}$,787_${ik}$,3577_${ik}$] mm(35_${ik}$,1_${ik}$:4_${ik}$)=[1267_${ik}$,234_${ik}$,2125_${ik}$,77_${ik}$] mm(36_${ik}$,1_${ik}$:4_${ik}$)=[1828_${ik}$,641_${ik}$,2364_${ik}$,3761_${ik}$] mm(37_${ik}$,1_${ik}$:4_${ik}$)=[164_${ik}$,4005_${ik}$,2460_${ik}$,2149_${ik}$] mm(38_${ik}$,1_${ik}$:4_${ik}$)=[3798_${ik}$,1122_${ik}$,257_${ik}$,1449_${ik}$] mm(39_${ik}$,1_${ik}$:4_${ik}$)=[3087_${ik}$,3135_${ik}$,1574_${ik}$,3005_${ik}$] mm(40_${ik}$,1_${ik}$:4_${ik}$)=[2400_${ik}$,2640_${ik}$,3912_${ik}$,225_${ik}$] mm(41_${ik}$,1_${ik}$:4_${ik}$)=[2870_${ik}$,2302_${ik}$,1216_${ik}$,85_${ik}$] mm(42_${ik}$,1_${ik}$:4_${ik}$)=[3876_${ik}$,40_${ik}$,3248_${ik}$,3673_${ik}$] mm(43_${ik}$,1_${ik}$:4_${ik}$)=[1905_${ik}$,1832_${ik}$,3401_${ik}$,3117_${ik}$] mm(44_${ik}$,1_${ik}$:4_${ik}$)=[1593_${ik}$,2247_${ik}$,2124_${ik}$,3089_${ik}$] mm(45_${ik}$,1_${ik}$:4_${ik}$)=[1797_${ik}$,2034_${ik}$,2762_${ik}$,1349_${ik}$] mm(46_${ik}$,1_${ik}$:4_${ik}$)=[1234_${ik}$,2637_${ik}$,149_${ik}$,2057_${ik}$] mm(47_${ik}$,1_${ik}$:4_${ik}$)=[3460_${ik}$,1287_${ik}$,2245_${ik}$,413_${ik}$] mm(48_${ik}$,1_${ik}$:4_${ik}$)=[328_${ik}$,1691_${ik}$,166_${ik}$,65_${ik}$] mm(49_${ik}$,1_${ik}$:4_${ik}$)=[2861_${ik}$,496_${ik}$,466_${ik}$,1845_${ik}$] mm(50_${ik}$,1_${ik}$:4_${ik}$)=[1950_${ik}$,1597_${ik}$,4018_${ik}$,697_${ik}$] mm(51_${ik}$,1_${ik}$:4_${ik}$)=[617_${ik}$,2394_${ik}$,1399_${ik}$,3085_${ik}$] mm(52_${ik}$,1_${ik}$:4_${ik}$)=[2070_${ik}$,2584_${ik}$,190_${ik}$,3441_${ik}$] mm(53_${ik}$,1_${ik}$:4_${ik}$)=[3331_${ik}$,1843_${ik}$,2879_${ik}$,1573_${ik}$] mm(54_${ik}$,1_${ik}$:4_${ik}$)=[769_${ik}$,336_${ik}$,153_${ik}$,3689_${ik}$] mm(55_${ik}$,1_${ik}$:4_${ik}$)=[1558_${ik}$,1472_${ik}$,2320_${ik}$,2941_${ik}$] mm(56_${ik}$,1_${ik}$:4_${ik}$)=[2412_${ik}$,2407_${ik}$,18_${ik}$,929_${ik}$] mm(57_${ik}$,1_${ik}$:4_${ik}$)=[2800_${ik}$,433_${ik}$,712_${ik}$,533_${ik}$] mm(58_${ik}$,1_${ik}$:4_${ik}$)=[189_${ik}$,2096_${ik}$,2159_${ik}$,2841_${ik}$] mm(59_${ik}$,1_${ik}$:4_${ik}$)=[287_${ik}$,1761_${ik}$,2318_${ik}$,4077_${ik}$] mm(60_${ik}$,1_${ik}$:4_${ik}$)=[2045_${ik}$,2810_${ik}$,2091_${ik}$,721_${ik}$] mm(61_${ik}$,1_${ik}$:4_${ik}$)=[1227_${ik}$,566_${ik}$,3443_${ik}$,2821_${ik}$] mm(62_${ik}$,1_${ik}$:4_${ik}$)=[2838_${ik}$,442_${ik}$,1510_${ik}$,2249_${ik}$] mm(63_${ik}$,1_${ik}$:4_${ik}$)=[209_${ik}$,41_${ik}$,449_${ik}$,2397_${ik}$] mm(64_${ik}$,1_${ik}$:4_${ik}$)=[2770_${ik}$,1238_${ik}$,1956_${ik}$,2817_${ik}$] mm(65_${ik}$,1_${ik}$:4_${ik}$)=[3654_${ik}$,1086_${ik}$,2201_${ik}$,245_${ik}$] mm(66_${ik}$,1_${ik}$:4_${ik}$)=[3993_${ik}$,603_${ik}$,3137_${ik}$,1913_${ik}$] mm(67_${ik}$,1_${ik}$:4_${ik}$)=[192_${ik}$,840_${ik}$,3399_${ik}$,1997_${ik}$] mm(68_${ik}$,1_${ik}$:4_${ik}$)=[2253_${ik}$,3168_${ik}$,1321_${ik}$,3121_${ik}$] mm(69_${ik}$,1_${ik}$:4_${ik}$)=[3491_${ik}$,1499_${ik}$,2271_${ik}$,997_${ik}$] mm(70_${ik}$,1_${ik}$:4_${ik}$)=[2889_${ik}$,1084_${ik}$,3667_${ik}$,1833_${ik}$] mm(71_${ik}$,1_${ik}$:4_${ik}$)=[2857_${ik}$,3438_${ik}$,2703_${ik}$,2877_${ik}$] mm(72_${ik}$,1_${ik}$:4_${ik}$)=[2094_${ik}$,2408_${ik}$,629_${ik}$,1633_${ik}$] mm(73_${ik}$,1_${ik}$:4_${ik}$)=[1818_${ik}$,1589_${ik}$,2365_${ik}$,981_${ik}$] mm(74_${ik}$,1_${ik}$:4_${ik}$)=[688_${ik}$,2391_${ik}$,2431_${ik}$,2009_${ik}$] mm(75_${ik}$,1_${ik}$:4_${ik}$)=[1407_${ik}$,288_${ik}$,1113_${ik}$,941_${ik}$] mm(76_${ik}$,1_${ik}$:4_${ik}$)=[634_${ik}$,26_${ik}$,3922_${ik}$,2449_${ik}$] mm(77_${ik}$,1_${ik}$:4_${ik}$)=[3231_${ik}$,512_${ik}$,2554_${ik}$,197_${ik}$] mm(78_${ik}$,1_${ik}$:4_${ik}$)=[815_${ik}$,1456_${ik}$,184_${ik}$,2441_${ik}$] mm(79_${ik}$,1_${ik}$:4_${ik}$)=[3524_${ik}$,171_${ik}$,2099_${ik}$,285_${ik}$] mm(80_${ik}$,1_${ik}$:4_${ik}$)=[1914_${ik}$,1677_${ik}$,3228_${ik}$,1473_${ik}$] mm(81_${ik}$,1_${ik}$:4_${ik}$)=[516_${ik}$,2657_${ik}$,4012_${ik}$,2741_${ik}$] mm(82_${ik}$,1_${ik}$:4_${ik}$)=[164_${ik}$,2270_${ik}$,1921_${ik}$,3129_${ik}$] mm(83_${ik}$,1_${ik}$:4_${ik}$)=[303_${ik}$,2587_${ik}$,3452_${ik}$,909_${ik}$] mm(84_${ik}$,1_${ik}$:4_${ik}$)=[2144_${ik}$,2961_${ik}$,3901_${ik}$,2801_${ik}$] mm(85_${ik}$,1_${ik}$:4_${ik}$)=[3480_${ik}$,1970_${ik}$,572_${ik}$,421_${ik}$] mm(86_${ik}$,1_${ik}$:4_${ik}$)=[119_${ik}$,1817_${ik}$,3309_${ik}$,4073_${ik}$] mm(87_${ik}$,1_${ik}$:4_${ik}$)=[3357_${ik}$,676_${ik}$,3171_${ik}$,2813_${ik}$] mm(88_${ik}$,1_${ik}$:4_${ik}$)=[837_${ik}$,1410_${ik}$,817_${ik}$,2337_${ik}$] mm(89_${ik}$,1_${ik}$:4_${ik}$)=[2826_${ik}$,3723_${ik}$,3039_${ik}$,1429_${ik}$] mm(90_${ik}$,1_${ik}$:4_${ik}$)=[2332_${ik}$,2803_${ik}$,1696_${ik}$,1177_${ik}$] mm(91_${ik}$,1_${ik}$:4_${ik}$)=[2089_${ik}$,3185_${ik}$,1256_${ik}$,1901_${ik}$] mm(92_${ik}$,1_${ik}$:4_${ik}$)=[3780_${ik}$,184_${ik}$,3715_${ik}$,81_${ik}$] mm(93_${ik}$,1_${ik}$:4_${ik}$)=[1700_${ik}$,663_${ik}$,2077_${ik}$,1669_${ik}$] mm(94_${ik}$,1_${ik}$:4_${ik}$)=[3712_${ik}$,499_${ik}$,3019_${ik}$,2633_${ik}$] mm(95_${ik}$,1_${ik}$:4_${ik}$)=[150_${ik}$,3784_${ik}$,1497_${ik}$,2269_${ik}$] mm(96_${ik}$,1_${ik}$:4_${ik}$)=[2000_${ik}$,1631_${ik}$,1101_${ik}$,129_${ik}$] mm(97_${ik}$,1_${ik}$:4_${ik}$)=[3375_${ik}$,1925_${ik}$,717_${ik}$,1141_${ik}$] mm(98_${ik}$,1_${ik}$:4_${ik}$)=[1621_${ik}$,3912_${ik}$,51_${ik}$,249_${ik}$] mm(99_${ik}$,1_${ik}$:4_${ik}$)=[3090_${ik}$,1398_${ik}$,981_${ik}$,3917_${ik}$] mm(100_${ik}$,1_${ik}$:4_${ik}$)=[3765_${ik}$,1349_${ik}$,1978_${ik}$,2481_${ik}$] mm(101_${ik}$,1_${ik}$:4_${ik}$)=[1149_${ik}$,1441_${ik}$,1813_${ik}$,3941_${ik}$] mm(102_${ik}$,1_${ik}$:4_${ik}$)=[3146_${ik}$,2224_${ik}$,3881_${ik}$,2217_${ik}$] mm(103_${ik}$,1_${ik}$:4_${ik}$)=[33_${ik}$,2411_${ik}$,76_${ik}$,2749_${ik}$] mm(104_${ik}$,1_${ik}$:4_${ik}$)=[3082_${ik}$,1907_${ik}$,3846_${ik}$,3041_${ik}$] mm(105_${ik}$,1_${ik}$:4_${ik}$)=[2741_${ik}$,3192_${ik}$,3694_${ik}$,1877_${ik}$] mm(106_${ik}$,1_${ik}$:4_${ik}$)=[359_${ik}$,2786_${ik}$,1682_${ik}$,345_${ik}$] mm(107_${ik}$,1_${ik}$:4_${ik}$)=[3316_${ik}$,382_${ik}$,124_${ik}$,2861_${ik}$] mm(108_${ik}$,1_${ik}$:4_${ik}$)=[1749_${ik}$,37_${ik}$,1660_${ik}$,1809_${ik}$] mm(109_${ik}$,1_${ik}$:4_${ik}$)=[185_${ik}$,759_${ik}$,3997_${ik}$,3141_${ik}$] mm(110_${ik}$,1_${ik}$:4_${ik}$)=[2784_${ik}$,2948_${ik}$,479_${ik}$,2825_${ik}$] mm(111_${ik}$,1_${ik}$:4_${ik}$)=[2202_${ik}$,1862_${ik}$,1141_${ik}$,157_${ik}$] mm(112_${ik}$,1_${ik}$:4_${ik}$)=[2199_${ik}$,3802_${ik}$,886_${ik}$,2881_${ik}$] mm(113_${ik}$,1_${ik}$:4_${ik}$)=[1364_${ik}$,2423_${ik}$,3514_${ik}$,3637_${ik}$] mm(114_${ik}$,1_${ik}$:4_${ik}$)=[1244_${ik}$,2051_${ik}$,1301_${ik}$,1465_${ik}$] mm(115_${ik}$,1_${ik}$:4_${ik}$)=[2020_${ik}$,2295_${ik}$,3604_${ik}$,2829_${ik}$] mm(116_${ik}$,1_${ik}$:4_${ik}$)=[3160_${ik}$,1332_${ik}$,1888_${ik}$,2161_${ik}$] mm(117_${ik}$,1_${ik}$:4_${ik}$)=[2785_${ik}$,1832_${ik}$,1836_${ik}$,3365_${ik}$] mm(118_${ik}$,1_${ik}$:4_${ik}$)=[2772_${ik}$,2405_${ik}$,1990_${ik}$,361_${ik}$] mm(119_${ik}$,1_${ik}$:4_${ik}$)=[1217_${ik}$,3638_${ik}$,2058_${ik}$,2685_${ik}$] mm(120_${ik}$,1_${ik}$:4_${ik}$)=[1822_${ik}$,3661_${ik}$,692_${ik}$,3745_${ik}$] mm(121_${ik}$,1_${ik}$:4_${ik}$)=[1245_${ik}$,327_${ik}$,1194_${ik}$,2325_${ik}$] mm(122_${ik}$,1_${ik}$:4_${ik}$)=[2252_${ik}$,3660_${ik}$,20_${ik}$,3609_${ik}$] mm(123_${ik}$,1_${ik}$:4_${ik}$)=[3904_${ik}$,716_${ik}$,3285_${ik}$,3821_${ik}$] mm(124_${ik}$,1_${ik}$:4_${ik}$)=[2774_${ik}$,1842_${ik}$,2046_${ik}$,3537_${ik}$] mm(125_${ik}$,1_${ik}$:4_${ik}$)=[997_${ik}$,3987_${ik}$,2107_${ik}$,517_${ik}$] mm(126_${ik}$,1_${ik}$:4_${ik}$)=[2573_${ik}$,1368_${ik}$,3508_${ik}$,3017_${ik}$] mm(127_${ik}$,1_${ik}$:4_${ik}$)=[1148_${ik}$,1848_${ik}$,3525_${ik}$,2141_${ik}$] mm(128_${ik}$,1_${ik}$:4_${ik}$)=[545_${ik}$,2366_${ik}$,3801_${ik}$,1537_${ik}$] ! Executable Statements i1 = iseed( 1_${ik}$ ) i2 = iseed( 2_${ik}$ ) i3 = iseed( 3_${ik}$ ) i4 = iseed( 4_${ik}$ ) loop_10: do i = 1, min( n, lv ) 20 continue ! multiply the seed by i-th power of the multiplier modulo 2**48 it4 = i4*mm( i, 4_${ik}$ ) it3 = it4 / ipw2 it4 = it4 - ipw2*it3 it3 = it3 + i3*mm( i, 4_${ik}$ ) + i4*mm( i, 3_${ik}$ ) it2 = it3 / ipw2 it3 = it3 - ipw2*it2 it2 = it2 + i2*mm( i, 4_${ik}$ ) + i3*mm( i, 3_${ik}$ ) + i4*mm( i, 2_${ik}$ ) it1 = it2 / ipw2 it2 = it2 - ipw2*it1 it1 = it1 + i1*mm( i, 4_${ik}$ ) + i2*mm( i, 3_${ik}$ ) + i3*mm( i, 2_${ik}$ ) +i4*mm( i, 1_${ik}$ ) it1 = mod( it1, ipw2 ) ! convert 48-bit integer to a realnumber in the interval (0,1,KIND=dp) x( i ) = r*( real( it1,KIND=dp)+r*( real( it2,KIND=dp)+r*( real( it3,KIND=dp)+& r*real( it4,KIND=dp) ) ) ) if (x( i )==one) then ! if a real number has n bits of precision, and the first ! n bits of the 48-bit integer above happen to be all 1 (which ! will occur about once every 2**n calls), then x( i ) will ! be rounded to exactly one. ! since x( i ) is not supposed to return exactly 0.0_dp or one, ! the statistically correct thing to do in this situation is ! simply to iterate again. ! n.b. the case x( i ) = 0.0_dp should not be possible. i1 = i1 + 2_${ik}$ i2 = i2 + 2_${ik}$ i3 = i3 + 2_${ik}$ i4 = i4 + 2_${ik}$ goto 20 end if end do loop_10 ! return final value of seed iseed( 1_${ik}$ ) = it1 iseed( 2_${ik}$ ) = it2 iseed( 3_${ik}$ ) = it3 iseed( 4_${ik}$ ) = it4 return end subroutine stdlib${ii}$_dlaruv #:for rk,rt,ri in REAL_KINDS_TYPES #:if not rk in ["sp","dp"] pure module subroutine stdlib${ii}$_${ri}$laruv( iseed, n, x ) !! DLARUV: returns a vector of n random real numbers from a uniform (0,1) !! distribution (n <= 128). !! This is an auxiliary routine called by DLARNV and ZLARNV. ! -- lapack auxiliary routine -- ! -- lapack is a software package provided by univ. of tennessee, -- ! -- univ. of california berkeley, univ. of colorado denver and nag ltd..-- use stdlib_blas_constants_${rk}$, only: one ! Scalar Arguments integer(${ik}$), intent(in) :: n ! Array Arguments integer(${ik}$), intent(inout) :: iseed(4_${ik}$) real(${rk}$), intent(out) :: x(n) ! ===================================================================== ! Parameters integer(${ik}$), parameter :: lv = 128_${ik}$ integer(${ik}$), parameter :: ipw2 = 4096_${ik}$ real(${rk}$), parameter :: r = one/ipw2 ! Local Scalars integer(${ik}$) :: i, i1, i2, i3, i4, it1, it2, it3, it4, j ! Local Arrays integer(${ik}$) :: mm(lv,4_${ik}$) ! Intrinsic Functions ! Data Statements mm(1_${ik}$,1_${ik}$:4_${ik}$)=[494_${ik}$,322_${ik}$,2508_${ik}$,2549_${ik}$] mm(2_${ik}$,1_${ik}$:4_${ik}$)=[2637_${ik}$,789_${ik}$,3754_${ik}$,1145_${ik}$] mm(3_${ik}$,1_${ik}$:4_${ik}$)=[255_${ik}$,1440_${ik}$,1766_${ik}$,2253_${ik}$] mm(4_${ik}$,1_${ik}$:4_${ik}$)=[2008_${ik}$,752_${ik}$,3572_${ik}$,305_${ik}$] mm(5_${ik}$,1_${ik}$:4_${ik}$)=[1253_${ik}$,2859_${ik}$,2893_${ik}$,3301_${ik}$] mm(6_${ik}$,1_${ik}$:4_${ik}$)=[3344_${ik}$,123_${ik}$,307_${ik}$,1065_${ik}$] mm(7_${ik}$,1_${ik}$:4_${ik}$)=[4084_${ik}$,1848_${ik}$,1297_${ik}$,3133_${ik}$] mm(8_${ik}$,1_${ik}$:4_${ik}$)=[1739_${ik}$,643_${ik}$,3966_${ik}$,2913_${ik}$] mm(9_${ik}$,1_${ik}$:4_${ik}$)=[3143_${ik}$,2405_${ik}$,758_${ik}$,3285_${ik}$] mm(10_${ik}$,1_${ik}$:4_${ik}$)=[3468_${ik}$,2638_${ik}$,2598_${ik}$,1241_${ik}$] mm(11_${ik}$,1_${ik}$:4_${ik}$)=[688_${ik}$,2344_${ik}$,3406_${ik}$,1197_${ik}$] mm(12_${ik}$,1_${ik}$:4_${ik}$)=[1657_${ik}$,46_${ik}$,2922_${ik}$,3729_${ik}$] mm(13_${ik}$,1_${ik}$:4_${ik}$)=[1238_${ik}$,3814_${ik}$,1038_${ik}$,2501_${ik}$] mm(14_${ik}$,1_${ik}$:4_${ik}$)=[3166_${ik}$,913_${ik}$,2934_${ik}$,1673_${ik}$] mm(15_${ik}$,1_${ik}$:4_${ik}$)=[1292_${ik}$,3649_${ik}$,2091_${ik}$,541_${ik}$] mm(16_${ik}$,1_${ik}$:4_${ik}$)=[3422_${ik}$,339_${ik}$,2451_${ik}$,2753_${ik}$] mm(17_${ik}$,1_${ik}$:4_${ik}$)=[1270_${ik}$,3808_${ik}$,1580_${ik}$,949_${ik}$] mm(18_${ik}$,1_${ik}$:4_${ik}$)=[2016_${ik}$,822_${ik}$,1958_${ik}$,2361_${ik}$] mm(19_${ik}$,1_${ik}$:4_${ik}$)=[154_${ik}$,2832_${ik}$,2055_${ik}$,1165_${ik}$] mm(20_${ik}$,1_${ik}$:4_${ik}$)=[2862_${ik}$,3078_${ik}$,1507_${ik}$,4081_${ik}$] mm(21_${ik}$,1_${ik}$:4_${ik}$)=[697_${ik}$,3633_${ik}$,1078_${ik}$,2725_${ik}$] mm(22_${ik}$,1_${ik}$:4_${ik}$)=[1706_${ik}$,2970_${ik}$,3273_${ik}$,3305_${ik}$] mm(23_${ik}$,1_${ik}$:4_${ik}$)=[491_${ik}$,637_${ik}$,17_${ik}$,3069_${ik}$] mm(24_${ik}$,1_${ik}$:4_${ik}$)=[931_${ik}$,2249_${ik}$,854_${ik}$,3617_${ik}$] mm(25_${ik}$,1_${ik}$:4_${ik}$)=[1444_${ik}$,2081_${ik}$,2916_${ik}$,3733_${ik}$] mm(26_${ik}$,1_${ik}$:4_${ik}$)=[444_${ik}$,4019_${ik}$,3971_${ik}$,409_${ik}$] mm(27_${ik}$,1_${ik}$:4_${ik}$)=[3577_${ik}$,1478_${ik}$,2889_${ik}$,2157_${ik}$] mm(28_${ik}$,1_${ik}$:4_${ik}$)=[3944_${ik}$,242_${ik}$,3831_${ik}$,1361_${ik}$] mm(29_${ik}$,1_${ik}$:4_${ik}$)=[2184_${ik}$,481_${ik}$,2621_${ik}$,3973_${ik}$] mm(30_${ik}$,1_${ik}$:4_${ik}$)=[1661_${ik}$,2075_${ik}$,1541_${ik}$,1865_${ik}$] mm(31_${ik}$,1_${ik}$:4_${ik}$)=[3482_${ik}$,4058_${ik}$,893_${ik}$,2525_${ik}$] mm(32_${ik}$,1_${ik}$:4_${ik}$)=[657_${ik}$,622_${ik}$,736_${ik}$,1409_${ik}$] mm(33_${ik}$,1_${ik}$:4_${ik}$)=[3023_${ik}$,3376_${ik}$,3992_${ik}$,3445_${ik}$] mm(34_${ik}$,1_${ik}$:4_${ik}$)=[3618_${ik}$,812_${ik}$,787_${ik}$,3577_${ik}$] mm(35_${ik}$,1_${ik}$:4_${ik}$)=[1267_${ik}$,234_${ik}$,2125_${ik}$,77_${ik}$] mm(36_${ik}$,1_${ik}$:4_${ik}$)=[1828_${ik}$,641_${ik}$,2364_${ik}$,3761_${ik}$] mm(37_${ik}$,1_${ik}$:4_${ik}$)=[164_${ik}$,4005_${ik}$,2460_${ik}$,2149_${ik}$] mm(38_${ik}$,1_${ik}$:4_${ik}$)=[3798_${ik}$,1122_${ik}$,257_${ik}$,1449_${ik}$] mm(39_${ik}$,1_${ik}$:4_${ik}$)=[3087_${ik}$,3135_${ik}$,1574_${ik}$,3005_${ik}$] mm(40_${ik}$,1_${ik}$:4_${ik}$)=[2400_${ik}$,2640_${ik}$,3912_${ik}$,225_${ik}$] mm(41_${ik}$,1_${ik}$:4_${ik}$)=[2870_${ik}$,2302_${ik}$,1216_${ik}$,85_${ik}$] mm(42_${ik}$,1_${ik}$:4_${ik}$)=[3876_${ik}$,40_${ik}$,3248_${ik}$,3673_${ik}$] mm(43_${ik}$,1_${ik}$:4_${ik}$)=[1905_${ik}$,1832_${ik}$,3401_${ik}$,3117_${ik}$] mm(44_${ik}$,1_${ik}$:4_${ik}$)=[1593_${ik}$,2247_${ik}$,2124_${ik}$,3089_${ik}$] mm(45_${ik}$,1_${ik}$:4_${ik}$)=[1797_${ik}$,2034_${ik}$,2762_${ik}$,1349_${ik}$] mm(46_${ik}$,1_${ik}$:4_${ik}$)=[1234_${ik}$,2637_${ik}$,149_${ik}$,2057_${ik}$] mm(47_${ik}$,1_${ik}$:4_${ik}$)=[3460_${ik}$,1287_${ik}$,2245_${ik}$,413_${ik}$] mm(48_${ik}$,1_${ik}$:4_${ik}$)=[328_${ik}$,1691_${ik}$,166_${ik}$,65_${ik}$] mm(49_${ik}$,1_${ik}$:4_${ik}$)=[2861_${ik}$,496_${ik}$,466_${ik}$,1845_${ik}$] mm(50_${ik}$,1_${ik}$:4_${ik}$)=[1950_${ik}$,1597_${ik}$,4018_${ik}$,697_${ik}$] mm(51_${ik}$,1_${ik}$:4_${ik}$)=[617_${ik}$,2394_${ik}$,1399_${ik}$,3085_${ik}$] mm(52_${ik}$,1_${ik}$:4_${ik}$)=[2070_${ik}$,2584_${ik}$,190_${ik}$,3441_${ik}$] mm(53_${ik}$,1_${ik}$:4_${ik}$)=[3331_${ik}$,1843_${ik}$,2879_${ik}$,1573_${ik}$] mm(54_${ik}$,1_${ik}$:4_${ik}$)=[769_${ik}$,336_${ik}$,153_${ik}$,3689_${ik}$] mm(55_${ik}$,1_${ik}$:4_${ik}$)=[1558_${ik}$,1472_${ik}$,2320_${ik}$,2941_${ik}$] mm(56_${ik}$,1_${ik}$:4_${ik}$)=[2412_${ik}$,2407_${ik}$,18_${ik}$,929_${ik}$] mm(57_${ik}$,1_${ik}$:4_${ik}$)=[2800_${ik}$,433_${ik}$,712_${ik}$,533_${ik}$] mm(58_${ik}$,1_${ik}$:4_${ik}$)=[189_${ik}$,2096_${ik}$,2159_${ik}$,2841_${ik}$] mm(59_${ik}$,1_${ik}$:4_${ik}$)=[287_${ik}$,1761_${ik}$,2318_${ik}$,4077_${ik}$] mm(60_${ik}$,1_${ik}$:4_${ik}$)=[2045_${ik}$,2810_${ik}$,2091_${ik}$,721_${ik}$] mm(61_${ik}$,1_${ik}$:4_${ik}$)=[1227_${ik}$,566_${ik}$,3443_${ik}$,2821_${ik}$] mm(62_${ik}$,1_${ik}$:4_${ik}$)=[2838_${ik}$,442_${ik}$,1510_${ik}$,2249_${ik}$] mm(63_${ik}$,1_${ik}$:4_${ik}$)=[209_${ik}$,41_${ik}$,449_${ik}$,2397_${ik}$] mm(64_${ik}$,1_${ik}$:4_${ik}$)=[2770_${ik}$,1238_${ik}$,1956_${ik}$,2817_${ik}$] mm(65_${ik}$,1_${ik}$:4_${ik}$)=[3654_${ik}$,1086_${ik}$,2201_${ik}$,245_${ik}$] mm(66_${ik}$,1_${ik}$:4_${ik}$)=[3993_${ik}$,603_${ik}$,3137_${ik}$,1913_${ik}$] mm(67_${ik}$,1_${ik}$:4_${ik}$)=[192_${ik}$,840_${ik}$,3399_${ik}$,1997_${ik}$] mm(68_${ik}$,1_${ik}$:4_${ik}$)=[2253_${ik}$,3168_${ik}$,1321_${ik}$,3121_${ik}$] mm(69_${ik}$,1_${ik}$:4_${ik}$)=[3491_${ik}$,1499_${ik}$,2271_${ik}$,997_${ik}$] mm(70_${ik}$,1_${ik}$:4_${ik}$)=[2889_${ik}$,1084_${ik}$,3667_${ik}$,1833_${ik}$] mm(71_${ik}$,1_${ik}$:4_${ik}$)=[2857_${ik}$,3438_${ik}$,2703_${ik}$,2877_${ik}$] mm(72_${ik}$,1_${ik}$:4_${ik}$)=[2094_${ik}$,2408_${ik}$,629_${ik}$,1633_${ik}$] mm(73_${ik}$,1_${ik}$:4_${ik}$)=[1818_${ik}$,1589_${ik}$,2365_${ik}$,981_${ik}$] mm(74_${ik}$,1_${ik}$:4_${ik}$)=[688_${ik}$,2391_${ik}$,2431_${ik}$,2009_${ik}$] mm(75_${ik}$,1_${ik}$:4_${ik}$)=[1407_${ik}$,288_${ik}$,1113_${ik}$,941_${ik}$] mm(76_${ik}$,1_${ik}$:4_${ik}$)=[634_${ik}$,26_${ik}$,3922_${ik}$,2449_${ik}$] mm(77_${ik}$,1_${ik}$:4_${ik}$)=[3231_${ik}$,512_${ik}$,2554_${ik}$,197_${ik}$] mm(78_${ik}$,1_${ik}$:4_${ik}$)=[815_${ik}$,1456_${ik}$,184_${ik}$,2441_${ik}$] mm(79_${ik}$,1_${ik}$:4_${ik}$)=[3524_${ik}$,171_${ik}$,2099_${ik}$,285_${ik}$] mm(80_${ik}$,1_${ik}$:4_${ik}$)=[1914_${ik}$,1677_${ik}$,3228_${ik}$,1473_${ik}$] mm(81_${ik}$,1_${ik}$:4_${ik}$)=[516_${ik}$,2657_${ik}$,4012_${ik}$,2741_${ik}$] mm(82_${ik}$,1_${ik}$:4_${ik}$)=[164_${ik}$,2270_${ik}$,1921_${ik}$,3129_${ik}$] mm(83_${ik}$,1_${ik}$:4_${ik}$)=[303_${ik}$,2587_${ik}$,3452_${ik}$,909_${ik}$] mm(84_${ik}$,1_${ik}$:4_${ik}$)=[2144_${ik}$,2961_${ik}$,3901_${ik}$,2801_${ik}$] mm(85_${ik}$,1_${ik}$:4_${ik}$)=[3480_${ik}$,1970_${ik}$,572_${ik}$,421_${ik}$] mm(86_${ik}$,1_${ik}$:4_${ik}$)=[119_${ik}$,1817_${ik}$,3309_${ik}$,4073_${ik}$] mm(87_${ik}$,1_${ik}$:4_${ik}$)=[3357_${ik}$,676_${ik}$,3171_${ik}$,2813_${ik}$] mm(88_${ik}$,1_${ik}$:4_${ik}$)=[837_${ik}$,1410_${ik}$,817_${ik}$,2337_${ik}$] mm(89_${ik}$,1_${ik}$:4_${ik}$)=[2826_${ik}$,3723_${ik}$,3039_${ik}$,1429_${ik}$] mm(90_${ik}$,1_${ik}$:4_${ik}$)=[2332_${ik}$,2803_${ik}$,1696_${ik}$,1177_${ik}$] mm(91_${ik}$,1_${ik}$:4_${ik}$)=[2089_${ik}$,3185_${ik}$,1256_${ik}$,1901_${ik}$] mm(92_${ik}$,1_${ik}$:4_${ik}$)=[3780_${ik}$,184_${ik}$,3715_${ik}$,81_${ik}$] mm(93_${ik}$,1_${ik}$:4_${ik}$)=[1700_${ik}$,663_${ik}$,2077_${ik}$,1669_${ik}$] mm(94_${ik}$,1_${ik}$:4_${ik}$)=[3712_${ik}$,499_${ik}$,3019_${ik}$,2633_${ik}$] mm(95_${ik}$,1_${ik}$:4_${ik}$)=[150_${ik}$,3784_${ik}$,1497_${ik}$,2269_${ik}$] mm(96_${ik}$,1_${ik}$:4_${ik}$)=[2000_${ik}$,1631_${ik}$,1101_${ik}$,129_${ik}$] mm(97_${ik}$,1_${ik}$:4_${ik}$)=[3375_${ik}$,1925_${ik}$,717_${ik}$,1141_${ik}$] mm(98_${ik}$,1_${ik}$:4_${ik}$)=[1621_${ik}$,3912_${ik}$,51_${ik}$,249_${ik}$] mm(99_${ik}$,1_${ik}$:4_${ik}$)=[3090_${ik}$,1398_${ik}$,981_${ik}$,3917_${ik}$] mm(100_${ik}$,1_${ik}$:4_${ik}$)=[3765_${ik}$,1349_${ik}$,1978_${ik}$,2481_${ik}$] mm(101_${ik}$,1_${ik}$:4_${ik}$)=[1149_${ik}$,1441_${ik}$,1813_${ik}$,3941_${ik}$] mm(102_${ik}$,1_${ik}$:4_${ik}$)=[3146_${ik}$,2224_${ik}$,3881_${ik}$,2217_${ik}$] mm(103_${ik}$,1_${ik}$:4_${ik}$)=[33_${ik}$,2411_${ik}$,76_${ik}$,2749_${ik}$] mm(104_${ik}$,1_${ik}$:4_${ik}$)=[3082_${ik}$,1907_${ik}$,3846_${ik}$,3041_${ik}$] mm(105_${ik}$,1_${ik}$:4_${ik}$)=[2741_${ik}$,3192_${ik}$,3694_${ik}$,1877_${ik}$] mm(106_${ik}$,1_${ik}$:4_${ik}$)=[359_${ik}$,2786_${ik}$,1682_${ik}$,345_${ik}$] mm(107_${ik}$,1_${ik}$:4_${ik}$)=[3316_${ik}$,382_${ik}$,124_${ik}$,2861_${ik}$] mm(108_${ik}$,1_${ik}$:4_${ik}$)=[1749_${ik}$,37_${ik}$,1660_${ik}$,1809_${ik}$] mm(109_${ik}$,1_${ik}$:4_${ik}$)=[185_${ik}$,759_${ik}$,3997_${ik}$,3141_${ik}$] mm(110_${ik}$,1_${ik}$:4_${ik}$)=[2784_${ik}$,2948_${ik}$,479_${ik}$,2825_${ik}$] mm(111_${ik}$,1_${ik}$:4_${ik}$)=[2202_${ik}$,1862_${ik}$,1141_${ik}$,157_${ik}$] mm(112_${ik}$,1_${ik}$:4_${ik}$)=[2199_${ik}$,3802_${ik}$,886_${ik}$,2881_${ik}$] mm(113_${ik}$,1_${ik}$:4_${ik}$)=[1364_${ik}$,2423_${ik}$,3514_${ik}$,3637_${ik}$] mm(114_${ik}$,1_${ik}$:4_${ik}$)=[1244_${ik}$,2051_${ik}$,1301_${ik}$,1465_${ik}$] mm(115_${ik}$,1_${ik}$:4_${ik}$)=[2020_${ik}$,2295_${ik}$,3604_${ik}$,2829_${ik}$] mm(116_${ik}$,1_${ik}$:4_${ik}$)=[3160_${ik}$,1332_${ik}$,1888_${ik}$,2161_${ik}$] mm(117_${ik}$,1_${ik}$:4_${ik}$)=[2785_${ik}$,1832_${ik}$,1836_${ik}$,3365_${ik}$] mm(118_${ik}$,1_${ik}$:4_${ik}$)=[2772_${ik}$,2405_${ik}$,1990_${ik}$,361_${ik}$] mm(119_${ik}$,1_${ik}$:4_${ik}$)=[1217_${ik}$,3638_${ik}$,2058_${ik}$,2685_${ik}$] mm(120_${ik}$,1_${ik}$:4_${ik}$)=[1822_${ik}$,3661_${ik}$,692_${ik}$,3745_${ik}$] mm(121_${ik}$,1_${ik}$:4_${ik}$)=[1245_${ik}$,327_${ik}$,1194_${ik}$,2325_${ik}$] mm(122_${ik}$,1_${ik}$:4_${ik}$)=[2252_${ik}$,3660_${ik}$,20_${ik}$,3609_${ik}$] mm(123_${ik}$,1_${ik}$:4_${ik}$)=[3904_${ik}$,716_${ik}$,3285_${ik}$,3821_${ik}$] mm(124_${ik}$,1_${ik}$:4_${ik}$)=[2774_${ik}$,1842_${ik}$,2046_${ik}$,3537_${ik}$] mm(125_${ik}$,1_${ik}$:4_${ik}$)=[997_${ik}$,3987_${ik}$,2107_${ik}$,517_${ik}$] mm(126_${ik}$,1_${ik}$:4_${ik}$)=[2573_${ik}$,1368_${ik}$,3508_${ik}$,3017_${ik}$] mm(127_${ik}$,1_${ik}$:4_${ik}$)=[1148_${ik}$,1848_${ik}$,3525_${ik}$,2141_${ik}$] mm(128_${ik}$,1_${ik}$:4_${ik}$)=[545_${ik}$,2366_${ik}$,3801_${ik}$,1537_${ik}$] ! Executable Statements i1 = iseed( 1_${ik}$ ) i2 = iseed( 2_${ik}$ ) i3 = iseed( 3_${ik}$ ) i4 = iseed( 4_${ik}$ ) loop_10: do i = 1, min( n, lv ) 20 continue ! multiply the seed by i-th power of the multiplier modulo 2**48 it4 = i4*mm( i, 4_${ik}$ ) it3 = it4 / ipw2 it4 = it4 - ipw2*it3 it3 = it3 + i3*mm( i, 4_${ik}$ ) + i4*mm( i, 3_${ik}$ ) it2 = it3 / ipw2 it3 = it3 - ipw2*it2 it2 = it2 + i2*mm( i, 4_${ik}$ ) + i3*mm( i, 3_${ik}$ ) + i4*mm( i, 2_${ik}$ ) it1 = it2 / ipw2 it2 = it2 - ipw2*it1 it1 = it1 + i1*mm( i, 4_${ik}$ ) + i2*mm( i, 3_${ik}$ ) + i3*mm( i, 2_${ik}$ ) +i4*mm( i, 1_${ik}$ ) it1 = mod( it1, ipw2 ) ! convert 48-bit integer to a realnumber in the interval (0,1,KIND=${rk}$) x( i ) = r*( real( it1,KIND=${rk}$)+r*( real( it2,KIND=${rk}$)+r*( real( it3,KIND=${rk}$)+& r*real( it4,KIND=${rk}$) ) ) ) if (x( i )==one) then ! if a real number has n bits of precision, and the first ! n bits of the 48-bit integer above happen to be all 1 (which ! will occur about once every 2**n calls), then x( i ) will ! be rounded to exactly one. ! since x( i ) is not supposed to return exactly 0.0_${rk}$ or one, ! the statistically correct thing to do in this situation is ! simply to iterate again. ! n.b. the case x( i ) = 0.0_${rk}$ should not be possible. i1 = i1 + 2_${ik}$ i2 = i2 + 2_${ik}$ i3 = i3 + 2_${ik}$ i4 = i4 + 2_${ik}$ goto 20 end if end do loop_10 ! return final value of seed iseed( 1_${ik}$ ) = it1 iseed( 2_${ik}$ ) = it2 iseed( 3_${ik}$ ) = it3 iseed( 4_${ik}$ ) = it4 return end subroutine stdlib${ii}$_${ri}$laruv #:endif #:endfor pure module subroutine stdlib${ii}$_slacpy( uplo, m, n, a, lda, b, ldb ) !! SLACPY copies all or part of a two-dimensional matrix A to another !! matrix B. ! -- lapack auxiliary routine -- ! -- lapack is a software package provided by univ. of tennessee, -- ! -- univ. of california berkeley, univ. of colorado denver and nag ltd..-- ! Scalar Arguments character, intent(in) :: uplo integer(${ik}$), intent(in) :: lda, ldb, m, n ! Array Arguments real(sp), intent(in) :: a(lda,*) real(sp), intent(out) :: b(ldb,*) ! ===================================================================== ! Local Scalars integer(${ik}$) :: i, j ! Intrinsic Functions ! Executable Statements if( stdlib_lsame( uplo, 'U' ) ) then do j = 1, n do i = 1, min( j, m ) b( i, j ) = a( i, j ) end do end do else if( stdlib_lsame( uplo, 'L' ) ) then do j = 1, n do i = j, m b( i, j ) = a( i, j ) end do end do else do j = 1, n do i = 1, m b( i, j ) = a( i, j ) end do end do end if return end subroutine stdlib${ii}$_slacpy pure module subroutine stdlib${ii}$_dlacpy( uplo, m, n, a, lda, b, ldb ) !! DLACPY copies all or part of a two-dimensional matrix A to another !! matrix B. ! -- lapack auxiliary routine -- ! -- lapack is a software package provided by univ. of tennessee, -- ! -- univ. of california berkeley, univ. of colorado denver and nag ltd..-- ! Scalar Arguments character, intent(in) :: uplo integer(${ik}$), intent(in) :: lda, ldb, m, n ! Array Arguments real(dp), intent(in) :: a(lda,*) real(dp), intent(out) :: b(ldb,*) ! ===================================================================== ! Local Scalars integer(${ik}$) :: i, j ! Intrinsic Functions ! Executable Statements if( stdlib_lsame( uplo, 'U' ) ) then do j = 1, n do i = 1, min( j, m ) b( i, j ) = a( i, j ) end do end do else if( stdlib_lsame( uplo, 'L' ) ) then do j = 1, n do i = j, m b( i, j ) = a( i, j ) end do end do else do j = 1, n do i = 1, m b( i, j ) = a( i, j ) end do end do end if return end subroutine stdlib${ii}$_dlacpy #:for rk,rt,ri in REAL_KINDS_TYPES #:if not rk in ["sp","dp"] pure module subroutine stdlib${ii}$_${ri}$lacpy( uplo, m, n, a, lda, b, ldb ) !! DLACPY: copies all or part of a two-dimensional matrix A to another !! matrix B. ! -- lapack auxiliary routine -- ! -- lapack is a software package provided by univ. of tennessee, -- ! -- univ. of california berkeley, univ. of colorado denver and nag ltd..-- ! Scalar Arguments character, intent(in) :: uplo integer(${ik}$), intent(in) :: lda, ldb, m, n ! Array Arguments real(${rk}$), intent(in) :: a(lda,*) real(${rk}$), intent(out) :: b(ldb,*) ! ===================================================================== ! Local Scalars integer(${ik}$) :: i, j ! Intrinsic Functions ! Executable Statements if( stdlib_lsame( uplo, 'U' ) ) then do j = 1, n do i = 1, min( j, m ) b( i, j ) = a( i, j ) end do end do else if( stdlib_lsame( uplo, 'L' ) ) then do j = 1, n do i = j, m b( i, j ) = a( i, j ) end do end do else do j = 1, n do i = 1, m b( i, j ) = a( i, j ) end do end do end if return end subroutine stdlib${ii}$_${ri}$lacpy #:endif #:endfor pure module subroutine stdlib${ii}$_clacpy( uplo, m, n, a, lda, b, ldb ) !! CLACPY copies all or part of a two-dimensional matrix A to another !! matrix B. ! -- lapack auxiliary routine -- ! -- lapack is a software package provided by univ. of tennessee, -- ! -- univ. of california berkeley, univ. of colorado denver and nag ltd..-- ! Scalar Arguments character, intent(in) :: uplo integer(${ik}$), intent(in) :: lda, ldb, m, n ! Array Arguments complex(sp), intent(in) :: a(lda,*) complex(sp), intent(out) :: b(ldb,*) ! ===================================================================== ! Local Scalars integer(${ik}$) :: i, j ! Intrinsic Functions ! Executable Statements if( stdlib_lsame( uplo, 'U' ) ) then do j = 1, n do i = 1, min( j, m ) b( i, j ) = a( i, j ) end do end do else if( stdlib_lsame( uplo, 'L' ) ) then do j = 1, n do i = j, m b( i, j ) = a( i, j ) end do end do else do j = 1, n do i = 1, m b( i, j ) = a( i, j ) end do end do end if return end subroutine stdlib${ii}$_clacpy pure module subroutine stdlib${ii}$_zlacpy( uplo, m, n, a, lda, b, ldb ) !! ZLACPY copies all or part of a two-dimensional matrix A to another !! matrix B. ! -- lapack auxiliary routine -- ! -- lapack is a software package provided by univ. of tennessee, -- ! -- univ. of california berkeley, univ. of colorado denver and nag ltd..-- ! Scalar Arguments character, intent(in) :: uplo integer(${ik}$), intent(in) :: lda, ldb, m, n ! Array Arguments complex(dp), intent(in) :: a(lda,*) complex(dp), intent(out) :: b(ldb,*) ! ===================================================================== ! Local Scalars integer(${ik}$) :: i, j ! Intrinsic Functions ! Executable Statements if( stdlib_lsame( uplo, 'U' ) ) then do j = 1, n do i = 1, min( j, m ) b( i, j ) = a( i, j ) end do end do else if( stdlib_lsame( uplo, 'L' ) ) then do j = 1, n do i = j, m b( i, j ) = a( i, j ) end do end do else do j = 1, n do i = 1, m b( i, j ) = a( i, j ) end do end do end if return end subroutine stdlib${ii}$_zlacpy #:for ck,ct,ci in CMPLX_KINDS_TYPES #:if not ck in ["sp","dp"] pure module subroutine stdlib${ii}$_${ci}$lacpy( uplo, m, n, a, lda, b, ldb ) !! ZLACPY: copies all or part of a two-dimensional matrix A to another !! matrix B. ! -- lapack auxiliary routine -- ! -- lapack is a software package provided by univ. of tennessee, -- ! -- univ. of california berkeley, univ. of colorado denver and nag ltd..-- ! Scalar Arguments character, intent(in) :: uplo integer(${ik}$), intent(in) :: lda, ldb, m, n ! Array Arguments complex(${ck}$), intent(in) :: a(lda,*) complex(${ck}$), intent(out) :: b(ldb,*) ! ===================================================================== ! Local Scalars integer(${ik}$) :: i, j ! Intrinsic Functions ! Executable Statements if( stdlib_lsame( uplo, 'U' ) ) then do j = 1, n do i = 1, min( j, m ) b( i, j ) = a( i, j ) end do end do else if( stdlib_lsame( uplo, 'L' ) ) then do j = 1, n do i = j, m b( i, j ) = a( i, j ) end do end do else do j = 1, n do i = 1, m b( i, j ) = a( i, j ) end do end do end if return end subroutine stdlib${ii}$_${ci}$lacpy #:endif #:endfor pure module subroutine stdlib${ii}$_clacp2( uplo, m, n, a, lda, b, ldb ) !! CLACP2 copies all or part of a real two-dimensional matrix A to a !! complex matrix B. ! -- lapack auxiliary routine -- ! -- lapack is a software package provided by univ. of tennessee, -- ! -- univ. of california berkeley, univ. of colorado denver and nag ltd..-- ! Scalar Arguments character, intent(in) :: uplo integer(${ik}$), intent(in) :: lda, ldb, m, n ! Array Arguments real(sp), intent(in) :: a(lda,*) complex(sp), intent(out) :: b(ldb,*) ! ===================================================================== ! Local Scalars integer(${ik}$) :: i, j ! Intrinsic Functions ! Executable Statements if( stdlib_lsame( uplo, 'U' ) ) then do j = 1, n do i = 1, min( j, m ) b( i, j ) = a( i, j ) end do end do else if( stdlib_lsame( uplo, 'L' ) ) then do j = 1, n do i = j, m b( i, j ) = a( i, j ) end do end do else do j = 1, n do i = 1, m b( i, j ) = a( i, j ) end do end do end if return end subroutine stdlib${ii}$_clacp2 pure module subroutine stdlib${ii}$_zlacp2( uplo, m, n, a, lda, b, ldb ) !! ZLACP2 copies all or part of a real two-dimensional matrix A to a !! complex matrix B. ! -- lapack auxiliary routine -- ! -- lapack is a software package provided by univ. of tennessee, -- ! -- univ. of california berkeley, univ. of colorado denver and nag ltd..-- ! Scalar Arguments character, intent(in) :: uplo integer(${ik}$), intent(in) :: lda, ldb, m, n ! Array Arguments real(dp), intent(in) :: a(lda,*) complex(dp), intent(out) :: b(ldb,*) ! ===================================================================== ! Local Scalars integer(${ik}$) :: i, j ! Intrinsic Functions ! Executable Statements if( stdlib_lsame( uplo, 'U' ) ) then do j = 1, n do i = 1, min( j, m ) b( i, j ) = a( i, j ) end do end do else if( stdlib_lsame( uplo, 'L' ) ) then do j = 1, n do i = j, m b( i, j ) = a( i, j ) end do end do else do j = 1, n do i = 1, m b( i, j ) = a( i, j ) end do end do end if return end subroutine stdlib${ii}$_zlacp2 #:for ck,ct,ci in CMPLX_KINDS_TYPES #:if not ck in ["sp","dp"] pure module subroutine stdlib${ii}$_${ci}$lacp2( uplo, m, n, a, lda, b, ldb ) !! ZLACP2: copies all or part of a real two-dimensional matrix A to a !! complex matrix B. ! -- lapack auxiliary routine -- ! -- lapack is a software package provided by univ. of tennessee, -- ! -- univ. of california berkeley, univ. of colorado denver and nag ltd..-- ! Scalar Arguments character, intent(in) :: uplo integer(${ik}$), intent(in) :: lda, ldb, m, n ! Array Arguments real(${ck}$), intent(in) :: a(lda,*) complex(${ck}$), intent(out) :: b(ldb,*) ! ===================================================================== ! Local Scalars integer(${ik}$) :: i, j ! Intrinsic Functions ! Executable Statements if( stdlib_lsame( uplo, 'U' ) ) then do j = 1, n do i = 1, min( j, m ) b( i, j ) = a( i, j ) end do end do else if( stdlib_lsame( uplo, 'L' ) ) then do j = 1, n do i = j, m b( i, j ) = a( i, j ) end do end do else do j = 1, n do i = 1, m b( i, j ) = a( i, j ) end do end do end if return end subroutine stdlib${ii}$_${ci}$lacp2 #:endif #:endfor pure module subroutine stdlib${ii}$_stfttp( transr, uplo, n, arf, ap, info ) !! STFTTP copies a triangular matrix A from rectangular full packed !! format (TF) to standard packed format (TP). ! -- lapack computational routine -- ! -- lapack is a software package provided by univ. of tennessee, -- ! -- univ. of california berkeley, univ. of colorado denver and nag ltd..-- ! Scalar Arguments character, intent(in) :: transr, uplo integer(${ik}$), intent(out) :: info integer(${ik}$), intent(in) :: n ! Array Arguments real(sp), intent(out) :: ap(0_${ik}$:*) real(sp), intent(in) :: arf(0_${ik}$:*) ! ===================================================================== ! Parameters ! Local Scalars logical(lk) :: lower, nisodd, normaltransr integer(${ik}$) :: n1, n2, k, nt integer(${ik}$) :: i, j, ij integer(${ik}$) :: ijp, jp, lda, js ! Executable Statements ! test the input parameters. info = 0_${ik}$ normaltransr = stdlib_lsame( transr, 'N' ) lower = stdlib_lsame( uplo, 'L' ) if( .not.normaltransr .and. .not.stdlib_lsame( transr, 'T' ) ) then info = -1_${ik}$ else if( .not.lower .and. .not.stdlib_lsame( uplo, 'U' ) ) then info = -2_${ik}$ else if( n<0_${ik}$ ) then info = -3_${ik}$ end if if( info/=0_${ik}$ ) then call stdlib${ii}$_xerbla( 'STFTTP', -info ) return end if ! quick return if possible if( n==0 )return if( n==1_${ik}$ ) then if( normaltransr ) then ap( 0_${ik}$ ) = arf( 0_${ik}$ ) else ap( 0_${ik}$ ) = arf( 0_${ik}$ ) end if return end if ! size of array arf(0:nt-1) nt = n*( n+1 ) / 2_${ik}$ ! set n1 and n2 depending on lower if( lower ) then n2 = n / 2_${ik}$ n1 = n - n2 else n1 = n / 2_${ik}$ n2 = n - n1 end if ! if n is odd, set nisodd = .true. ! if n is even, set k = n/2 and nisodd = .false. ! set lda of arf^c; arf^c is (0:(n+1)/2-1,0:n-noe) ! where noe = 0 if n is even, noe = 1 if n is odd if( mod( n, 2_${ik}$ )==0_${ik}$ ) then k = n / 2_${ik}$ nisodd = .false. lda = n + 1_${ik}$ else nisodd = .true. lda = n end if ! arf^c has lda rows and n+1-noe cols if( .not.normaltransr )lda = ( n+1 ) / 2_${ik}$ ! start execution: there are eight cases if( nisodd ) then ! n is odd if( normaltransr ) then ! n is odd and transr = 'n' if( lower ) then ! srpa for lower, normal and n is odd ( a(0:n-1,0:n1-1) ) ! t1 -> a(0,0), t2 -> a(0,1), s -> a(n1,0) ! t1 -> a(0), t2 -> a(n), s -> a(n1); lda = n ijp = 0_${ik}$ jp = 0_${ik}$ do j = 0, n2 do i = j, n - 1 ij = i + jp ap( ijp ) = arf( ij ) ijp = ijp + 1_${ik}$ end do jp = jp + lda end do do i = 0, n2 - 1 do j = 1 + i, n2 ij = i + j*lda ap( ijp ) = arf( ij ) ijp = ijp + 1_${ik}$ end do end do else ! srpa for upper, normal and n is odd ( a(0:n-1,0:n2-1) ! t1 -> a(n1+1,0), t2 -> a(n1,0), s -> a(0,0) ! t1 -> a(n2), t2 -> a(n1), s -> a(0) ijp = 0_${ik}$ do j = 0, n1 - 1 ij = n2 + j do i = 0, j ap( ijp ) = arf( ij ) ijp = ijp + 1_${ik}$ ij = ij + lda end do end do js = 0_${ik}$ do j = n1, n - 1 ij = js do ij = js, js + j ap( ijp ) = arf( ij ) ijp = ijp + 1_${ik}$ end do js = js + lda end do end if else ! n is odd and transr = 't' if( lower ) then ! srpa for lower, transpose and n is odd ! t1 -> a(0,0) , t2 -> a(1,0) , s -> a(0,n1) ! t1 -> a(0+0) , t2 -> a(1+0) , s -> a(0+n1*n1); lda=n1 ijp = 0_${ik}$ do i = 0, n2 do ij = i*( lda+1 ), n*lda - 1, lda ap( ijp ) = arf( ij ) ijp = ijp + 1_${ik}$ end do end do js = 1_${ik}$ do j = 0, n2 - 1 do ij = js, js + n2 - j - 1 ap( ijp ) = arf( ij ) ijp = ijp + 1_${ik}$ end do js = js + lda + 1_${ik}$ end do else ! srpa for upper, transpose and n is odd ! t1 -> a(0,n1+1), t2 -> a(0,n1), s -> a(0,0) ! t1 -> a(n2*n2), t2 -> a(n1*n2), s -> a(0); lda = n2 ijp = 0_${ik}$ js = n2*lda do j = 0, n1 - 1 do ij = js, js + j ap( ijp ) = arf( ij ) ijp = ijp + 1_${ik}$ end do js = js + lda end do do i = 0, n1 do ij = i, i + ( n1+i )*lda, lda ap( ijp ) = arf( ij ) ijp = ijp + 1_${ik}$ end do end do end if end if else ! n is even if( normaltransr ) then ! n is even and transr = 'n' if( lower ) then ! srpa for lower, normal, and n is even ( a(0:n,0:k-1) ) ! t1 -> a(1,0), t2 -> a(0,0), s -> a(k+1,0) ! t1 -> a(1), t2 -> a(0), s -> a(k+1) ijp = 0_${ik}$ jp = 0_${ik}$ do j = 0, k - 1 do i = j, n - 1 ij = 1_${ik}$ + i + jp ap( ijp ) = arf( ij ) ijp = ijp + 1_${ik}$ end do jp = jp + lda end do do i = 0, k - 1 do j = i, k - 1 ij = i + j*lda ap( ijp ) = arf( ij ) ijp = ijp + 1_${ik}$ end do end do else ! srpa for upper, normal, and n is even ( a(0:n,0:k-1) ) ! t1 -> a(k+1,0) , t2 -> a(k,0), s -> a(0,0) ! t1 -> a(k+1), t2 -> a(k), s -> a(0) ijp = 0_${ik}$ do j = 0, k - 1 ij = k + 1_${ik}$ + j do i = 0, j ap( ijp ) = arf( ij ) ijp = ijp + 1_${ik}$ ij = ij + lda end do end do js = 0_${ik}$ do j = k, n - 1 ij = js do ij = js, js + j ap( ijp ) = arf( ij ) ijp = ijp + 1_${ik}$ end do js = js + lda end do end if else ! n is even and transr = 't' if( lower ) then ! srpa for lower, transpose and n is even (see paper) ! t1 -> b(0,1), t2 -> b(0,0), s -> b(0,k+1) ! t1 -> a(0+k), t2 -> a(0+0), s -> a(0+k*(k+1)); lda=k ijp = 0_${ik}$ do i = 0, k - 1 do ij = i + ( i+1 )*lda, ( n+1 )*lda - 1, lda ap( ijp ) = arf( ij ) ijp = ijp + 1_${ik}$ end do end do js = 0_${ik}$ do j = 0, k - 1 do ij = js, js + k - j - 1 ap( ijp ) = arf( ij ) ijp = ijp + 1_${ik}$ end do js = js + lda + 1_${ik}$ end do else ! srpa for upper, transpose and n is even (see paper) ! t1 -> b(0,k+1), t2 -> b(0,k), s -> b(0,0) ! t1 -> a(0+k*(k+1)), t2 -> a(0+k*k), s -> a(0+0)); lda=k ijp = 0_${ik}$ js = ( k+1 )*lda do j = 0, k - 1 do ij = js, js + j ap( ijp ) = arf( ij ) ijp = ijp + 1_${ik}$ end do js = js + lda end do do i = 0, k - 1 do ij = i, i + ( k+i )*lda, lda ap( ijp ) = arf( ij ) ijp = ijp + 1_${ik}$ end do end do end if end if end if return end subroutine stdlib${ii}$_stfttp pure module subroutine stdlib${ii}$_dtfttp( transr, uplo, n, arf, ap, info ) !! DTFTTP copies a triangular matrix A from rectangular full packed !! format (TF) to standard packed format (TP). ! -- lapack computational routine -- ! -- lapack is a software package provided by univ. of tennessee, -- ! -- univ. of california berkeley, univ. of colorado denver and nag ltd..-- ! Scalar Arguments character, intent(in) :: transr, uplo integer(${ik}$), intent(out) :: info integer(${ik}$), intent(in) :: n ! Array Arguments real(dp), intent(out) :: ap(0_${ik}$:*) real(dp), intent(in) :: arf(0_${ik}$:*) ! ===================================================================== ! Parameters ! Local Scalars logical(lk) :: lower, nisodd, normaltransr integer(${ik}$) :: n1, n2, k, nt integer(${ik}$) :: i, j, ij integer(${ik}$) :: ijp, jp, lda, js ! Executable Statements ! test the input parameters. info = 0_${ik}$ normaltransr = stdlib_lsame( transr, 'N' ) lower = stdlib_lsame( uplo, 'L' ) if( .not.normaltransr .and. .not.stdlib_lsame( transr, 'T' ) ) then info = -1_${ik}$ else if( .not.lower .and. .not.stdlib_lsame( uplo, 'U' ) ) then info = -2_${ik}$ else if( n<0_${ik}$ ) then info = -3_${ik}$ end if if( info/=0_${ik}$ ) then call stdlib${ii}$_xerbla( 'DTFTTP', -info ) return end if ! quick return if possible if( n==0 )return if( n==1_${ik}$ ) then if( normaltransr ) then ap( 0_${ik}$ ) = arf( 0_${ik}$ ) else ap( 0_${ik}$ ) = arf( 0_${ik}$ ) end if return end if ! size of array arf(0:nt-1) nt = n*( n+1 ) / 2_${ik}$ ! set n1 and n2 depending on lower if( lower ) then n2 = n / 2_${ik}$ n1 = n - n2 else n1 = n / 2_${ik}$ n2 = n - n1 end if ! if n is odd, set nisodd = .true. ! if n is even, set k = n/2 and nisodd = .false. ! set lda of arf^c; arf^c is (0:(n+1)/2-1,0:n-noe) ! where noe = 0 if n is even, noe = 1 if n is odd if( mod( n, 2_${ik}$ )==0_${ik}$ ) then k = n / 2_${ik}$ nisodd = .false. lda = n + 1_${ik}$ else nisodd = .true. lda = n end if ! arf^c has lda rows and n+1-noe cols if( .not.normaltransr )lda = ( n+1 ) / 2_${ik}$ ! start execution: there are eight cases if( nisodd ) then ! n is odd if( normaltransr ) then ! n is odd and transr = 'n' if( lower ) then ! srpa for lower, normal and n is odd ( a(0:n-1,0:n1-1) ) ! t1 -> a(0,0), t2 -> a(0,1), s -> a(n1,0) ! t1 -> a(0), t2 -> a(n), s -> a(n1); lda = n ijp = 0_${ik}$ jp = 0_${ik}$ do j = 0, n2 do i = j, n - 1 ij = i + jp ap( ijp ) = arf( ij ) ijp = ijp + 1_${ik}$ end do jp = jp + lda end do do i = 0, n2 - 1 do j = 1 + i, n2 ij = i + j*lda ap( ijp ) = arf( ij ) ijp = ijp + 1_${ik}$ end do end do else ! srpa for upper, normal and n is odd ( a(0:n-1,0:n2-1) ! t1 -> a(n1+1,0), t2 -> a(n1,0), s -> a(0,0) ! t1 -> a(n2), t2 -> a(n1), s -> a(0) ijp = 0_${ik}$ do j = 0, n1 - 1 ij = n2 + j do i = 0, j ap( ijp ) = arf( ij ) ijp = ijp + 1_${ik}$ ij = ij + lda end do end do js = 0_${ik}$ do j = n1, n - 1 ij = js do ij = js, js + j ap( ijp ) = arf( ij ) ijp = ijp + 1_${ik}$ end do js = js + lda end do end if else ! n is odd and transr = 't' if( lower ) then ! srpa for lower, transpose and n is odd ! t1 -> a(0,0) , t2 -> a(1,0) , s -> a(0,n1) ! t1 -> a(0+0) , t2 -> a(1+0) , s -> a(0+n1*n1); lda=n1 ijp = 0_${ik}$ do i = 0, n2 do ij = i*( lda+1 ), n*lda - 1, lda ap( ijp ) = arf( ij ) ijp = ijp + 1_${ik}$ end do end do js = 1_${ik}$ do j = 0, n2 - 1 do ij = js, js + n2 - j - 1 ap( ijp ) = arf( ij ) ijp = ijp + 1_${ik}$ end do js = js + lda + 1_${ik}$ end do else ! srpa for upper, transpose and n is odd ! t1 -> a(0,n1+1), t2 -> a(0,n1), s -> a(0,0) ! t1 -> a(n2*n2), t2 -> a(n1*n2), s -> a(0); lda = n2 ijp = 0_${ik}$ js = n2*lda do j = 0, n1 - 1 do ij = js, js + j ap( ijp ) = arf( ij ) ijp = ijp + 1_${ik}$ end do js = js + lda end do do i = 0, n1 do ij = i, i + ( n1+i )*lda, lda ap( ijp ) = arf( ij ) ijp = ijp + 1_${ik}$ end do end do end if end if else ! n is even if( normaltransr ) then ! n is even and transr = 'n' if( lower ) then ! srpa for lower, normal, and n is even ( a(0:n,0:k-1) ) ! t1 -> a(1,0), t2 -> a(0,0), s -> a(k+1,0) ! t1 -> a(1), t2 -> a(0), s -> a(k+1) ijp = 0_${ik}$ jp = 0_${ik}$ do j = 0, k - 1 do i = j, n - 1 ij = 1_${ik}$ + i + jp ap( ijp ) = arf( ij ) ijp = ijp + 1_${ik}$ end do jp = jp + lda end do do i = 0, k - 1 do j = i, k - 1 ij = i + j*lda ap( ijp ) = arf( ij ) ijp = ijp + 1_${ik}$ end do end do else ! srpa for upper, normal, and n is even ( a(0:n,0:k-1) ) ! t1 -> a(k+1,0) , t2 -> a(k,0), s -> a(0,0) ! t1 -> a(k+1), t2 -> a(k), s -> a(0) ijp = 0_${ik}$ do j = 0, k - 1 ij = k + 1_${ik}$ + j do i = 0, j ap( ijp ) = arf( ij ) ijp = ijp + 1_${ik}$ ij = ij + lda end do end do js = 0_${ik}$ do j = k, n - 1 ij = js do ij = js, js + j ap( ijp ) = arf( ij ) ijp = ijp + 1_${ik}$ end do js = js + lda end do end if else ! n is even and transr = 't' if( lower ) then ! srpa for lower, transpose and n is even (see paper) ! t1 -> b(0,1), t2 -> b(0,0), s -> b(0,k+1) ! t1 -> a(0+k), t2 -> a(0+0), s -> a(0+k*(k+1)); lda=k ijp = 0_${ik}$ do i = 0, k - 1 do ij = i + ( i+1 )*lda, ( n+1 )*lda - 1, lda ap( ijp ) = arf( ij ) ijp = ijp + 1_${ik}$ end do end do js = 0_${ik}$ do j = 0, k - 1 do ij = js, js + k - j - 1 ap( ijp ) = arf( ij ) ijp = ijp + 1_${ik}$ end do js = js + lda + 1_${ik}$ end do else ! srpa for upper, transpose and n is even (see paper) ! t1 -> b(0,k+1), t2 -> b(0,k), s -> b(0,0) ! t1 -> a(0+k*(k+1)), t2 -> a(0+k*k), s -> a(0+0)); lda=k ijp = 0_${ik}$ js = ( k+1 )*lda do j = 0, k - 1 do ij = js, js + j ap( ijp ) = arf( ij ) ijp = ijp + 1_${ik}$ end do js = js + lda end do do i = 0, k - 1 do ij = i, i + ( k+i )*lda, lda ap( ijp ) = arf( ij ) ijp = ijp + 1_${ik}$ end do end do end if end if end if return end subroutine stdlib${ii}$_dtfttp #:for rk,rt,ri in REAL_KINDS_TYPES #:if not rk in ["sp","dp"] pure module subroutine stdlib${ii}$_${ri}$tfttp( transr, uplo, n, arf, ap, info ) !! DTFTTP: copies a triangular matrix A from rectangular full packed !! format (TF) to standard packed format (TP). ! -- lapack computational routine -- ! -- lapack is a software package provided by univ. of tennessee, -- ! -- univ. of california berkeley, univ. of colorado denver and nag ltd..-- ! Scalar Arguments character, intent(in) :: transr, uplo integer(${ik}$), intent(out) :: info integer(${ik}$), intent(in) :: n ! Array Arguments real(${rk}$), intent(out) :: ap(0_${ik}$:*) real(${rk}$), intent(in) :: arf(0_${ik}$:*) ! ===================================================================== ! Parameters ! Local Scalars logical(lk) :: lower, nisodd, normaltransr integer(${ik}$) :: n1, n2, k, nt integer(${ik}$) :: i, j, ij integer(${ik}$) :: ijp, jp, lda, js ! Executable Statements ! test the input parameters. info = 0_${ik}$ normaltransr = stdlib_lsame( transr, 'N' ) lower = stdlib_lsame( uplo, 'L' ) if( .not.normaltransr .and. .not.stdlib_lsame( transr, 'T' ) ) then info = -1_${ik}$ else if( .not.lower .and. .not.stdlib_lsame( uplo, 'U' ) ) then info = -2_${ik}$ else if( n<0_${ik}$ ) then info = -3_${ik}$ end if if( info/=0_${ik}$ ) then call stdlib${ii}$_xerbla( 'DTFTTP', -info ) return end if ! quick return if possible if( n==0 )return if( n==1_${ik}$ ) then if( normaltransr ) then ap( 0_${ik}$ ) = arf( 0_${ik}$ ) else ap( 0_${ik}$ ) = arf( 0_${ik}$ ) end if return end if ! size of array arf(0:nt-1) nt = n*( n+1 ) / 2_${ik}$ ! set n1 and n2 depending on lower if( lower ) then n2 = n / 2_${ik}$ n1 = n - n2 else n1 = n / 2_${ik}$ n2 = n - n1 end if ! if n is odd, set nisodd = .true. ! if n is even, set k = n/2 and nisodd = .false. ! set lda of arf^c; arf^c is (0:(n+1)/2-1,0:n-noe) ! where noe = 0 if n is even, noe = 1 if n is odd if( mod( n, 2_${ik}$ )==0_${ik}$ ) then k = n / 2_${ik}$ nisodd = .false. lda = n + 1_${ik}$ else nisodd = .true. lda = n end if ! arf^c has lda rows and n+1-noe cols if( .not.normaltransr )lda = ( n+1 ) / 2_${ik}$ ! start execution: there are eight cases if( nisodd ) then ! n is odd if( normaltransr ) then ! n is odd and transr = 'n' if( lower ) then ! srpa for lower, normal and n is odd ( a(0:n-1,0:n1-1) ) ! t1 -> a(0,0), t2 -> a(0,1), s -> a(n1,0) ! t1 -> a(0), t2 -> a(n), s -> a(n1); lda = n ijp = 0_${ik}$ jp = 0_${ik}$ do j = 0, n2 do i = j, n - 1 ij = i + jp ap( ijp ) = arf( ij ) ijp = ijp + 1_${ik}$ end do jp = jp + lda end do do i = 0, n2 - 1 do j = 1 + i, n2 ij = i + j*lda ap( ijp ) = arf( ij ) ijp = ijp + 1_${ik}$ end do end do else ! srpa for upper, normal and n is odd ( a(0:n-1,0:n2-1) ! t1 -> a(n1+1,0), t2 -> a(n1,0), s -> a(0,0) ! t1 -> a(n2), t2 -> a(n1), s -> a(0) ijp = 0_${ik}$ do j = 0, n1 - 1 ij = n2 + j do i = 0, j ap( ijp ) = arf( ij ) ijp = ijp + 1_${ik}$ ij = ij + lda end do end do js = 0_${ik}$ do j = n1, n - 1 ij = js do ij = js, js + j ap( ijp ) = arf( ij ) ijp = ijp + 1_${ik}$ end do js = js + lda end do end if else ! n is odd and transr = 't' if( lower ) then ! srpa for lower, transpose and n is odd ! t1 -> a(0,0) , t2 -> a(1,0) , s -> a(0,n1) ! t1 -> a(0+0) , t2 -> a(1+0) , s -> a(0+n1*n1); lda=n1 ijp = 0_${ik}$ do i = 0, n2 do ij = i*( lda+1 ), n*lda - 1, lda ap( ijp ) = arf( ij ) ijp = ijp + 1_${ik}$ end do end do js = 1_${ik}$ do j = 0, n2 - 1 do ij = js, js + n2 - j - 1 ap( ijp ) = arf( ij ) ijp = ijp + 1_${ik}$ end do js = js + lda + 1_${ik}$ end do else ! srpa for upper, transpose and n is odd ! t1 -> a(0,n1+1), t2 -> a(0,n1), s -> a(0,0) ! t1 -> a(n2*n2), t2 -> a(n1*n2), s -> a(0); lda = n2 ijp = 0_${ik}$ js = n2*lda do j = 0, n1 - 1 do ij = js, js + j ap( ijp ) = arf( ij ) ijp = ijp + 1_${ik}$ end do js = js + lda end do do i = 0, n1 do ij = i, i + ( n1+i )*lda, lda ap( ijp ) = arf( ij ) ijp = ijp + 1_${ik}$ end do end do end if end if else ! n is even if( normaltransr ) then ! n is even and transr = 'n' if( lower ) then ! srpa for lower, normal, and n is even ( a(0:n,0:k-1) ) ! t1 -> a(1,0), t2 -> a(0,0), s -> a(k+1,0) ! t1 -> a(1), t2 -> a(0), s -> a(k+1) ijp = 0_${ik}$ jp = 0_${ik}$ do j = 0, k - 1 do i = j, n - 1 ij = 1_${ik}$ + i + jp ap( ijp ) = arf( ij ) ijp = ijp + 1_${ik}$ end do jp = jp + lda end do do i = 0, k - 1 do j = i, k - 1 ij = i + j*lda ap( ijp ) = arf( ij ) ijp = ijp + 1_${ik}$ end do end do else ! srpa for upper, normal, and n is even ( a(0:n,0:k-1) ) ! t1 -> a(k+1,0) , t2 -> a(k,0), s -> a(0,0) ! t1 -> a(k+1), t2 -> a(k), s -> a(0) ijp = 0_${ik}$ do j = 0, k - 1 ij = k + 1_${ik}$ + j do i = 0, j ap( ijp ) = arf( ij ) ijp = ijp + 1_${ik}$ ij = ij + lda end do end do js = 0_${ik}$ do j = k, n - 1 ij = js do ij = js, js + j ap( ijp ) = arf( ij ) ijp = ijp + 1_${ik}$ end do js = js + lda end do end if else ! n is even and transr = 't' if( lower ) then ! srpa for lower, transpose and n is even (see paper) ! t1 -> b(0,1), t2 -> b(0,0), s -> b(0,k+1) ! t1 -> a(0+k), t2 -> a(0+0), s -> a(0+k*(k+1)); lda=k ijp = 0_${ik}$ do i = 0, k - 1 do ij = i + ( i+1 )*lda, ( n+1 )*lda - 1, lda ap( ijp ) = arf( ij ) ijp = ijp + 1_${ik}$ end do end do js = 0_${ik}$ do j = 0, k - 1 do ij = js, js + k - j - 1 ap( ijp ) = arf( ij ) ijp = ijp + 1_${ik}$ end do js = js + lda + 1_${ik}$ end do else ! srpa for upper, transpose and n is even (see paper) ! t1 -> b(0,k+1), t2 -> b(0,k), s -> b(0,0) ! t1 -> a(0+k*(k+1)), t2 -> a(0+k*k), s -> a(0+0)); lda=k ijp = 0_${ik}$ js = ( k+1 )*lda do j = 0, k - 1 do ij = js, js + j ap( ijp ) = arf( ij ) ijp = ijp + 1_${ik}$ end do js = js + lda end do do i = 0, k - 1 do ij = i, i + ( k+i )*lda, lda ap( ijp ) = arf( ij ) ijp = ijp + 1_${ik}$ end do end do end if end if end if return end subroutine stdlib${ii}$_${ri}$tfttp #:endif #:endfor pure module subroutine stdlib${ii}$_ctfttp( transr, uplo, n, arf, ap, info ) !! CTFTTP copies a triangular matrix A from rectangular full packed !! format (TF) to standard packed format (TP). ! -- lapack computational routine -- ! -- lapack is a software package provided by univ. of tennessee, -- ! -- univ. of california berkeley, univ. of colorado denver and nag ltd..-- ! Scalar Arguments character, intent(in) :: transr, uplo integer(${ik}$), intent(out) :: info integer(${ik}$), intent(in) :: n ! Array Arguments complex(sp), intent(out) :: ap(0_${ik}$:*) complex(sp), intent(in) :: arf(0_${ik}$:*) ! ===================================================================== ! Parameters ! Local Scalars logical(lk) :: lower, nisodd, normaltransr integer(${ik}$) :: n1, n2, k, nt integer(${ik}$) :: i, j, ij integer(${ik}$) :: ijp, jp, lda, js ! Intrinsic Functions ! Intrinsic Functions ! Executable Statements ! test the input parameters. info = 0_${ik}$ normaltransr = stdlib_lsame( transr, 'N' ) lower = stdlib_lsame( uplo, 'L' ) if( .not.normaltransr .and. .not.stdlib_lsame( transr, 'C' ) ) then info = -1_${ik}$ else if( .not.lower .and. .not.stdlib_lsame( uplo, 'U' ) ) then info = -2_${ik}$ else if( n<0_${ik}$ ) then info = -3_${ik}$ end if if( info/=0_${ik}$ ) then call stdlib${ii}$_xerbla( 'CTFTTP', -info ) return end if ! quick return if possible if( n==0 )return if( n==1_${ik}$ ) then if( normaltransr ) then ap( 0_${ik}$ ) = arf( 0_${ik}$ ) else ap( 0_${ik}$ ) = conjg( arf( 0_${ik}$ ) ) end if return end if ! size of array arf(0:nt-1) nt = n*( n+1 ) / 2_${ik}$ ! set n1 and n2 depending on lower if( lower ) then n2 = n / 2_${ik}$ n1 = n - n2 else n1 = n / 2_${ik}$ n2 = n - n1 end if ! if n is odd, set nisodd = .true. ! if n is even, set k = n/2 and nisodd = .false. ! set lda of arf^c; arf^c is (0:(n+1)/2-1,0:n-noe)