diff --git a/src/lapack/stdlib_linalg_lapack_aux.fypp b/src/lapack/stdlib_linalg_lapack_aux.fypp
index ba2ec8241..442999ca3 100644
--- a/src/lapack/stdlib_linalg_lapack_aux.fypp
+++ b/src/lapack/stdlib_linalg_lapack_aux.fypp
@@ -3,6 +3,8 @@
 module stdlib_linalg_lapack_aux
      use stdlib_linalg_constants
      use stdlib_linalg_blas
+     use stdlib_linalg_state, only: linalg_state_type, linalg_error_handling, LINALG_ERROR, &
+         LINALG_INTERNAL_ERROR, LINALG_VALUE_ERROR, LINALG_SUCCESS
      use ieee_arithmetic, only: ieee_support_inf, ieee_support_nan
      implicit none
      private
@@ -38,6 +40,17 @@ module stdlib_linalg_lapack_aux
      public :: stdlib_select_${ri}$
      public :: stdlib_selctg_${ri}$     
      #:endfor   
+     public :: handle_potrf_info
+     public :: handle_getri_info
+     public :: handle_gesdd_info
+     public :: handle_gesv_info
+     public :: handle_gees_info
+     public :: handle_geqrf_info
+     public :: handle_orgqr_info
+     public :: handle_gelsd_info
+     public :: handle_geev_info
+     public :: handle_ggev_info
+     public :: handle_heev_info
      
      ! SELCTG is a LOGICAL FUNCTION of three DOUBLE PRECISION arguments 
      ! used to select eigenvalues to sort to the top left of the Schur form. 
@@ -1275,4 +1288,323 @@ module stdlib_linalg_lapack_aux
 
      #:endfor
 
+   !----------------------------------------------------------------------------
+   !-----                                                                  -----
+   !-----     AUXILIARY INFO HANDLING FUNCTIONS FOR LAPACK SUBROUTINES     -----
+   !-----                                                                  -----
+   !----------------------------------------------------------------------------
+   
+   ! Cholesky factorization
+   elemental subroutine handle_potrf_info(this,info,triangle,lda,n,err)
+      character(len=*), intent(in) :: this
+      character, intent(in) :: triangle
+      integer(ilp), intent(in) :: info,lda,n
+      type(linalg_state_type), intent(out) :: err
+
+      ! Process output
+      select case (info)
+      case (0)
+         ! Success
+      case (-1)
+         err = linalg_state_type(this,LINALG_INTERNAL_ERROR,'invalid triangle selection: ', &
+                  triangle,'. should be U/L')
+      case (-2)
+         err = linalg_state_type(this,LINALG_VALUE_ERROR,'invalid matrix size n=',n)
+      case (-4)
+         err = linalg_state_type(this,LINALG_VALUE_ERROR,'invalid lda=',lda,': is < n = ',n)
+      case (1:)
+         err = linalg_state_type(this,LINALG_ERROR,'cannot complete factorization:',info, &
+               '-th order leading minor is not positive definite')
+      case default
+         err = linalg_state_type(this,LINALG_INTERNAL_ERROR,'catastrophic error')
+      end select
+
+   end subroutine handle_potrf_info
+
+   elemental subroutine handle_getri_info(this,info,lda,n,err)
+      character(len=*), intent(in) :: this
+      integer(ilp), intent(in) :: info,lda,n
+      type(linalg_state_type), intent(out) :: err
+
+      ! Process output
+      select case (info)
+         case (0)
+             ! Success
+         case (:-1)
+             err = linalg_state_type(this,LINALG_ERROR,'invalid matrix size a=',[lda,n])
+         case (1:)
+             ! Matrix is singular
+             err = linalg_state_type(this,LINALG_ERROR,'singular matrix')
+         case default
+             err = linalg_state_type(this,LINALG_INTERNAL_ERROR,'catastrophic error')
+      end select
+   end subroutine handle_getri_info
+
+   elemental subroutine handle_gesdd_info(this,err,info,m,n)
+        character(len=*), intent(in) :: this
+        !> Error handler
+        type(linalg_state_type), intent(inout) :: err
+        !> GESDD return flag
+        integer(ilp), intent(in) :: info
+        !> Input matrix size
+        integer(ilp), intent(in) :: m,n
+
+        select case (info)
+           case (0)
+               ! Success!
+               err%state = LINALG_SUCCESS
+           case (-1)
+               err = linalg_state_type(this,LINALG_INTERNAL_ERROR,'Invalid task ID on input to GESDD.')
+           case (-5,-3:-2)
+               err = linalg_state_type(this,LINALG_VALUE_ERROR,'invalid matrix size: a=',[m,n])
+           case (-8)
+               err = linalg_state_type(this,LINALG_VALUE_ERROR,'invalid matrix U size, with a=',[m,n])
+           case (-10)
+               err = linalg_state_type(this,LINALG_VALUE_ERROR,'invalid matrix V size, with a=',[m,n])
+           case (-4)
+               err = linalg_state_type(this,LINALG_VALUE_ERROR,'A contains invalid/NaN values.')
+           case (1:)
+               err = linalg_state_type(this,LINALG_ERROR,'SVD computation did not converge.')
+           case default
+               err = linalg_state_type(this,LINALG_INTERNAL_ERROR,'Unknown error returned by GESDD.')
+        end select
+
+     end subroutine handle_gesdd_info
+
+     elemental subroutine handle_gesv_info(this,info,lda,n,nrhs,err)
+         character(len=*), intent(in) :: this
+         integer(ilp), intent(in) :: info,lda,n,nrhs
+         type(linalg_state_type), intent(out) :: err
+
+         ! Process output
+         select case (info)
+            case (0)
+                ! Success
+            case (-1)
+                err = linalg_state_type(this,LINALG_VALUE_ERROR,'invalid problem size n=',n)
+            case (-2)
+                err = linalg_state_type(this,LINALG_VALUE_ERROR,'invalid rhs size n=',nrhs)
+            case (-4)
+                err = linalg_state_type(this,LINALG_VALUE_ERROR,'invalid matrix size a=',[lda,n])
+            case (-7)
+                err = linalg_state_type(this,LINALG_ERROR,'invalid matrix size a=',[lda,n])
+            case (1:)
+                err = linalg_state_type(this,LINALG_ERROR,'singular matrix')
+            case default
+                err = linalg_state_type(this,LINALG_INTERNAL_ERROR,'catastrophic error')
+         end select
+         
+     end subroutine handle_gesv_info
+
+   !> Wrapper function to handle GEES error codes
+    elemental subroutine handle_gees_info(this, info, m, n, ldvs, err)
+        character(len=*), intent(in) :: this
+        integer(ilp), intent(in) :: info, m, n, ldvs
+        type(linalg_state_type), intent(out) :: err
+
+        ! Process GEES output
+        select case (info)
+        case (0_ilp)
+            ! Success
+        case (-1_ilp)
+            ! Vector not wanted, but task is wrong
+            err = linalg_state_type(this, LINALG_INTERNAL_ERROR,'Invalid Schur vector task request')
+        case (-2_ilp)
+            ! Vector not wanted, but task is wrong
+            err = linalg_state_type(this, LINALG_INTERNAL_ERROR,'Invalid sorting task request')            
+        case (-4_ilp,-6_ilp)
+            ! Vector not wanted, but task is wrong
+            err = linalg_state_type(this, LINALG_VALUE_ERROR,'Invalid/non-square input matrix size:',[m,n])  
+        case (-11_ilp)
+            err = linalg_state_type(this, LINALG_VALUE_ERROR,'Schur vector matrix has insufficient size',[ldvs,n])
+        case (-13_ilp)
+            err = linalg_state_type(this, LINALG_INTERNAL_ERROR,'Insufficient working storage size')
+        case (1_ilp:)
+            
+            if (info==n+2) then 
+                err = linalg_state_type(this, LINALG_ERROR, 'Ill-conditioned problem: could not sort eigenvalues')
+            elseif (info==n+1) then 
+                err = linalg_state_type(this, LINALG_ERROR, 'Some selected eigenvalues lost property due to sorting')                
+            elseif (info==n) then 
+                err = linalg_state_type(this, LINALG_ERROR, 'Convergence failure: no converged eigenvalues')
+            else
+                err = linalg_state_type(this, LINALG_ERROR, 'Convergence failure; converged range is',[info,n])                            
+            end if
+            
+        case default
+            
+            err = linalg_state_type(this, LINALG_INTERNAL_ERROR, 'GEES catastrophic error: info=', info)
+
+        end select
+        
+    end subroutine handle_gees_info
+
+   elemental subroutine handle_geqrf_info(this,info,m,n,lwork,err)
+         character(len=*), intent(in) :: this
+         integer(ilp), intent(in) :: info,m,n,lwork
+         type(linalg_state_type), intent(out) :: err
+
+         ! Process output
+         select case (info)
+            case (0)
+                ! Success
+            case (-1)
+                err = linalg_state_type(this,LINALG_VALUE_ERROR,'invalid matrix size m=',m)
+            case (-2)
+                err = linalg_state_type(this,LINALG_VALUE_ERROR,'invalid matrix size n=',n)
+            case (-4)
+                err = linalg_state_type(this,LINALG_VALUE_ERROR,'invalid matrix size a=',[m,n])
+            case (-7)
+                err = linalg_state_type(this,LINALG_ERROR,'invalid input for lwork=',lwork)
+            case default
+                err = linalg_state_type(this,LINALG_INTERNAL_ERROR,'catastrophic error')
+         end select
+
+     end subroutine handle_geqrf_info
+
+     elemental subroutine handle_orgqr_info(this,info,m,n,k,lwork,err)
+         character(len=*), intent(in) :: this
+         integer(ilp), intent(in) :: info,m,n,k,lwork
+         type(linalg_state_type), intent(out) :: err
+
+         ! Process output
+         select case (info)
+            case (0)
+                ! Success
+            case (-1)
+                err = linalg_state_type(this,LINALG_VALUE_ERROR,'invalid matrix size m=',m)
+            case (-2)
+                err = linalg_state_type(this,LINALG_VALUE_ERROR,'invalid matrix size n=',n)
+            case (-4)
+                err = linalg_state_type(this,LINALG_VALUE_ERROR,'invalid k=min(m,n)=',k)
+            case (-5)
+                err = linalg_state_type(this,LINALG_VALUE_ERROR,'invalid matrix size a=',[m,n])                
+            case (-8)
+                err = linalg_state_type(this,LINALG_ERROR,'invalid input for lwork=',lwork)
+            case default
+                err = linalg_state_type(this,LINALG_INTERNAL_ERROR,'catastrophic error')
+         end select
+
+     end subroutine handle_orgqr_info     
+
+     elemental subroutine handle_gelsd_info(this,info,lda,n,ldb,nrhs,err)
+          character(len=*), intent(in) :: this
+          integer(ilp), intent(in) :: info,lda,n,ldb,nrhs
+          type(linalg_state_type), intent(out) :: err
+
+          ! Process output
+          select case (info)
+             case (0)
+                 ! Success
+             case (:-1)
+                 err = linalg_state_type(this,LINALG_VALUE_ERROR,'invalid problem size a=',[lda,n], &
+                                                                                    ', b=',[ldb,nrhs])
+             case (1:)
+                 err = linalg_state_type(this,LINALG_ERROR,'SVD did not converge.')
+             case default
+                 err = linalg_state_type(this,LINALG_INTERNAL_ERROR,'catastrophic error')
+          end select
+
+      end subroutine handle_gelsd_info
+
+     !> Process GEEV output flags
+     pure subroutine handle_geev_info(this,err,info,shapea)
+        character(len=*), intent(in) :: this
+        !> Error handler
+        type(linalg_state_type), intent(inout) :: err
+        !> GEEV return flag
+        integer(ilp), intent(in) :: info
+        !> Input matrix size
+        integer(ilp), intent(in) :: shapea(2)
+
+        select case (info)
+           case (0)
+               ! Success!
+               err%state = LINALG_SUCCESS
+           case (-1)
+               err = linalg_state_type(this,LINALG_INTERNAL_ERROR,'Invalid task ID: left eigenvectors.')
+           case (-2)
+               err = linalg_state_type(this,LINALG_INTERNAL_ERROR,'Invalid task ID: right eigenvectors.')
+           case (-5,-3)
+               err = linalg_state_type(this,LINALG_VALUE_ERROR,'invalid matrix size: a=',shapea)
+           case (-9)
+               err = linalg_state_type(this,LINALG_VALUE_ERROR,'insufficient left vector matrix size.')
+           case (-11)
+               err = linalg_state_type(this,LINALG_VALUE_ERROR,'insufficient right vector matrix size.')
+           case (-13)
+               err = linalg_state_type(this,LINALG_INTERNAL_ERROR,'Insufficient work array size.')
+           case (1:)
+               err = linalg_state_type(this,LINALG_ERROR,'Eigenvalue computation did not converge.')
+           case default
+               err = linalg_state_type(this,LINALG_INTERNAL_ERROR,'Unknown error returned by geev.')
+        end select
+
+     end subroutine handle_geev_info
+
+     !> Process GGEV output flags
+     pure subroutine handle_ggev_info(this,err,info,shapea,shapeb)
+        character(len=*), intent(in) :: this
+        !> Error handler
+        type(linalg_state_type), intent(inout) :: err
+        !> GEEV return flag
+        integer(ilp), intent(in) :: info
+        !> Input matrix size
+        integer(ilp), intent(in) :: shapea(2),shapeb(2)
+
+        select case (info)
+           case (0)
+               ! Success!
+               err%state = LINALG_SUCCESS
+           case (-1)
+               err = linalg_state_type(this,LINALG_INTERNAL_ERROR,'Invalid task ID: left eigenvectors.')
+           case (-2)
+               err = linalg_state_type(this,LINALG_INTERNAL_ERROR,'Invalid task ID: right eigenvectors.')
+           case (-5,-3)
+               err = linalg_state_type(this,LINALG_VALUE_ERROR,'invalid matrix size: a=',shapea)
+           case (-7)
+               err = linalg_state_type(this,LINALG_VALUE_ERROR,'invalid matrix size: b=',shapeb)               
+           case (-12)
+               err = linalg_state_type(this,LINALG_VALUE_ERROR,'insufficient left vector matrix size.')
+           case (-14)
+               err = linalg_state_type(this,LINALG_VALUE_ERROR,'insufficient right vector matrix size.')
+           case (-16)
+               err = linalg_state_type(this,LINALG_INTERNAL_ERROR,'Insufficient work array size.')
+           case (1:)
+               err = linalg_state_type(this,LINALG_ERROR,'Eigenvalue computation did not converge.')
+           case default
+               err = linalg_state_type(this,LINALG_INTERNAL_ERROR,'Unknown error returned by ggev.')
+        end select
+
+     end subroutine handle_ggev_info
+
+     !> Process SYEV/HEEV output flags
+     elemental subroutine handle_heev_info(this,err,info,m,n)
+        character(len=*), intent(in) :: this
+        !> Error handler
+        type(linalg_state_type), intent(inout) :: err
+        !> SYEV/HEEV return flag
+        integer(ilp), intent(in) :: info
+        !> Input matrix size
+        integer(ilp), intent(in) :: m,n
+
+        select case (info)
+           case (0)
+               ! Success!
+               err%state = LINALG_SUCCESS
+           case (-1)
+               err = linalg_state_type(this,LINALG_INTERNAL_ERROR,'Invalid eigenvector request.')
+           case (-2)
+               err = linalg_state_type(this,LINALG_INTERNAL_ERROR,'Invalid triangular section request.')
+           case (-5,-3)
+               err = linalg_state_type(this,LINALG_VALUE_ERROR,'invalid matrix size: a=',[m,n])
+           case (-8)
+               err = linalg_state_type(this,LINALG_INTERNAL_ERROR,'insufficient workspace size.')
+           case (1:)
+               err = linalg_state_type(this,LINALG_ERROR,'Eigenvalue computation did not converge.')
+           case default
+               err = linalg_state_type(this,LINALG_INTERNAL_ERROR,'Unknown error returned by syev/heev.')
+        end select
+
+     end subroutine handle_heev_info
+
 end module stdlib_linalg_lapack_aux
diff --git a/src/stdlib_linalg_cholesky.fypp b/src/stdlib_linalg_cholesky.fypp
index 73d79845b..9a17c9815 100644
--- a/src/stdlib_linalg_cholesky.fypp
+++ b/src/stdlib_linalg_cholesky.fypp
@@ -4,6 +4,7 @@
 submodule (stdlib_linalg) stdlib_linalg_cholesky
      use stdlib_linalg_constants
      use stdlib_linalg_lapack, only: potrf
+     use stdlib_linalg_lapack_aux, only: handle_potrf_info
      use stdlib_linalg_state, only: linalg_state_type, linalg_error_handling, LINALG_ERROR, &
          LINALG_INTERNAL_ERROR, LINALG_VALUE_ERROR
      implicit none
@@ -13,31 +14,6 @@ submodule (stdlib_linalg) stdlib_linalg_cholesky
 
      contains
 
-     elemental subroutine handle_potrf_info(info,triangle,lda,n,err)
-         character, intent(in) :: triangle
-         integer(ilp), intent(in) :: info,lda,n
-         type(linalg_state_type), intent(out) :: err
-
-         ! Process output
-         select case (info)
-            case (0)
-                ! Success
-            case (-1)
-                err = linalg_state_type(this,LINALG_INTERNAL_ERROR,'invalid triangle selection: ', &
-                                   triangle,'. should be U/L')
-            case (-2)
-                err = linalg_state_type(this,LINALG_VALUE_ERROR,'invalid matrix size n=',n)
-            case (-4)
-                err = linalg_state_type(this,LINALG_VALUE_ERROR,'invalid lda=',lda,': is < n = ',n)
-            case (1:)
-                err = linalg_state_type(this,LINALG_ERROR,'cannot complete factorization:',info, &
-                                   '-th order leading minor is not positive definite')
-            case default
-                err = linalg_state_type(this,LINALG_INTERNAL_ERROR,'catastrophic error')
-         end select
-
-     end subroutine handle_potrf_info
-     
      #:for rk,rt,ri in RC_KINDS_TYPES
      
      ! Compute the Cholesky factorization of a symmetric / Hermitian matrix, A = L*L^T = U^T*U. 
@@ -84,7 +60,7 @@ submodule (stdlib_linalg) stdlib_linalg_cholesky
          
              ! Compute factorization
              call potrf(triangle,n,a,lda,info)
-             call handle_potrf_info(info,triangle,lda,n,err0)
+             call handle_potrf_info(this, info,triangle,lda,n,err0)
                       
              ! Zero-out the unused part of matrix A
              clean_unused: if (other_zeroed_ .and. err0%ok()) then    
diff --git a/src/stdlib_linalg_eigenvalues.fypp b/src/stdlib_linalg_eigenvalues.fypp
index 4123c9ef7..3a9bce178 100644
--- a/src/stdlib_linalg_eigenvalues.fypp
+++ b/src/stdlib_linalg_eigenvalues.fypp
@@ -7,6 +7,7 @@ submodule (stdlib_linalg) stdlib_linalg_eigenvalues
 !! Compute eigenvalues and eigenvectors    
      use stdlib_linalg_constants
      use stdlib_linalg_lapack, only: geev, ggev, heev, syev
+     use stdlib_linalg_lapack_aux, only: handle_geev_info, handle_ggev_info, handle_heev_info
      use stdlib_linalg_state, only: linalg_state_type, linalg_error_handling, LINALG_ERROR, &
           LINALG_INTERNAL_ERROR, LINALG_VALUE_ERROR, LINALG_SUCCESS     
      use, intrinsic:: ieee_arithmetic, only: ieee_value, ieee_positive_inf, ieee_quiet_nan
@@ -36,103 +37,6 @@ submodule (stdlib_linalg) stdlib_linalg_eigenvalues
         if (present(upper)) symmetric_triangle_task = merge('U','L',upper)
      end function symmetric_triangle_task
 
-     !> Process GEEV output flags
-     pure subroutine handle_geev_info(err,info,shapea)
-        !> Error handler
-        type(linalg_state_type), intent(inout) :: err
-        !> GEEV return flag
-        integer(ilp), intent(in) :: info
-        !> Input matrix size
-        integer(ilp), intent(in) :: shapea(2)
-
-        select case (info)
-           case (0)
-               ! Success!
-               err%state = LINALG_SUCCESS
-           case (-1)
-               err = linalg_state_type(this,LINALG_INTERNAL_ERROR,'Invalid task ID: left eigenvectors.')
-           case (-2)
-               err = linalg_state_type(this,LINALG_INTERNAL_ERROR,'Invalid task ID: right eigenvectors.')
-           case (-5,-3)
-               err = linalg_state_type(this,LINALG_VALUE_ERROR,'invalid matrix size: a=',shapea)
-           case (-9)
-               err = linalg_state_type(this,LINALG_VALUE_ERROR,'insufficient left vector matrix size.')
-           case (-11)
-               err = linalg_state_type(this,LINALG_VALUE_ERROR,'insufficient right vector matrix size.')
-           case (-13)
-               err = linalg_state_type(this,LINALG_INTERNAL_ERROR,'Insufficient work array size.')
-           case (1:)
-               err = linalg_state_type(this,LINALG_ERROR,'Eigenvalue computation did not converge.')
-           case default
-               err = linalg_state_type(this,LINALG_INTERNAL_ERROR,'Unknown error returned by geev.')
-        end select
-
-     end subroutine handle_geev_info
-
-     !> Process GGEV output flags
-     pure subroutine handle_ggev_info(err,info,shapea,shapeb)
-        !> Error handler
-        type(linalg_state_type), intent(inout) :: err
-        !> GEEV return flag
-        integer(ilp), intent(in) :: info
-        !> Input matrix size
-        integer(ilp), intent(in) :: shapea(2),shapeb(2)
-
-        select case (info)
-           case (0)
-               ! Success!
-               err%state = LINALG_SUCCESS
-           case (-1)
-               err = linalg_state_type(this,LINALG_INTERNAL_ERROR,'Invalid task ID: left eigenvectors.')
-           case (-2)
-               err = linalg_state_type(this,LINALG_INTERNAL_ERROR,'Invalid task ID: right eigenvectors.')
-           case (-5,-3)
-               err = linalg_state_type(this,LINALG_VALUE_ERROR,'invalid matrix size: a=',shapea)
-           case (-7)
-               err = linalg_state_type(this,LINALG_VALUE_ERROR,'invalid matrix size: b=',shapeb)               
-           case (-12)
-               err = linalg_state_type(this,LINALG_VALUE_ERROR,'insufficient left vector matrix size.')
-           case (-14)
-               err = linalg_state_type(this,LINALG_VALUE_ERROR,'insufficient right vector matrix size.')
-           case (-16)
-               err = linalg_state_type(this,LINALG_INTERNAL_ERROR,'Insufficient work array size.')
-           case (1:)
-               err = linalg_state_type(this,LINALG_ERROR,'Eigenvalue computation did not converge.')
-           case default
-               err = linalg_state_type(this,LINALG_INTERNAL_ERROR,'Unknown error returned by ggev.')
-        end select
-
-     end subroutine handle_ggev_info
-
-     !> Process SYEV/HEEV output flags
-     elemental subroutine handle_heev_info(err,info,m,n)
-        !> Error handler
-        type(linalg_state_type), intent(inout) :: err
-        !> SYEV/HEEV return flag
-        integer(ilp), intent(in) :: info
-        !> Input matrix size
-        integer(ilp), intent(in) :: m,n
-
-        select case (info)
-           case (0)
-               ! Success!
-               err%state = LINALG_SUCCESS
-           case (-1)
-               err = linalg_state_type(this,LINALG_INTERNAL_ERROR,'Invalid eigenvector request.')
-           case (-2)
-               err = linalg_state_type(this,LINALG_INTERNAL_ERROR,'Invalid triangular section request.')
-           case (-5,-3)
-               err = linalg_state_type(this,LINALG_VALUE_ERROR,'invalid matrix size: a=',[m,n])
-           case (-8)
-               err = linalg_state_type(this,LINALG_INTERNAL_ERROR,'insufficient workspace size.')
-           case (1:)
-               err = linalg_state_type(this,LINALG_ERROR,'Eigenvalue computation did not converge.')
-           case default
-               err = linalg_state_type(this,LINALG_INTERNAL_ERROR,'Unknown error returned by syev/heev.')
-        end select
-
-     end subroutine handle_heev_info
-
      #:for rk,rt,ri in RC_KINDS_TYPES
      #:for ep,ei in EIG_PROBLEM_LIST
 
@@ -370,7 +274,7 @@ submodule (stdlib_linalg) stdlib_linalg_eigenvalues
                        #:endif                       
                        umat,ldu,vmat,ldv,&
                        work_dummy,lwork,#{if rt.startswith('complex')}#rwork,#{endif}#info)
-             call handle_${ei}$_info(err0,info,shape(amat)#{if ei=='ggev'}#,shape(bmat)#{endif}#)
+             call handle_${ei}$_info(this,err0,info,shape(amat)#{if ei=='ggev'}#,shape(bmat)#{endif}#)
 
              ! Compute eigenvalues
              if (info==0) then
@@ -390,7 +294,7 @@ submodule (stdlib_linalg) stdlib_linalg_eigenvalues
                           #:endif                                                 
                           umat,ldu,vmat,ldv,&            
                           work,lwork,#{if rt.startswith('complex')}#rwork,#{endif}#info)
-                call handle_${ei}$_info(err0,info,shape(amat)#{if ei=='ggev'}#,shape(bmat)#{endif}#)
+                call handle_${ei}$_info(this,err0,info,shape(amat)#{if ei=='ggev'}#,shape(bmat)#{endif}#)
 
              endif
              
@@ -584,7 +488,7 @@ submodule (stdlib_linalg) stdlib_linalg_eigenvalues
          #:else
          call syev(task,triangle,n,amat,lda,lambda,work_dummy,lwork,info)
          #:endif
-         call handle_heev_info(err0,info,m,n)
+         call handle_heev_info(this,err0,info,m,n)
 
          ! Compute eigenvalues
          if (info==0) then
@@ -599,7 +503,7 @@ submodule (stdlib_linalg) stdlib_linalg_eigenvalues
             #:else
             call syev(task,triangle,n,amat,lda,lambda,work,lwork,info)
             #:endif
-            call handle_heev_info(err0,info,m,n)
+            call handle_heev_info(this,err0,info,m,n)
 
          endif
          
diff --git a/src/stdlib_linalg_inverse.fypp b/src/stdlib_linalg_inverse.fypp
index b38fe97ea..ef3de9579 100644
--- a/src/stdlib_linalg_inverse.fypp
+++ b/src/stdlib_linalg_inverse.fypp
@@ -4,6 +4,7 @@ submodule (stdlib_linalg) stdlib_linalg_inverse
 !! Compute inverse of a square matrix    
      use stdlib_linalg_constants
      use stdlib_linalg_lapack, only: getri,getrf,stdlib_ilaenv
+     use stdlib_linalg_lapack_aux, only: handle_getri_info
      use stdlib_linalg_state, only: linalg_state_type, linalg_error_handling, LINALG_ERROR, &
          LINALG_INTERNAL_ERROR, LINALG_VALUE_ERROR
      use ieee_arithmetic, only: ieee_value, ieee_quiet_nan
@@ -13,25 +14,6 @@ submodule (stdlib_linalg) stdlib_linalg_inverse
 
      contains
 
-     elemental subroutine handle_getri_info(info,lda,n,err)
-         integer(ilp), intent(in) :: info,lda,n
-         type(linalg_state_type), intent(out) :: err
-
-         ! Process output
-         select case (info)
-            case (0)
-                ! Success
-            case (:-1)
-                err = linalg_state_type(this,LINALG_ERROR,'invalid matrix size a=',[lda,n])
-            case (1:)
-                ! Matrix is singular
-                err = linalg_state_type(this,LINALG_ERROR,'singular matrix')
-            case default
-                err = linalg_state_type(this,LINALG_INTERNAL_ERROR,'catastrophic error')
-         end select
-         
-     end subroutine handle_getri_info
-
      #:for rk,rt,ri in RC_KINDS_TYPES
      ! Compute the in-place square matrix inverse of a
      module subroutine stdlib_linalg_invert_inplace_${ri}$(a,pivot,err)
@@ -86,7 +68,7 @@ submodule (stdlib_linalg) stdlib_linalg_inverse
          endif
 
          ! Process output
-         call handle_getri_info(info,lda,n,err0)
+         call handle_getri_info(this,info,lda,n,err0)
 
          ! Process output and return
          if (.not.present(pivot)) deallocate(ipiv)
diff --git a/src/stdlib_linalg_least_squares.fypp b/src/stdlib_linalg_least_squares.fypp
index 14b281360..8cb374689 100644
--- a/src/stdlib_linalg_least_squares.fypp
+++ b/src/stdlib_linalg_least_squares.fypp
@@ -8,6 +8,7 @@ submodule (stdlib_linalg) stdlib_linalg_least_squares
 !! Least-squares solution to Ax=b
      use stdlib_linalg_constants
      use stdlib_linalg_lapack, only: gelsd, stdlib_ilaenv
+     use stdlib_linalg_lapack_aux, only: handle_gelsd_info
      use stdlib_linalg_state, only: linalg_state_type, linalg_error_handling, LINALG_ERROR, &
          LINALG_INTERNAL_ERROR, LINALG_VALUE_ERROR
      implicit none
@@ -16,25 +17,6 @@ submodule (stdlib_linalg) stdlib_linalg_least_squares
 
      contains
 
-     elemental subroutine handle_gelsd_info(info,lda,n,ldb,nrhs,err)
-          integer(ilp), intent(in) :: info,lda,n,ldb,nrhs
-          type(linalg_state_type), intent(out) :: err
-
-          ! Process output
-          select case (info)
-             case (0)
-                 ! Success
-             case (:-1)
-                 err = linalg_state_type(this,LINALG_VALUE_ERROR,'invalid problem size a=',[lda,n], &
-                                                                                    ', b=',[ldb,nrhs])
-             case (1:)
-                 err = linalg_state_type(this,LINALG_ERROR,'SVD did not converge.')
-             case default
-                 err = linalg_state_type(this,LINALG_INTERNAL_ERROR,'catastrophic error')
-          end select
-
-      end subroutine handle_gelsd_info
-
      #:for rk,rt,ri in RC_KINDS_TYPES
      ! Workspace needed by gelsd
      elemental subroutine ${ri}$gelsd_space(m,n,nrhs,lrwork,liwork,lcwork)
@@ -334,7 +316,7 @@ submodule (stdlib_linalg) stdlib_linalg_least_squares
             acond = singular(1)/singular(mnmin)
 
             ! Process output
-            call handle_gelsd_info(info,lda,n,ldb,nrhs,err0)
+            call handle_gelsd_info(this,info,lda,n,ldb,nrhs,err0)
          
          endif           
                   
diff --git a/src/stdlib_linalg_qr.fypp b/src/stdlib_linalg_qr.fypp
index d2656ca34..2212c39c3 100644
--- a/src/stdlib_linalg_qr.fypp
+++ b/src/stdlib_linalg_qr.fypp
@@ -3,6 +3,7 @@
 submodule (stdlib_linalg) stdlib_linalg_qr
      use stdlib_linalg_constants
      use stdlib_linalg_lapack, only: geqrf, orgqr, ungqr
+     use stdlib_linalg_lapack_aux, only: handle_geqrf_info, handle_orgqr_info
      use stdlib_linalg_state, only: linalg_state_type, linalg_error_handling, LINALG_ERROR, &
          LINALG_INTERNAL_ERROR, LINALG_VALUE_ERROR     
      implicit none
@@ -42,51 +43,6 @@ submodule (stdlib_linalg) stdlib_linalg_qr
          
      end subroutine check_problem_size
      
-     elemental subroutine handle_orgqr_info(info,m,n,k,lwork,err)
-         integer(ilp), intent(in) :: info,m,n,k,lwork
-         type(linalg_state_type), intent(out) :: err
-
-         ! Process output
-         select case (info)
-            case (0)
-                ! Success
-            case (-1)
-                err = linalg_state_type(this,LINALG_VALUE_ERROR,'invalid matrix size m=',m)
-            case (-2)
-                err = linalg_state_type(this,LINALG_VALUE_ERROR,'invalid matrix size n=',n)
-            case (-4)
-                err = linalg_state_type(this,LINALG_VALUE_ERROR,'invalid k=min(m,n)=',k)
-            case (-5)
-                err = linalg_state_type(this,LINALG_VALUE_ERROR,'invalid matrix size a=',[m,n])                
-            case (-8)
-                err = linalg_state_type(this,LINALG_ERROR,'invalid input for lwork=',lwork)
-            case default
-                err = linalg_state_type(this,LINALG_INTERNAL_ERROR,'catastrophic error')
-         end select
-
-     end subroutine handle_orgqr_info
-     
-     elemental subroutine handle_geqrf_info(info,m,n,lwork,err)
-         integer(ilp), intent(in) :: info,m,n,lwork
-         type(linalg_state_type), intent(out) :: err
-
-         ! Process output
-         select case (info)
-            case (0)
-                ! Success
-            case (-1)
-                err = linalg_state_type(this,LINALG_VALUE_ERROR,'invalid matrix size m=',m)
-            case (-2)
-                err = linalg_state_type(this,LINALG_VALUE_ERROR,'invalid matrix size n=',n)
-            case (-4)
-                err = linalg_state_type(this,LINALG_VALUE_ERROR,'invalid matrix size a=',[m,n])
-            case (-7)
-                err = linalg_state_type(this,LINALG_ERROR,'invalid input for lwork=',lwork)
-            case default
-                err = linalg_state_type(this,LINALG_INTERNAL_ERROR,'catastrophic error')
-         end select
-
-     end subroutine handle_geqrf_info
 
      #:for rk,rt,ri in RC_KINDS_TYPES
      
@@ -113,7 +69,7 @@ submodule (stdlib_linalg) stdlib_linalg_qr
          ! QR space
          lwork_qr = -1_ilp
          call geqrf(m,n,a_dummy,m,tau_dummy,work_dummy,lwork_qr,info)    
-         call handle_geqrf_info(info,m,n,lwork_qr,err0)
+         call handle_geqrf_info(this,info,m,n,lwork_qr,err0)
          if (err0%error()) then 
             call linalg_error_handling(err0,err)
             return
@@ -124,7 +80,7 @@ submodule (stdlib_linalg) stdlib_linalg_qr
          lwork_ord = -1_ilp
          call #{if rt.startswith('complex')}# ungqr #{else}# orgqr #{endif}#  &
               (m,m,k,a_dummy,m,tau_dummy,work_dummy,lwork_ord,info)
-         call handle_orgqr_info(info,m,n,k,lwork_ord,err0)   
+         call handle_orgqr_info(this,info,m,n,k,lwork_ord,err0)   
          if (err0%error()) then 
             call linalg_error_handling(err0,err)
             return
@@ -224,7 +180,7 @@ submodule (stdlib_linalg) stdlib_linalg_qr
              
              ! Compute factorization. 
              call geqrf(m,n,amat,m,tau,work,lwork,info) 
-             call handle_geqrf_info(info,m,n,lwork,err0)
+             call handle_geqrf_info(this,info,m,n,lwork,err0)
              
              if (err0%ok()) then      
                 
@@ -236,7 +192,7 @@ submodule (stdlib_linalg) stdlib_linalg_qr
                  ! Convert K elementary reflectors tau(1:k) -> orthogonal matrix Q
                  call #{if rt.startswith('complex')}# ungqr #{else}# orgqr #{endif}#  &
                       (q1,q2,k,amat,lda,tau,work,lwork,info)
-                 call handle_orgqr_info(info,m,n,k,lwork,err0)      
+                 call handle_orgqr_info(this,info,m,n,k,lwork,err0)      
                       
                  ! Copy result back to Q
                  if (.not.use_q_matrix) q = amat(:q1,:q2) 
diff --git a/src/stdlib_linalg_schur.fypp b/src/stdlib_linalg_schur.fypp
index 68b526249..7d28200cb 100644
--- a/src/stdlib_linalg_schur.fypp
+++ b/src/stdlib_linalg_schur.fypp
@@ -3,6 +3,7 @@
 submodule (stdlib_linalg) stdlib_linalg_schur
     use stdlib_linalg_constants
     use stdlib_linalg_lapack, only: gees
+    use stdlib_linalg_lapack_aux, only: handle_gees_info
     use stdlib_linalg_state, only: linalg_state_type, linalg_error_handling, LINALG_ERROR, &
         LINALG_INTERNAL_ERROR, LINALG_VALUE_ERROR
     implicit none
@@ -35,48 +36,6 @@ submodule (stdlib_linalg) stdlib_linalg_schur
         gees_sort_eigs = merge(GEES_SORTED_VECTORS,GEES_NOT,sorted)
     end function gees_sort_eigs    
 
-    !> Wrapper function to handle GEES error codes
-    elemental subroutine handle_gees_info(info, m, n, ldvs, err)
-        integer(ilp), intent(in) :: info, m, n, ldvs
-        type(linalg_state_type), intent(out) :: err
-
-        ! Process GEES output
-        select case (info)
-        case (0_ilp)
-            ! Success
-        case (-1_ilp)
-            ! Vector not wanted, but task is wrong
-            err = linalg_state_type(this, LINALG_INTERNAL_ERROR,'Invalid Schur vector task request')
-        case (-2_ilp)
-            ! Vector not wanted, but task is wrong
-            err = linalg_state_type(this, LINALG_INTERNAL_ERROR,'Invalid sorting task request')            
-        case (-4_ilp,-6_ilp)
-            ! Vector not wanted, but task is wrong
-            err = linalg_state_type(this, LINALG_VALUE_ERROR,'Invalid/non-square input matrix size:',[m,n])  
-        case (-11_ilp)
-            err = linalg_state_type(this, LINALG_VALUE_ERROR,'Schur vector matrix has insufficient size',[ldvs,n])
-        case (-13_ilp)
-            err = linalg_state_type(this, LINALG_INTERNAL_ERROR,'Insufficient working storage size')
-        case (1_ilp:)
-            
-            if (info==n+2) then 
-                err = linalg_state_type(this, LINALG_ERROR, 'Ill-conditioned problem: could not sort eigenvalues')
-            elseif (info==n+1) then 
-                err = linalg_state_type(this, LINALG_ERROR, 'Some selected eigenvalues lost property due to sorting')                
-            elseif (info==n) then 
-                err = linalg_state_type(this, LINALG_ERROR, 'Convergence failure: no converged eigenvalues')
-            else
-                err = linalg_state_type(this, LINALG_ERROR, 'Convergence failure; converged range is',[info,n])                            
-            end if
-            
-        case default
-            
-            err = linalg_state_type(this, LINALG_INTERNAL_ERROR, 'GEES catastrophic error: info=', info)
-
-        end select
-        
-    end subroutine handle_gees_info
-
     #:for rk, rt, ri in RC_KINDS_TYPES    
     !> Workspace query
     module subroutine get_schur_${ri}$_workspace(a,lwork,err)
@@ -112,7 +71,7 @@ submodule (stdlib_linalg) stdlib_linalg_schur
         call gees(jobvs,sort,do_not_select,n,amat,m,sdim,wr_dummy,#{if rt.startswith('r')}#wi_dummy, #{endif}#&
                   vs_dummy,m,work_dummy,lwork,#{if rt.startswith('c')}#rwork_dummy,#{endif}#bwork_dummy,info)
         if (info==0) lwork = nint(real(work_dummy(1),kind=${rk}$),kind=ilp)
-        call handle_gees_info(info,m,n,m,err0)
+        call handle_gees_info(this,info,m,n,m,err0)
         call linalg_error_handling(err0,err)
         
         contains
@@ -275,7 +234,7 @@ submodule (stdlib_linalg) stdlib_linalg_schur
             ! Compute Schur decomposition
             call gees(jobvs,sort,eig_select,nt,t,mt,sdim,eigs,#{if rt.startswith('r')}#eigi,#{endif}# &
                       vs,ldvs,work,lwork,#{if rt.startswith('c')}#rwork,#{endif}#bwork,info)
-            call handle_gees_info(info,m,n,m,err0)            
+            call handle_gees_info(this,info,m,n,m,err0)            
             
         
         end if  
diff --git a/src/stdlib_linalg_solve.fypp b/src/stdlib_linalg_solve.fypp
index 1c700683b..b4df33133 100644
--- a/src/stdlib_linalg_solve.fypp
+++ b/src/stdlib_linalg_solve.fypp
@@ -8,6 +8,7 @@ submodule (stdlib_linalg) stdlib_linalg_solve
 !! Solve linear system Ax=b
      use stdlib_linalg_constants
      use stdlib_linalg_lapack, only: gesv
+     use stdlib_linalg_lapack_aux, only: handle_gesv_info
      use stdlib_linalg_state, only: linalg_state_type, linalg_error_handling, LINALG_ERROR, &
           LINALG_INTERNAL_ERROR, LINALG_VALUE_ERROR     
      implicit none
@@ -16,29 +17,6 @@ submodule (stdlib_linalg) stdlib_linalg_solve
 
      contains
      
-     elemental subroutine handle_gesv_info(info,lda,n,nrhs,err)
-         integer(ilp), intent(in) :: info,lda,n,nrhs
-         type(linalg_state_type), intent(out) :: err
-
-         ! Process output
-         select case (info)
-            case (0)
-                ! Success
-            case (-1)
-                err = linalg_state_type(this,LINALG_VALUE_ERROR,'invalid problem size n=',n)
-            case (-2)
-                err = linalg_state_type(this,LINALG_VALUE_ERROR,'invalid rhs size n=',nrhs)
-            case (-4)
-                err = linalg_state_type(this,LINALG_VALUE_ERROR,'invalid matrix size a=',[lda,n])
-            case (-7)
-                err = linalg_state_type(this,LINALG_ERROR,'invalid matrix size a=',[lda,n])
-            case (1:)
-                err = linalg_state_type(this,LINALG_ERROR,'singular matrix')
-            case default
-                err = linalg_state_type(this,LINALG_INTERNAL_ERROR,'catastrophic error')
-         end select
-         
-     end subroutine handle_gesv_info
 
      #:for nd,ndsuf,nde in ALL_RHS
      #:for rk,rt,ri in RC_KINDS_TYPES
@@ -152,7 +130,7 @@ submodule (stdlib_linalg) stdlib_linalg_solve
          call gesv(n,nrhs,amat,lda,ipiv,xmat,ldb,info)
 
          ! Process output
-         call handle_gesv_info(info,lda,n,nrhs,err0)
+         call handle_gesv_info(this,info,lda,n,nrhs,err0)
 
          if (copy_a) deallocate(amat)
          if (.not.present(pivot)) deallocate(ipiv)
diff --git a/src/stdlib_linalg_svd.fypp b/src/stdlib_linalg_svd.fypp
index 5c1befe79..3d9e730b8 100644
--- a/src/stdlib_linalg_svd.fypp
+++ b/src/stdlib_linalg_svd.fypp
@@ -4,6 +4,7 @@ submodule(stdlib_linalg) stdlib_linalg_svd
 !! Singular-Value Decomposition    
      use stdlib_linalg_constants
      use stdlib_linalg_lapack, only: gesdd
+     use stdlib_linalg_lapack_aux, only: handle_gesdd_info
      use stdlib_linalg_state, only: linalg_state_type, linalg_error_handling, LINALG_ERROR, &
          LINALG_INTERNAL_ERROR, LINALG_VALUE_ERROR, LINALG_SUCCESS
      implicit none
@@ -26,38 +27,6 @@ submodule(stdlib_linalg) stdlib_linalg_svd
 
      contains
 
-     !> Process GESDD output flag
-     elemental subroutine handle_gesdd_info(err,info,m,n)
-        !> Error handler
-        type(linalg_state_type), intent(inout) :: err
-        !> GESDD return flag
-        integer(ilp), intent(in) :: info
-        !> Input matrix size
-        integer(ilp), intent(in) :: m,n
-
-        select case (info)
-           case (0)
-               ! Success!
-               err%state = LINALG_SUCCESS
-           case (-1)
-               err = linalg_state_type(this,LINALG_INTERNAL_ERROR,'Invalid task ID on input to GESDD.')
-           case (-5,-3:-2)
-               err = linalg_state_type(this,LINALG_VALUE_ERROR,'invalid matrix size: a=',[m,n])
-           case (-8)
-               err = linalg_state_type(this,LINALG_VALUE_ERROR,'invalid matrix U size, with a=',[m,n])
-           case (-10)
-               err = linalg_state_type(this,LINALG_VALUE_ERROR,'invalid matrix V size, with a=',[m,n])
-           case (-4)
-               err = linalg_state_type(this,LINALG_VALUE_ERROR,'A contains invalid/NaN values.')
-           case (1:)
-               err = linalg_state_type(this,LINALG_ERROR,'SVD computation did not converge.')
-           case default
-               err = linalg_state_type(this,LINALG_INTERNAL_ERROR,'Unknown error returned by GESDD.')
-        end select
-
-     end subroutine handle_gesdd_info
-
-
      #:for rk,rt,ri in RC_KINDS_TYPES
 
      !> Singular values of matrix A
@@ -265,7 +234,7 @@ submodule(stdlib_linalg) stdlib_linalg_svd
          lwork = -1_ilp
          call gesdd(task,m,n,amat,lda,s,umat,ldu,vtmat,ldvt,&
                     work_dummy,lwork,#{if rt.startswith('complex')}#rwork,#{endif}#iwork,info)
-         call handle_gesdd_info(err0,info,m,n)
+         call handle_gesdd_info(this,err0,info,m,n)
 
          ! Compute SVD
          if (info==0) then
@@ -281,7 +250,7 @@ submodule(stdlib_linalg) stdlib_linalg_svd
             !> Compute SVD
             call gesdd(task,m,n,amat,lda,s,umat,ldu,vtmat,ldvt,&
                        work,lwork,#{if rt.startswith('complex')}#rwork,#{endif}#iwork,info)
-            call handle_gesdd_info(err0,info,m,n)
+            call handle_gesdd_info(this,err0,info,m,n)
 
          endif