From 9f6c2c10a23f93edf749812a8f710e350a9605be Mon Sep 17 00:00:00 2001
From: Sudhanshu Thakur <sudhanshukit2021@gmail.com>
Date: Wed, 25 Mar 2026 17:50:39 +0530
Subject: [PATCH 01/10] feat(linalg): add permutation matrix type with
 efficient vector storage (#891)

---
 src/linalg/CMakeLists.txt                 |   1 +
 src/linalg/stdlib_linalg_permutation.fypp | 254 ++++++++++++++++++++
 test/linalg/CMakeLists.txt                |   2 +
 test/linalg/test_linalg_permutation.fypp  | 278 ++++++++++++++++++++++
 4 files changed, 535 insertions(+)
 create mode 100644 src/linalg/stdlib_linalg_permutation.fypp
 create mode 100644 test/linalg/test_linalg_permutation.fypp

diff --git a/src/linalg/CMakeLists.txt b/src/linalg/CMakeLists.txt
index e0a6c10dd..778ed3353 100644
--- a/src/linalg/CMakeLists.txt
+++ b/src/linalg/CMakeLists.txt
@@ -11,6 +11,7 @@ set(linalg_fppFiles
     stdlib_linalg_matrix_functions.fypp
     stdlib_linalg_norms.fypp
     stdlib_linalg_outer_product.fypp
+    stdlib_linalg_permutation.fypp
     stdlib_linalg_pinv.fypp
     stdlib_linalg_qr.fypp
     stdlib_linalg_schur.fypp
diff --git a/src/linalg/stdlib_linalg_permutation.fypp b/src/linalg/stdlib_linalg_permutation.fypp
new file mode 100644
index 000000000..98e5e982e
--- /dev/null
+++ b/src/linalg/stdlib_linalg_permutation.fypp
@@ -0,0 +1,254 @@
+!> Permutation matrix operations using efficient vector storage.
+!>
+!> Instead of storing permutations as dense n×n matrices (O(n²) memory),
+!> we represent them as integer vectors of length n (O(n) memory).
+!> All operations (inversion, left/right multiplication, composition)
+!> run in O(n) time per column/row rather than O(n²).
+!>
+!> Two formats are supported:
+!>   - Finalized: perm(i) = j means row i maps to row j (direct mapping)
+!>   - LAPACK pivot: ipiv(i) = j means swap rows i and j at step i
+!>
+!> Reference: Issue #891 (https://github.com/fortran-lang/stdlib/issues/891)
+
+#:include "common.fypp"
+#:set REAL_KINDS_TYPES = [("sp","real(sp)"), ("dp","real(dp)")]
+#:set INT_KINDS_TYPES = [("int32","integer(int32)"), ("int64","integer(int64)")]
+
+module stdlib_linalg_permutation
+  use stdlib_kinds, only: sp, dp, int32, int64
+  implicit none
+  private
+
+  !> Permutation type: stores permutation as integer vector
+  type, public :: permutation_t
+    integer, allocatable :: perm(:)
+    integer :: n = 0
+    logical :: finalized = .true.
+  end type
+
+  ! Public procedures
+  public :: permutation_identity
+  public :: permutation_init
+  public :: permutation_is_valid
+  public :: permutation_invert
+  public :: permutation_compose
+  public :: permutation_from_lapack
+  public :: permutation_to_matrix
+  public :: permutation_apply_left
+  public :: permutation_apply_right
+  public :: permutation_apply_vector
+
+  ! Generic interfaces for multi-type support
+  interface permutation_apply_left
+#:for rk, rt in REAL_KINDS_TYPES
+    module procedure permutation_apply_left_${rk}$
+#:endfor
+  end interface
+
+  interface permutation_apply_right
+#:for rk, rt in REAL_KINDS_TYPES
+    module procedure permutation_apply_right_${rk}$
+#:endfor
+  end interface
+
+  interface permutation_apply_vector
+#:for rk, rt in REAL_KINDS_TYPES
+    module procedure permutation_apply_vector_${rk}$
+#:endfor
+  end interface
+
+  interface permutation_to_matrix
+#:for rk, rt in REAL_KINDS_TYPES
+    module procedure permutation_to_matrix_${rk}$
+#:endfor
+  end interface
+
+contains
+
+  !> Create an identity permutation of size n
+  function permutation_identity(n) result(p)
+    integer, intent(in) :: n
+    type(permutation_t) :: p
+    integer :: i
+
+    p%n = n
+    p%finalized = .true.
+    allocate(p%perm(n))
+    do i = 1, n
+      p%perm(i) = i
+    end do
+  end function
+
+  !> Create a permutation from a given integer vector
+  function permutation_init(perm_vec, finalized) result(p)
+    integer, intent(in) :: perm_vec(:)
+    logical, intent(in), optional :: finalized
+    type(permutation_t) :: p
+
+    p%n = size(perm_vec)
+    allocate(p%perm(p%n), source=perm_vec)
+    if (present(finalized)) then
+      p%finalized = finalized
+    else
+      p%finalized = .true.
+    end if
+  end function
+
+  !> Check if a permutation is valid (each integer 1..n appears exactly once)
+  function permutation_is_valid(p) result(valid)
+    type(permutation_t), intent(in) :: p
+    logical :: valid
+    logical, allocatable :: seen(:)
+    integer :: i
+
+    valid = .false.
+    if (.not. allocated(p%perm)) return
+    if (p%n < 1) return
+    if (size(p%perm) /= p%n) return
+
+    allocate(seen(p%n), source=.false.)
+    do i = 1, p%n
+      if (p%perm(i) < 1 .or. p%perm(i) > p%n) return
+      if (seen(p%perm(i))) return
+      seen(p%perm(i)) = .true.
+    end do
+    valid = .true.
+  end function
+
+  !> Invert a finalized permutation
+  !> If perm(i) = j, then inv(j) = i
+  recursive function permutation_invert(p) result(inv)
+    type(permutation_t), intent(in) :: p
+    type(permutation_t) :: inv
+    integer :: i
+
+    inv%n = p%n
+    inv%finalized = .true.
+    allocate(inv%perm(p%n))
+
+    if (p%finalized) then
+      ! Direct inversion: swap keys and values
+      do i = 1, p%n
+        inv%perm(p%perm(i)) = i
+      end do
+    else
+      ! For LAPACK format: first finalize, then invert
+      inv = permutation_invert(permutation_finalize(p))
+    end if
+  end function
+
+  !> Convert LAPACK pivot vector (swap sequence) to finalized permutation
+  !> ipiv(i) = j means "swap rows i and j" applied at step i
+  function permutation_from_lapack(ipiv, n) result(p)
+    integer, intent(in) :: ipiv(:)
+    integer, intent(in) :: n
+    type(permutation_t) :: p
+    integer :: i, tmp
+
+    p%n = n
+    p%finalized = .true.
+    allocate(p%perm(n))
+
+    ! Start with identity
+    do i = 1, n
+      p%perm(i) = i
+    end do
+
+    ! Replay swap sequence forward
+    do i = 1, min(n, size(ipiv))
+      if (ipiv(i) /= i .and. ipiv(i) >= 1 .and. ipiv(i) <= n) then
+        tmp = p%perm(i)
+        p%perm(i) = p%perm(ipiv(i))
+        p%perm(ipiv(i)) = tmp
+      end if
+    end do
+  end function
+
+  !> Convert nonfinalized (swap) to finalized (direct) permutation
+  function permutation_finalize(p) result(fp)
+    type(permutation_t), intent(in) :: p
+    type(permutation_t) :: fp
+
+    if (p%finalized) then
+      fp = p
+      return
+    end if
+    ! Treat swap sequence as LAPACK pivots
+    fp = permutation_from_lapack(p%perm, p%n)
+  end function
+
+  !> Compose two finalized permutations: result = p1 * p2
+  !> Apply p2 first, then p1
+  function permutation_compose(p1, p2) result(p)
+    type(permutation_t), intent(in) :: p1, p2
+    type(permutation_t) :: p
+    integer :: i
+
+    p%n = p1%n
+    p%finalized = .true.
+    allocate(p%perm(p%n))
+
+    do i = 1, p%n
+      p%perm(i) = p1%perm(p2%perm(i))
+    end do
+  end function
+
+  !> ----- Templated operations for multiple real kinds -----
+
+#:for rk, rt in REAL_KINDS_TYPES
+
+  !> Apply permutation from the left: B = P * A (row permutation)
+  !> Row i of B comes from row perm(i) of A
+  subroutine permutation_apply_left_${rk}$ (p, A, B)
+    type(permutation_t), intent(in) :: p
+    ${rt}$, intent(in) :: A(:,:)
+    ${rt}$, intent(out) :: B(:,:)
+    integer :: i
+
+    do i = 1, p%n
+      B(i, :) = A(p%perm(i), :)
+    end do
+  end subroutine
+
+  !> Apply permutation from the right: B = A * P (column permutation)
+  !> Column perm(j) of B comes from column j of A
+  subroutine permutation_apply_right_${rk}$ (p, A, B)
+    type(permutation_t), intent(in) :: p
+    ${rt}$, intent(in) :: A(:,:)
+    ${rt}$, intent(out) :: B(:,:)
+    integer :: j
+
+    do j = 1, p%n
+      B(:, p%perm(j)) = A(:, j)
+    end do
+  end subroutine
+
+  !> Apply permutation to a vector: b = P * a
+  !> b(i) = a(perm(i))
+  subroutine permutation_apply_vector_${rk}$ (p, a, b)
+    type(permutation_t), intent(in) :: p
+    ${rt}$, intent(in) :: a(:)
+    ${rt}$, intent(out) :: b(:)
+    integer :: i
+
+    do i = 1, p%n
+      b(i) = a(p%perm(i))
+    end do
+  end subroutine
+
+  !> Convert permutation to dense matrix (for testing / small cases)
+  subroutine permutation_to_matrix_${rk}$ (p, A)
+    type(permutation_t), intent(in) :: p
+    ${rt}$, intent(out) :: A(:,:)
+    integer :: i
+
+    A = 0.0_${rk}$
+    do i = 1, p%n
+      A(i, p%perm(i)) = 1.0_${rk}$
+    end do
+  end subroutine
+
+#:endfor
+
+end module stdlib_linalg_permutation
diff --git a/test/linalg/CMakeLists.txt b/test/linalg/CMakeLists.txt
index 097591893..f2e5fdf50 100644
--- a/test/linalg/CMakeLists.txt
+++ b/test/linalg/CMakeLists.txt
@@ -21,6 +21,7 @@ set(
   "test_linalg_cholesky.fypp"
   "test_linalg_solve_chol.fypp"
   "test_linalg_expm.fypp"
+  "test_linalg_permutation.fypp"
 )
 
 # Preprocessed files to contain preprocessor directives -> .F90
@@ -57,4 +58,5 @@ ADDTEST(linalg_sparse)
 if (STDLIB_SPECIALMATRICES)
   ADDTEST(linalg_specialmatrices)
 endif()
+ADDTEST(linalg_permutation)
 ADDTESTPP(blas_lapack)
diff --git a/test/linalg/test_linalg_permutation.fypp b/test/linalg/test_linalg_permutation.fypp
new file mode 100644
index 000000000..918292e70
--- /dev/null
+++ b/test/linalg/test_linalg_permutation.fypp
@@ -0,0 +1,278 @@
+!> Tests for stdlib_linalg_permutation module
+#:include "common.fypp"
+
+module test_linalg_permutation
+  use testdrive, only: new_unittest, unittest_type, error_type, check
+  use stdlib_kinds, only: sp, dp, int32
+  use stdlib_linalg_permutation
+  implicit none
+  private
+
+  public :: collect_permutation_tests
+
+contains
+
+  subroutine collect_permutation_tests(testsuite)
+    type(unittest_type), allocatable, intent(out) :: testsuite(:)
+    testsuite = [ &
+      new_unittest("identity", test_identity), &
+      new_unittest("is_valid", test_is_valid), &
+      new_unittest("invalid_rejects", test_invalid_rejects), &
+      new_unittest("invert", test_invert), &
+      new_unittest("invert_involution", test_invert_involution), &
+      new_unittest("apply_left", test_apply_left), &
+      new_unittest("apply_right", test_apply_right), &
+      new_unittest("apply_vector", test_apply_vector), &
+      new_unittest("compose", test_compose), &
+      new_unittest("from_lapack", test_from_lapack), &
+      new_unittest("to_matrix", test_to_matrix), &
+      new_unittest("identity_noop", test_identity_noop) &
+    ]
+  end subroutine
+
+  !> Test 1: Identity permutation
+  subroutine test_identity(error)
+    type(error_type), allocatable, intent(out) :: error
+    type(permutation_t) :: p
+    integer :: i
+
+    p = permutation_identity(5)
+    call check(error, p%n == 5, "identity size")
+    if (allocated(error)) return
+    do i = 1, 5
+      call check(error, p%perm(i) == i, "identity element")
+      if (allocated(error)) return
+    end do
+    call check(error, permutation_is_valid(p), "identity valid")
+  end subroutine
+
+  !> Test 2: Valid permutation accepted
+  subroutine test_is_valid(error)
+    type(error_type), allocatable, intent(out) :: error
+    type(permutation_t) :: p
+
+    p = permutation_init([3, 1, 4, 2])
+    call check(error, permutation_is_valid(p), "valid permutation")
+  end subroutine
+
+  !> Test 3: Invalid permutations rejected
+  subroutine test_invalid_rejects(error)
+    type(error_type), allocatable, intent(out) :: error
+    type(permutation_t) :: p
+
+    ! Duplicate entry
+    p = permutation_init([1, 1, 3])
+    call check(error, .not. permutation_is_valid(p), "reject duplicate")
+    if (allocated(error)) return
+
+    ! Out of range
+    p = permutation_init([1, 5, 3])
+    call check(error, .not. permutation_is_valid(p), "reject out of range")
+  end subroutine
+
+  !> Test 4: Inversion correctness
+  subroutine test_invert(error)
+    type(error_type), allocatable, intent(out) :: error
+    type(permutation_t) :: p, inv
+
+    ! perm = [3,1,2] means: pos1->3, pos2->1, pos3->2
+    ! inverse: pos3->1, pos1->2, pos2->3 => inv = [2,3,1]
+    p = permutation_init([3, 1, 2])
+    inv = permutation_invert(p)
+    call check(error, inv%perm(1) == 2, "inv(1)")
+    if (allocated(error)) return
+    call check(error, inv%perm(2) == 3, "inv(2)")
+    if (allocated(error)) return
+    call check(error, inv%perm(3) == 1, "inv(3)")
+  end subroutine
+
+  !> Test 5: P * P^-1 = Identity (involution property)
+  subroutine test_invert_involution(error)
+    type(error_type), allocatable, intent(out) :: error
+    type(permutation_t) :: p, inv, composed
+    integer :: i
+
+    p = permutation_init([3, 1, 4, 2])
+    inv = permutation_invert(p)
+    composed = permutation_compose(p, inv)
+    do i = 1, 4
+      call check(error, composed%perm(i) == i, "P*P_inv = I")
+      if (allocated(error)) return
+    end do
+  end subroutine
+
+  !> Test 6: Left multiplication (row permutation)
+  subroutine test_apply_left(error)
+    type(error_type), allocatable, intent(out) :: error
+    type(permutation_t) :: p
+    real(dp) :: A(3,2), B(3,2)
+
+    ! A = [[1,2],[3,4],[5,6]]
+    A(1,:) = [1.0_dp, 2.0_dp]
+    A(2,:) = [3.0_dp, 4.0_dp]
+    A(3,:) = [5.0_dp, 6.0_dp]
+
+    ! perm = [3,1,2]: row1 <- row3, row2 <- row1, row3 <- row2
+    p = permutation_init([3, 1, 2])
+    call permutation_apply_left(p, A, B)
+
+    ! B should be [[5,6],[1,2],[3,4]]
+    call check(error, abs(B(1,1) - 5.0_dp) < 1.0e-10_dp, "left(1,1)")
+    if (allocated(error)) return
+    call check(error, abs(B(2,1) - 1.0_dp) < 1.0e-10_dp, "left(2,1)")
+    if (allocated(error)) return
+    call check(error, abs(B(3,1) - 3.0_dp) < 1.0e-10_dp, "left(3,1)")
+  end subroutine
+
+  !> Test 7: Right multiplication (column permutation)
+  subroutine test_apply_right(error)
+    type(error_type), allocatable, intent(out) :: error
+    type(permutation_t) :: p
+    real(dp) :: A(2,3), B(2,3)
+
+    A(1,:) = [1.0_dp, 2.0_dp, 3.0_dp]
+    A(2,:) = [4.0_dp, 5.0_dp, 6.0_dp]
+
+    ! perm = [2,3,1]: col_perm(2) <- col1, col_perm(3) <- col2, col_perm(1) <- col3
+    p = permutation_init([2, 3, 1])
+    call permutation_apply_right(p, A, B)
+
+    ! B(:, perm(j)) = A(:, j)
+    ! B(:,2) = A(:,1) = [1,4], B(:,3) = A(:,2) = [2,5], B(:,1) = A(:,3) = [3,6]
+    call check(error, abs(B(1,1) - 3.0_dp) < 1.0e-10_dp, "right(1,1)")
+    if (allocated(error)) return
+    call check(error, abs(B(1,2) - 1.0_dp) < 1.0e-10_dp, "right(1,2)")
+    if (allocated(error)) return
+    call check(error, abs(B(1,3) - 2.0_dp) < 1.0e-10_dp, "right(1,3)")
+  end subroutine
+
+  !> Test 8: Vector permutation
+  subroutine test_apply_vector(error)
+    type(error_type), allocatable, intent(out) :: error
+    type(permutation_t) :: p
+    real(dp) :: a(4), b(4)
+
+    a = [10.0_dp, 20.0_dp, 30.0_dp, 40.0_dp]
+    p = permutation_init([3, 1, 4, 2])
+    call permutation_apply_vector(p, a, b)
+
+    ! b(1) = a(3) = 30, b(2) = a(1) = 10, b(3) = a(4) = 40, b(4) = a(2) = 20
+    call check(error, abs(b(1) - 30.0_dp) < 1.0e-10_dp, "vec(1)")
+    if (allocated(error)) return
+    call check(error, abs(b(2) - 10.0_dp) < 1.0e-10_dp, "vec(2)")
+    if (allocated(error)) return
+    call check(error, abs(b(3) - 40.0_dp) < 1.0e-10_dp, "vec(3)")
+    if (allocated(error)) return
+    call check(error, abs(b(4) - 20.0_dp) < 1.0e-10_dp, "vec(4)")
+  end subroutine
+
+  !> Test 9: Composition P1 * P2
+  subroutine test_compose(error)
+    type(error_type), allocatable, intent(out) :: error
+    type(permutation_t) :: p1, p2, composed
+    real(dp) :: a(3), b1(3), b2(3), b_composed(3)
+
+    p1 = permutation_init([2, 3, 1])
+    p2 = permutation_init([3, 1, 2])
+    composed = permutation_compose(p1, p2)
+
+    ! Verify: composed applied to vector = p1 applied to (p2 applied to vector)
+    a = [10.0_dp, 20.0_dp, 30.0_dp]
+    call permutation_apply_vector(p2, a, b1)       ! b1 = P2 * a
+    call permutation_apply_vector(p1, b1, b2)      ! b2 = P1 * P2 * a
+    call permutation_apply_vector(composed, a, b_composed) ! should equal b2
+
+    call check(error, abs(b2(1) - b_composed(1)) < 1.0e-10_dp, "compose(1)")
+    if (allocated(error)) return
+    call check(error, abs(b2(2) - b_composed(2)) < 1.0e-10_dp, "compose(2)")
+    if (allocated(error)) return
+    call check(error, abs(b2(3) - b_composed(3)) < 1.0e-10_dp, "compose(3)")
+  end subroutine
+
+  !> Test 10: LAPACK pivot vector conversion
+  subroutine test_from_lapack(error)
+    type(error_type), allocatable, intent(out) :: error
+    type(permutation_t) :: p
+
+    ! LAPACK pivot: [2, 3, 3]
+    ! Step 1: swap rows 1,2 -> [2,1,3]
+    ! Step 2: swap rows 2,3 -> [2,3,1]
+    ! Step 3: swap rows 3,3 -> [2,3,1] (no-op)
+    p = permutation_from_lapack([2, 3, 3], 3)
+
+    call check(error, p%perm(1) == 2, "lapack(1)")
+    if (allocated(error)) return
+    call check(error, p%perm(2) == 3, "lapack(2)")
+    if (allocated(error)) return
+    call check(error, p%perm(3) == 1, "lapack(3)")
+  end subroutine
+
+  !> Test 11: Dense matrix matches permutation operations
+  subroutine test_to_matrix(error)
+    type(error_type), allocatable, intent(out) :: error
+    type(permutation_t) :: p
+    real(dp) :: P_dense(3,3)
+
+    p = permutation_init([3, 1, 2])
+    call permutation_to_matrix(p, P_dense)
+
+    ! Row 1, col 3 should be 1 (perm(1)=3)
+    call check(error, abs(P_dense(1,3) - 1.0_dp) < 1.0e-10_dp, "dense(1,3)")
+    if (allocated(error)) return
+    ! Row 2, col 1 should be 1 (perm(2)=1)
+    call check(error, abs(P_dense(2,1) - 1.0_dp) < 1.0e-10_dp, "dense(2,1)")
+    if (allocated(error)) return
+    ! Row 1, col 1 should be 0
+    call check(error, abs(P_dense(1,1)) < 1.0e-10_dp, "dense(1,1)=0")
+  end subroutine
+
+  !> Test 12: Identity permutation is a no-op
+  subroutine test_identity_noop(error)
+    type(error_type), allocatable, intent(out) :: error
+    type(permutation_t) :: p
+    real(dp) :: A(3,3), B(3,3)
+    integer :: i, j
+
+    p = permutation_identity(3)
+    do i = 1, 3
+      do j = 1, 3
+        A(i,j) = real(i*10 + j, dp)
+      end do
+    end do
+
+    call permutation_apply_left(p, A, B)
+    do i = 1, 3
+      do j = 1, 3
+        call check(error, abs(A(i,j) - B(i,j)) < 1.0e-10_dp, "identity noop")
+        if (allocated(error)) return
+      end do
+    end do
+  end subroutine
+
+end module test_linalg_permutation
+
+program test_permutation
+    use, intrinsic :: iso_fortran_env, only : error_unit
+    use testdrive, only : run_testsuite, new_testsuite, testsuite_type
+    use test_linalg_permutation, only : collect_permutation_tests
+    implicit none
+    integer :: stat, is
+    type(testsuite_type), allocatable :: testsuites(:)
+    character(len=*), parameter :: fmt = '("#", *(1x, a))'
+
+    stat = 0
+
+    testsuites = [ &
+        new_testsuite("linalg_permutation", collect_permutation_tests) &
+        ]
+
+    do is = 1, size(testsuites)
+        write(error_unit, fmt) "Testing:", testsuites(is)%name
+        call run_testsuite(testsuites(is)%collect, error_unit, stat)
+    end do
+
+    if (stat > 0) then
+        write(error_unit, '(i0, 1x, a)') stat, "test(s) failed!"
+        error stop
+    end if
+end program test_permutation

From e1fb158e73e25e41c2480073272a0376777e6dd2 Mon Sep 17 00:00:00 2001
From: Sudhanshu Thakur <sudhanshukit2021@gmail.com>
Date: Wed, 25 Mar 2026 18:15:04 +0530
Subject: [PATCH 02/10] refactor: remove recursive keyword from
 permutation_invert

---
 src/linalg/stdlib_linalg_permutation.fypp | 25 ++++++++++++-----------
 1 file changed, 13 insertions(+), 12 deletions(-)

diff --git a/src/linalg/stdlib_linalg_permutation.fypp b/src/linalg/stdlib_linalg_permutation.fypp
index 98e5e982e..af1c9398d 100644
--- a/src/linalg/stdlib_linalg_permutation.fypp
+++ b/src/linalg/stdlib_linalg_permutation.fypp
@@ -118,24 +118,25 @@ contains
 
   !> Invert a finalized permutation
   !> If perm(i) = j, then inv(j) = i
-  recursive function permutation_invert(p) result(inv)
+  function permutation_invert(p) result(inv)
     type(permutation_t), intent(in) :: p
-    type(permutation_t) :: inv
+    type(permutation_t) :: inv, fp
     integer :: i
 
-    inv%n = p%n
-    inv%finalized = .true.
-    allocate(inv%perm(p%n))
-
     if (p%finalized) then
-      ! Direct inversion: swap keys and values
-      do i = 1, p%n
-        inv%perm(p%perm(i)) = i
-      end do
+      fp = p
     else
-      ! For LAPACK format: first finalize, then invert
-      inv = permutation_invert(permutation_finalize(p))
+      fp = permutation_finalize(p)
     end if
+
+    inv%n = fp%n
+    inv%finalized = .true.
+    allocate(inv%perm(fp%n))
+
+    ! Direct inversion: swap keys and values
+    do i = 1, fp%n
+      inv%perm(fp%perm(i)) = i
+    end do
   end function
 
   !> Convert LAPACK pivot vector (swap sequence) to finalized permutation

From 64068116f7275d29d837357fde9850d748985acf Mon Sep 17 00:00:00 2001
From: Sudhanshu Thakur <sudhanshukit2021@gmail.com>
Date: Wed, 25 Mar 2026 22:56:48 +0530
Subject: [PATCH 03/10] refactor(linalg): adressing the PR review

---
 src/linalg/stdlib_linalg_permutation.fypp | 111 +++++++++++-----------
 test/linalg/test_linalg_permutation.fypp  |  39 ++++----
 2 files changed, 74 insertions(+), 76 deletions(-)

diff --git a/src/linalg/stdlib_linalg_permutation.fypp b/src/linalg/stdlib_linalg_permutation.fypp
index af1c9398d..92825e8ea 100644
--- a/src/linalg/stdlib_linalg_permutation.fypp
+++ b/src/linalg/stdlib_linalg_permutation.fypp
@@ -12,18 +12,19 @@
 !> Reference: Issue #891 (https://github.com/fortran-lang/stdlib/issues/891)
 
 #:include "common.fypp"
-#:set REAL_KINDS_TYPES = [("sp","real(sp)"), ("dp","real(dp)")]
-#:set INT_KINDS_TYPES = [("int32","integer(int32)"), ("int64","integer(int64)")]
+#:set RC_KINDS_TYPES = REAL_KINDS_TYPES + CMPLX_KINDS_TYPES
+#:set RCI_KINDS_TYPES = RC_KINDS_TYPES + INT_KINDS_TYPES
 
 module stdlib_linalg_permutation
-  use stdlib_kinds, only: sp, dp, int32, int64
+  use stdlib_kinds, only: xdp, int8, int16, int32, int64
+  use stdlib_linalg_constants, only: sp, dp, qp, lk, ilp
   implicit none
   private
 
   !> Permutation type: stores permutation as integer vector
-  type, public :: permutation_t
-    integer, allocatable :: perm(:)
-    integer :: n = 0
+  type, public :: permutation_type
+    integer(ilp), allocatable :: perm(:)
+    integer(ilp) :: n = 0
     logical :: finalized = .true.
   end type
 
@@ -41,26 +42,26 @@ module stdlib_linalg_permutation
 
   ! Generic interfaces for multi-type support
   interface permutation_apply_left
-#:for rk, rt in REAL_KINDS_TYPES
-    module procedure permutation_apply_left_${rk}$
+#:for k1, t1 in RCI_KINDS_TYPES
+      module procedure permutation_apply_left_${t1[0]}$${k1}$
 #:endfor
   end interface
 
   interface permutation_apply_right
-#:for rk, rt in REAL_KINDS_TYPES
-    module procedure permutation_apply_right_${rk}$
+#:for k1, t1 in RCI_KINDS_TYPES
+      module procedure permutation_apply_right_${t1[0]}$${k1}$
 #:endfor
   end interface
 
   interface permutation_apply_vector
-#:for rk, rt in REAL_KINDS_TYPES
-    module procedure permutation_apply_vector_${rk}$
+#:for k1, t1 in RCI_KINDS_TYPES
+      module procedure permutation_apply_vector_${t1[0]}$${k1}$
 #:endfor
   end interface
 
   interface permutation_to_matrix
-#:for rk, rt in REAL_KINDS_TYPES
-    module procedure permutation_to_matrix_${rk}$
+#:for k1, t1 in RCI_KINDS_TYPES
+      module procedure permutation_to_matrix_${t1[0]}$${k1}$
 #:endfor
   end interface
 
@@ -68,9 +69,9 @@ contains
 
   !> Create an identity permutation of size n
   function permutation_identity(n) result(p)
-    integer, intent(in) :: n
-    type(permutation_t) :: p
-    integer :: i
+    integer(ilp), intent(in) :: n
+    type(permutation_type) :: p
+    integer(ilp) :: i
 
     p%n = n
     p%finalized = .true.
@@ -82,9 +83,9 @@ contains
 
   !> Create a permutation from a given integer vector
   function permutation_init(perm_vec, finalized) result(p)
-    integer, intent(in) :: perm_vec(:)
+    integer(ilp), intent(in) :: perm_vec(:)
     logical, intent(in), optional :: finalized
-    type(permutation_t) :: p
+    type(permutation_type) :: p
 
     p%n = size(perm_vec)
     allocate(p%perm(p%n), source=perm_vec)
@@ -97,10 +98,10 @@ contains
 
   !> Check if a permutation is valid (each integer 1..n appears exactly once)
   function permutation_is_valid(p) result(valid)
-    type(permutation_t), intent(in) :: p
+    type(permutation_type), intent(in) :: p
     logical :: valid
     logical, allocatable :: seen(:)
-    integer :: i
+    integer(ilp) :: i
 
     valid = .false.
     if (.not. allocated(p%perm)) return
@@ -119,9 +120,9 @@ contains
   !> Invert a finalized permutation
   !> If perm(i) = j, then inv(j) = i
   function permutation_invert(p) result(inv)
-    type(permutation_t), intent(in) :: p
-    type(permutation_t) :: inv, fp
-    integer :: i
+    type(permutation_type), intent(in) :: p
+    type(permutation_type) :: inv, fp
+    integer(ilp) :: i
 
     if (p%finalized) then
       fp = p
@@ -142,10 +143,10 @@ contains
   !> Convert LAPACK pivot vector (swap sequence) to finalized permutation
   !> ipiv(i) = j means "swap rows i and j" applied at step i
   function permutation_from_lapack(ipiv, n) result(p)
-    integer, intent(in) :: ipiv(:)
-    integer, intent(in) :: n
-    type(permutation_t) :: p
-    integer :: i, tmp
+    integer(ilp), intent(in) :: ipiv(:)
+    integer(ilp), intent(in) :: n
+    type(permutation_type) :: p
+    integer(ilp) :: i, tmp
 
     p%n = n
     p%finalized = .true.
@@ -168,8 +169,8 @@ contains
 
   !> Convert nonfinalized (swap) to finalized (direct) permutation
   function permutation_finalize(p) result(fp)
-    type(permutation_t), intent(in) :: p
-    type(permutation_t) :: fp
+    type(permutation_type), intent(in) :: p
+    type(permutation_type) :: fp
 
     if (p%finalized) then
       fp = p
@@ -182,9 +183,9 @@ contains
   !> Compose two finalized permutations: result = p1 * p2
   !> Apply p2 first, then p1
   function permutation_compose(p1, p2) result(p)
-    type(permutation_t), intent(in) :: p1, p2
-    type(permutation_t) :: p
-    integer :: i
+    type(permutation_type), intent(in) :: p1, p2
+    type(permutation_type) :: p
+    integer(ilp) :: i
 
     p%n = p1%n
     p%finalized = .true.
@@ -197,15 +198,15 @@ contains
 
   !> ----- Templated operations for multiple real kinds -----
 
-#:for rk, rt in REAL_KINDS_TYPES
+#:for k1, t1 in RCI_KINDS_TYPES
 
   !> Apply permutation from the left: B = P * A (row permutation)
   !> Row i of B comes from row perm(i) of A
-  subroutine permutation_apply_left_${rk}$ (p, A, B)
-    type(permutation_t), intent(in) :: p
-    ${rt}$, intent(in) :: A(:,:)
-    ${rt}$, intent(out) :: B(:,:)
-    integer :: i
+  subroutine permutation_apply_left_${t1[0]}$${k1}$ (p, A, B)
+    type(permutation_type), intent(in) :: p
+    ${t1}$, intent(in) :: A(:,:)
+    ${t1}$, intent(out) :: B(:,:)
+    integer(ilp) :: i
 
     do i = 1, p%n
       B(i, :) = A(p%perm(i), :)
@@ -214,11 +215,11 @@ contains
 
   !> Apply permutation from the right: B = A * P (column permutation)
   !> Column perm(j) of B comes from column j of A
-  subroutine permutation_apply_right_${rk}$ (p, A, B)
-    type(permutation_t), intent(in) :: p
-    ${rt}$, intent(in) :: A(:,:)
-    ${rt}$, intent(out) :: B(:,:)
-    integer :: j
+  subroutine permutation_apply_right_${t1[0]}$${k1}$ (p, A, B)
+    type(permutation_type), intent(in) :: p
+    ${t1}$, intent(in) :: A(:,:)
+    ${t1}$, intent(out) :: B(:,:)
+    integer(ilp) :: j
 
     do j = 1, p%n
       B(:, p%perm(j)) = A(:, j)
@@ -227,11 +228,11 @@ contains
 
   !> Apply permutation to a vector: b = P * a
   !> b(i) = a(perm(i))
-  subroutine permutation_apply_vector_${rk}$ (p, a, b)
-    type(permutation_t), intent(in) :: p
-    ${rt}$, intent(in) :: a(:)
-    ${rt}$, intent(out) :: b(:)
-    integer :: i
+  subroutine permutation_apply_vector_${t1[0]}$${k1}$ (p, a, b)
+    type(permutation_type), intent(in) :: p
+    ${t1}$, intent(in) :: a(:)
+    ${t1}$, intent(out) :: b(:)
+    integer(ilp) :: i
 
     do i = 1, p%n
       b(i) = a(p%perm(i))
@@ -239,14 +240,14 @@ contains
   end subroutine
 
   !> Convert permutation to dense matrix (for testing / small cases)
-  subroutine permutation_to_matrix_${rk}$ (p, A)
-    type(permutation_t), intent(in) :: p
-    ${rt}$, intent(out) :: A(:,:)
-    integer :: i
+  subroutine permutation_to_matrix_${t1[0]}$${k1}$ (p, A)
+    type(permutation_type), intent(in) :: p
+    ${t1}$, intent(out) :: A(:,:)
+    integer(ilp) :: i
 
-    A = 0.0_${rk}$
+    A = 0
     do i = 1, p%n
-      A(i, p%perm(i)) = 1.0_${rk}$
+      A(i, p%perm(i)) = 1
     end do
   end subroutine
 
diff --git a/test/linalg/test_linalg_permutation.fypp b/test/linalg/test_linalg_permutation.fypp
index 918292e70..177818350 100644
--- a/test/linalg/test_linalg_permutation.fypp
+++ b/test/linalg/test_linalg_permutation.fypp
@@ -4,6 +4,7 @@
 module test_linalg_permutation
   use testdrive, only: new_unittest, unittest_type, error_type, check
   use stdlib_kinds, only: sp, dp, int32
+  use stdlib_linalg_constants, only: ilp
   use stdlib_linalg_permutation
   implicit none
   private
@@ -33,8 +34,8 @@ contains
   !> Test 1: Identity permutation
   subroutine test_identity(error)
     type(error_type), allocatable, intent(out) :: error
-    type(permutation_t) :: p
-    integer :: i
+    type(permutation_type) :: p
+    integer(ilp) :: i
 
     p = permutation_identity(5)
     call check(error, p%n == 5, "identity size")
@@ -49,7 +50,7 @@ contains
   !> Test 2: Valid permutation accepted
   subroutine test_is_valid(error)
     type(error_type), allocatable, intent(out) :: error
-    type(permutation_t) :: p
+    type(permutation_type) :: p
 
     p = permutation_init([3, 1, 4, 2])
     call check(error, permutation_is_valid(p), "valid permutation")
@@ -58,7 +59,7 @@ contains
   !> Test 3: Invalid permutations rejected
   subroutine test_invalid_rejects(error)
     type(error_type), allocatable, intent(out) :: error
-    type(permutation_t) :: p
+    type(permutation_type) :: p
 
     ! Duplicate entry
     p = permutation_init([1, 1, 3])
@@ -73,7 +74,7 @@ contains
   !> Test 4: Inversion correctness
   subroutine test_invert(error)
     type(error_type), allocatable, intent(out) :: error
-    type(permutation_t) :: p, inv
+    type(permutation_type) :: p, inv
 
     ! perm = [3,1,2] means: pos1->3, pos2->1, pos3->2
     ! inverse: pos3->1, pos1->2, pos2->3 => inv = [2,3,1]
@@ -89,8 +90,8 @@ contains
   !> Test 5: P * P^-1 = Identity (involution property)
   subroutine test_invert_involution(error)
     type(error_type), allocatable, intent(out) :: error
-    type(permutation_t) :: p, inv, composed
-    integer :: i
+    type(permutation_type) :: p, inv, composed
+    integer(ilp) :: i
 
     p = permutation_init([3, 1, 4, 2])
     inv = permutation_invert(p)
@@ -104,7 +105,7 @@ contains
   !> Test 6: Left multiplication (row permutation)
   subroutine test_apply_left(error)
     type(error_type), allocatable, intent(out) :: error
-    type(permutation_t) :: p
+    type(permutation_type) :: p
     real(dp) :: A(3,2), B(3,2)
 
     ! A = [[1,2],[3,4],[5,6]]
@@ -127,7 +128,7 @@ contains
   !> Test 7: Right multiplication (column permutation)
   subroutine test_apply_right(error)
     type(error_type), allocatable, intent(out) :: error
-    type(permutation_t) :: p
+    type(permutation_type) :: p
     real(dp) :: A(2,3), B(2,3)
 
     A(1,:) = [1.0_dp, 2.0_dp, 3.0_dp]
@@ -149,7 +150,7 @@ contains
   !> Test 8: Vector permutation
   subroutine test_apply_vector(error)
     type(error_type), allocatable, intent(out) :: error
-    type(permutation_t) :: p
+    type(permutation_type) :: p
     real(dp) :: a(4), b(4)
 
     a = [10.0_dp, 20.0_dp, 30.0_dp, 40.0_dp]
@@ -169,7 +170,7 @@ contains
   !> Test 9: Composition P1 * P2
   subroutine test_compose(error)
     type(error_type), allocatable, intent(out) :: error
-    type(permutation_t) :: p1, p2, composed
+    type(permutation_type) :: p1, p2, composed
     real(dp) :: a(3), b1(3), b2(3), b_composed(3)
 
     p1 = permutation_init([2, 3, 1])
@@ -192,7 +193,7 @@ contains
   !> Test 10: LAPACK pivot vector conversion
   subroutine test_from_lapack(error)
     type(error_type), allocatable, intent(out) :: error
-    type(permutation_t) :: p
+    type(permutation_type) :: p
 
     ! LAPACK pivot: [2, 3, 3]
     ! Step 1: swap rows 1,2 -> [2,1,3]
@@ -210,7 +211,7 @@ contains
   !> Test 11: Dense matrix matches permutation operations
   subroutine test_to_matrix(error)
     type(error_type), allocatable, intent(out) :: error
-    type(permutation_t) :: p
+    type(permutation_type) :: p
     real(dp) :: P_dense(3,3)
 
     p = permutation_init([3, 1, 2])
@@ -229,9 +230,9 @@ contains
   !> Test 12: Identity permutation is a no-op
   subroutine test_identity_noop(error)
     type(error_type), allocatable, intent(out) :: error
-    type(permutation_t) :: p
+    type(permutation_type) :: p
     real(dp) :: A(3,3), B(3,3)
-    integer :: i, j
+    integer(ilp) :: i, j
 
     p = permutation_identity(3)
     do i = 1, 3
@@ -241,12 +242,8 @@ contains
     end do
 
     call permutation_apply_left(p, A, B)
-    do i = 1, 3
-      do j = 1, 3
-        call check(error, abs(A(i,j) - B(i,j)) < 1.0e-10_dp, "identity noop")
-        if (allocated(error)) return
-      end do
-    end do
+    call check(error, all(abs(A - B) < 1.0e-10_dp), "identity noop")
+    if (allocated(error)) return
   end subroutine
 
 end module test_linalg_permutation

From 5dd11ea7e58a2d4a5c0cb66da189db9fbcbe4da6 Mon Sep 17 00:00:00 2001
From: Sudhanshu Thakur <sudhanshukit2021@gmail.com>
Date: Wed, 25 Mar 2026 23:51:04 +0530
Subject: [PATCH 04/10] fix(linalg): address Copilot AI semantic and
 mathematical review feedback

---
 src/linalg/stdlib_linalg_permutation.fypp | 35 ++++++++++++++++-------
 test/linalg/test_linalg_permutation.fypp  |  4 +--
 2 files changed, 26 insertions(+), 13 deletions(-)

diff --git a/src/linalg/stdlib_linalg_permutation.fypp b/src/linalg/stdlib_linalg_permutation.fypp
index 92825e8ea..8c641a362 100644
--- a/src/linalg/stdlib_linalg_permutation.fypp
+++ b/src/linalg/stdlib_linalg_permutation.fypp
@@ -6,8 +6,9 @@
 !> run in O(n) time per column/row rather than O(n²).
 !>
 !> Two formats are supported:
-!>   - Finalized: perm(i) = j means row i maps to row j (direct mapping)
-!>   - LAPACK pivot: ipiv(i) = j means swap rows i and j at step i
+!>   - Finalized: perm(i) = j means that output row i is taken from input row j
+!>     (i.e., after applying the permutation on the left, B(i,:) = A(perm(i),:))
+!>   - LAPACK pivot: ipiv(i) = j means swap rows i and j at elimination step i
 !>
 !> Reference: Issue #891 (https://github.com/fortran-lang/stdlib/issues/891)
 
@@ -108,12 +109,19 @@ contains
     if (p%n < 1) return
     if (size(p%perm) /= p%n) return
 
-    allocate(seen(p%n), source=.false.)
-    do i = 1, p%n
-      if (p%perm(i) < 1 .or. p%perm(i) > p%n) return
-      if (seen(p%perm(i))) return
-      seen(p%perm(i)) = .true.
-    end do
+    if (p%finalized) then
+      allocate(seen(p%n), source=.false.)
+      do i = 1, p%n
+        if (p%perm(i) < 1 .or. p%perm(i) > p%n) return
+        if (seen(p%perm(i))) return
+        seen(p%perm(i)) = .true.
+      end do
+    else
+      ! For LAPACK pivots, just check bounds (duplicates are valid for swapping)
+      do i = 1, size(p%perm)
+        if (p%perm(i) < 1 .or. p%perm(i) > p%n) return
+      end do
+    end if
     valid = .true.
   end function
 
@@ -158,8 +166,9 @@ contains
     end do
 
     ! Replay swap sequence forward
-    do i = 1, min(n, size(ipiv))
-      if (ipiv(i) /= i .and. ipiv(i) >= 1 .and. ipiv(i) <= n) then
+    do i = 1, min(n, int(size(ipiv), ilp))
+      if (ipiv(i) < 1 .or. ipiv(i) > n) error stop "permutation_from_lapack: Invalid LAPACK pivot index out of bounds"
+      if (ipiv(i) /= i) then
         tmp = p%perm(i)
         p%perm(i) = p%perm(ipiv(i))
         p%perm(ipiv(i)) = tmp
@@ -192,7 +201,7 @@ contains
     allocate(p%perm(p%n))
 
     do i = 1, p%n
-      p%perm(i) = p1%perm(p2%perm(i))
+      p%perm(i) = p2%perm(p1%perm(i))
     end do
   end function
 
@@ -208,6 +217,7 @@ contains
     ${t1}$, intent(out) :: B(:,:)
     integer(ilp) :: i
 
+    if (.not. p%finalized) error stop "Permutation must be finalized before applying."
     do i = 1, p%n
       B(i, :) = A(p%perm(i), :)
     end do
@@ -221,6 +231,7 @@ contains
     ${t1}$, intent(out) :: B(:,:)
     integer(ilp) :: j
 
+    if (.not. p%finalized) error stop "Permutation must be finalized before applying."
     do j = 1, p%n
       B(:, p%perm(j)) = A(:, j)
     end do
@@ -234,6 +245,7 @@ contains
     ${t1}$, intent(out) :: b(:)
     integer(ilp) :: i
 
+    if (.not. p%finalized) error stop "Permutation must be finalized before applying."
     do i = 1, p%n
       b(i) = a(p%perm(i))
     end do
@@ -245,6 +257,7 @@ contains
     ${t1}$, intent(out) :: A(:,:)
     integer(ilp) :: i
 
+    if (.not. p%finalized) error stop "Permutation must be finalized to convert to matrix."
     A = 0
     do i = 1, p%n
       A(i, p%perm(i)) = 1
diff --git a/test/linalg/test_linalg_permutation.fypp b/test/linalg/test_linalg_permutation.fypp
index 177818350..8953dd6c3 100644
--- a/test/linalg/test_linalg_permutation.fypp
+++ b/test/linalg/test_linalg_permutation.fypp
@@ -76,8 +76,8 @@ contains
     type(error_type), allocatable, intent(out) :: error
     type(permutation_type) :: p, inv
 
-    ! perm = [3,1,2] means: pos1->3, pos2->1, pos3->2
-    ! inverse: pos3->1, pos1->2, pos2->3 => inv = [2,3,1]
+    ! perm = [3,1,2] means: dest1 <- src3, dest2 <- src1, dest3 <- src2
+    ! inverse: dest3 <- src1, dest1 <- src2, dest2 <- src3 => inv = [2,3,1]
     p = permutation_init([3, 1, 2])
     inv = permutation_invert(p)
     call check(error, inv%perm(1) == 2, "inv(1)")

From 02d9f773158b71e90ef19ef16bdca7342337d931 Mon Sep 17 00:00:00 2001
From: Sudhanshu Thakur <sudhanshukit2021@gmail.com>
Date: Wed, 25 Mar 2026 23:52:54 +0530
Subject: [PATCH 05/10] test(linalg): add non-inverse composition test and test
 complex/integer support

---
 test/linalg/test_linalg_permutation.fypp | 34 +++++++++++++++++++-----
 1 file changed, 27 insertions(+), 7 deletions(-)

diff --git a/test/linalg/test_linalg_permutation.fypp b/test/linalg/test_linalg_permutation.fypp
index 8953dd6c3..ae1651bd2 100644
--- a/test/linalg/test_linalg_permutation.fypp
+++ b/test/linalg/test_linalg_permutation.fypp
@@ -171,23 +171,19 @@ contains
   subroutine test_compose(error)
     type(error_type), allocatable, intent(out) :: error
     type(permutation_type) :: p1, p2, composed
-    real(dp) :: a(3), b1(3), b2(3), b_composed(3)
+    real(dp) :: a(4), b1(4), b2(4), b_composed(4)
 
     p1 = permutation_init([2, 3, 1])
     p2 = permutation_init([3, 1, 2])
     composed = permutation_compose(p1, p2)
 
     ! Verify: composed applied to vector = p1 applied to (p2 applied to vector)
-    a = [10.0_dp, 20.0_dp, 30.0_dp]
+    a = [10.0_dp, 20.0_dp, 30.0_dp, 40.0_dp]
     call permutation_apply_vector(p2, a, b1)       ! b1 = P2 * a
     call permutation_apply_vector(p1, b1, b2)      ! b2 = P1 * P2 * a
     call permutation_apply_vector(composed, a, b_composed) ! should equal b2
 
-    call check(error, abs(b2(1) - b_composed(1)) < 1.0e-10_dp, "compose(1)")
-    if (allocated(error)) return
-    call check(error, abs(b2(2) - b_composed(2)) < 1.0e-10_dp, "compose(2)")
-    if (allocated(error)) return
-    call check(error, abs(b2(3) - b_composed(3)) < 1.0e-10_dp, "compose(3)")
+    call check(error, all(abs(b2 - b_composed) < 1.0e-10_dp), "compose_matches_explicit")
   end subroutine
 
   !> Test 10: LAPACK pivot vector conversion
@@ -246,6 +242,30 @@ contains
     if (allocated(error)) return
   end subroutine
 
+
+  !> Test permutation with different generic data types (complex and integer)
+  subroutine test_types_support(error)
+    type(error_type), allocatable, intent(out) :: error
+    type(permutation_type) :: p
+    complex(dp) :: c_in(3), c_out(3)
+    integer(int32) :: i_in(3), i_out(3)
+
+    p = permutation_init([3_ilp, 1_ilp, 2_ilp])
+
+    ! Complex support
+    c_in = [(1.0_dp, -1.0_dp), (2.0_dp, -2.0_dp), (3.0_dp, -3.0_dp)]
+    call permutation_apply_vector(p, c_in, c_out)
+    call check(error, abs(c_out(1) - c_in(3)) < 1.0e-10_dp, "complex")
+    if (allocated(error)) return
+
+    ! Integer support
+    i_in = [10_int32, 20_int32, 30_int32]
+    call permutation_apply_vector(p, i_in, i_out)
+    call check(error, i_out(1) == i_in(3), "integer")
+    if (allocated(error)) return
+
+  end subroutine test_types_support
+
 end module test_linalg_permutation
 
 program test_permutation

From 894733c9c5f11d9544f63183ffe6787c567b974d Mon Sep 17 00:00:00 2001
From: Sudhanshu Thakur <sudhanshukit2021@gmail.com>
Date: Thu, 26 Mar 2026 15:25:05 +0530
Subject: [PATCH 06/10] test(linalg): fix uninitialized array in compose test

---
 test/linalg/test_linalg_permutation.fypp | 6 ++++--
 1 file changed, 4 insertions(+), 2 deletions(-)

diff --git a/test/linalg/test_linalg_permutation.fypp b/test/linalg/test_linalg_permutation.fypp
index ae1651bd2..9f841b824 100644
--- a/test/linalg/test_linalg_permutation.fypp
+++ b/test/linalg/test_linalg_permutation.fypp
@@ -173,12 +173,14 @@ contains
     type(permutation_type) :: p1, p2, composed
     real(dp) :: a(4), b1(4), b2(4), b_composed(4)
 
-    p1 = permutation_init([2, 3, 1])
-    p2 = permutation_init([3, 1, 2])
+    p1 = permutation_init([2_ilp, 3_ilp, 4_ilp, 1_ilp])
+    p2 = permutation_init([2_ilp, 1_ilp, 4_ilp, 3_ilp])
     composed = permutation_compose(p1, p2)
 
     ! Verify: composed applied to vector = p1 applied to (p2 applied to vector)
     a = [10.0_dp, 20.0_dp, 30.0_dp, 40.0_dp]
+    b1 = 0.0_dp; b2 = 0.0_dp; b_composed = 0.0_dp ! Initialize to avoid garbage
+
     call permutation_apply_vector(p2, a, b1)       ! b1 = P2 * a
     call permutation_apply_vector(p1, b1, b2)      ! b2 = P1 * P2 * a
     call permutation_apply_vector(composed, a, b_composed) ! should equal b2

From 83334b0a63fdc8fd1bdf60a4148bb5c323046618 Mon Sep 17 00:00:00 2001
From: Sudhanshu Thakur <sudhanshukit2021@gmail.com>
Date: Sun, 29 Mar 2026 22:49:33 +0530
Subject: [PATCH 07/10] fix fypp macro indents and use array checks in tests

---
 src/linalg/stdlib_linalg_permutation.fypp | 28 ++++----
 test/linalg/test_linalg_permutation.fypp  | 82 ++++++++---------------
 2 files changed, 42 insertions(+), 68 deletions(-)

diff --git a/src/linalg/stdlib_linalg_permutation.fypp b/src/linalg/stdlib_linalg_permutation.fypp
index 8c641a362..b600736cd 100644
--- a/src/linalg/stdlib_linalg_permutation.fypp
+++ b/src/linalg/stdlib_linalg_permutation.fypp
@@ -43,27 +43,27 @@ module stdlib_linalg_permutation
 
   ! Generic interfaces for multi-type support
   interface permutation_apply_left
-#:for k1, t1 in RCI_KINDS_TYPES
-      module procedure permutation_apply_left_${t1[0]}$${k1}$
-#:endfor
+    #:for k1, t1 in RCI_KINDS_TYPES
+    module procedure permutation_apply_left_${t1[0]}$${k1}$
+    #:endfor
   end interface
 
   interface permutation_apply_right
-#:for k1, t1 in RCI_KINDS_TYPES
-      module procedure permutation_apply_right_${t1[0]}$${k1}$
-#:endfor
+    #:for k1, t1 in RCI_KINDS_TYPES
+    module procedure permutation_apply_right_${t1[0]}$${k1}$
+    #:endfor
   end interface
 
   interface permutation_apply_vector
-#:for k1, t1 in RCI_KINDS_TYPES
-      module procedure permutation_apply_vector_${t1[0]}$${k1}$
-#:endfor
+    #:for k1, t1 in RCI_KINDS_TYPES
+    module procedure permutation_apply_vector_${t1[0]}$${k1}$
+    #:endfor
   end interface
 
   interface permutation_to_matrix
-#:for k1, t1 in RCI_KINDS_TYPES
-      module procedure permutation_to_matrix_${t1[0]}$${k1}$
-#:endfor
+    #:for k1, t1 in RCI_KINDS_TYPES
+    module procedure permutation_to_matrix_${t1[0]}$${k1}$
+    #:endfor
   end interface
 
 contains
@@ -207,7 +207,7 @@ contains
 
   !> ----- Templated operations for multiple real kinds -----
 
-#:for k1, t1 in RCI_KINDS_TYPES
+  #:for k1, t1 in RCI_KINDS_TYPES
 
   !> Apply permutation from the left: B = P * A (row permutation)
   !> Row i of B comes from row perm(i) of A
@@ -264,6 +264,6 @@ contains
     end do
   end subroutine
 
-#:endfor
+  #:endfor
 
 end module stdlib_linalg_permutation
diff --git a/test/linalg/test_linalg_permutation.fypp b/test/linalg/test_linalg_permutation.fypp
index 9f841b824..eefad5271 100644
--- a/test/linalg/test_linalg_permutation.fypp
+++ b/test/linalg/test_linalg_permutation.fypp
@@ -35,15 +35,11 @@ contains
   subroutine test_identity(error)
     type(error_type), allocatable, intent(out) :: error
     type(permutation_type) :: p
-    integer(ilp) :: i
-
     p = permutation_identity(5)
     call check(error, p%n == 5, "identity size")
     if (allocated(error)) return
-    do i = 1, 5
-      call check(error, p%perm(i) == i, "identity element")
-      if (allocated(error)) return
-    end do
+    call check(error, all(p%perm == [1_ilp, 2_ilp, 3_ilp, 4_ilp, 5_ilp]), "identity elements")
+    if (allocated(error)) return
     call check(error, permutation_is_valid(p), "identity valid")
   end subroutine
 
@@ -80,71 +76,61 @@ contains
     ! inverse: dest3 <- src1, dest1 <- src2, dest2 <- src3 => inv = [2,3,1]
     p = permutation_init([3, 1, 2])
     inv = permutation_invert(p)
-    call check(error, inv%perm(1) == 2, "inv(1)")
-    if (allocated(error)) return
-    call check(error, inv%perm(2) == 3, "inv(2)")
-    if (allocated(error)) return
-    call check(error, inv%perm(3) == 1, "inv(3)")
+    call check(error, all(inv%perm == [2_ilp, 3_ilp, 1_ilp]), "inv(1:3)")
   end subroutine
 
   !> Test 5: P * P^-1 = Identity (involution property)
   subroutine test_invert_involution(error)
     type(error_type), allocatable, intent(out) :: error
     type(permutation_type) :: p, inv, composed
-    integer(ilp) :: i
-
-    p = permutation_init([3, 1, 4, 2])
+    p = permutation_init([3_ilp, 1_ilp, 4_ilp, 2_ilp])
     inv = permutation_invert(p)
     composed = permutation_compose(p, inv)
-    do i = 1, 4
-      call check(error, composed%perm(i) == i, "P*P_inv = I")
-      if (allocated(error)) return
-    end do
+    call check(error, all(composed%perm == [1_ilp, 2_ilp, 3_ilp, 4_ilp]), "P*P_inv = I")
   end subroutine
 
   !> Test 6: Left multiplication (row permutation)
   subroutine test_apply_left(error)
     type(error_type), allocatable, intent(out) :: error
     type(permutation_type) :: p
-    real(dp) :: A(3,2), B(3,2)
+    real(dp) :: A(3,2), B(3,2), B_ref(3,2)
 
     ! A = [[1,2],[3,4],[5,6]]
     A(1,:) = [1.0_dp, 2.0_dp]
     A(2,:) = [3.0_dp, 4.0_dp]
     A(3,:) = [5.0_dp, 6.0_dp]
 
+    ! B should be [[5,6],[1,2],[3,4]]
+    B_ref(1,:) = [5.0_dp, 6.0_dp]
+    B_ref(2,:) = [1.0_dp, 2.0_dp]
+    B_ref(3,:) = [3.0_dp, 4.0_dp]
+
     ! perm = [3,1,2]: row1 <- row3, row2 <- row1, row3 <- row2
     p = permutation_init([3, 1, 2])
     call permutation_apply_left(p, A, B)
 
-    ! B should be [[5,6],[1,2],[3,4]]
-    call check(error, abs(B(1,1) - 5.0_dp) < 1.0e-10_dp, "left(1,1)")
-    if (allocated(error)) return
-    call check(error, abs(B(2,1) - 1.0_dp) < 1.0e-10_dp, "left(2,1)")
-    if (allocated(error)) return
-    call check(error, abs(B(3,1) - 3.0_dp) < 1.0e-10_dp, "left(3,1)")
+    call check(error, all(abs(B - B_ref) < 1.0e-10_dp), "apply_left whole array")
   end subroutine
 
   !> Test 7: Right multiplication (column permutation)
   subroutine test_apply_right(error)
     type(error_type), allocatable, intent(out) :: error
     type(permutation_type) :: p
-    real(dp) :: A(2,3), B(2,3)
+    real(dp) :: A(2,3), B(2,3), B_ref(2,3)
 
     A(1,:) = [1.0_dp, 2.0_dp, 3.0_dp]
     A(2,:) = [4.0_dp, 5.0_dp, 6.0_dp]
 
+    ! B(:,2) = A(:,1) = [1,4], B(:,3) = A(:,2) = [2,5], B(:,1) = A(:,3) = [3,6]
+    B_ref(:,1) = [3.0_dp, 6.0_dp]
+    B_ref(:,2) = [1.0_dp, 4.0_dp]
+    B_ref(:,3) = [2.0_dp, 5.0_dp]
+
     ! perm = [2,3,1]: col_perm(2) <- col1, col_perm(3) <- col2, col_perm(1) <- col3
     p = permutation_init([2, 3, 1])
     call permutation_apply_right(p, A, B)
 
-    ! B(:, perm(j)) = A(:, j)
-    ! B(:,2) = A(:,1) = [1,4], B(:,3) = A(:,2) = [2,5], B(:,1) = A(:,3) = [3,6]
-    call check(error, abs(B(1,1) - 3.0_dp) < 1.0e-10_dp, "right(1,1)")
-    if (allocated(error)) return
-    call check(error, abs(B(1,2) - 1.0_dp) < 1.0e-10_dp, "right(1,2)")
-    if (allocated(error)) return
-    call check(error, abs(B(1,3) - 2.0_dp) < 1.0e-10_dp, "right(1,3)")
+    call check(error, all(abs(B - B_ref) < 1.0e-10_dp), "apply_right whole array")
   end subroutine
 
   !> Test 8: Vector permutation
@@ -158,13 +144,7 @@ contains
     call permutation_apply_vector(p, a, b)
 
     ! b(1) = a(3) = 30, b(2) = a(1) = 10, b(3) = a(4) = 40, b(4) = a(2) = 20
-    call check(error, abs(b(1) - 30.0_dp) < 1.0e-10_dp, "vec(1)")
-    if (allocated(error)) return
-    call check(error, abs(b(2) - 10.0_dp) < 1.0e-10_dp, "vec(2)")
-    if (allocated(error)) return
-    call check(error, abs(b(3) - 40.0_dp) < 1.0e-10_dp, "vec(3)")
-    if (allocated(error)) return
-    call check(error, abs(b(4) - 20.0_dp) < 1.0e-10_dp, "vec(4)")
+    call check(error, all(abs(b - [30.0_dp, 10.0_dp, 40.0_dp, 20.0_dp]) < 1.0e-10_dp), "apply_vector whole array")
   end subroutine
 
   !> Test 9: Composition P1 * P2
@@ -199,30 +179,24 @@ contains
     ! Step 3: swap rows 3,3 -> [2,3,1] (no-op)
     p = permutation_from_lapack([2, 3, 3], 3)
 
-    call check(error, p%perm(1) == 2, "lapack(1)")
-    if (allocated(error)) return
-    call check(error, p%perm(2) == 3, "lapack(2)")
-    if (allocated(error)) return
-    call check(error, p%perm(3) == 1, "lapack(3)")
+    call check(error, all(p%perm == [2_ilp, 3_ilp, 1_ilp]), "from_lapack whole array")
   end subroutine
 
   !> Test 11: Dense matrix matches permutation operations
   subroutine test_to_matrix(error)
     type(error_type), allocatable, intent(out) :: error
     type(permutation_type) :: p
-    real(dp) :: P_dense(3,3)
+    real(dp) :: P_dense(3,3), P_ref(3,3)
 
     p = permutation_init([3, 1, 2])
     call permutation_to_matrix(p, P_dense)
 
-    ! Row 1, col 3 should be 1 (perm(1)=3)
-    call check(error, abs(P_dense(1,3) - 1.0_dp) < 1.0e-10_dp, "dense(1,3)")
-    if (allocated(error)) return
-    ! Row 2, col 1 should be 1 (perm(2)=1)
-    call check(error, abs(P_dense(2,1) - 1.0_dp) < 1.0e-10_dp, "dense(2,1)")
-    if (allocated(error)) return
-    ! Row 1, col 1 should be 0
-    call check(error, abs(P_dense(1,1)) < 1.0e-10_dp, "dense(1,1)=0")
+    P_ref = 0.0_dp
+    P_ref(1,3) = 1.0_dp
+    P_ref(2,1) = 1.0_dp
+    P_ref(3,2) = 1.0_dp
+
+    call check(error, all(abs(P_dense - P_ref) < 1.0e-10_dp), "to_matrix whole array")
   end subroutine
 
   !> Test 12: Identity permutation is a no-op

From 4780ae2ec8048d3789e136a88c5f160cc7b119e1 Mon Sep 17 00:00:00 2001
From: Sudhanshu Thakur <sudhanshukit2021@gmail.com>
Date: Mon, 30 Mar 2026 23:01:29 +0530
Subject: [PATCH 08/10] convert permutation to submodule of stdlib_linalg

---
 src/linalg/CMakeLists.txt                 |   2 +-
 src/linalg/stdlib_linalg.fypp             | 139 ++++++++++++++++++++
 src/linalg/stdlib_linalg_permutation.fypp | 153 +++++++++-------------
 test/linalg/test_linalg_permutation.fypp  |   5 +-
 4 files changed, 208 insertions(+), 91 deletions(-)

diff --git a/src/linalg/CMakeLists.txt b/src/linalg/CMakeLists.txt
index 778ed3353..31a76ef94 100644
--- a/src/linalg/CMakeLists.txt
+++ b/src/linalg/CMakeLists.txt
@@ -27,4 +27,4 @@ set(linalg_f90Files
 
 configure_stdlib_target(${PROJECT_NAME}_linalg linalg_f90Files linalg_fppFiles linalg_cppFiles)
 
-target_link_libraries(${PROJECT_NAME}_linalg PUBLIC ${PROJECT_NAME}_core ${PROJECT_NAME}_linalg_core ${PROJECT_NAME}_constants ${PROJECT_NAME}_lapack ${PROJECT_NAME}_lapack_extended ${PROJECT_NAME}_sorting ${PROJECT_NAME}_intrinsics)
+target_link_libraries(${PROJECT_NAME}_linalg PUBLIC ${PROJECT_NAME}_core ${PROJECT_NAME}_linalg_core ${PROJECT_NAME}_constants ${PROJECT_NAME}_lapack ${PROJECT_NAME}_lapack_extended ${PROJECT_NAME}_sorting ${PROJECT_NAME}_intrinsics ${PROJECT_NAME}_math)
diff --git a/src/linalg/stdlib_linalg.fypp b/src/linalg/stdlib_linalg.fypp
index 5c1b3b7e1..30a9f98ef 100644
--- a/src/linalg/stdlib_linalg.fypp
+++ b/src/linalg/stdlib_linalg.fypp
@@ -70,6 +70,19 @@ module stdlib_linalg
   public :: is_hermitian
   public :: is_triangular
   public :: is_hessenberg
+
+  ! Permutation matrix operations
+  public :: permutation_type
+  public :: permutation_identity
+  public :: permutation_init
+  public :: permutation_is_valid
+  public :: permutation_invert
+  public :: permutation_compose
+  public :: permutation_from_lapack
+  public :: permutation_to_matrix
+  public :: permutation_apply_left
+  public :: permutation_apply_right
+  public :: permutation_apply_vector
   
   ! Export linalg error handling
   public :: linalg_state_type, linalg_error_handling
@@ -2103,6 +2116,132 @@ module stdlib_linalg
     end subroutine stdlib_linalg_${ri}$_expm
     #:endfor
   end interface matrix_exp
+
+  !> Permutation type: stores permutation as integer vector
+  type :: permutation_type
+    integer(ilp), allocatable :: perm(:)
+    integer(ilp) :: n = 0
+    logical :: finalized = .true.
+  end type
+
+  interface permutation_identity
+    !! version: experimental
+    !!
+    !! Creates an identity permutation of size n
+    module function permutation_identity(n) result(p)
+      integer(ilp), intent(in) :: n
+      type(permutation_type) :: p
+    end function
+  end interface
+
+  interface permutation_init
+    !! version: experimental
+    !!
+    !! Creates a permutation from a given integer vector
+    module function permutation_init(perm_vec, finalized) result(p)
+      integer(ilp), intent(in) :: perm_vec(:)
+      logical, intent(in), optional :: finalized
+      type(permutation_type) :: p
+    end function
+  end interface
+
+  interface permutation_is_valid
+    !! version: experimental
+    !!
+    !! Checks if a permutation is valid (each integer 1..n appears exactly once)
+    module function permutation_is_valid(p) result(valid)
+      type(permutation_type), intent(in) :: p
+      logical :: valid
+    end function
+  end interface
+
+  interface permutation_invert
+    !! version: experimental
+    !!
+    !! Inverts a finalized permutation
+    module function permutation_invert(p) result(inv)
+      type(permutation_type), intent(in) :: p
+      type(permutation_type) :: inv
+    end function
+  end interface
+
+  interface permutation_compose
+    !! version: experimental
+    !!
+    !! Composes two finalized permutations: result = p1 * p2
+    module function permutation_compose(p1, p2) result(p)
+      type(permutation_type), intent(in) :: p1, p2
+      type(permutation_type) :: p
+    end function
+  end interface
+
+  interface permutation_from_lapack
+    !! version: experimental
+    !!
+    !! Converts LAPACK pivot vector to finalized permutation
+    module function permutation_from_lapack(ipiv, n, err) result(p)
+      integer(ilp), intent(in) :: ipiv(:)
+      integer(ilp), intent(in) :: n
+      type(linalg_state_type), intent(out), optional :: err
+      type(permutation_type) :: p
+    end function
+  end interface
+
+  interface permutation_apply_left
+    !! version: experimental
+    !!
+    !! Applies permutation from the left: B = P * A (row permutation)
+    #:for k1, t1 in RCI_KINDS_TYPES
+    module subroutine permutation_apply_left_${t1[0]}$${k1}$ (p, A, B, err)
+      type(permutation_type), intent(in) :: p
+      ${t1}$, intent(in) :: A(:,:)
+      ${t1}$, intent(out) :: B(:,:)
+      type(linalg_state_type), intent(out), optional :: err
+    end subroutine
+    #:endfor
+  end interface
+
+  interface permutation_apply_right
+    !! version: experimental
+    !!
+    !! Applies permutation from the right: B = A * P (column permutation)
+    #:for k1, t1 in RCI_KINDS_TYPES
+    module subroutine permutation_apply_right_${t1[0]}$${k1}$ (p, A, B, err)
+      type(permutation_type), intent(in) :: p
+      ${t1}$, intent(in) :: A(:,:)
+      ${t1}$, intent(out) :: B(:,:)
+      type(linalg_state_type), intent(out), optional :: err
+    end subroutine
+    #:endfor
+  end interface
+
+  interface permutation_apply_vector
+    !! version: experimental
+    !!
+    !! Applies permutation to a vector: b = P * a
+    #:for k1, t1 in RCI_KINDS_TYPES
+    module subroutine permutation_apply_vector_${t1[0]}$${k1}$ (p, a, b, err)
+      type(permutation_type), intent(in) :: p
+      ${t1}$, intent(in) :: a(:)
+      ${t1}$, intent(out) :: b(:)
+      type(linalg_state_type), intent(out), optional :: err
+    end subroutine
+    #:endfor
+  end interface
+
+  interface permutation_to_matrix
+    !! version: experimental
+    !!
+    !! Converts permutation to dense matrix
+    #:for k1, t1 in RCI_KINDS_TYPES
+    module subroutine permutation_to_matrix_${t1[0]}$${k1}$ (p, A, err)
+      type(permutation_type), intent(in) :: p
+      ${t1}$, intent(out) :: A(:,:)
+      type(linalg_state_type), intent(out), optional :: err
+    end subroutine
+    #:endfor
+  end interface
+
 contains
 
 
diff --git a/src/linalg/stdlib_linalg_permutation.fypp b/src/linalg/stdlib_linalg_permutation.fypp
index b600736cd..890213cc6 100644
--- a/src/linalg/stdlib_linalg_permutation.fypp
+++ b/src/linalg/stdlib_linalg_permutation.fypp
@@ -16,89 +16,37 @@
 #:set RC_KINDS_TYPES = REAL_KINDS_TYPES + CMPLX_KINDS_TYPES
 #:set RCI_KINDS_TYPES = RC_KINDS_TYPES + INT_KINDS_TYPES
 
-module stdlib_linalg_permutation
-  use stdlib_kinds, only: xdp, int8, int16, int32, int64
-  use stdlib_linalg_constants, only: sp, dp, qp, lk, ilp
+submodule (stdlib_linalg) stdlib_linalg_permutation
+  use stdlib_math, only: arange
+  use stdlib_linalg_state, only: LINALG_VALUE_ERROR
   implicit none
-  private
-
-  !> Permutation type: stores permutation as integer vector
-  type, public :: permutation_type
-    integer(ilp), allocatable :: perm(:)
-    integer(ilp) :: n = 0
-    logical :: finalized = .true.
-  end type
-
-  ! Public procedures
-  public :: permutation_identity
-  public :: permutation_init
-  public :: permutation_is_valid
-  public :: permutation_invert
-  public :: permutation_compose
-  public :: permutation_from_lapack
-  public :: permutation_to_matrix
-  public :: permutation_apply_left
-  public :: permutation_apply_right
-  public :: permutation_apply_vector
-
-  ! Generic interfaces for multi-type support
-  interface permutation_apply_left
-    #:for k1, t1 in RCI_KINDS_TYPES
-    module procedure permutation_apply_left_${t1[0]}$${k1}$
-    #:endfor
-  end interface
-
-  interface permutation_apply_right
-    #:for k1, t1 in RCI_KINDS_TYPES
-    module procedure permutation_apply_right_${t1[0]}$${k1}$
-    #:endfor
-  end interface
-
-  interface permutation_apply_vector
-    #:for k1, t1 in RCI_KINDS_TYPES
-    module procedure permutation_apply_vector_${t1[0]}$${k1}$
-    #:endfor
-  end interface
-
-  interface permutation_to_matrix
-    #:for k1, t1 in RCI_KINDS_TYPES
-    module procedure permutation_to_matrix_${t1[0]}$${k1}$
-    #:endfor
-  end interface
 
 contains
 
   !> Create an identity permutation of size n
-  function permutation_identity(n) result(p)
+  module function permutation_identity(n) result(p)
     integer(ilp), intent(in) :: n
     type(permutation_type) :: p
-    integer(ilp) :: i
 
     p%n = n
     p%finalized = .true.
     allocate(p%perm(n))
-    do i = 1, n
-      p%perm(i) = i
-    end do
+    p%perm = arange(1_ilp, n)
   end function
 
   !> Create a permutation from a given integer vector
-  function permutation_init(perm_vec, finalized) result(p)
+  module function permutation_init(perm_vec, finalized) result(p)
     integer(ilp), intent(in) :: perm_vec(:)
     logical, intent(in), optional :: finalized
     type(permutation_type) :: p
 
     p%n = size(perm_vec)
     allocate(p%perm(p%n), source=perm_vec)
-    if (present(finalized)) then
-      p%finalized = finalized
-    else
-      p%finalized = .true.
-    end if
+    p%finalized = optval(finalized, .true.)
   end function
 
   !> Check if a permutation is valid (each integer 1..n appears exactly once)
-  function permutation_is_valid(p) result(valid)
+  module function permutation_is_valid(p) result(valid)
     type(permutation_type), intent(in) :: p
     logical :: valid
     logical, allocatable :: seen(:)
@@ -127,7 +75,7 @@ contains
 
   !> Invert a finalized permutation
   !> If perm(i) = j, then inv(j) = i
-  function permutation_invert(p) result(inv)
+  module function permutation_invert(p) result(inv)
     type(permutation_type), intent(in) :: p
     type(permutation_type) :: inv, fp
     integer(ilp) :: i
@@ -150,10 +98,12 @@ contains
 
   !> Convert LAPACK pivot vector (swap sequence) to finalized permutation
   !> ipiv(i) = j means "swap rows i and j" applied at step i
-  function permutation_from_lapack(ipiv, n) result(p)
+  module function permutation_from_lapack(ipiv, n, err) result(p)
     integer(ilp), intent(in) :: ipiv(:)
     integer(ilp), intent(in) :: n
+    type(linalg_state_type), intent(out), optional :: err
     type(permutation_type) :: p
+    type(linalg_state_type) :: err0
     integer(ilp) :: i, tmp
 
     p%n = n
@@ -161,19 +111,22 @@ contains
     allocate(p%perm(n))
 
     ! Start with identity
-    do i = 1, n
-      p%perm(i) = i
-    end do
+    p%perm = arange(1_ilp, n)
 
     ! Replay swap sequence forward
     do i = 1, min(n, int(size(ipiv), ilp))
-      if (ipiv(i) < 1 .or. ipiv(i) > n) error stop "permutation_from_lapack: Invalid LAPACK pivot index out of bounds"
+      if (ipiv(i) < 1 .or. ipiv(i) > n) then
+        err0 = linalg_state_type("permutation_from_lapack", LINALG_VALUE_ERROR, "Invalid LAPACK pivot index out of bounds")
+        call linalg_error_handling(err0, err)
+        return
+      end if
       if (ipiv(i) /= i) then
         tmp = p%perm(i)
         p%perm(i) = p%perm(ipiv(i))
         p%perm(ipiv(i)) = tmp
       end if
     end do
+    call linalg_error_handling(err0, err)
   end function
 
   !> Convert nonfinalized (swap) to finalized (direct) permutation
@@ -183,26 +136,23 @@ contains
 
     if (p%finalized) then
       fp = p
-      return
+    else
+      ! Treat swap sequence as LAPACK pivots
+      fp = permutation_from_lapack(p%perm, p%n)
     end if
-    ! Treat swap sequence as LAPACK pivots
-    fp = permutation_from_lapack(p%perm, p%n)
   end function
 
   !> Compose two finalized permutations: result = p1 * p2
   !> Apply p2 first, then p1
-  function permutation_compose(p1, p2) result(p)
+  module function permutation_compose(p1, p2) result(p)
     type(permutation_type), intent(in) :: p1, p2
     type(permutation_type) :: p
-    integer(ilp) :: i
 
     p%n = p1%n
     p%finalized = .true.
     allocate(p%perm(p%n))
 
-    do i = 1, p%n
-      p%perm(i) = p2%perm(p1%perm(i))
-    end do
+    p%perm = p2%perm(p1%perm)
   end function
 
   !> ----- Templated operations for multiple real kinds -----
@@ -211,59 +161,84 @@ contains
 
   !> Apply permutation from the left: B = P * A (row permutation)
   !> Row i of B comes from row perm(i) of A
-  subroutine permutation_apply_left_${t1[0]}$${k1}$ (p, A, B)
+  module subroutine permutation_apply_left_${t1[0]}$${k1}$ (p, A, B, err)
     type(permutation_type), intent(in) :: p
     ${t1}$, intent(in) :: A(:,:)
     ${t1}$, intent(out) :: B(:,:)
+    type(linalg_state_type), intent(out), optional :: err
+    type(linalg_state_type) :: err0
     integer(ilp) :: i
 
-    if (.not. p%finalized) error stop "Permutation must be finalized before applying."
-    do i = 1, p%n
+    if (.not. p%finalized) then
+      err0 = linalg_state_type("permutation_apply_left", LINALG_VALUE_ERROR, "Permutation must be finalized before applying.")
+      call linalg_error_handling(err0, err)
+      return
+    end if
+    do concurrent (i = 1:p%n)
       B(i, :) = A(p%perm(i), :)
     end do
+    call linalg_error_handling(err0, err)
   end subroutine
 
   !> Apply permutation from the right: B = A * P (column permutation)
   !> Column perm(j) of B comes from column j of A
-  subroutine permutation_apply_right_${t1[0]}$${k1}$ (p, A, B)
+  module subroutine permutation_apply_right_${t1[0]}$${k1}$ (p, A, B, err)
     type(permutation_type), intent(in) :: p
     ${t1}$, intent(in) :: A(:,:)
     ${t1}$, intent(out) :: B(:,:)
+    type(linalg_state_type), intent(out), optional :: err
+    type(linalg_state_type) :: err0
     integer(ilp) :: j
 
-    if (.not. p%finalized) error stop "Permutation must be finalized before applying."
-    do j = 1, p%n
+    if (.not. p%finalized) then
+      err0 = linalg_state_type("permutation_apply_right", LINALG_VALUE_ERROR, "Permutation must be finalized before applying.")
+      call linalg_error_handling(err0, err)
+      return
+    end if
+    do concurrent (j = 1:p%n)
       B(:, p%perm(j)) = A(:, j)
     end do
+    call linalg_error_handling(err0, err)
   end subroutine
 
   !> Apply permutation to a vector: b = P * a
   !> b(i) = a(perm(i))
-  subroutine permutation_apply_vector_${t1[0]}$${k1}$ (p, a, b)
+  module subroutine permutation_apply_vector_${t1[0]}$${k1}$ (p, a, b, err)
     type(permutation_type), intent(in) :: p
     ${t1}$, intent(in) :: a(:)
     ${t1}$, intent(out) :: b(:)
-    integer(ilp) :: i
+    type(linalg_state_type), intent(out), optional :: err
+    type(linalg_state_type) :: err0
 
-    if (.not. p%finalized) error stop "Permutation must be finalized before applying."
-    do i = 1, p%n
-      b(i) = a(p%perm(i))
-    end do
+    if (.not. p%finalized) then
+      err0 = linalg_state_type("permutation_apply_vector", LINALG_VALUE_ERROR, "Permutation must be finalized before applying.")
+      call linalg_error_handling(err0, err)
+      return
+    end if
+    b = a(p%perm)
+    call linalg_error_handling(err0, err)
   end subroutine
 
   !> Convert permutation to dense matrix (for testing / small cases)
-  subroutine permutation_to_matrix_${t1[0]}$${k1}$ (p, A)
+  module subroutine permutation_to_matrix_${t1[0]}$${k1}$ (p, A, err)
     type(permutation_type), intent(in) :: p
     ${t1}$, intent(out) :: A(:,:)
+    type(linalg_state_type), intent(out), optional :: err
+    type(linalg_state_type) :: err0
     integer(ilp) :: i
 
-    if (.not. p%finalized) error stop "Permutation must be finalized to convert to matrix."
+    if (.not. p%finalized) then
+      err0 = linalg_state_type("permutation_to_matrix", LINALG_VALUE_ERROR, "Permutation must be finalized to convert to matrix.")
+      call linalg_error_handling(err0, err)
+      return
+    end if
     A = 0
     do i = 1, p%n
       A(i, p%perm(i)) = 1
     end do
+    call linalg_error_handling(err0, err)
   end subroutine
 
   #:endfor
 
-end module stdlib_linalg_permutation
+end submodule stdlib_linalg_permutation
diff --git a/test/linalg/test_linalg_permutation.fypp b/test/linalg/test_linalg_permutation.fypp
index eefad5271..8d2aaa447 100644
--- a/test/linalg/test_linalg_permutation.fypp
+++ b/test/linalg/test_linalg_permutation.fypp
@@ -5,7 +5,10 @@ module test_linalg_permutation
   use testdrive, only: new_unittest, unittest_type, error_type, check
   use stdlib_kinds, only: sp, dp, int32
   use stdlib_linalg_constants, only: ilp
-  use stdlib_linalg_permutation
+  use stdlib_linalg, only: permutation_type, permutation_identity, permutation_init, &
+    permutation_is_valid, permutation_invert, permutation_compose, &
+    permutation_from_lapack, permutation_to_matrix, &
+    permutation_apply_left, permutation_apply_right, permutation_apply_vector
   implicit none
   private
 

From 540c57c582e2971431d614c3f8e9054257e2e594 Mon Sep 17 00:00:00 2001
From: edish-github <sudhanshukit2021@gmail.com>
Date: Wed, 29 Apr 2026 03:51:20 +0530
Subject: [PATCH 09/10] Update src/linalg/stdlib_linalg_permutation.fypp

Co-authored-by: Jeremie Vandenplas <jeremie.vandenplas@gmail.com>
---
 src/linalg/stdlib_linalg_permutation.fypp | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/src/linalg/stdlib_linalg_permutation.fypp b/src/linalg/stdlib_linalg_permutation.fypp
index 890213cc6..bd3c40ae1 100644
--- a/src/linalg/stdlib_linalg_permutation.fypp
+++ b/src/linalg/stdlib_linalg_permutation.fypp
@@ -41,7 +41,7 @@ contains
     type(permutation_type) :: p
 
     p%n = size(perm_vec)
-    allocate(p%perm(p%n), source=perm_vec)
+    allocate(p%perm, source=perm_vec)
     p%finalized = optval(finalized, .true.)
   end function
 

From 50c9bebbd828a181f829a524fab90306aad48eea Mon Sep 17 00:00:00 2001
From: Sudhanshu Thakur <sudhanshukit2021@gmail.com>
Date: Wed, 29 Apr 2026 04:06:26 +0530
Subject: [PATCH 10/10] refactor(linalg): address PR review (pure, do
 concurrent, naming)

---
 src/linalg/stdlib_linalg.fypp             | 10 ++++----
 src/linalg/stdlib_linalg_permutation.fypp |  8 +++----
 test/linalg/test_linalg_permutation.fypp  | 28 +++++++++++------------
 3 files changed, 23 insertions(+), 23 deletions(-)

diff --git a/src/linalg/stdlib_linalg.fypp b/src/linalg/stdlib_linalg.fypp
index 30a9f98ef..0f42e5fc5 100644
--- a/src/linalg/stdlib_linalg.fypp
+++ b/src/linalg/stdlib_linalg.fypp
@@ -75,7 +75,7 @@ module stdlib_linalg
   public :: permutation_type
   public :: permutation_identity
   public :: permutation_init
-  public :: permutation_is_valid
+  public :: is_valid_permutation
   public :: permutation_invert
   public :: permutation_compose
   public :: permutation_from_lapack
@@ -2128,7 +2128,7 @@ module stdlib_linalg
     !! version: experimental
     !!
     !! Creates an identity permutation of size n
-    module function permutation_identity(n) result(p)
+    pure module function permutation_identity(n) result(p)
       integer(ilp), intent(in) :: n
       type(permutation_type) :: p
     end function
@@ -2138,18 +2138,18 @@ module stdlib_linalg
     !! version: experimental
     !!
     !! Creates a permutation from a given integer vector
-    module function permutation_init(perm_vec, finalized) result(p)
+    pure module function permutation_init(perm_vec, finalized) result(p)
       integer(ilp), intent(in) :: perm_vec(:)
       logical, intent(in), optional :: finalized
       type(permutation_type) :: p
     end function
   end interface
 
-  interface permutation_is_valid
+  interface is_valid_permutation
     !! version: experimental
     !!
     !! Checks if a permutation is valid (each integer 1..n appears exactly once)
-    module function permutation_is_valid(p) result(valid)
+    pure module function is_valid_permutation(p) result(valid)
       type(permutation_type), intent(in) :: p
       logical :: valid
     end function
diff --git a/src/linalg/stdlib_linalg_permutation.fypp b/src/linalg/stdlib_linalg_permutation.fypp
index bd3c40ae1..bb1e695c2 100644
--- a/src/linalg/stdlib_linalg_permutation.fypp
+++ b/src/linalg/stdlib_linalg_permutation.fypp
@@ -24,7 +24,7 @@ submodule (stdlib_linalg) stdlib_linalg_permutation
 contains
 
   !> Create an identity permutation of size n
-  module function permutation_identity(n) result(p)
+  pure module function permutation_identity(n) result(p)
     integer(ilp), intent(in) :: n
     type(permutation_type) :: p
 
@@ -35,7 +35,7 @@ contains
   end function
 
   !> Create a permutation from a given integer vector
-  module function permutation_init(perm_vec, finalized) result(p)
+  pure module function permutation_init(perm_vec, finalized) result(p)
     integer(ilp), intent(in) :: perm_vec(:)
     logical, intent(in), optional :: finalized
     type(permutation_type) :: p
@@ -46,7 +46,7 @@ contains
   end function
 
   !> Check if a permutation is valid (each integer 1..n appears exactly once)
-  module function permutation_is_valid(p) result(valid)
+  pure module function is_valid_permutation(p) result(valid)
     type(permutation_type), intent(in) :: p
     logical :: valid
     logical, allocatable :: seen(:)
@@ -91,7 +91,7 @@ contains
     allocate(inv%perm(fp%n))
 
     ! Direct inversion: swap keys and values
-    do i = 1, fp%n
+    do concurrent (i = 1:fp%n)
       inv%perm(fp%perm(i)) = i
     end do
   end function
diff --git a/test/linalg/test_linalg_permutation.fypp b/test/linalg/test_linalg_permutation.fypp
index 8d2aaa447..52fd441cf 100644
--- a/test/linalg/test_linalg_permutation.fypp
+++ b/test/linalg/test_linalg_permutation.fypp
@@ -6,7 +6,7 @@ module test_linalg_permutation
   use stdlib_kinds, only: sp, dp, int32
   use stdlib_linalg_constants, only: ilp
   use stdlib_linalg, only: permutation_type, permutation_identity, permutation_init, &
-    permutation_is_valid, permutation_invert, permutation_compose, &
+    is_valid_permutation, permutation_invert, permutation_compose, &
     permutation_from_lapack, permutation_to_matrix, &
     permutation_apply_left, permutation_apply_right, permutation_apply_vector
   implicit none
@@ -43,7 +43,7 @@ contains
     if (allocated(error)) return
     call check(error, all(p%perm == [1_ilp, 2_ilp, 3_ilp, 4_ilp, 5_ilp]), "identity elements")
     if (allocated(error)) return
-    call check(error, permutation_is_valid(p), "identity valid")
+    call check(error, is_valid_permutation(p), "identity valid")
   end subroutine
 
   !> Test 2: Valid permutation accepted
@@ -52,7 +52,7 @@ contains
     type(permutation_type) :: p
 
     p = permutation_init([3, 1, 4, 2])
-    call check(error, permutation_is_valid(p), "valid permutation")
+    call check(error, is_valid_permutation(p), "valid permutation")
   end subroutine
 
   !> Test 3: Invalid permutations rejected
@@ -62,12 +62,12 @@ contains
 
     ! Duplicate entry
     p = permutation_init([1, 1, 3])
-    call check(error, .not. permutation_is_valid(p), "reject duplicate")
+    call check(error, .not. is_valid_permutation(p), "reject duplicate")
     if (allocated(error)) return
 
     ! Out of range
     p = permutation_init([1, 5, 3])
-    call check(error, .not. permutation_is_valid(p), "reject out of range")
+    call check(error, .not. is_valid_permutation(p), "reject out of range")
   end subroutine
 
   !> Test 4: Inversion correctness
@@ -140,9 +140,9 @@ contains
   subroutine test_apply_vector(error)
     type(error_type), allocatable, intent(out) :: error
     type(permutation_type) :: p
-    real(dp) :: a(4), b(4)
+    real(dp), parameter :: a(4) = [10.0_dp, 20.0_dp, 30.0_dp, 40.0_dp]
+    real(dp) :: b(4)
 
-    a = [10.0_dp, 20.0_dp, 30.0_dp, 40.0_dp]
     p = permutation_init([3, 1, 4, 2])
     call permutation_apply_vector(p, a, b)
 
@@ -162,7 +162,9 @@ contains
 
     ! Verify: composed applied to vector = p1 applied to (p2 applied to vector)
     a = [10.0_dp, 20.0_dp, 30.0_dp, 40.0_dp]
-    b1 = 0.0_dp; b2 = 0.0_dp; b_composed = 0.0_dp ! Initialize to avoid garbage
+    b1 = 0.0_dp
+    b2 = 0.0_dp
+    b_composed = 0.0_dp
 
     call permutation_apply_vector(p2, a, b1)       ! b1 = P2 * a
     call permutation_apply_vector(p1, b1, b2)      ! b2 = P1 * P2 * a
@@ -210,10 +212,8 @@ contains
     integer(ilp) :: i, j
 
     p = permutation_identity(3)
-    do i = 1, 3
-      do j = 1, 3
-        A(i,j) = real(i*10 + j, dp)
-      end do
+    do concurrent (i = 1:3, j = 1:3)
+      A(i,j) = real(i*10 + j, dp)
     end do
 
     call permutation_apply_left(p, A, B)
@@ -226,13 +226,13 @@ contains
   subroutine test_types_support(error)
     type(error_type), allocatable, intent(out) :: error
     type(permutation_type) :: p
-    complex(dp) :: c_in(3), c_out(3)
+    complex(dp), parameter :: c_in(3) = [(1.0_dp, -1.0_dp), (2.0_dp, -2.0_dp), (3.0_dp, -3.0_dp)]
+    complex(dp) :: c_out(3)
     integer(int32) :: i_in(3), i_out(3)
 
     p = permutation_init([3_ilp, 1_ilp, 2_ilp])
 
     ! Complex support
-    c_in = [(1.0_dp, -1.0_dp), (2.0_dp, -2.0_dp), (3.0_dp, -3.0_dp)]
     call permutation_apply_vector(p, c_in, c_out)
     call check(error, abs(c_out(1) - c_in(3)) < 1.0e-10_dp, "complex")
     if (allocated(error)) return