#:include "common.fypp" submodule(stdlib_lapack_solve) stdlib_lapack_solve_lu_comp implicit none contains #:for ik,it,ii in LINALG_INT_KINDS_TYPES pure module subroutine stdlib${ii}$_sgecon( norm, n, a, lda, anorm, rcond, work, iwork,info ) !! SGECON estimates the reciprocal of the condition number of a general !! real matrix A, in either the 1-norm or the infinity-norm, using !! the LU factorization computed by SGETRF. !! An estimate is obtained for norm(inv(A)), and the reciprocal of the !! condition number is computed as !! RCOND = 1 / ( norm(A) * norm(inv(A)) ). ! -- lapack computational 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: negone, zero, half, one, two, three, four, eight, ten, czero, chalf, cone, cnegone ! Scalar Arguments character, intent(in) :: norm integer(${ik}$), intent(out) :: info integer(${ik}$), intent(in) :: lda, n real(sp), intent(in) :: anorm real(sp), intent(out) :: rcond ! Array Arguments integer(${ik}$), intent(out) :: iwork(*) real(sp), intent(inout) :: a(lda,*) real(sp), intent(out) :: work(*) ! ===================================================================== ! Local Scalars logical(lk) :: onenrm character :: normin integer(${ik}$) :: ix, kase, kase1 real(sp) :: ainvnm, scale, sl, smlnum, su ! Local Arrays integer(${ik}$) :: isave(3_${ik}$) ! Intrinsic Functions ! Executable Statements ! test the input parameters. info = 0_${ik}$ onenrm = norm=='1' .or. stdlib_lsame( norm, 'O' ) if( .not.onenrm .and. .not.stdlib_lsame( norm, 'I' ) ) then info = -1_${ik}$ else if( n<0_${ik}$ ) then info = -2_${ik}$ else if( lda<max( 1_${ik}$, n ) ) then info = -4_${ik}$ else if( anorm<zero ) then info = -5_${ik}$ end if if( info/=0_${ik}$ ) then call stdlib${ii}$_xerbla( 'SGECON', -info ) return end if ! quick return if possible rcond = zero if( n==0_${ik}$ ) then rcond = one return else if( anorm==zero ) then return end if smlnum = stdlib${ii}$_slamch( 'SAFE MINIMUM' ) ! estimate the norm of inv(a). ainvnm = zero normin = 'N' if( onenrm ) then kase1 = 1_${ik}$ else kase1 = 2_${ik}$ end if kase = 0_${ik}$ 10 continue call stdlib${ii}$_slacn2( n, work( n+1 ), work, iwork, ainvnm, kase, isave ) if( kase/=0_${ik}$ ) then if( kase==kase1 ) then ! multiply by inv(l). call stdlib${ii}$_slatrs( 'LOWER', 'NO TRANSPOSE', 'UNIT', normin, n, a,lda, work, sl, & work( 2_${ik}$*n+1 ), info ) ! multiply by inv(u). call stdlib${ii}$_slatrs( 'UPPER', 'NO TRANSPOSE', 'NON-UNIT', normin, n,a, lda, work, & su, work( 3_${ik}$*n+1 ), info ) else ! multiply by inv(u**t). call stdlib${ii}$_slatrs( 'UPPER', 'TRANSPOSE', 'NON-UNIT', normin, n, a,lda, work, su,& work( 3_${ik}$*n+1 ), info ) ! multiply by inv(l**t). call stdlib${ii}$_slatrs( 'LOWER', 'TRANSPOSE', 'UNIT', normin, n, a,lda, work, sl, & work( 2_${ik}$*n+1 ), info ) end if ! divide x by 1/(sl*su) if doing so will not cause overflow. scale = sl*su normin = 'Y' if( scale/=one ) then ix = stdlib${ii}$_isamax( n, work, 1_${ik}$ ) if( scale<abs( work( ix ) )*smlnum .or. scale==zero )go to 20 call stdlib${ii}$_srscl( n, scale, work, 1_${ik}$ ) end if go to 10 end if ! compute the estimate of the reciprocal condition number. if( ainvnm/=zero )rcond = ( one / ainvnm ) / anorm 20 continue return end subroutine stdlib${ii}$_sgecon pure module subroutine stdlib${ii}$_dgecon( norm, n, a, lda, anorm, rcond, work, iwork,info ) !! DGECON estimates the reciprocal of the condition number of a general !! real matrix A, in either the 1-norm or the infinity-norm, using !! the LU factorization computed by DGETRF. !! An estimate is obtained for norm(inv(A)), and the reciprocal of the !! condition number is computed as !! RCOND = 1 / ( norm(A) * norm(inv(A)) ). ! -- lapack computational 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: negone, zero, half, one, two, three, four, eight, ten, czero, chalf, cone, cnegone ! Scalar Arguments character, intent(in) :: norm integer(${ik}$), intent(out) :: info integer(${ik}$), intent(in) :: lda, n real(dp), intent(in) :: anorm real(dp), intent(out) :: rcond ! Array Arguments integer(${ik}$), intent(out) :: iwork(*) real(dp), intent(inout) :: a(lda,*) real(dp), intent(out) :: work(*) ! ===================================================================== ! Local Scalars logical(lk) :: onenrm character :: normin integer(${ik}$) :: ix, kase, kase1 real(dp) :: ainvnm, scale, sl, smlnum, su ! Local Arrays integer(${ik}$) :: isave(3_${ik}$) ! Intrinsic Functions ! Executable Statements ! test the input parameters. info = 0_${ik}$ onenrm = norm=='1' .or. stdlib_lsame( norm, 'O' ) if( .not.onenrm .and. .not.stdlib_lsame( norm, 'I' ) ) then info = -1_${ik}$ else if( n<0_${ik}$ ) then info = -2_${ik}$ else if( lda<max( 1_${ik}$, n ) ) then info = -4_${ik}$ else if( anorm<zero ) then info = -5_${ik}$ end if if( info/=0_${ik}$ ) then call stdlib${ii}$_xerbla( 'DGECON', -info ) return end if ! quick return if possible rcond = zero if( n==0_${ik}$ ) then rcond = one return else if( anorm==zero ) then return end if smlnum = stdlib${ii}$_dlamch( 'SAFE MINIMUM' ) ! estimate the norm of inv(a). ainvnm = zero normin = 'N' if( onenrm ) then kase1 = 1_${ik}$ else kase1 = 2_${ik}$ end if kase = 0_${ik}$ 10 continue call stdlib${ii}$_dlacn2( n, work( n+1 ), work, iwork, ainvnm, kase, isave ) if( kase/=0_${ik}$ ) then if( kase==kase1 ) then ! multiply by inv(l). call stdlib${ii}$_dlatrs( 'LOWER', 'NO TRANSPOSE', 'UNIT', normin, n, a,lda, work, sl, & work( 2_${ik}$*n+1 ), info ) ! multiply by inv(u). call stdlib${ii}$_dlatrs( 'UPPER', 'NO TRANSPOSE', 'NON-UNIT', normin, n,a, lda, work, & su, work( 3_${ik}$*n+1 ), info ) else ! multiply by inv(u**t). call stdlib${ii}$_dlatrs( 'UPPER', 'TRANSPOSE', 'NON-UNIT', normin, n, a,lda, work, su,& work( 3_${ik}$*n+1 ), info ) ! multiply by inv(l**t). call stdlib${ii}$_dlatrs( 'LOWER', 'TRANSPOSE', 'UNIT', normin, n, a,lda, work, sl, & work( 2_${ik}$*n+1 ), info ) end if ! divide x by 1/(sl*su) if doing so will not cause overflow. scale = sl*su normin = 'Y' if( scale/=one ) then ix = stdlib${ii}$_idamax( n, work, 1_${ik}$ ) if( scale<abs( work( ix ) )*smlnum .or. scale==zero )go to 20 call stdlib${ii}$_drscl( n, scale, work, 1_${ik}$ ) end if go to 10 end if ! compute the estimate of the reciprocal condition number. if( ainvnm/=zero )rcond = ( one / ainvnm ) / anorm 20 continue return end subroutine stdlib${ii}$_dgecon #:for rk,rt,ri in REAL_KINDS_TYPES #:if not rk in ["sp","dp"] pure module subroutine stdlib${ii}$_${ri}$gecon( norm, n, a, lda, anorm, rcond, work, iwork,info ) !! DGECON: estimates the reciprocal of the condition number of a general !! real matrix A, in either the 1-norm or the infinity-norm, using !! the LU factorization computed by DGETRF. !! An estimate is obtained for norm(inv(A)), and the reciprocal of the !! condition number is computed as !! RCOND = 1 / ( norm(A) * norm(inv(A)) ). ! -- lapack computational 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: negone, zero, half, one, two, three, four, eight, ten, czero, chalf, cone, cnegone ! Scalar Arguments character, intent(in) :: norm integer(${ik}$), intent(out) :: info integer(${ik}$), intent(in) :: lda, n real(${rk}$), intent(in) :: anorm real(${rk}$), intent(out) :: rcond ! Array Arguments integer(${ik}$), intent(out) :: iwork(*) real(${rk}$), intent(inout) :: a(lda,*) real(${rk}$), intent(out) :: work(*) ! ===================================================================== ! Local Scalars logical(lk) :: onenrm character :: normin integer(${ik}$) :: ix, kase, kase1 real(${rk}$) :: ainvnm, scale, sl, smlnum, su ! Local Arrays integer(${ik}$) :: isave(3_${ik}$) ! Intrinsic Functions ! Executable Statements ! test the input parameters. info = 0_${ik}$ onenrm = norm=='1' .or. stdlib_lsame( norm, 'O' ) if( .not.onenrm .and. .not.stdlib_lsame( norm, 'I' ) ) then info = -1_${ik}$ else if( n<0_${ik}$ ) then info = -2_${ik}$ else if( lda<max( 1_${ik}$, n ) ) then info = -4_${ik}$ else if( anorm<zero ) then info = -5_${ik}$ end if if( info/=0_${ik}$ ) then call stdlib${ii}$_xerbla( 'DGECON', -info ) return end if ! quick return if possible rcond = zero if( n==0_${ik}$ ) then rcond = one return else if( anorm==zero ) then return end if smlnum = stdlib${ii}$_${ri}$lamch( 'SAFE MINIMUM' ) ! estimate the norm of inv(a). ainvnm = zero normin = 'N' if( onenrm ) then kase1 = 1_${ik}$ else kase1 = 2_${ik}$ end if kase = 0_${ik}$ 10 continue call stdlib${ii}$_${ri}$lacn2( n, work( n+1 ), work, iwork, ainvnm, kase, isave ) if( kase/=0_${ik}$ ) then if( kase==kase1 ) then ! multiply by inv(l). call stdlib${ii}$_${ri}$latrs( 'LOWER', 'NO TRANSPOSE', 'UNIT', normin, n, a,lda, work, sl, & work( 2_${ik}$*n+1 ), info ) ! multiply by inv(u). call stdlib${ii}$_${ri}$latrs( 'UPPER', 'NO TRANSPOSE', 'NON-UNIT', normin, n,a, lda, work, & su, work( 3_${ik}$*n+1 ), info ) else ! multiply by inv(u**t). call stdlib${ii}$_${ri}$latrs( 'UPPER', 'TRANSPOSE', 'NON-UNIT', normin, n, a,lda, work, su,& work( 3_${ik}$*n+1 ), info ) ! multiply by inv(l**t). call stdlib${ii}$_${ri}$latrs( 'LOWER', 'TRANSPOSE', 'UNIT', normin, n, a,lda, work, sl, & work( 2_${ik}$*n+1 ), info ) end if ! divide x by 1/(sl*su) if doing so will not cause overflow. scale = sl*su normin = 'Y' if( scale/=one ) then ix = stdlib${ii}$_i${ri}$amax( n, work, 1_${ik}$ ) if( scale<abs( work( ix ) )*smlnum .or. scale==zero )go to 20 call stdlib${ii}$_${ri}$rscl( n, scale, work, 1_${ik}$ ) end if go to 10 end if ! compute the estimate of the reciprocal condition number. if( ainvnm/=zero )rcond = ( one / ainvnm ) / anorm 20 continue return end subroutine stdlib${ii}$_${ri}$gecon #:endif #:endfor pure module subroutine stdlib${ii}$_cgecon( norm, n, a, lda, anorm, rcond, work, rwork,info ) !! CGECON estimates the reciprocal of the condition number of a general !! complex matrix A, in either the 1-norm or the infinity-norm, using !! the LU factorization computed by CGETRF. !! An estimate is obtained for norm(inv(A)), and the reciprocal of the !! condition number is computed as !! RCOND = 1 / ( norm(A) * norm(inv(A)) ). ! -- lapack computational 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: negone, zero, half, one, two, three, four, eight, ten, czero, chalf, cone, cnegone ! Scalar Arguments character, intent(in) :: norm integer(${ik}$), intent(out) :: info integer(${ik}$), intent(in) :: lda, n real(sp), intent(in) :: anorm real(sp), intent(out) :: rcond ! Array Arguments real(sp), intent(out) :: rwork(*) complex(sp), intent(inout) :: a(lda,*) complex(sp), intent(out) :: work(*) ! ===================================================================== ! Local Scalars logical(lk) :: onenrm character :: normin integer(${ik}$) :: ix, kase, kase1 real(sp) :: ainvnm, scale, sl, smlnum, su complex(sp) :: zdum ! Local Arrays integer(${ik}$) :: isave(3_${ik}$) ! Intrinsic Functions ! Statement Functions real(sp) :: cabs1 ! Statement Function Definitions cabs1( zdum ) = abs( real( zdum,KIND=sp) ) + abs( aimag( zdum ) ) ! Executable Statements ! test the input parameters. info = 0_${ik}$ onenrm = norm=='1' .or. stdlib_lsame( norm, 'O' ) if( .not.onenrm .and. .not.stdlib_lsame( norm, 'I' ) ) then info = -1_${ik}$ else if( n<0_${ik}$ ) then info = -2_${ik}$ else if( lda<max( 1_${ik}$, n ) ) then info = -4_${ik}$ else if( anorm<zero ) then info = -5_${ik}$ end if if( info/=0_${ik}$ ) then call stdlib${ii}$_xerbla( 'CGECON', -info ) return end if ! quick return if possible rcond = zero if( n==0_${ik}$ ) then rcond = one return else if( anorm==zero ) then return end if smlnum = stdlib${ii}$_slamch( 'SAFE MINIMUM' ) ! estimate the norm of inv(a). ainvnm = zero normin = 'N' if( onenrm ) then kase1 = 1_${ik}$ else kase1 = 2_${ik}$ end if kase = 0_${ik}$ 10 continue call stdlib${ii}$_clacn2( n, work( n+1 ), work, ainvnm, kase, isave ) if( kase/=0_${ik}$ ) then if( kase==kase1 ) then ! multiply by inv(l). call stdlib${ii}$_clatrs( 'LOWER', 'NO TRANSPOSE', 'UNIT', normin, n, a,lda, work, sl, & rwork, info ) ! multiply by inv(u). call stdlib${ii}$_clatrs( 'UPPER', 'NO TRANSPOSE', 'NON-UNIT', normin, n,a, lda, work, & su, rwork( n+1 ), info ) else ! multiply by inv(u**h). call stdlib${ii}$_clatrs( 'UPPER', 'CONJUGATE TRANSPOSE', 'NON-UNIT',normin, n, a, lda,& work, su, rwork( n+1 ),info ) ! multiply by inv(l**h). call stdlib${ii}$_clatrs( 'LOWER', 'CONJUGATE TRANSPOSE', 'UNIT', normin,n, a, lda, & work, sl, rwork, info ) end if ! divide x by 1/(sl*su) if doing so will not cause overflow. scale = sl*su normin = 'Y' if( scale/=one ) then ix = stdlib${ii}$_icamax( n, work, 1_${ik}$ ) if( scale<cabs1( work( ix ) )*smlnum .or. scale==zero )go to 20 call stdlib${ii}$_csrscl( n, scale, work, 1_${ik}$ ) end if go to 10 end if ! compute the estimate of the reciprocal condition number. if( ainvnm/=zero )rcond = ( one / ainvnm ) / anorm 20 continue return end subroutine stdlib${ii}$_cgecon pure module subroutine stdlib${ii}$_zgecon( norm, n, a, lda, anorm, rcond, work, rwork,info ) !! ZGECON estimates the reciprocal of the condition number of a general !! complex matrix A, in either the 1-norm or the infinity-norm, using !! the LU factorization computed by ZGETRF. !! An estimate is obtained for norm(inv(A)), and the reciprocal of the !! condition number is computed as !! RCOND = 1 / ( norm(A) * norm(inv(A)) ). ! -- lapack computational 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: negone, zero, half, one, two, three, four, eight, ten, czero, chalf, cone, cnegone ! Scalar Arguments character, intent(in) :: norm integer(${ik}$), intent(out) :: info integer(${ik}$), intent(in) :: lda, n real(dp), intent(in) :: anorm real(dp), intent(out) :: rcond ! Array Arguments real(dp), intent(out) :: rwork(*) complex(dp), intent(inout) :: a(lda,*) complex(dp), intent(out) :: work(*) ! ===================================================================== ! Local Scalars logical(lk) :: onenrm character :: normin integer(${ik}$) :: ix, kase, kase1 real(dp) :: ainvnm, scale, sl, smlnum, su complex(dp) :: zdum ! Local Arrays integer(${ik}$) :: isave(3_${ik}$) ! Intrinsic Functions ! Statement Functions real(dp) :: cabs1 ! Statement Function Definitions cabs1( zdum ) = abs( real( zdum,KIND=dp) ) + abs( aimag( zdum ) ) ! Executable Statements ! test the input parameters. info = 0_${ik}$ onenrm = norm=='1' .or. stdlib_lsame( norm, 'O' ) if( .not.onenrm .and. .not.stdlib_lsame( norm, 'I' ) ) then info = -1_${ik}$ else if( n<0_${ik}$ ) then info = -2_${ik}$ else if( lda<max( 1_${ik}$, n ) ) then info = -4_${ik}$ else if( anorm<zero ) then info = -5_${ik}$ end if if( info/=0_${ik}$ ) then call stdlib${ii}$_xerbla( 'ZGECON', -info ) return end if ! quick return if possible rcond = zero if( n==0_${ik}$ ) then rcond = one return else if( anorm==zero ) then return end if smlnum = stdlib${ii}$_dlamch( 'SAFE MINIMUM' ) ! estimate the norm of inv(a). ainvnm = zero normin = 'N' if( onenrm ) then kase1 = 1_${ik}$ else kase1 = 2_${ik}$ end if kase = 0_${ik}$ 10 continue call stdlib${ii}$_zlacn2( n, work( n+1 ), work, ainvnm, kase, isave ) if( kase/=0_${ik}$ ) then if( kase==kase1 ) then ! multiply by inv(l). call stdlib${ii}$_zlatrs( 'LOWER', 'NO TRANSPOSE', 'UNIT', normin, n, a,lda, work, sl, & rwork, info ) ! multiply by inv(u). call stdlib${ii}$_zlatrs( 'UPPER', 'NO TRANSPOSE', 'NON-UNIT', normin, n,a, lda, work, & su, rwork( n+1 ), info ) else ! multiply by inv(u**h). call stdlib${ii}$_zlatrs( 'UPPER', 'CONJUGATE TRANSPOSE', 'NON-UNIT',normin, n, a, lda,& work, su, rwork( n+1 ),info ) ! multiply by inv(l**h). call stdlib${ii}$_zlatrs( 'LOWER', 'CONJUGATE TRANSPOSE', 'UNIT', normin,n, a, lda, & work, sl, rwork, info ) end if ! divide x by 1/(sl*su) if doing so will not cause overflow. scale = sl*su normin = 'Y' if( scale/=one ) then ix = stdlib${ii}$_izamax( n, work, 1_${ik}$ ) if( scale<cabs1( work( ix ) )*smlnum .or. scale==zero )go to 20 call stdlib${ii}$_zdrscl( n, scale, work, 1_${ik}$ ) end if go to 10 end if ! compute the estimate of the reciprocal condition number. if( ainvnm/=zero )rcond = ( one / ainvnm ) / anorm 20 continue return end subroutine stdlib${ii}$_zgecon #:for ck,ct,ci in CMPLX_KINDS_TYPES #:if not ck in ["sp","dp"] pure module subroutine stdlib${ii}$_${ci}$gecon( norm, n, a, lda, anorm, rcond, work, rwork,info ) !! ZGECON: estimates the reciprocal of the condition number of a general !! complex matrix A, in either the 1-norm or the infinity-norm, using !! the LU factorization computed by ZGETRF. !! An estimate is obtained for norm(inv(A)), and the reciprocal of the !! condition number is computed as !! RCOND = 1 / ( norm(A) * norm(inv(A)) ). ! -- lapack computational 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: negone, zero, half, one, two, three, four, eight, ten, czero, chalf, cone, cnegone ! Scalar Arguments character, intent(in) :: norm integer(${ik}$), intent(out) :: info integer(${ik}$), intent(in) :: lda, n real(${ck}$), intent(in) :: anorm real(${ck}$), intent(out) :: rcond ! Array Arguments real(${ck}$), intent(out) :: rwork(*) complex(${ck}$), intent(inout) :: a(lda,*) complex(${ck}$), intent(out) :: work(*) ! ===================================================================== ! Local Scalars logical(lk) :: onenrm character :: normin integer(${ik}$) :: ix, kase, kase1 real(${ck}$) :: ainvnm, scale, sl, smlnum, su complex(${ck}$) :: zdum ! Local Arrays integer(${ik}$) :: isave(3_${ik}$) ! Intrinsic Functions ! Statement Functions real(${ck}$) :: cabs1 ! Statement Function Definitions cabs1( zdum ) = abs( real( zdum,KIND=${ck}$) ) + abs( aimag( zdum ) ) ! Executable Statements ! test the input parameters. info = 0_${ik}$ onenrm = norm=='1' .or. stdlib_lsame( norm, 'O' ) if( .not.onenrm .and. .not.stdlib_lsame( norm, 'I' ) ) then info = -1_${ik}$ else if( n<0_${ik}$ ) then info = -2_${ik}$ else if( lda<max( 1_${ik}$, n ) ) then info = -4_${ik}$ else if( anorm<zero ) then info = -5_${ik}$ end if if( info/=0_${ik}$ ) then call stdlib${ii}$_xerbla( 'ZGECON', -info ) return end if ! quick return if possible rcond = zero if( n==0_${ik}$ ) then rcond = one return else if( anorm==zero ) then return end if smlnum = stdlib${ii}$_${c2ri(ci)}$lamch( 'SAFE MINIMUM' ) ! estimate the norm of inv(a). ainvnm = zero normin = 'N' if( onenrm ) then kase1 = 1_${ik}$ else kase1 = 2_${ik}$ end if kase = 0_${ik}$ 10 continue call stdlib${ii}$_${ci}$lacn2( n, work( n+1 ), work, ainvnm, kase, isave ) if( kase/=0_${ik}$ ) then if( kase==kase1 ) then ! multiply by inv(l). call stdlib${ii}$_${ci}$latrs( 'LOWER', 'NO TRANSPOSE', 'UNIT', normin, n, a,lda, work, sl, & rwork, info ) ! multiply by inv(u). call stdlib${ii}$_${ci}$latrs( 'UPPER', 'NO TRANSPOSE', 'NON-UNIT', normin, n,a, lda, work, & su, rwork( n+1 ), info ) else ! multiply by inv(u**h). call stdlib${ii}$_${ci}$latrs( 'UPPER', 'CONJUGATE TRANSPOSE', 'NON-UNIT',normin, n, a, lda,& work, su, rwork( n+1 ),info ) ! multiply by inv(l**h). call stdlib${ii}$_${ci}$latrs( 'LOWER', 'CONJUGATE TRANSPOSE', 'UNIT', normin,n, a, lda, & work, sl, rwork, info ) end if ! divide x by 1/(sl*su) if doing so will not cause overflow. scale = sl*su normin = 'Y' if( scale/=one ) then ix = stdlib${ii}$_i${ci}$amax( n, work, 1_${ik}$ ) if( scale<cabs1( work( ix ) )*smlnum .or. scale==zero )go to 20 call stdlib${ii}$_${ci}$drscl( n, scale, work, 1_${ik}$ ) end if go to 10 end if ! compute the estimate of the reciprocal condition number. if( ainvnm/=zero )rcond = ( one / ainvnm ) / anorm 20 continue return end subroutine stdlib${ii}$_${ci}$gecon #:endif #:endfor pure module subroutine stdlib${ii}$_sgetrf( m, n, a, lda, ipiv, info ) !! SGETRF computes an LU factorization of a general M-by-N matrix A !! using partial pivoting with row interchanges. !! The factorization has the form !! A = P * L * U !! where P is a permutation matrix, L is lower triangular with unit !! diagonal elements (lower trapezoidal if m > n), and U is upper !! triangular (upper trapezoidal if m < n). !! This is the right-looking Level 3 BLAS version of the algorithm. ! -- lapack computational 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: negone, zero, half, one, two, three, four, eight, ten, czero, chalf, cone, cnegone ! Scalar Arguments integer(${ik}$), intent(out) :: info integer(${ik}$), intent(in) :: lda, m, n ! Array Arguments integer(${ik}$), intent(out) :: ipiv(*) real(sp), intent(inout) :: a(lda,*) ! ===================================================================== ! Local Scalars integer(${ik}$) :: i, iinfo, j, jb, nb ! Intrinsic Functions ! Executable Statements ! test the input parameters. info = 0_${ik}$ if( m<0_${ik}$ ) then info = -1_${ik}$ else if( n<0_${ik}$ ) then info = -2_${ik}$ else if( lda<max( 1_${ik}$, m ) ) then info = -4_${ik}$ end if if( info/=0_${ik}$ ) then call stdlib${ii}$_xerbla( 'SGETRF', -info ) return end if ! quick return if possible if( m==0 .or. n==0 )return ! determine the block size for this environment. nb = stdlib${ii}$_ilaenv( 1_${ik}$, 'SGETRF', ' ', m, n, -1_${ik}$, -1_${ik}$ ) if( nb<=1_${ik}$ .or. nb>=min( m, n ) ) then ! use unblocked code. call stdlib${ii}$_sgetrf2( m, n, a, lda, ipiv, info ) else ! use blocked code. do j = 1, min( m, n ), nb jb = min( min( m, n )-j+1, nb ) ! factor diagonal and subdiagonal blocks and test for exact ! singularity. call stdlib${ii}$_sgetrf2( m-j+1, jb, a( j, j ), lda, ipiv( j ), iinfo ) ! adjust info and the pivot indices. if( info==0_${ik}$ .and. iinfo>0_${ik}$ )info = iinfo + j - 1_${ik}$ do i = j, min( m, j+jb-1 ) ipiv( i ) = j - 1_${ik}$ + ipiv( i ) end do ! apply interchanges to columns 1:j-1. call stdlib${ii}$_slaswp( j-1, a, lda, j, j+jb-1, ipiv, 1_${ik}$ ) if( j+jb<=n ) then ! apply interchanges to columns j+jb:n. call stdlib${ii}$_slaswp( n-j-jb+1, a( 1_${ik}$, j+jb ), lda, j, j+jb-1,ipiv, 1_${ik}$ ) ! compute block row of u. call stdlib${ii}$_strsm( 'LEFT', 'LOWER', 'NO TRANSPOSE', 'UNIT', jb,n-j-jb+1, one, & a( j, j ), lda, a( j, j+jb ),lda ) if( j+jb<=m ) then ! update trailing submatrix. call stdlib${ii}$_sgemm( 'NO TRANSPOSE', 'NO TRANSPOSE', m-j-jb+1,n-j-jb+1, jb, -& one, a( j+jb, j ), lda,a( j, j+jb ), lda, one, a( j+jb, j+jb ),lda ) end if end if end do end if return end subroutine stdlib${ii}$_sgetrf pure module subroutine stdlib${ii}$_dgetrf( m, n, a, lda, ipiv, info ) !! DGETRF computes an LU factorization of a general M-by-N matrix A !! using partial pivoting with row interchanges. !! The factorization has the form !! A = P * L * U !! where P is a permutation matrix, L is lower triangular with unit !! diagonal elements (lower trapezoidal if m > n), and U is upper !! triangular (upper trapezoidal if m < n). !! This is the right-looking Level 3 BLAS version of the algorithm. ! -- lapack computational 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: negone, zero, half, one, two, three, four, eight, ten, czero, chalf, cone, cnegone ! Scalar Arguments integer(${ik}$), intent(out) :: info integer(${ik}$), intent(in) :: lda, m, n ! Array Arguments integer(${ik}$), intent(out) :: ipiv(*) real(dp), intent(inout) :: a(lda,*) ! ===================================================================== ! Local Scalars integer(${ik}$) :: i, iinfo, j, jb, nb ! Intrinsic Functions ! Executable Statements ! test the input parameters. info = 0_${ik}$ if( m<0_${ik}$ ) then info = -1_${ik}$ else if( n<0_${ik}$ ) then info = -2_${ik}$ else if( lda<max( 1_${ik}$, m ) ) then info = -4_${ik}$ end if if( info/=0_${ik}$ ) then call stdlib${ii}$_xerbla( 'DGETRF', -info ) return end if ! quick return if possible if( m==0 .or. n==0 )return ! determine the block size for this environment. nb = stdlib${ii}$_ilaenv( 1_${ik}$, 'DGETRF', ' ', m, n, -1_${ik}$, -1_${ik}$ ) if( nb<=1_${ik}$ .or. nb>=min( m, n ) ) then ! use unblocked code. call stdlib${ii}$_dgetrf2( m, n, a, lda, ipiv, info ) else ! use blocked code. do j = 1, min( m, n ), nb jb = min( min( m, n )-j+1, nb ) ! factor diagonal and subdiagonal blocks and test for exact ! singularity. call stdlib${ii}$_dgetrf2( m-j+1, jb, a( j, j ), lda, ipiv( j ), iinfo ) ! adjust info and the pivot indices. if( info==0_${ik}$ .and. iinfo>0_${ik}$ )info = iinfo + j - 1_${ik}$ do i = j, min( m, j+jb-1 ) ipiv( i ) = j - 1_${ik}$ + ipiv( i ) end do ! apply interchanges to columns 1:j-1. call stdlib${ii}$_dlaswp( j-1, a, lda, j, j+jb-1, ipiv, 1_${ik}$ ) if( j+jb<=n ) then ! apply interchanges to columns j+jb:n. call stdlib${ii}$_dlaswp( n-j-jb+1, a( 1_${ik}$, j+jb ), lda, j, j+jb-1,ipiv, 1_${ik}$ ) ! compute block row of u. call stdlib${ii}$_dtrsm( 'LEFT', 'LOWER', 'NO TRANSPOSE', 'UNIT', jb,n-j-jb+1, one, & a( j, j ), lda, a( j, j+jb ),lda ) if( j+jb<=m ) then ! update trailing submatrix. call stdlib${ii}$_dgemm( 'NO TRANSPOSE', 'NO TRANSPOSE', m-j-jb+1,n-j-jb+1, jb, -& one, a( j+jb, j ), lda,a( j, j+jb ), lda, one, a( j+jb, j+jb ),lda ) end if end if end do end if return end subroutine stdlib${ii}$_dgetrf #:for rk,rt,ri in REAL_KINDS_TYPES #:if not rk in ["sp","dp"] pure module subroutine stdlib${ii}$_${ri}$getrf( m, n, a, lda, ipiv, info ) !! DGETRF: computes an LU factorization of a general M-by-N matrix A !! using partial pivoting with row interchanges. !! The factorization has the form !! A = P * L * U !! where P is a permutation matrix, L is lower triangular with unit !! diagonal elements (lower trapezoidal if m > n), and U is upper !! triangular (upper trapezoidal if m < n). !! This is the right-looking Level 3 BLAS version of the algorithm. ! -- lapack computational 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: negone, zero, half, one, two, three, four, eight, ten, czero, chalf, cone, cnegone ! Scalar Arguments integer(${ik}$), intent(out) :: info integer(${ik}$), intent(in) :: lda, m, n ! Array Arguments integer(${ik}$), intent(out) :: ipiv(*) real(${rk}$), intent(inout) :: a(lda,*) ! ===================================================================== ! Local Scalars integer(${ik}$) :: i, iinfo, j, jb, nb ! Intrinsic Functions ! Executable Statements ! test the input parameters. info = 0_${ik}$ if( m<0_${ik}$ ) then info = -1_${ik}$ else if( n<0_${ik}$ ) then info = -2_${ik}$ else if( lda<max( 1_${ik}$, m ) ) then info = -4_${ik}$ end if if( info/=0_${ik}$ ) then call stdlib${ii}$_xerbla( 'DGETRF', -info ) return end if ! quick return if possible if( m==0 .or. n==0 )return ! determine the block size for this environment. nb = stdlib${ii}$_ilaenv( 1_${ik}$, 'DGETRF', ' ', m, n, -1_${ik}$, -1_${ik}$ ) if( nb<=1_${ik}$ .or. nb>=min( m, n ) ) then ! use unblocked code. call stdlib${ii}$_${ri}$getrf2( m, n, a, lda, ipiv, info ) else ! use blocked code. do j = 1, min( m, n ), nb jb = min( min( m, n )-j+1, nb ) ! factor diagonal and subdiagonal blocks and test for exact ! singularity. call stdlib${ii}$_${ri}$getrf2( m-j+1, jb, a( j, j ), lda, ipiv( j ), iinfo ) ! adjust info and the pivot indices. if( info==0_${ik}$ .and. iinfo>0_${ik}$ )info = iinfo + j - 1_${ik}$ do i = j, min( m, j+jb-1 ) ipiv( i ) = j - 1_${ik}$ + ipiv( i ) end do ! apply interchanges to columns 1:j-1. call stdlib${ii}$_${ri}$laswp( j-1, a, lda, j, j+jb-1, ipiv, 1_${ik}$ ) if( j+jb<=n ) then ! apply interchanges to columns j+jb:n. call stdlib${ii}$_${ri}$laswp( n-j-jb+1, a( 1_${ik}$, j+jb ), lda, j, j+jb-1,ipiv, 1_${ik}$ ) ! compute block row of u. call stdlib${ii}$_${ri}$trsm( 'LEFT', 'LOWER', 'NO TRANSPOSE', 'UNIT', jb,n-j-jb+1, one, & a( j, j ), lda, a( j, j+jb ),lda ) if( j+jb<=m ) then ! update trailing submatrix. call stdlib${ii}$_${ri}$gemm( 'NO TRANSPOSE', 'NO TRANSPOSE', m-j-jb+1,n-j-jb+1, jb, -& one, a( j+jb, j ), lda,a( j, j+jb ), lda, one, a( j+jb, j+jb ),lda ) end if end if end do end if return end subroutine stdlib${ii}$_${ri}$getrf #:endif #:endfor pure module subroutine stdlib${ii}$_cgetrf( m, n, a, lda, ipiv, info ) !! CGETRF computes an LU factorization of a general M-by-N matrix A !! using partial pivoting with row interchanges. !! The factorization has the form !! A = P * L * U !! where P is a permutation matrix, L is lower triangular with unit !! diagonal elements (lower trapezoidal if m > n), and U is upper !! triangular (upper trapezoidal if m < n). !! This is the right-looking Level 3 BLAS version of the algorithm. ! -- lapack computational 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: negone, zero, half, one, two, three, four, eight, ten, czero, chalf, cone, cnegone ! Scalar Arguments integer(${ik}$), intent(out) :: info integer(${ik}$), intent(in) :: lda, m, n ! Array Arguments integer(${ik}$), intent(out) :: ipiv(*) complex(sp), intent(inout) :: a(lda,*) ! ===================================================================== ! Local Scalars integer(${ik}$) :: i, iinfo, j, jb, nb ! Intrinsic Functions ! Executable Statements ! test the input parameters. info = 0_${ik}$ if( m<0_${ik}$ ) then info = -1_${ik}$ else if( n<0_${ik}$ ) then info = -2_${ik}$ else if( lda<max( 1_${ik}$, m ) ) then info = -4_${ik}$ end if if( info/=0_${ik}$ ) then call stdlib${ii}$_xerbla( 'CGETRF', -info ) return end if ! quick return if possible if( m==0 .or. n==0 )return ! determine the block size for this environment. nb = stdlib${ii}$_ilaenv( 1_${ik}$, 'CGETRF', ' ', m, n, -1_${ik}$, -1_${ik}$ ) if( nb<=1_${ik}$ .or. nb>=min( m, n ) ) then ! use unblocked code. call stdlib${ii}$_cgetrf2( m, n, a, lda, ipiv, info ) else ! use blocked code. do j = 1, min( m, n ), nb jb = min( min( m, n )-j+1, nb ) ! factor diagonal and subdiagonal blocks and test for exact ! singularity. call stdlib${ii}$_cgetrf2( m-j+1, jb, a( j, j ), lda, ipiv( j ), iinfo ) ! adjust info and the pivot indices. if( info==0_${ik}$ .and. iinfo>0_${ik}$ )info = iinfo + j - 1_${ik}$ do i = j, min( m, j+jb-1 ) ipiv( i ) = j - 1_${ik}$ + ipiv( i ) end do ! apply interchanges to columns 1:j-1. call stdlib${ii}$_claswp( j-1, a, lda, j, j+jb-1, ipiv, 1_${ik}$ ) if( j+jb<=n ) then ! apply interchanges to columns j+jb:n. call stdlib${ii}$_claswp( n-j-jb+1, a( 1_${ik}$, j+jb ), lda, j, j+jb-1,ipiv, 1_${ik}$ ) ! compute block row of u. call stdlib${ii}$_ctrsm( 'LEFT', 'LOWER', 'NO TRANSPOSE', 'UNIT', jb,n-j-jb+1, cone,& a( j, j ), lda, a( j, j+jb ),lda ) if( j+jb<=m ) then ! update trailing submatrix. call stdlib${ii}$_cgemm( 'NO TRANSPOSE', 'NO TRANSPOSE', m-j-jb+1,n-j-jb+1, jb, -& cone, a( j+jb, j ), lda,a( j, j+jb ), lda, cone, a( j+jb, j+jb ),lda ) end if end if end do end if return end subroutine stdlib${ii}$_cgetrf pure module subroutine stdlib${ii}$_zgetrf( m, n, a, lda, ipiv, info ) !! ZGETRF computes an LU factorization of a general M-by-N matrix A !! using partial pivoting with row interchanges. !! The factorization has the form !! A = P * L * U !! where P is a permutation matrix, L is lower triangular with unit !! diagonal elements (lower trapezoidal if m > n), and U is upper !! triangular (upper trapezoidal if m < n). !! This is the right-looking Level 3 BLAS version of the algorithm. ! -- lapack computational 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: negone, zero, half, one, two, three, four, eight, ten, czero, chalf, cone, cnegone ! Scalar Arguments integer(${ik}$), intent(out) :: info integer(${ik}$), intent(in) :: lda, m, n ! Array Arguments integer(${ik}$), intent(out) :: ipiv(*) complex(dp), intent(inout) :: a(lda,*) ! ===================================================================== ! Local Scalars integer(${ik}$) :: i, iinfo, j, jb, nb ! Intrinsic Functions ! Executable Statements ! test the input parameters. info = 0_${ik}$ if( m<0_${ik}$ ) then info = -1_${ik}$ else if( n<0_${ik}$ ) then info = -2_${ik}$ else if( lda<max( 1_${ik}$, m ) ) then info = -4_${ik}$ end if if( info/=0_${ik}$ ) then call stdlib${ii}$_xerbla( 'ZGETRF', -info ) return end if ! quick return if possible if( m==0 .or. n==0 )return ! determine the block size for this environment. nb = stdlib${ii}$_ilaenv( 1_${ik}$, 'ZGETRF', ' ', m, n, -1_${ik}$, -1_${ik}$ ) if( nb<=1_${ik}$ .or. nb>=min( m, n ) ) then ! use unblocked code. call stdlib${ii}$_zgetrf2( m, n, a, lda, ipiv, info ) else ! use blocked code. do j = 1, min( m, n ), nb jb = min( min( m, n )-j+1, nb ) ! factor diagonal and subdiagonal blocks and test for exact ! singularity. call stdlib${ii}$_zgetrf2( m-j+1, jb, a( j, j ), lda, ipiv( j ), iinfo ) ! adjust info and the pivot indices. if( info==0_${ik}$ .and. iinfo>0_${ik}$ )info = iinfo + j - 1_${ik}$ do i = j, min( m, j+jb-1 ) ipiv( i ) = j - 1_${ik}$ + ipiv( i ) end do ! apply interchanges to columns 1:j-1. call stdlib${ii}$_zlaswp( j-1, a, lda, j, j+jb-1, ipiv, 1_${ik}$ ) if( j+jb<=n ) then ! apply interchanges to columns j+jb:n. call stdlib${ii}$_zlaswp( n-j-jb+1, a( 1_${ik}$, j+jb ), lda, j, j+jb-1,ipiv, 1_${ik}$ ) ! compute block row of u. call stdlib${ii}$_ztrsm( 'LEFT', 'LOWER', 'NO TRANSPOSE', 'UNIT', jb,n-j-jb+1, cone,& a( j, j ), lda, a( j, j+jb ),lda ) if( j+jb<=m ) then ! update trailing submatrix. call stdlib${ii}$_zgemm( 'NO TRANSPOSE', 'NO TRANSPOSE', m-j-jb+1,n-j-jb+1, jb, -& cone, a( j+jb, j ), lda,a( j, j+jb ), lda, cone, a( j+jb, j+jb ),lda ) end if end if end do end if return end subroutine stdlib${ii}$_zgetrf #:for ck,ct,ci in CMPLX_KINDS_TYPES #:if not ck in ["sp","dp"] pure module subroutine stdlib${ii}$_${ci}$getrf( m, n, a, lda, ipiv, info ) !! ZGETRF: computes an LU factorization of a general M-by-N matrix A !! using partial pivoting with row interchanges. !! The factorization has the form !! A = P * L * U !! where P is a permutation matrix, L is lower triangular with unit !! diagonal elements (lower trapezoidal if m > n), and U is upper !! triangular (upper trapezoidal if m < n). !! This is the right-looking Level 3 BLAS version of the algorithm. ! -- lapack computational 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: negone, zero, half, one, two, three, four, eight, ten, czero, chalf, cone, cnegone ! Scalar Arguments integer(${ik}$), intent(out) :: info integer(${ik}$), intent(in) :: lda, m, n ! Array Arguments integer(${ik}$), intent(out) :: ipiv(*) complex(${ck}$), intent(inout) :: a(lda,*) ! ===================================================================== ! Local Scalars integer(${ik}$) :: i, iinfo, j, jb, nb ! Intrinsic Functions ! Executable Statements ! test the input parameters. info = 0_${ik}$ if( m<0_${ik}$ ) then info = -1_${ik}$ else if( n<0_${ik}$ ) then info = -2_${ik}$ else if( lda<max( 1_${ik}$, m ) ) then info = -4_${ik}$ end if if( info/=0_${ik}$ ) then call stdlib${ii}$_xerbla( 'ZGETRF', -info ) return end if ! quick return if possible if( m==0 .or. n==0 )return ! determine the block size for this environment. nb = stdlib${ii}$_ilaenv( 1_${ik}$, 'ZGETRF', ' ', m, n, -1_${ik}$, -1_${ik}$ ) if( nb<=1_${ik}$ .or. nb>=min( m, n ) ) then ! use unblocked code. call stdlib${ii}$_${ci}$getrf2( m, n, a, lda, ipiv, info ) else ! use blocked code. do j = 1, min( m, n ), nb jb = min( min( m, n )-j+1, nb ) ! factor diagonal and subdiagonal blocks and test for exact ! singularity. call stdlib${ii}$_${ci}$getrf2( m-j+1, jb, a( j, j ), lda, ipiv( j ), iinfo ) ! adjust info and the pivot indices. if( info==0_${ik}$ .and. iinfo>0_${ik}$ )info = iinfo + j - 1_${ik}$ do i = j, min( m, j+jb-1 ) ipiv( i ) = j - 1_${ik}$ + ipiv( i ) end do ! apply interchanges to columns 1:j-1. call stdlib${ii}$_${ci}$laswp( j-1, a, lda, j, j+jb-1, ipiv, 1_${ik}$ ) if( j+jb<=n ) then ! apply interchanges to columns j+jb:n. call stdlib${ii}$_${ci}$laswp( n-j-jb+1, a( 1_${ik}$, j+jb ), lda, j, j+jb-1,ipiv, 1_${ik}$ ) ! compute block row of u. call stdlib${ii}$_${ci}$trsm( 'LEFT', 'LOWER', 'NO TRANSPOSE', 'UNIT', jb,n-j-jb+1, cone,& a( j, j ), lda, a( j, j+jb ),lda ) if( j+jb<=m ) then ! update trailing submatrix. call stdlib${ii}$_${ci}$gemm( 'NO TRANSPOSE', 'NO TRANSPOSE', m-j-jb+1,n-j-jb+1, jb, -& cone, a( j+jb, j ), lda,a( j, j+jb ), lda, cone, a( j+jb, j+jb ),lda ) end if end if end do end if return end subroutine stdlib${ii}$_${ci}$getrf #:endif #:endfor pure recursive module subroutine stdlib${ii}$_sgetrf2( m, n, a, lda, ipiv, info ) !! SGETRF2 computes an LU factorization of a general M-by-N matrix A !! using partial pivoting with row interchanges. !! The factorization has the form !! A = P * L * U !! where P is a permutation matrix, L is lower triangular with unit !! diagonal elements (lower trapezoidal if m > n), and U is upper !! triangular (upper trapezoidal if m < n). !! This is the recursive version of the algorithm. It divides !! the matrix into four submatrices: !! [ A11 | A12 ] where A11 is n1 by n1 and A22 is n2 by n2 !! A = [ -----|----- ] with n1 = min(m,n)/2 !! [ A21 | A22 ] n2 = n-n1 !! [ A11 ] !! The subroutine calls itself to factor [ --- ], !! [ A12 ] !! [ A12 ] !! do the swaps on [ --- ], solve A12, update A22, !! [ A22 ] !! then calls itself to factor A22 and do the swaps on A21. ! -- lapack computational 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: negone, zero, half, one, two, three, four, eight, ten, czero, chalf, cone, cnegone ! Scalar Arguments integer(${ik}$), intent(out) :: info integer(${ik}$), intent(in) :: lda, m, n ! Array Arguments integer(${ik}$), intent(out) :: ipiv(*) real(sp), intent(inout) :: a(lda,*) ! ===================================================================== ! Local Scalars real(sp) :: sfmin, temp integer(${ik}$) :: i, iinfo, n1, n2 ! Intrinsic Functions ! Executable Statements ! test the input parameters info = 0_${ik}$ if( m<0_${ik}$ ) then info = -1_${ik}$ else if( n<0_${ik}$ ) then info = -2_${ik}$ else if( lda<max( 1_${ik}$, m ) ) then info = -4_${ik}$ end if if( info/=0_${ik}$ ) then call stdlib${ii}$_xerbla( 'SGETRF2', -info ) return end if ! quick return if possible if( m==0 .or. n==0 )return if ( m==1_${ik}$ ) then ! use unblocked code for one row case ! just need to handle ipiv and info ipiv( 1_${ik}$ ) = 1_${ik}$ if ( a(1_${ik}$,1_${ik}$)==zero )info = 1_${ik}$ else if( n==1_${ik}$ ) then ! use unblocked code for one column case ! compute machine safe minimum sfmin = stdlib${ii}$_slamch('S') ! find pivot and test for singularity i = stdlib${ii}$_isamax( m, a( 1_${ik}$, 1_${ik}$ ), 1_${ik}$ ) ipiv( 1_${ik}$ ) = i if( a( i, 1_${ik}$ )/=zero ) then ! apply the interchange if( i/=1_${ik}$ ) then temp = a( 1_${ik}$, 1_${ik}$ ) a( 1_${ik}$, 1_${ik}$ ) = a( i, 1_${ik}$ ) a( i, 1_${ik}$ ) = temp end if ! compute elements 2:m of the column if( abs(a( 1_${ik}$, 1_${ik}$ )) >= sfmin ) then call stdlib${ii}$_sscal( m-1, one / a( 1_${ik}$, 1_${ik}$ ), a( 2_${ik}$, 1_${ik}$ ), 1_${ik}$ ) else do i = 1, m-1 a( 1_${ik}$+i, 1_${ik}$ ) = a( 1_${ik}$+i, 1_${ik}$ ) / a( 1_${ik}$, 1_${ik}$ ) end do end if else info = 1_${ik}$ end if else ! use recursive code n1 = min( m, n ) / 2_${ik}$ n2 = n-n1 ! [ a11 ] ! factor [ --- ] ! [ a21 ] call stdlib${ii}$_sgetrf2( m, n1, a, lda, ipiv, iinfo ) if ( info==0_${ik}$ .and. iinfo>0_${ik}$ )info = iinfo ! [ a12 ] ! apply interchanges to [ --- ] ! [ a22 ] call stdlib${ii}$_slaswp( n2, a( 1_${ik}$, n1+1 ), lda, 1_${ik}$, n1, ipiv, 1_${ik}$ ) ! solve a12 call stdlib${ii}$_strsm( 'L', 'L', 'N', 'U', n1, n2, one, a, lda,a( 1_${ik}$, n1+1 ), lda ) ! update a22 call stdlib${ii}$_sgemm( 'N', 'N', m-n1, n2, n1, -one, a( n1+1, 1_${ik}$ ), lda,a( 1_${ik}$, n1+1 ), & lda, one, a( n1+1, n1+1 ), lda ) ! factor a22 call stdlib${ii}$_sgetrf2( m-n1, n2, a( n1+1, n1+1 ), lda, ipiv( n1+1 ),iinfo ) ! adjust info and the pivot indices if ( info==0_${ik}$ .and. iinfo>0_${ik}$ )info = iinfo + n1 do i = n1+1, min( m, n ) ipiv( i ) = ipiv( i ) + n1 end do ! apply interchanges to a21 call stdlib${ii}$_slaswp( n1, a( 1_${ik}$, 1_${ik}$ ), lda, n1+1, min( m, n), ipiv, 1_${ik}$ ) end if return end subroutine stdlib${ii}$_sgetrf2 pure recursive module subroutine stdlib${ii}$_dgetrf2( m, n, a, lda, ipiv, info ) !! DGETRF2 computes an LU factorization of a general M-by-N matrix A !! using partial pivoting with row interchanges. !! The factorization has the form !! A = P * L * U !! where P is a permutation matrix, L is lower triangular with unit !! diagonal elements (lower trapezoidal if m > n), and U is upper !! triangular (upper trapezoidal if m < n). !! This is the recursive version of the algorithm. It divides !! the matrix into four submatrices: !! [ A11 | A12 ] where A11 is n1 by n1 and A22 is n2 by n2 !! A = [ -----|----- ] with n1 = min(m,n)/2 !! [ A21 | A22 ] n2 = n-n1 !! [ A11 ] !! The subroutine calls itself to factor [ --- ], !! [ A12 ] !! [ A12 ] !! do the swaps on [ --- ], solve A12, update A22, !! [ A22 ] !! then calls itself to factor A22 and do the swaps on A21. ! -- lapack computational 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: negone, zero, half, one, two, three, four, eight, ten, czero, chalf, cone, cnegone ! Scalar Arguments integer(${ik}$), intent(out) :: info integer(${ik}$), intent(in) :: lda, m, n ! Array Arguments integer(${ik}$), intent(out) :: ipiv(*) real(dp), intent(inout) :: a(lda,*) ! ===================================================================== ! Local Scalars real(dp) :: sfmin, temp integer(${ik}$) :: i, iinfo, n1, n2 ! Intrinsic Functions ! Executable Statements ! test the input parameters info = 0_${ik}$ if( m<0_${ik}$ ) then info = -1_${ik}$ else if( n<0_${ik}$ ) then info = -2_${ik}$ else if( lda<max( 1_${ik}$, m ) ) then info = -4_${ik}$ end if if( info/=0_${ik}$ ) then call stdlib${ii}$_xerbla( 'DGETRF2', -info ) return end if ! quick return if possible if( m==0 .or. n==0 )return if ( m==1_${ik}$ ) then ! use unblocked code for one row case ! just need to handle ipiv and info ipiv( 1_${ik}$ ) = 1_${ik}$ if ( a(1_${ik}$,1_${ik}$)==zero )info = 1_${ik}$ else if( n==1_${ik}$ ) then ! use unblocked code for one column case ! compute machine safe minimum sfmin = stdlib${ii}$_dlamch('S') ! find pivot and test for singularity i = stdlib${ii}$_idamax( m, a( 1_${ik}$, 1_${ik}$ ), 1_${ik}$ ) ipiv( 1_${ik}$ ) = i if( a( i, 1_${ik}$ )/=zero ) then ! apply the interchange if( i/=1_${ik}$ ) then temp = a( 1_${ik}$, 1_${ik}$ ) a( 1_${ik}$, 1_${ik}$ ) = a( i, 1_${ik}$ ) a( i, 1_${ik}$ ) = temp end if ! compute elements 2:m of the column if( abs(a( 1_${ik}$, 1_${ik}$ )) >= sfmin ) then call stdlib${ii}$_dscal( m-1, one / a( 1_${ik}$, 1_${ik}$ ), a( 2_${ik}$, 1_${ik}$ ), 1_${ik}$ ) else do i = 1, m-1 a( 1_${ik}$+i, 1_${ik}$ ) = a( 1_${ik}$+i, 1_${ik}$ ) / a( 1_${ik}$, 1_${ik}$ ) end do end if else info = 1_${ik}$ end if else ! use recursive code n1 = min( m, n ) / 2_${ik}$ n2 = n-n1 ! [ a11 ] ! factor [ --- ] ! [ a21 ] call stdlib${ii}$_dgetrf2( m, n1, a, lda, ipiv, iinfo ) if ( info==0_${ik}$ .and. iinfo>0_${ik}$ )info = iinfo ! [ a12 ] ! apply interchanges to [ --- ] ! [ a22 ] call stdlib${ii}$_dlaswp( n2, a( 1_${ik}$, n1+1 ), lda, 1_${ik}$, n1, ipiv, 1_${ik}$ ) ! solve a12 call stdlib${ii}$_dtrsm( 'L', 'L', 'N', 'U', n1, n2, one, a, lda,a( 1_${ik}$, n1+1 ), lda ) ! update a22 call stdlib${ii}$_dgemm( 'N', 'N', m-n1, n2, n1, -one, a( n1+1, 1_${ik}$ ), lda,a( 1_${ik}$, n1+1 ), & lda, one, a( n1+1, n1+1 ), lda ) ! factor a22 call stdlib${ii}$_dgetrf2( m-n1, n2, a( n1+1, n1+1 ), lda, ipiv( n1+1 ),iinfo ) ! adjust info and the pivot indices if ( info==0_${ik}$ .and. iinfo>0_${ik}$ )info = iinfo + n1 do i = n1+1, min( m, n ) ipiv( i ) = ipiv( i ) + n1 end do ! apply interchanges to a21 call stdlib${ii}$_dlaswp( n1, a( 1_${ik}$, 1_${ik}$ ), lda, n1+1, min( m, n), ipiv, 1_${ik}$ ) end if return end subroutine stdlib${ii}$_dgetrf2 #:for rk,rt,ri in REAL_KINDS_TYPES #:if not rk in ["sp","dp"] pure recursive module subroutine stdlib${ii}$_${ri}$getrf2( m, n, a, lda, ipiv, info ) !! DGETRF2: computes an LU factorization of a general M-by-N matrix A !! using partial pivoting with row interchanges. !! The factorization has the form !! A = P * L * U !! where P is a permutation matrix, L is lower triangular with unit !! diagonal elements (lower trapezoidal if m > n), and U is upper !! triangular (upper trapezoidal if m < n). !! This is the recursive version of the algorithm. It divides !! the matrix into four submatrices: !! [ A11 | A12 ] where A11 is n1 by n1 and A22 is n2 by n2 !! A = [ -----|----- ] with n1 = min(m,n)/2 !! [ A21 | A22 ] n2 = n-n1 !! [ A11 ] !! The subroutine calls itself to factor [ --- ], !! [ A12 ] !! [ A12 ] !! do the swaps on [ --- ], solve A12, update A22, !! [ A22 ] !! then calls itself to factor A22 and do the swaps on A21. ! -- lapack computational 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: negone, zero, half, one, two, three, four, eight, ten, czero, chalf, cone, cnegone ! Scalar Arguments integer(${ik}$), intent(out) :: info integer(${ik}$), intent(in) :: lda, m, n ! Array Arguments integer(${ik}$), intent(out) :: ipiv(*) real(${rk}$), intent(inout) :: a(lda,*) ! ===================================================================== ! Local Scalars real(${rk}$) :: sfmin, temp integer(${ik}$) :: i, iinfo, n1, n2 ! Intrinsic Functions ! Executable Statements ! test the input parameters info = 0_${ik}$ if( m<0_${ik}$ ) then info = -1_${ik}$ else if( n<0_${ik}$ ) then info = -2_${ik}$ else if( lda<max( 1_${ik}$, m ) ) then info = -4_${ik}$ end if if( info/=0_${ik}$ ) then call stdlib${ii}$_xerbla( 'DGETRF2', -info ) return end if ! quick return if possible if( m==0 .or. n==0 )return if ( m==1_${ik}$ ) then ! use unblocked code for one row case ! just need to handle ipiv and info ipiv( 1_${ik}$ ) = 1_${ik}$ if ( a(1_${ik}$,1_${ik}$)==zero )info = 1_${ik}$ else if( n==1_${ik}$ ) then ! use unblocked code for one column case ! compute machine safe minimum sfmin = stdlib${ii}$_${ri}$lamch('S') ! find pivot and test for singularity i = stdlib${ii}$_i${ri}$amax( m, a( 1_${ik}$, 1_${ik}$ ), 1_${ik}$ ) ipiv( 1_${ik}$ ) = i if( a( i, 1_${ik}$ )/=zero ) then ! apply the interchange if( i/=1_${ik}$ ) then temp = a( 1_${ik}$, 1_${ik}$ ) a( 1_${ik}$, 1_${ik}$ ) = a( i, 1_${ik}$ ) a( i, 1_${ik}$ ) = temp end if ! compute elements 2:m of the column if( abs(a( 1_${ik}$, 1_${ik}$ )) >= sfmin ) then call stdlib${ii}$_${ri}$scal( m-1, one / a( 1_${ik}$, 1_${ik}$ ), a( 2_${ik}$, 1_${ik}$ ), 1_${ik}$ ) else do i = 1, m-1 a( 1_${ik}$+i, 1_${ik}$ ) = a( 1_${ik}$+i, 1_${ik}$ ) / a( 1_${ik}$, 1_${ik}$ ) end do end if else info = 1_${ik}$ end if else ! use recursive code n1 = min( m, n ) / 2_${ik}$ n2 = n-n1 ! [ a11 ] ! factor [ --- ] ! [ a21 ] call stdlib${ii}$_${ri}$getrf2( m, n1, a, lda, ipiv, iinfo ) if ( info==0_${ik}$ .and. iinfo>0_${ik}$ )info = iinfo ! [ a12 ] ! apply interchanges to [ --- ] ! [ a22 ] call stdlib${ii}$_${ri}$laswp( n2, a( 1_${ik}$, n1+1 ), lda, 1_${ik}$, n1, ipiv, 1_${ik}$ ) ! solve a12 call stdlib${ii}$_${ri}$trsm( 'L', 'L', 'N', 'U', n1, n2, one, a, lda,a( 1_${ik}$, n1+1 ), lda ) ! update a22 call stdlib${ii}$_${ri}$gemm( 'N', 'N', m-n1, n2, n1, -one, a( n1+1, 1_${ik}$ ), lda,a( 1_${ik}$, n1+1 ), & lda, one, a( n1+1, n1+1 ), lda ) ! factor a22 call stdlib${ii}$_${ri}$getrf2( m-n1, n2, a( n1+1, n1+1 ), lda, ipiv( n1+1 ),iinfo ) ! adjust info and the pivot indices if ( info==0_${ik}$ .and. iinfo>0_${ik}$ )info = iinfo + n1 do i = n1+1, min( m, n ) ipiv( i ) = ipiv( i ) + n1 end do ! apply interchanges to a21 call stdlib${ii}$_${ri}$laswp( n1, a( 1_${ik}$, 1_${ik}$ ), lda, n1+1, min( m, n), ipiv, 1_${ik}$ ) end if return end subroutine stdlib${ii}$_${ri}$getrf2 #:endif #:endfor pure recursive module subroutine stdlib${ii}$_cgetrf2( m, n, a, lda, ipiv, info ) !! CGETRF2 computes an LU factorization of a general M-by-N matrix A !! using partial pivoting with row interchanges. !! The factorization has the form !! A = P * L * U !! where P is a permutation matrix, L is lower triangular with unit !! diagonal elements (lower trapezoidal if m > n), and U is upper !! triangular (upper trapezoidal if m < n). !! This is the recursive version of the algorithm. It divides !! the matrix into four submatrices: !! [ A11 | A12 ] where A11 is n1 by n1 and A22 is n2 by n2 !! A = [ -----|----- ] with n1 = min(m,n)/2 !! [ A21 | A22 ] n2 = n-n1 !! [ A11 ] !! The subroutine calls itself to factor [ --- ], !! [ A12 ] !! [ A12 ] !! do the swaps on [ --- ], solve A12, update A22, !! [ A22 ] !! then calls itself to factor A22 and do the swaps on A21. ! -- lapack computational 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: negone, zero, half, one, two, three, four, eight, ten, czero, chalf, cone, cnegone ! Scalar Arguments integer(${ik}$), intent(out) :: info integer(${ik}$), intent(in) :: lda, m, n ! Array Arguments integer(${ik}$), intent(out) :: ipiv(*) complex(sp), intent(inout) :: a(lda,*) ! ===================================================================== ! Local Scalars real(sp) :: sfmin complex(sp) :: temp integer(${ik}$) :: i, iinfo, n1, n2 ! Intrinsic Functions ! Executable Statements ! test the input parameters info = 0_${ik}$ if( m<0_${ik}$ ) then info = -1_${ik}$ else if( n<0_${ik}$ ) then info = -2_${ik}$ else if( lda<max( 1_${ik}$, m ) ) then info = -4_${ik}$ end if if( info/=0_${ik}$ ) then call stdlib${ii}$_xerbla( 'CGETRF2', -info ) return end if ! quick return if possible if( m==0 .or. n==0 )return if ( m==1_${ik}$ ) then ! use unblocked code for cone row case ! just need to handle ipiv and info ipiv( 1_${ik}$ ) = 1_${ik}$ if ( a(1_${ik}$,1_${ik}$)==czero )info = 1_${ik}$ else if( n==1_${ik}$ ) then ! use unblocked code for cone column case ! compute machine safe minimum sfmin = stdlib${ii}$_slamch('S') ! find pivot and test for singularity i = stdlib${ii}$_icamax( m, a( 1_${ik}$, 1_${ik}$ ), 1_${ik}$ ) ipiv( 1_${ik}$ ) = i if( a( i, 1_${ik}$ )/=czero ) then ! apply the interchange if( i/=1_${ik}$ ) then temp = a( 1_${ik}$, 1_${ik}$ ) a( 1_${ik}$, 1_${ik}$ ) = a( i, 1_${ik}$ ) a( i, 1_${ik}$ ) = temp end if ! compute elements 2:m of the column if( abs(a( 1_${ik}$, 1_${ik}$ )) >= sfmin ) then call stdlib${ii}$_cscal( m-1, cone / a( 1_${ik}$, 1_${ik}$ ), a( 2_${ik}$, 1_${ik}$ ), 1_${ik}$ ) else do i = 1, m-1 a( 1_${ik}$+i, 1_${ik}$ ) = a( 1_${ik}$+i, 1_${ik}$ ) / a( 1_${ik}$, 1_${ik}$ ) end do end if else info = 1_${ik}$ end if else ! use recursive code n1 = min( m, n ) / 2_${ik}$ n2 = n-n1 ! [ a11 ] ! factor [ --- ] ! [ a21 ] call stdlib${ii}$_cgetrf2( m, n1, a, lda, ipiv, iinfo ) if ( info==0_${ik}$ .and. iinfo>0_${ik}$ )info = iinfo ! [ a12 ] ! apply interchanges to [ --- ] ! [ a22 ] call stdlib${ii}$_claswp( n2, a( 1_${ik}$, n1+1 ), lda, 1_${ik}$, n1, ipiv, 1_${ik}$ ) ! solve a12 call stdlib${ii}$_ctrsm( 'L', 'L', 'N', 'U', n1, n2, cone, a, lda,a( 1_${ik}$, n1+1 ), lda ) ! update a22 call stdlib${ii}$_cgemm( 'N', 'N', m-n1, n2, n1, -cone, a( n1+1, 1_${ik}$ ), lda,a( 1_${ik}$, n1+1 ), & lda, cone, a( n1+1, n1+1 ), lda ) ! factor a22 call stdlib${ii}$_cgetrf2( m-n1, n2, a( n1+1, n1+1 ), lda, ipiv( n1+1 ),iinfo ) ! adjust info and the pivot indices if ( info==0_${ik}$ .and. iinfo>0_${ik}$ )info = iinfo + n1 do i = n1+1, min( m, n ) ipiv( i ) = ipiv( i ) + n1 end do ! apply interchanges to a21 call stdlib${ii}$_claswp( n1, a( 1_${ik}$, 1_${ik}$ ), lda, n1+1, min( m, n), ipiv, 1_${ik}$ ) end if return end subroutine stdlib${ii}$_cgetrf2 pure recursive module subroutine stdlib${ii}$_zgetrf2( m, n, a, lda, ipiv, info ) !! ZGETRF2 computes an LU factorization of a general M-by-N matrix A !! using partial pivoting with row interchanges. !! The factorization has the form !! A = P * L * U !! where P is a permutation matrix, L is lower triangular with unit !! diagonal elements (lower trapezoidal if m > n), and U is upper !! triangular (upper trapezoidal if m < n). !! This is the recursive version of the algorithm. It divides !! the matrix into four submatrices: !! [ A11 | A12 ] where A11 is n1 by n1 and A22 is n2 by n2 !! A = [ -----|----- ] with n1 = min(m,n)/2 !! [ A21 | A22 ] n2 = n-n1 !! [ A11 ] !! The subroutine calls itself to factor [ --- ], !! [ A12 ] !! [ A12 ] !! do the swaps on [ --- ], solve A12, update A22, !! [ A22 ] !! then calls itself to factor A22 and do the swaps on A21. ! -- lapack computational 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: negone, zero, half, one, two, three, four, eight, ten, czero, chalf, cone, cnegone ! Scalar Arguments integer(${ik}$), intent(out) :: info integer(${ik}$), intent(in) :: lda, m, n ! Array Arguments integer(${ik}$), intent(out) :: ipiv(*) complex(dp), intent(inout) :: a(lda,*) ! ===================================================================== ! Local Scalars real(dp) :: sfmin complex(dp) :: temp integer(${ik}$) :: i, iinfo, n1, n2 ! Intrinsic Functions ! Executable Statements ! test the input parameters info = 0_${ik}$ if( m<0_${ik}$ ) then info = -1_${ik}$ else if( n<0_${ik}$ ) then info = -2_${ik}$ else if( lda<max( 1_${ik}$, m ) ) then info = -4_${ik}$ end if if( info/=0_${ik}$ ) then call stdlib${ii}$_xerbla( 'ZGETRF2', -info ) return end if ! quick return if possible if( m==0 .or. n==0 )return if ( m==1_${ik}$ ) then ! use unblocked code for cone row case ! just need to handle ipiv and info ipiv( 1_${ik}$ ) = 1_${ik}$ if ( a(1_${ik}$,1_${ik}$)==czero )info = 1_${ik}$ else if( n==1_${ik}$ ) then ! use unblocked code for cone column case ! compute machine safe minimum sfmin = stdlib${ii}$_dlamch('S') ! find pivot and test for singularity i = stdlib${ii}$_izamax( m, a( 1_${ik}$, 1_${ik}$ ), 1_${ik}$ ) ipiv( 1_${ik}$ ) = i if( a( i, 1_${ik}$ )/=czero ) then ! apply the interchange if( i/=1_${ik}$ ) then temp = a( 1_${ik}$, 1_${ik}$ ) a( 1_${ik}$, 1_${ik}$ ) = a( i, 1_${ik}$ ) a( i, 1_${ik}$ ) = temp end if ! compute elements 2:m of the column if( abs(a( 1_${ik}$, 1_${ik}$ )) >= sfmin ) then call stdlib${ii}$_zscal( m-1, cone / a( 1_${ik}$, 1_${ik}$ ), a( 2_${ik}$, 1_${ik}$ ), 1_${ik}$ ) else do i = 1, m-1 a( 1_${ik}$+i, 1_${ik}$ ) = a( 1_${ik}$+i, 1_${ik}$ ) / a( 1_${ik}$, 1_${ik}$ ) end do end if else info = 1_${ik}$ end if else ! use recursive code n1 = min( m, n ) / 2_${ik}$ n2 = n-n1 ! [ a11 ] ! factor [ --- ] ! [ a21 ] call stdlib${ii}$_zgetrf2( m, n1, a, lda, ipiv, iinfo ) if ( info==0_${ik}$ .and. iinfo>0_${ik}$ )info = iinfo ! [ a12 ] ! apply interchanges to [ --- ] ! [ a22 ] call stdlib${ii}$_zlaswp( n2, a( 1_${ik}$, n1+1 ), lda, 1_${ik}$, n1, ipiv, 1_${ik}$ ) ! solve a12 call stdlib${ii}$_ztrsm( 'L', 'L', 'N', 'U', n1, n2, cone, a, lda,a( 1_${ik}$, n1+1 ), lda ) ! update a22 call stdlib${ii}$_zgemm( 'N', 'N', m-n1, n2, n1, -cone, a( n1+1, 1_${ik}$ ), lda,a( 1_${ik}$, n1+1 ), & lda, cone, a( n1+1, n1+1 ), lda ) ! factor a22 call stdlib${ii}$_zgetrf2( m-n1, n2, a( n1+1, n1+1 ), lda, ipiv( n1+1 ),iinfo ) ! adjust info and the pivot indices if ( info==0_${ik}$ .and. iinfo>0_${ik}$ )info = iinfo + n1 do i = n1+1, min( m, n ) ipiv( i ) = ipiv( i ) + n1 end do ! apply interchanges to a21 call stdlib${ii}$_zlaswp( n1, a( 1_${ik}$, 1_${ik}$ ), lda, n1+1, min( m, n), ipiv, 1_${ik}$ ) end if return end subroutine stdlib${ii}$_zgetrf2 #:for ck,ct,ci in CMPLX_KINDS_TYPES #:if not ck in ["sp","dp"] pure recursive module subroutine stdlib${ii}$_${ci}$getrf2( m, n, a, lda, ipiv, info ) !! ZGETRF2: computes an LU factorization of a general M-by-N matrix A !! using partial pivoting with row interchanges. !! The factorization has the form !! A = P * L * U !! where P is a permutation matrix, L is lower triangular with unit !! diagonal elements (lower trapezoidal if m > n), and U is upper !! triangular (upper trapezoidal if m < n). !! This is the recursive version of the algorithm. It divides !! the matrix into four submatrices: !! [ A11 | A12 ] where A11 is n1 by n1 and A22 is n2 by n2 !! A = [ -----|----- ] with n1 = min(m,n)/2 !! [ A21 | A22 ] n2 = n-n1 !! [ A11 ] !! The subroutine calls itself to factor [ --- ], !! [ A12 ] !! [ A12 ] !! do the swaps on [ --- ], solve A12, update A22, !! [ A22 ] !! then calls itself to factor A22 and do the swaps on A21. ! -- lapack computational 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: negone, zero, half, one, two, three, four, eight, ten, czero, chalf, cone, cnegone ! Scalar Arguments integer(${ik}$), intent(out) :: info integer(${ik}$), intent(in) :: lda, m, n ! Array Arguments integer(${ik}$), intent(out) :: ipiv(*) complex(${ck}$), intent(inout) :: a(lda,*) ! ===================================================================== ! Local Scalars real(${ck}$) :: sfmin complex(${ck}$) :: temp integer(${ik}$) :: i, iinfo, n1, n2 ! Intrinsic Functions ! Executable Statements ! test the input parameters info = 0_${ik}$ if( m<0_${ik}$ ) then info = -1_${ik}$ else if( n<0_${ik}$ ) then info = -2_${ik}$ else if( lda<max( 1_${ik}$, m ) ) then info = -4_${ik}$ end if if( info/=0_${ik}$ ) then call stdlib${ii}$_xerbla( 'ZGETRF2', -info ) return end if ! quick return if possible if( m==0 .or. n==0 )return if ( m==1_${ik}$ ) then ! use unblocked code for cone row case ! just need to handle ipiv and info ipiv( 1_${ik}$ ) = 1_${ik}$ if ( a(1_${ik}$,1_${ik}$)==czero )info = 1_${ik}$ else if( n==1_${ik}$ ) then ! use unblocked code for cone column case ! compute machine safe minimum sfmin = stdlib${ii}$_${c2ri(ci)}$lamch('S') ! find pivot and test for singularity i = stdlib${ii}$_i${ci}$amax( m, a( 1_${ik}$, 1_${ik}$ ), 1_${ik}$ ) ipiv( 1_${ik}$ ) = i if( a( i, 1_${ik}$ )/=czero ) then ! apply the interchange if( i/=1_${ik}$ ) then temp = a( 1_${ik}$, 1_${ik}$ ) a( 1_${ik}$, 1_${ik}$ ) = a( i, 1_${ik}$ ) a( i, 1_${ik}$ ) = temp end if ! compute elements 2:m of the column if( abs(a( 1_${ik}$, 1_${ik}$ )) >= sfmin ) then call stdlib${ii}$_${ci}$scal( m-1, cone / a( 1_${ik}$, 1_${ik}$ ), a( 2_${ik}$, 1_${ik}$ ), 1_${ik}$ ) else do i = 1, m-1 a( 1_${ik}$+i, 1_${ik}$ ) = a( 1_${ik}$+i, 1_${ik}$ ) / a( 1_${ik}$, 1_${ik}$ ) end do end if else info = 1_${ik}$ end if else ! use recursive code n1 = min( m, n ) / 2_${ik}$ n2 = n-n1 ! [ a11 ] ! factor [ --- ] ! [ a21 ] call stdlib${ii}$_${ci}$getrf2( m, n1, a, lda, ipiv, iinfo ) if ( info==0_${ik}$ .and. iinfo>0_${ik}$ )info = iinfo ! [ a12 ] ! apply interchanges to [ --- ] ! [ a22 ] call stdlib${ii}$_${ci}$laswp( n2, a( 1_${ik}$, n1+1 ), lda, 1_${ik}$, n1, ipiv, 1_${ik}$ ) ! solve a12 call stdlib${ii}$_${ci}$trsm( 'L', 'L', 'N', 'U', n1, n2, cone, a, lda,a( 1_${ik}$, n1+1 ), lda ) ! update a22 call stdlib${ii}$_${ci}$gemm( 'N', 'N', m-n1, n2, n1, -cone, a( n1+1, 1_${ik}$ ), lda,a( 1_${ik}$, n1+1 ), & lda, cone, a( n1+1, n1+1 ), lda ) ! factor a22 call stdlib${ii}$_${ci}$getrf2( m-n1, n2, a( n1+1, n1+1 ), lda, ipiv( n1+1 ),iinfo ) ! adjust info and the pivot indices if ( info==0_${ik}$ .and. iinfo>0_${ik}$ )info = iinfo + n1 do i = n1+1, min( m, n ) ipiv( i ) = ipiv( i ) + n1 end do ! apply interchanges to a21 call stdlib${ii}$_${ci}$laswp( n1, a( 1_${ik}$, 1_${ik}$ ), lda, n1+1, min( m, n), ipiv, 1_${ik}$ ) end if return end subroutine stdlib${ii}$_${ci}$getrf2 #:endif #:endfor pure module subroutine stdlib${ii}$_sgetf2( m, n, a, lda, ipiv, info ) !! SGETF2 computes an LU factorization of a general m-by-n matrix A !! using partial pivoting with row interchanges. !! The factorization has the form !! A = P * L * U !! where P is a permutation matrix, L is lower triangular with unit !! diagonal elements (lower trapezoidal if m > n), and U is upper !! triangular (upper trapezoidal if m < n). !! This is the right-looking Level 2 BLAS version of the algorithm. ! -- lapack computational 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: negone, zero, half, one, two, three, four, eight, ten, czero, chalf, cone, cnegone ! Scalar Arguments integer(${ik}$), intent(out) :: info integer(${ik}$), intent(in) :: lda, m, n ! Array Arguments integer(${ik}$), intent(out) :: ipiv(*) real(sp), intent(inout) :: a(lda,*) ! ===================================================================== ! Local Scalars real(sp) :: sfmin integer(${ik}$) :: i, j, jp ! Intrinsic Functions ! Executable Statements ! test the input parameters. info = 0_${ik}$ if( m<0_${ik}$ ) then info = -1_${ik}$ else if( n<0_${ik}$ ) then info = -2_${ik}$ else if( lda<max( 1_${ik}$, m ) ) then info = -4_${ik}$ end if if( info/=0_${ik}$ ) then call stdlib${ii}$_xerbla( 'SGETF2', -info ) return end if ! quick return if possible if( m==0 .or. n==0 )return ! compute machine safe minimum sfmin = stdlib${ii}$_slamch('S') do j = 1, min( m, n ) ! find pivot and test for singularity. jp = j - 1_${ik}$ + stdlib${ii}$_isamax( m-j+1, a( j, j ), 1_${ik}$ ) ipiv( j ) = jp if( a( jp, j )/=zero ) then ! apply the interchange to columns 1:n. if( jp/=j )call stdlib${ii}$_sswap( n, a( j, 1_${ik}$ ), lda, a( jp, 1_${ik}$ ), lda ) ! compute elements j+1:m of j-th column. if( j<m ) then if( abs(a( j, j )) >= sfmin ) then call stdlib${ii}$_sscal( m-j, one / a( j, j ), a( j+1, j ), 1_${ik}$ ) else do i = 1, m-j a( j+i, j ) = a( j+i, j ) / a( j, j ) end do end if end if else if( info==0_${ik}$ ) then info = j end if if( j<min( m, n ) ) then ! update trailing submatrix. call stdlib${ii}$_sger( m-j, n-j, -one, a( j+1, j ), 1_${ik}$, a( j, j+1 ), lda,a( j+1, j+1 ),& lda ) end if end do return end subroutine stdlib${ii}$_sgetf2 pure module subroutine stdlib${ii}$_dgetf2( m, n, a, lda, ipiv, info ) !! DGETF2 computes an LU factorization of a general m-by-n matrix A !! using partial pivoting with row interchanges. !! The factorization has the form !! A = P * L * U !! where P is a permutation matrix, L is lower triangular with unit !! diagonal elements (lower trapezoidal if m > n), and U is upper !! triangular (upper trapezoidal if m < n). !! This is the right-looking Level 2 BLAS version of the algorithm. ! -- lapack computational 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: negone, zero, half, one, two, three, four, eight, ten, czero, chalf, cone, cnegone ! Scalar Arguments integer(${ik}$), intent(out) :: info integer(${ik}$), intent(in) :: lda, m, n ! Array Arguments integer(${ik}$), intent(out) :: ipiv(*) real(dp), intent(inout) :: a(lda,*) ! ===================================================================== ! Local Scalars real(dp) :: sfmin integer(${ik}$) :: i, j, jp ! Intrinsic Functions ! Executable Statements ! test the input parameters. info = 0_${ik}$ if( m<0_${ik}$ ) then info = -1_${ik}$ else if( n<0_${ik}$ ) then info = -2_${ik}$ else if( lda<max( 1_${ik}$, m ) ) then info = -4_${ik}$ end if if( info/=0_${ik}$ ) then call stdlib${ii}$_xerbla( 'DGETF2', -info ) return end if ! quick return if possible if( m==0 .or. n==0 )return ! compute machine safe minimum sfmin = stdlib${ii}$_dlamch('S') do j = 1, min( m, n ) ! find pivot and test for singularity. jp = j - 1_${ik}$ + stdlib${ii}$_idamax( m-j+1, a( j, j ), 1_${ik}$ ) ipiv( j ) = jp if( a( jp, j )/=zero ) then ! apply the interchange to columns 1:n. if( jp/=j )call stdlib${ii}$_dswap( n, a( j, 1_${ik}$ ), lda, a( jp, 1_${ik}$ ), lda ) ! compute elements j+1:m of j-th column. if( j<m ) then if( abs(a( j, j )) >= sfmin ) then call stdlib${ii}$_dscal( m-j, one / a( j, j ), a( j+1, j ), 1_${ik}$ ) else do i = 1, m-j a( j+i, j ) = a( j+i, j ) / a( j, j ) end do end if end if else if( info==0_${ik}$ ) then info = j end if if( j<min( m, n ) ) then ! update trailing submatrix. call stdlib${ii}$_dger( m-j, n-j, -one, a( j+1, j ), 1_${ik}$, a( j, j+1 ), lda,a( j+1, j+1 ),& lda ) end if end do return end subroutine stdlib${ii}$_dgetf2 #:for rk,rt,ri in REAL_KINDS_TYPES #:if not rk in ["sp","dp"] pure module subroutine stdlib${ii}$_${ri}$getf2( m, n, a, lda, ipiv, info ) !! DGETF2: computes an LU factorization of a general m-by-n matrix A !! using partial pivoting with row interchanges. !! The factorization has the form !! A = P * L * U !! where P is a permutation matrix, L is lower triangular with unit !! diagonal elements (lower trapezoidal if m > n), and U is upper !! triangular (upper trapezoidal if m < n). !! This is the right-looking Level 2 BLAS version of the algorithm. ! -- lapack computational 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: negone, zero, half, one, two, three, four, eight, ten, czero, chalf, cone, cnegone ! Scalar Arguments integer(${ik}$), intent(out) :: info integer(${ik}$), intent(in) :: lda, m, n ! Array Arguments integer(${ik}$), intent(out) :: ipiv(*) real(${rk}$), intent(inout) :: a(lda,*) ! ===================================================================== ! Local Scalars real(${rk}$) :: sfmin integer(${ik}$) :: i, j, jp ! Intrinsic Functions ! Executable Statements ! test the input parameters. info = 0_${ik}$ if( m<0_${ik}$ ) then info = -1_${ik}$ else if( n<0_${ik}$ ) then info = -2_${ik}$ else if( lda<max( 1_${ik}$, m ) ) then info = -4_${ik}$ end if if( info/=0_${ik}$ ) then call stdlib${ii}$_xerbla( 'DGETF2', -info ) return end if ! quick return if possible if( m==0 .or. n==0 )return ! compute machine safe minimum sfmin = stdlib${ii}$_${ri}$lamch('S') do j = 1, min( m, n ) ! find pivot and test for singularity. jp = j - 1_${ik}$ + stdlib${ii}$_i${ri}$amax( m-j+1, a( j, j ), 1_${ik}$ ) ipiv( j ) = jp if( a( jp, j )/=zero ) then ! apply the interchange to columns 1:n. if( jp/=j )call stdlib${ii}$_${ri}$swap( n, a( j, 1_${ik}$ ), lda, a( jp, 1_${ik}$ ), lda ) ! compute elements j+1:m of j-th column. if( j<m ) then if( abs(a( j, j )) >= sfmin ) then call stdlib${ii}$_${ri}$scal( m-j, one / a( j, j ), a( j+1, j ), 1_${ik}$ ) else do i = 1, m-j a( j+i, j ) = a( j+i, j ) / a( j, j ) end do end if end if else if( info==0_${ik}$ ) then info = j end if if( j<min( m, n ) ) then ! update trailing submatrix. call stdlib${ii}$_${ri}$ger( m-j, n-j, -one, a( j+1, j ), 1_${ik}$, a( j, j+1 ), lda,a( j+1, j+1 ),& lda ) end if end do return end subroutine stdlib${ii}$_${ri}$getf2 #:endif #:endfor pure module subroutine stdlib${ii}$_cgetf2( m, n, a, lda, ipiv, info ) !! CGETF2 computes an LU factorization of a general m-by-n matrix A !! using partial pivoting with row interchanges. !! The factorization has the form !! A = P * L * U !! where P is a permutation matrix, L is lower triangular with unit !! diagonal elements (lower trapezoidal if m > n), and U is upper !! triangular (upper trapezoidal if m < n). !! This is the right-looking Level 2 BLAS version of the algorithm. ! -- lapack computational 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: negone, zero, half, one, two, three, four, eight, ten, czero, chalf, cone, cnegone ! Scalar Arguments integer(${ik}$), intent(out) :: info integer(${ik}$), intent(in) :: lda, m, n ! Array Arguments integer(${ik}$), intent(out) :: ipiv(*) complex(sp), intent(inout) :: a(lda,*) ! ===================================================================== ! Local Scalars real(sp) :: sfmin integer(${ik}$) :: i, j, jp ! Intrinsic Functions ! Executable Statements ! test the input parameters. info = 0_${ik}$ if( m<0_${ik}$ ) then info = -1_${ik}$ else if( n<0_${ik}$ ) then info = -2_${ik}$ else if( lda<max( 1_${ik}$, m ) ) then info = -4_${ik}$ end if if( info/=0_${ik}$ ) then call stdlib${ii}$_xerbla( 'CGETF2', -info ) return end if ! quick return if possible if( m==0 .or. n==0 )return ! compute machine safe minimum sfmin = stdlib${ii}$_slamch('S') do j = 1, min( m, n ) ! find pivot and test for singularity. jp = j - 1_${ik}$ + stdlib${ii}$_icamax( m-j+1, a( j, j ), 1_${ik}$ ) ipiv( j ) = jp if( a( jp, j )/=czero ) then ! apply the interchange to columns 1:n. if( jp/=j )call stdlib${ii}$_cswap( n, a( j, 1_${ik}$ ), lda, a( jp, 1_${ik}$ ), lda ) ! compute elements j+1:m of j-th column. if( j<m ) then if( abs(a( j, j )) >= sfmin ) then call stdlib${ii}$_cscal( m-j, cone / a( j, j ), a( j+1, j ), 1_${ik}$ ) else do i = 1, m-j a( j+i, j ) = a( j+i, j ) / a( j, j ) end do end if end if else if( info==0_${ik}$ ) then info = j end if if( j<min( m, n ) ) then ! update trailing submatrix. call stdlib${ii}$_cgeru( m-j, n-j, -cone, a( j+1, j ), 1_${ik}$, a( j, j+1 ),lda, a( j+1, j+1 & ), lda ) end if end do return end subroutine stdlib${ii}$_cgetf2 pure module subroutine stdlib${ii}$_zgetf2( m, n, a, lda, ipiv, info ) !! ZGETF2 computes an LU factorization of a general m-by-n matrix A !! using partial pivoting with row interchanges. !! The factorization has the form !! A = P * L * U !! where P is a permutation matrix, L is lower triangular with unit !! diagonal elements (lower trapezoidal if m > n), and U is upper !! triangular (upper trapezoidal if m < n). !! This is the right-looking Level 2 BLAS version of the algorithm. ! -- lapack computational 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: negone, zero, half, one, two, three, four, eight, ten, czero, chalf, cone, cnegone ! Scalar Arguments integer(${ik}$), intent(out) :: info integer(${ik}$), intent(in) :: lda, m, n ! Array Arguments integer(${ik}$), intent(out) :: ipiv(*) complex(dp), intent(inout) :: a(lda,*) ! ===================================================================== ! Local Scalars real(dp) :: sfmin integer(${ik}$) :: i, j, jp ! Intrinsic Functions ! Executable Statements ! test the input parameters. info = 0_${ik}$ if( m<0_${ik}$ ) then info = -1_${ik}$ else if( n<0_${ik}$ ) then info = -2_${ik}$ else if( lda<max( 1_${ik}$, m ) ) then info = -4_${ik}$ end if if( info/=0_${ik}$ ) then call stdlib${ii}$_xerbla( 'ZGETF2', -info ) return end if ! quick return if possible if( m==0 .or. n==0 )return ! compute machine safe minimum sfmin = stdlib${ii}$_dlamch('S') do j = 1, min( m, n ) ! find pivot and test for singularity. jp = j - 1_${ik}$ + stdlib${ii}$_izamax( m-j+1, a( j, j ), 1_${ik}$ ) ipiv( j ) = jp if( a( jp, j )/=czero ) then ! apply the interchange to columns 1:n. if( jp/=j )call stdlib${ii}$_zswap( n, a( j, 1_${ik}$ ), lda, a( jp, 1_${ik}$ ), lda ) ! compute elements j+1:m of j-th column. if( j<m ) then if( abs(a( j, j )) >= sfmin ) then call stdlib${ii}$_zscal( m-j, cone / a( j, j ), a( j+1, j ), 1_${ik}$ ) else do i = 1, m-j a( j+i, j ) = a( j+i, j ) / a( j, j ) end do end if end if else if( info==0_${ik}$ ) then info = j end if if( j<min( m, n ) ) then ! update trailing submatrix. call stdlib${ii}$_zgeru( m-j, n-j, -cone, a( j+1, j ), 1_${ik}$, a( j, j+1 ),lda, a( j+1, j+1 & ), lda ) end if end do return end subroutine stdlib${ii}$_zgetf2 #:for ck,ct,ci in CMPLX_KINDS_TYPES #:if not ck in ["sp","dp"] pure module subroutine stdlib${ii}$_${ci}$getf2( m, n, a, lda, ipiv, info ) !! ZGETF2: computes an LU factorization of a general m-by-n matrix A !! using partial pivoting with row interchanges. !! The factorization has the form !! A = P * L * U !! where P is a permutation matrix, L is lower triangular with unit !! diagonal elements (lower trapezoidal if m > n), and U is upper !! triangular (upper trapezoidal if m < n). !! This is the right-looking Level 2 BLAS version of the algorithm. ! -- lapack computational 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: negone, zero, half, one, two, three, four, eight, ten, czero, chalf, cone, cnegone ! Scalar Arguments integer(${ik}$), intent(out) :: info integer(${ik}$), intent(in) :: lda, m, n ! Array Arguments integer(${ik}$), intent(out) :: ipiv(*) complex(${ck}$), intent(inout) :: a(lda,*) ! ===================================================================== ! Local Scalars real(${ck}$) :: sfmin integer(${ik}$) :: i, j, jp ! Intrinsic Functions ! Executable Statements ! test the input parameters. info = 0_${ik}$ if( m<0_${ik}$ ) then info = -1_${ik}$ else if( n<0_${ik}$ ) then info = -2_${ik}$ else if( lda<max( 1_${ik}$, m ) ) then info = -4_${ik}$ end if if( info/=0_${ik}$ ) then call stdlib${ii}$_xerbla( 'ZGETF2', -info ) return end if ! quick return if possible if( m==0 .or. n==0 )return ! compute machine safe minimum sfmin = stdlib${ii}$_${c2ri(ci)}$lamch('S') do j = 1, min( m, n ) ! find pivot and test for singularity. jp = j - 1_${ik}$ + stdlib${ii}$_i${ci}$amax( m-j+1, a( j, j ), 1_${ik}$ ) ipiv( j ) = jp if( a( jp, j )/=czero ) then ! apply the interchange to columns 1:n. if( jp/=j )call stdlib${ii}$_${ci}$swap( n, a( j, 1_${ik}$ ), lda, a( jp, 1_${ik}$ ), lda ) ! compute elements j+1:m of j-th column. if( j<m ) then if( abs(a( j, j )) >= sfmin ) then call stdlib${ii}$_${ci}$scal( m-j, cone / a( j, j ), a( j+1, j ), 1_${ik}$ ) else do i = 1, m-j a( j+i, j ) = a( j+i, j ) / a( j, j ) end do end if end if else if( info==0_${ik}$ ) then info = j end if if( j<min( m, n ) ) then ! update trailing submatrix. call stdlib${ii}$_${ci}$geru( m-j, n-j, -cone, a( j+1, j ), 1_${ik}$, a( j, j+1 ),lda, a( j+1, j+1 & ), lda ) end if end do return end subroutine stdlib${ii}$_${ci}$getf2 #:endif #:endfor pure module subroutine stdlib${ii}$_sgetrs( trans, n, nrhs, a, lda, ipiv, b, ldb, info ) !! SGETRS solves a system of linear equations !! A * X = B or A**T * X = B !! with a general N-by-N matrix A using the LU factorization computed !! by SGETRF. ! -- lapack computational 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: negone, zero, half, one, two, three, four, eight, ten, czero, chalf, cone, cnegone ! Scalar Arguments character, intent(in) :: trans integer(${ik}$), intent(out) :: info integer(${ik}$), intent(in) :: lda, ldb, n, nrhs ! Array Arguments integer(${ik}$), intent(in) :: ipiv(*) real(sp), intent(in) :: a(lda,*) real(sp), intent(inout) :: b(ldb,*) ! ===================================================================== ! Local Scalars logical(lk) :: notran ! Intrinsic Functions ! Executable Statements ! test the input parameters. info = 0_${ik}$ notran = stdlib_lsame( trans, 'N' ) if( .not.notran .and. .not.stdlib_lsame( trans, 'T' ) .and. .not.stdlib_lsame( trans, & 'C' ) ) then info = -1_${ik}$ else if( n<0_${ik}$ ) then info = -2_${ik}$ else if( nrhs<0_${ik}$ ) then info = -3_${ik}$ else if( lda<max( 1_${ik}$, n ) ) then info = -5_${ik}$ else if( ldb<max( 1_${ik}$, n ) ) then info = -8_${ik}$ end if if( info/=0_${ik}$ ) then call stdlib${ii}$_xerbla( 'SGETRS', -info ) return end if ! quick return if possible if( n==0 .or. nrhs==0 )return if( notran ) then ! solve a * x = b. ! apply row interchanges to the right hand sides. call stdlib${ii}$_slaswp( nrhs, b, ldb, 1_${ik}$, n, ipiv, 1_${ik}$ ) ! solve l*x = b, overwriting b with x. call stdlib${ii}$_strsm( 'LEFT', 'LOWER', 'NO TRANSPOSE', 'UNIT', n, nrhs,one, a, lda, b, & ldb ) ! solve u*x = b, overwriting b with x. call stdlib${ii}$_strsm( 'LEFT', 'UPPER', 'NO TRANSPOSE', 'NON-UNIT', n,nrhs, one, a, lda,& b, ldb ) else ! solve a**t * x = b. ! solve u**t *x = b, overwriting b with x. call stdlib${ii}$_strsm( 'LEFT', 'UPPER', 'TRANSPOSE', 'NON-UNIT', n, nrhs,one, a, lda, b,& ldb ) ! solve l**t *x = b, overwriting b with x. call stdlib${ii}$_strsm( 'LEFT', 'LOWER', 'TRANSPOSE', 'UNIT', n, nrhs, one,a, lda, b, & ldb ) ! apply row interchanges to the solution vectors. call stdlib${ii}$_slaswp( nrhs, b, ldb, 1_${ik}$, n, ipiv, -1_${ik}$ ) end if return end subroutine stdlib${ii}$_sgetrs pure module subroutine stdlib${ii}$_dgetrs( trans, n, nrhs, a, lda, ipiv, b, ldb, info ) !! DGETRS solves a system of linear equations !! A * X = B or A**T * X = B !! with a general N-by-N matrix A using the LU factorization computed !! by DGETRF. ! -- lapack computational 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: negone, zero, half, one, two, three, four, eight, ten, czero, chalf, cone, cnegone ! Scalar Arguments character, intent(in) :: trans integer(${ik}$), intent(out) :: info integer(${ik}$), intent(in) :: lda, ldb, n, nrhs ! Array Arguments integer(${ik}$), intent(in) :: ipiv(*) real(dp), intent(in) :: a(lda,*) real(dp), intent(inout) :: b(ldb,*) ! ===================================================================== ! Local Scalars logical(lk) :: notran ! Intrinsic Functions ! Executable Statements ! test the input parameters. info = 0_${ik}$ notran = stdlib_lsame( trans, 'N' ) if( .not.notran .and. .not.stdlib_lsame( trans, 'T' ) .and. .not.stdlib_lsame( trans, & 'C' ) ) then info = -1_${ik}$ else if( n<0_${ik}$ ) then info = -2_${ik}$ else if( nrhs<0_${ik}$ ) then info = -3_${ik}$ else if( lda<max( 1_${ik}$, n ) ) then info = -5_${ik}$ else if( ldb<max( 1_${ik}$, n ) ) then info = -8_${ik}$ end if if( info/=0_${ik}$ ) then call stdlib${ii}$_xerbla( 'DGETRS', -info ) return end if ! quick return if possible if( n==0 .or. nrhs==0 )return if( notran ) then ! solve a * x = b. ! apply row interchanges to the right hand sides. call stdlib${ii}$_dlaswp( nrhs, b, ldb, 1_${ik}$, n, ipiv, 1_${ik}$ ) ! solve l*x = b, overwriting b with x. call stdlib${ii}$_dtrsm( 'LEFT', 'LOWER', 'NO TRANSPOSE', 'UNIT', n, nrhs,one, a, lda, b, & ldb ) ! solve u*x = b, overwriting b with x. call stdlib${ii}$_dtrsm( 'LEFT', 'UPPER', 'NO TRANSPOSE', 'NON-UNIT', n,nrhs, one, a, lda,& b, ldb ) else ! solve a**t * x = b. ! solve u**t *x = b, overwriting b with x. call stdlib${ii}$_dtrsm( 'LEFT', 'UPPER', 'TRANSPOSE', 'NON-UNIT', n, nrhs,one, a, lda, b,& ldb ) ! solve l**t *x = b, overwriting b with x. call stdlib${ii}$_dtrsm( 'LEFT', 'LOWER', 'TRANSPOSE', 'UNIT', n, nrhs, one,a, lda, b, & ldb ) ! apply row interchanges to the solution vectors. call stdlib${ii}$_dlaswp( nrhs, b, ldb, 1_${ik}$, n, ipiv, -1_${ik}$ ) end if return end subroutine stdlib${ii}$_dgetrs #:for rk,rt,ri in REAL_KINDS_TYPES #:if not rk in ["sp","dp"] pure module subroutine stdlib${ii}$_${ri}$getrs( trans, n, nrhs, a, lda, ipiv, b, ldb, info ) !! DGETRS: solves a system of linear equations !! A * X = B or A**T * X = B !! with a general N-by-N matrix A using the LU factorization computed !! by DGETRF. ! -- lapack computational 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: negone, zero, half, one, two, three, four, eight, ten, czero, chalf, cone, cnegone ! Scalar Arguments character, intent(in) :: trans integer(${ik}$), intent(out) :: info integer(${ik}$), intent(in) :: lda, ldb, n, nrhs ! Array Arguments integer(${ik}$), intent(in) :: ipiv(*) real(${rk}$), intent(in) :: a(lda,*) real(${rk}$), intent(inout) :: b(ldb,*) ! ===================================================================== ! Local Scalars logical(lk) :: notran ! Intrinsic Functions ! Executable Statements ! test the input parameters. info = 0_${ik}$ notran = stdlib_lsame( trans, 'N' ) if( .not.notran .and. .not.stdlib_lsame( trans, 'T' ) .and. .not.stdlib_lsame( trans, & 'C' ) ) then info = -1_${ik}$ else if( n<0_${ik}$ ) then info = -2_${ik}$ else if( nrhs<0_${ik}$ ) then info = -3_${ik}$ else if( lda<max( 1_${ik}$, n ) ) then info = -5_${ik}$ else if( ldb<max( 1_${ik}$, n ) ) then info = -8_${ik}$ end if if( info/=0_${ik}$ ) then call stdlib${ii}$_xerbla( 'DGETRS', -info ) return end if ! quick return if possible if( n==0 .or. nrhs==0 )return if( notran ) then ! solve a * x = b. ! apply row interchanges to the right hand sides. call stdlib${ii}$_${ri}$laswp( nrhs, b, ldb, 1_${ik}$, n, ipiv, 1_${ik}$ ) ! solve l*x = b, overwriting b with x. call stdlib${ii}$_${ri}$trsm( 'LEFT', 'LOWER', 'NO TRANSPOSE', 'UNIT', n, nrhs,one, a, lda, b, & ldb ) ! solve u*x = b, overwriting b with x. call stdlib${ii}$_${ri}$trsm( 'LEFT', 'UPPER', 'NO TRANSPOSE', 'NON-UNIT', n,nrhs, one, a, lda,& b, ldb ) else ! solve a**t * x = b. ! solve u**t *x = b, overwriting b with x. call stdlib${ii}$_${ri}$trsm( 'LEFT', 'UPPER', 'TRANSPOSE', 'NON-UNIT', n, nrhs,one, a, lda, b,& ldb ) ! solve l**t *x = b, overwriting b with x. call stdlib${ii}$_${ri}$trsm( 'LEFT', 'LOWER', 'TRANSPOSE', 'UNIT', n, nrhs, one,a, lda, b, & ldb ) ! apply row interchanges to the solution vectors. call stdlib${ii}$_${ri}$laswp( nrhs, b, ldb, 1_${ik}$, n, ipiv, -1_${ik}$ ) end if return end subroutine stdlib${ii}$_${ri}$getrs #:endif #:endfor pure module subroutine stdlib${ii}$_cgetrs( trans, n, nrhs, a, lda, ipiv, b, ldb, info ) !! CGETRS solves a system of linear equations !! A * X = B, A**T * X = B, or A**H * X = B !! with a general N-by-N matrix A using the LU factorization computed !! by CGETRF. ! -- lapack computational 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: negone, zero, half, one, two, three, four, eight, ten, czero, chalf, cone, cnegone ! Scalar Arguments character, intent(in) :: trans integer(${ik}$), intent(out) :: info integer(${ik}$), intent(in) :: lda, ldb, n, nrhs ! Array Arguments integer(${ik}$), intent(in) :: ipiv(*) complex(sp), intent(in) :: a(lda,*) complex(sp), intent(inout) :: b(ldb,*) ! ===================================================================== ! Local Scalars logical(lk) :: notran ! Intrinsic Functions ! Executable Statements ! test the input parameters. info = 0_${ik}$ notran = stdlib_lsame( trans, 'N' ) if( .not.notran .and. .not.stdlib_lsame( trans, 'T' ) .and. .not.stdlib_lsame( trans, & 'C' ) ) then info = -1_${ik}$ else if( n<0_${ik}$ ) then info = -2_${ik}$ else if( nrhs<0_${ik}$ ) then info = -3_${ik}$ else if( lda<max( 1_${ik}$, n ) ) then info = -5_${ik}$ else if( ldb<max( 1_${ik}$, n ) ) then info = -8_${ik}$ end if if( info/=0_${ik}$ ) then call stdlib${ii}$_xerbla( 'CGETRS', -info ) return end if ! quick return if possible if( n==0 .or. nrhs==0 )return if( notran ) then ! solve a * x = b. ! apply row interchanges to the right hand sides. call stdlib${ii}$_claswp( nrhs, b, ldb, 1_${ik}$, n, ipiv, 1_${ik}$ ) ! solve l*x = b, overwriting b with x. call stdlib${ii}$_ctrsm( 'LEFT', 'LOWER', 'NO TRANSPOSE', 'UNIT', n, nrhs,cone, a, lda, b,& ldb ) ! solve u*x = b, overwriting b with x. call stdlib${ii}$_ctrsm( 'LEFT', 'UPPER', 'NO TRANSPOSE', 'NON-UNIT', n,nrhs, cone, a, & lda, b, ldb ) else ! solve a**t * x = b or a**h * x = b. ! solve u**t *x = b or u**h *x = b, overwriting b with x. call stdlib${ii}$_ctrsm( 'LEFT', 'UPPER', trans, 'NON-UNIT', n, nrhs, cone,a, lda, b, ldb & ) ! solve l**t *x = b, or l**h *x = b overwriting b with x. call stdlib${ii}$_ctrsm( 'LEFT', 'LOWER', trans, 'UNIT', n, nrhs, cone, a,lda, b, ldb ) ! apply row interchanges to the solution vectors. call stdlib${ii}$_claswp( nrhs, b, ldb, 1_${ik}$, n, ipiv, -1_${ik}$ ) end if return end subroutine stdlib${ii}$_cgetrs pure module subroutine stdlib${ii}$_zgetrs( trans, n, nrhs, a, lda, ipiv, b, ldb, info ) !! ZGETRS solves a system of linear equations !! A * X = B, A**T * X = B, or A**H * X = B !! with a general N-by-N matrix A using the LU factorization computed !! by ZGETRF. ! -- lapack computational 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: negone, zero, half, one, two, three, four, eight, ten, czero, chalf, cone, cnegone ! Scalar Arguments character, intent(in) :: trans integer(${ik}$), intent(out) :: info integer(${ik}$), intent(in) :: lda, ldb, n, nrhs ! Array Arguments integer(${ik}$), intent(in) :: ipiv(*) complex(dp), intent(in) :: a(lda,*) complex(dp), intent(inout) :: b(ldb,*) ! ===================================================================== ! Local Scalars logical(lk) :: notran ! Intrinsic Functions ! Executable Statements ! test the input parameters. info = 0_${ik}$ notran = stdlib_lsame( trans, 'N' ) if( .not.notran .and. .not.stdlib_lsame( trans, 'T' ) .and. .not.stdlib_lsame( trans, & 'C' ) ) then info = -1_${ik}$ else if( n<0_${ik}$ ) then info = -2_${ik}$ else if( nrhs<0_${ik}$ ) then info = -3_${ik}$ else if( lda<max( 1_${ik}$, n ) ) then info = -5_${ik}$ else if( ldb<max( 1_${ik}$, n ) ) then info = -8_${ik}$ end if if( info/=0_${ik}$ ) then call stdlib${ii}$_xerbla( 'ZGETRS', -info ) return end if ! quick return if possible if( n==0 .or. nrhs==0 )return if( notran ) then ! solve a * x = b. ! apply row interchanges to the right hand sides. call stdlib${ii}$_zlaswp( nrhs, b, ldb, 1_${ik}$, n, ipiv, 1_${ik}$ ) ! solve l*x = b, overwriting b with x. call stdlib${ii}$_ztrsm( 'LEFT', 'LOWER', 'NO TRANSPOSE', 'UNIT', n, nrhs,cone, a, lda, b,& ldb ) ! solve u*x = b, overwriting b with x. call stdlib${ii}$_ztrsm( 'LEFT', 'UPPER', 'NO TRANSPOSE', 'NON-UNIT', n,nrhs, cone, a, & lda, b, ldb ) else ! solve a**t * x = b or a**h * x = b. ! solve u**t *x = b or u**h *x = b, overwriting b with x. call stdlib${ii}$_ztrsm( 'LEFT', 'UPPER', trans, 'NON-UNIT', n, nrhs, cone,a, lda, b, ldb & ) ! solve l**t *x = b, or l**h *x = b overwriting b with x. call stdlib${ii}$_ztrsm( 'LEFT', 'LOWER', trans, 'UNIT', n, nrhs, cone, a,lda, b, ldb ) ! apply row interchanges to the solution vectors. call stdlib${ii}$_zlaswp( nrhs, b, ldb, 1_${ik}$, n, ipiv, -1_${ik}$ ) end if return end subroutine stdlib${ii}$_zgetrs #:for ck,ct,ci in CMPLX_KINDS_TYPES #:if not ck in ["sp","dp"] pure module subroutine stdlib${ii}$_${ci}$getrs( trans, n, nrhs, a, lda, ipiv, b, ldb, info ) !! ZGETRS: solves a system of linear equations !! A * X = B, A**T * X = B, or A**H * X = B !! with a general N-by-N matrix A using the LU factorization computed !! by ZGETRF. ! -- lapack computational 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: negone, zero, half, one, two, three, four, eight, ten, czero, chalf, cone, cnegone ! Scalar Arguments character, intent(in) :: trans integer(${ik}$), intent(out) :: info integer(${ik}$), intent(in) :: lda, ldb, n, nrhs ! Array Arguments integer(${ik}$), intent(in) :: ipiv(*) complex(${ck}$), intent(in) :: a(lda,*) complex(${ck}$), intent(inout) :: b(ldb,*) ! ===================================================================== ! Local Scalars logical(lk) :: notran ! Intrinsic Functions ! Executable Statements ! test the input parameters. info = 0_${ik}$ notran = stdlib_lsame( trans, 'N' ) if( .not.notran .and. .not.stdlib_lsame( trans, 'T' ) .and. .not.stdlib_lsame( trans, & 'C' ) ) then info = -1_${ik}$ else if( n<0_${ik}$ ) then info = -2_${ik}$ else if( nrhs<0_${ik}$ ) then info = -3_${ik}$ else if( lda<max( 1_${ik}$, n ) ) then info = -5_${ik}$ else if( ldb<max( 1_${ik}$, n ) ) then info = -8_${ik}$ end if if( info/=0_${ik}$ ) then call stdlib${ii}$_xerbla( 'ZGETRS', -info ) return end if ! quick return if possible if( n==0 .or. nrhs==0 )return if( notran ) then ! solve a * x = b. ! apply row interchanges to the right hand sides. call stdlib${ii}$_${ci}$laswp( nrhs, b, ldb, 1_${ik}$, n, ipiv, 1_${ik}$ ) ! solve l*x = b, overwriting b with x. call stdlib${ii}$_${ci}$trsm( 'LEFT', 'LOWER', 'NO TRANSPOSE', 'UNIT', n, nrhs,cone, a, lda, b,& ldb ) ! solve u*x = b, overwriting b with x. call stdlib${ii}$_${ci}$trsm( 'LEFT', 'UPPER', 'NO TRANSPOSE', 'NON-UNIT', n,nrhs, cone, a, & lda, b, ldb ) else ! solve a**t * x = b or a**h * x = b. ! solve u**t *x = b or u**h *x = b, overwriting b with x. call stdlib${ii}$_${ci}$trsm( 'LEFT', 'UPPER', trans, 'NON-UNIT', n, nrhs, cone,a, lda, b, ldb & ) ! solve l**t *x = b, or l**h *x = b overwriting b with x. call stdlib${ii}$_${ci}$trsm( 'LEFT', 'LOWER', trans, 'UNIT', n, nrhs, cone, a,lda, b, ldb ) ! apply row interchanges to the solution vectors. call stdlib${ii}$_${ci}$laswp( nrhs, b, ldb, 1_${ik}$, n, ipiv, -1_${ik}$ ) end if return end subroutine stdlib${ii}$_${ci}$getrs #:endif #:endfor pure module subroutine stdlib${ii}$_sgetri( n, a, lda, ipiv, work, lwork, info ) !! SGETRI computes the inverse of a matrix using the LU factorization !! computed by SGETRF. !! This method inverts U and then computes inv(A) by solving the system !! inv(A)*L = inv(U) for inv(A). ! -- lapack computational 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: negone, zero, half, one, two, three, four, eight, ten, czero, chalf, cone, cnegone ! Scalar Arguments integer(${ik}$), intent(out) :: info integer(${ik}$), intent(in) :: lda, lwork, n ! Array Arguments integer(${ik}$), intent(in) :: ipiv(*) real(sp), intent(inout) :: a(lda,*) real(sp), intent(out) :: work(*) ! ===================================================================== ! Local Scalars logical(lk) :: lquery integer(${ik}$) :: i, iws, j, jb, jj, jp, ldwork, lwkopt, nb, nbmin, nn ! Intrinsic Functions ! Executable Statements ! test the input parameters. info = 0_${ik}$ nb = stdlib${ii}$_ilaenv( 1_${ik}$, 'SGETRI', ' ', n, -1_${ik}$, -1_${ik}$, -1_${ik}$ ) lwkopt = n*nb work( 1_${ik}$ ) = lwkopt lquery = ( lwork==-1_${ik}$ ) if( n<0_${ik}$ ) then info = -1_${ik}$ else if( lda<max( 1_${ik}$, n ) ) then info = -3_${ik}$ else if( lwork<max( 1_${ik}$, n ) .and. .not.lquery ) then info = -6_${ik}$ end if if( info/=0_${ik}$ ) then call stdlib${ii}$_xerbla( 'SGETRI', -info ) return else if( lquery ) then return end if ! quick return if possible if( n==0 )return ! form inv(u). if info > 0 from stdlib${ii}$_strtri, then u is singular, ! and the inverse is not computed. call stdlib${ii}$_strtri( 'UPPER', 'NON-UNIT', n, a, lda, info ) if( info>0 )return nbmin = 2_${ik}$ ldwork = n if( nb>1_${ik}$ .and. nb<n ) then iws = max( ldwork*nb, 1_${ik}$ ) if( lwork<iws ) then nb = lwork / ldwork nbmin = max( 2_${ik}$, stdlib${ii}$_ilaenv( 2_${ik}$, 'SGETRI', ' ', n, -1_${ik}$, -1_${ik}$, -1_${ik}$ ) ) end if else iws = n end if ! solve the equation inv(a)*l = inv(u) for inv(a). if( nb<nbmin .or. nb>=n ) then ! use unblocked code. do j = n, 1, -1 ! copy current column of l to work and replace with zeros. do i = j + 1, n work( i ) = a( i, j ) a( i, j ) = zero end do ! compute current column of inv(a). if( j<n )call stdlib${ii}$_sgemv( 'NO TRANSPOSE', n, n-j, -one, a( 1_${ik}$, j+1 ),lda, work( & j+1 ), 1_${ik}$, one, a( 1_${ik}$, j ), 1_${ik}$ ) end do else ! use blocked code. nn = ( ( n-1 ) / nb )*nb + 1_${ik}$ do j = nn, 1, -nb jb = min( nb, n-j+1 ) ! copy current block column of l to work and replace with ! zeros. do jj = j, j + jb - 1 do i = jj + 1, n work( i+( jj-j )*ldwork ) = a( i, jj ) a( i, jj ) = zero end do end do ! compute current block column of inv(a). if( j+jb<=n )call stdlib${ii}$_sgemm( 'NO TRANSPOSE', 'NO TRANSPOSE', n, jb,n-j-jb+1, -& one, a( 1_${ik}$, j+jb ), lda,work( j+jb ), ldwork, one, a( 1_${ik}$, j ), lda ) call stdlib${ii}$_strsm( 'RIGHT', 'LOWER', 'NO TRANSPOSE', 'UNIT', n, jb,one, work( j )& , ldwork, a( 1_${ik}$, j ), lda ) end do end if ! apply column interchanges. do j = n - 1, 1, -1 jp = ipiv( j ) if( jp/=j )call stdlib${ii}$_sswap( n, a( 1_${ik}$, j ), 1_${ik}$, a( 1_${ik}$, jp ), 1_${ik}$ ) end do work( 1_${ik}$ ) = iws return end subroutine stdlib${ii}$_sgetri pure module subroutine stdlib${ii}$_dgetri( n, a, lda, ipiv, work, lwork, info ) !! DGETRI computes the inverse of a matrix using the LU factorization !! computed by DGETRF. !! This method inverts U and then computes inv(A) by solving the system !! inv(A)*L = inv(U) for inv(A). ! -- lapack computational 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: negone, zero, half, one, two, three, four, eight, ten, czero, chalf, cone, cnegone ! Scalar Arguments integer(${ik}$), intent(out) :: info integer(${ik}$), intent(in) :: lda, lwork, n ! Array Arguments integer(${ik}$), intent(in) :: ipiv(*) real(dp), intent(inout) :: a(lda,*) real(dp), intent(out) :: work(*) ! ===================================================================== ! Local Scalars logical(lk) :: lquery integer(${ik}$) :: i, iws, j, jb, jj, jp, ldwork, lwkopt, nb, nbmin, nn ! Intrinsic Functions ! Executable Statements ! test the input parameters. info = 0_${ik}$ nb = stdlib${ii}$_ilaenv( 1_${ik}$, 'DGETRI', ' ', n, -1_${ik}$, -1_${ik}$, -1_${ik}$ ) lwkopt = n*nb work( 1_${ik}$ ) = lwkopt lquery = ( lwork==-1_${ik}$ ) if( n<0_${ik}$ ) then info = -1_${ik}$ else if( lda<max( 1_${ik}$, n ) ) then info = -3_${ik}$ else if( lwork<max( 1_${ik}$, n ) .and. .not.lquery ) then info = -6_${ik}$ end if if( info/=0_${ik}$ ) then call stdlib${ii}$_xerbla( 'DGETRI', -info ) return else if( lquery ) then return end if ! quick return if possible if( n==0 )return ! form inv(u). if info > 0 from stdlib${ii}$_dtrtri, then u is singular, ! and the inverse is not computed. call stdlib${ii}$_dtrtri( 'UPPER', 'NON-UNIT', n, a, lda, info ) if( info>0 )return nbmin = 2_${ik}$ ldwork = n if( nb>1_${ik}$ .and. nb<n ) then iws = max( ldwork*nb, 1_${ik}$ ) if( lwork<iws ) then nb = lwork / ldwork nbmin = max( 2_${ik}$, stdlib${ii}$_ilaenv( 2_${ik}$, 'DGETRI', ' ', n, -1_${ik}$, -1_${ik}$, -1_${ik}$ ) ) end if else iws = n end if ! solve the equation inv(a)*l = inv(u) for inv(a). if( nb<nbmin .or. nb>=n ) then ! use unblocked code. do j = n, 1, -1 ! copy current column of l to work and replace with zeros. do i = j + 1, n work( i ) = a( i, j ) a( i, j ) = zero end do ! compute current column of inv(a). if( j<n )call stdlib${ii}$_dgemv( 'NO TRANSPOSE', n, n-j, -one, a( 1_${ik}$, j+1 ),lda, work( & j+1 ), 1_${ik}$, one, a( 1_${ik}$, j ), 1_${ik}$ ) end do else ! use blocked code. nn = ( ( n-1 ) / nb )*nb + 1_${ik}$ do j = nn, 1, -nb jb = min( nb, n-j+1 ) ! copy current block column of l to work and replace with ! zeros. do jj = j, j + jb - 1 do i = jj + 1, n work( i+( jj-j )*ldwork ) = a( i, jj ) a( i, jj ) = zero end do end do ! compute current block column of inv(a). if( j+jb<=n )call stdlib${ii}$_dgemm( 'NO TRANSPOSE', 'NO TRANSPOSE', n, jb,n-j-jb+1, -& one, a( 1_${ik}$, j+jb ), lda,work( j+jb ), ldwork, one, a( 1_${ik}$, j ), lda ) call stdlib${ii}$_dtrsm( 'RIGHT', 'LOWER', 'NO TRANSPOSE', 'UNIT', n, jb,one, work( j )& , ldwork, a( 1_${ik}$, j ), lda ) end do end if ! apply column interchanges. do j = n - 1, 1, -1 jp = ipiv( j ) if( jp/=j )call stdlib${ii}$_dswap( n, a( 1_${ik}$, j ), 1_${ik}$, a( 1_${ik}$, jp ), 1_${ik}$ ) end do work( 1_${ik}$ ) = iws return end subroutine stdlib${ii}$_dgetri #:for rk,rt,ri in REAL_KINDS_TYPES #:if not rk in ["sp","dp"] pure module subroutine stdlib${ii}$_${ri}$getri( n, a, lda, ipiv, work, lwork, info ) !! DGETRI: computes the inverse of a matrix using the LU factorization !! computed by DGETRF. !! This method inverts U and then computes inv(A) by solving the system !! inv(A)*L = inv(U) for inv(A). ! -- lapack computational 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: negone, zero, half, one, two, three, four, eight, ten, czero, chalf, cone, cnegone ! Scalar Arguments integer(${ik}$), intent(out) :: info integer(${ik}$), intent(in) :: lda, lwork, n ! Array Arguments integer(${ik}$), intent(in) :: ipiv(*) real(${rk}$), intent(inout) :: a(lda,*) real(${rk}$), intent(out) :: work(*) ! ===================================================================== ! Local Scalars logical(lk) :: lquery integer(${ik}$) :: i, iws, j, jb, jj, jp, ldwork, lwkopt, nb, nbmin, nn ! Intrinsic Functions ! Executable Statements ! test the input parameters. info = 0_${ik}$ nb = stdlib${ii}$_ilaenv( 1_${ik}$, 'DGETRI', ' ', n, -1_${ik}$, -1_${ik}$, -1_${ik}$ ) lwkopt = n*nb work( 1_${ik}$ ) = lwkopt lquery = ( lwork==-1_${ik}$ ) if( n<0_${ik}$ ) then info = -1_${ik}$ else if( lda<max( 1_${ik}$, n ) ) then info = -3_${ik}$ else if( lwork<max( 1_${ik}$, n ) .and. .not.lquery ) then info = -6_${ik}$ end if if( info/=0_${ik}$ ) then call stdlib${ii}$_xerbla( 'DGETRI', -info ) return else if( lquery ) then return end if ! quick return if possible if( n==0 )return ! form inv(u). if info > 0 from stdlib${ii}$_${ri}$trtri, then u is singular, ! and the inverse is not computed. call stdlib${ii}$_${ri}$trtri( 'UPPER', 'NON-UNIT', n, a, lda, info ) if( info>0 )return nbmin = 2_${ik}$ ldwork = n if( nb>1_${ik}$ .and. nb<n ) then iws = max( ldwork*nb, 1_${ik}$ ) if( lwork<iws ) then nb = lwork / ldwork nbmin = max( 2_${ik}$, stdlib${ii}$_ilaenv( 2_${ik}$, 'DGETRI', ' ', n, -1_${ik}$, -1_${ik}$, -1_${ik}$ ) ) end if else iws = n end if ! solve the equation inv(a)*l = inv(u) for inv(a). if( nb<nbmin .or. nb>=n ) then ! use unblocked code. do j = n, 1, -1 ! copy current column of l to work and replace with zeros. do i = j + 1,