stdlib_lapack_solve_lu_comp.fypp Source File


Source Code

#: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,