feat (linalg): permutation matrices

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

src/linalg/stdlib_linalg_permutation.fypp

@@ -0,0 +1,269 @@
+!> 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 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)
+
+#:include "common.fypp"
+#: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
+  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)
+    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
+  end function
+
+  !> Create a permutation from a given integer vector
+  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
+  end function
+
+  !> Check if a permutation is valid (each integer 1..n appears exactly once)
+  function permutation_is_valid(p) result(valid)
+    type(permutation_type), intent(in) :: p
+    logical :: valid
+    logical, allocatable :: seen(:)
+    integer(ilp) :: i
+
+    valid = .false.
+    if (.not. allocated(p%perm)) return
+    if (p%n < 1) return
+    if (size(p%perm) /= p%n) return
+
+    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
+
+  !> Invert a finalized permutation
+  !> If perm(i) = j, then inv(j) = i
+  function permutation_invert(p) result(inv)
+    type(permutation_type), intent(in) :: p
+    type(permutation_type) :: inv, fp
+    integer(ilp) :: i
+
+    if (p%finalized) then
+      fp = p
+    else
+      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
+  !> ipiv(i) = j means "swap rows i and j" applied at step i
+  function permutation_from_lapack(ipiv, n) result(p)
+    integer(ilp), intent(in) :: ipiv(:)
+    integer(ilp), intent(in) :: n
+    type(permutation_type) :: p
+    integer(ilp) :: 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, 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
+      end if
+    end do
+  end function
+
+  !> Convert nonfinalized (swap) to finalized (direct) permutation
+  function permutation_finalize(p) result(fp)
+    type(permutation_type), intent(in) :: p
+    type(permutation_type) :: 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_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
+  end function
+
+  !> ----- Templated operations for multiple real kinds -----
+
+  #: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_${t1[0]}$${k1}$ (p, A, B)
+    type(permutation_type), intent(in) :: p
+    ${t1}$, intent(in) :: A(:,:)
+    ${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
+  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)
+    type(permutation_type), intent(in) :: p
+    ${t1}$, intent(in) :: A(:,:)
+    ${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
+  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)
+    type(permutation_type), intent(in) :: p
+    ${t1}$, intent(in) :: a(:)
+    ${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
+  end subroutine
+
+  !> Convert permutation to dense matrix (for testing / small cases)
+  subroutine permutation_to_matrix_${t1[0]}$${k1}$ (p, A)
+    type(permutation_type), intent(in) :: p
+    ${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
+    end do
+  end subroutine
+
+  #:endfor
+
+end module stdlib_linalg_permutation

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)