feat (linalg): permutation matrices by edish-github · Pull Request #1155 · fortran-lang/stdlib

edish-github · 2026-03-25T14:11:33Z

Summary

Fixes #891.

Currently, applying permutation matrices requires storing them as dense matrices, wasting $O(n^2)$ memory and costing $O(n^2 k)$ time during multiplication.

This PR introduces a permutation_type type that stores the permutation as an integer array, bringing storage down to $O(n)$ and row/column operations down to $O(n)$.

Changes

Added permutation_type supporting both direct mapping (finalized) and LAPACK swap sequences (ipiv from dgetrf).
Added standard operations: identity, invert, compose, and array conversions.
Added permutation_apply_left (row swaps) and permutation_apply_right (column swaps).
Added 12 unit tests using testdrive covering all real, complex, and integer kinds..

Benchmark

I ran a quick sanity check to compare dense matmul vs the new vector approach on a 5000x100 matrix using gfortran -O3:

n=5000, k=100
Dense approach (P_dense @ A):  0.155507 s
Vector approach (row copy):    0.000703 s
Speedup: ~ 221x

…ge (fortran-lang#891)

jalvesz · 2026-03-25T16:31:47Z

src/linalg/stdlib_linalg_permutation.fypp

+!> Reference: Issue #891 (https://github.com/fortran-lang/stdlib/issues/891)
+
+#:include "common.fypp"
+#:set REAL_KINDS_TYPES = [("sp","real(sp)"), ("dp","real(dp)")]


Support for all real and complex kinds should be added for such functionality. This can be done quite easily with stdlib's fypp templates. Please look at other modules in the library, specially those in the linalg folder, to see how that can be done. Do not filter kinds unless there is a very good reason for it, if the motivation was LAPACK, stdlib's modernized BLAS/LAPACK covers all kinds thanks to fypp.

jalvesz · 2026-03-25T16:32:26Z

src/linalg/stdlib_linalg_permutation.fypp

+
+#:include "common.fypp"
+#:set REAL_KINDS_TYPES = [("sp","real(sp)"), ("dp","real(dp)")]
+#:set INT_KINDS_TYPES = [("int32","integer(int32)"), ("int64","integer(int64)")]


What's the purpose of defining integer kinds here in this manner ? if the reason is to define the integer size for the matrix maximun storage size, you should use ilp from stdlib_linalg_constants and procedure to use integer(ilp) instead of integer.
Now, if the reason is to support integer matrices, then you can use something like #:set RCI_KINDS_TYPES = REAL_KINDS_TYPES + CMPLX_KINDS_TYPES + INT_KINDS_TYPES for the template.

jalvesz · 2026-03-25T16:36:05Z

test/linalg/test_linalg_permutation.fypp

+    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")


errors checks should be done on full vectors or full matrices, not individual entries. you can use either all or a check on the L2 norm of the difference. This is much more robust.

jalvesz · 2026-03-25T16:56:17Z

src/linalg/stdlib_linalg_permutation.fypp

+  private
+
+  !> Permutation type: stores permutation as integer vector
+  type, public :: permutation_t


stdlib's convention is to use <>_type for public derived types.

jalvesz · 2026-03-25T16:57:31Z

src/linalg/stdlib_linalg_permutation.fypp

+
+  ! Generic interfaces for multi-type support
+  interface permutation_apply_left
+#:for rk, rt in REAL_KINDS_TYPES


please indent at the internal interface scope level

edish-github · 2026-03-25T17:49:20Z

Thanks for the detailed review @jalvesz, I've pushed the commit that addresses the review:

Swapped to permutation_type
Replaced the hardcoded kinds with RCI_KINDS_TYPES
Updated matrix sizing to use ilp
Switched the test assertions to array-level all(abs(...))
Fixed the interface formatting

Copilot

Pull request overview

Adds a dedicated linear-algebra permutation facility that represents permutation matrices as an index vector (instead of dense matrices) and wires it into the build + test suite.

Changes:

Introduced stdlib_linalg_permutation with permutation_type and basic permutation operations (identity/init/validate/invert/compose, LAPACK pivot conversion, apply left/right/vector, and dense-matrix conversion).
Added a new testdrive-based test program for permutation behavior and registered it in the linalg test CMake.
Added the new linalg source file to the linalg build CMake list.

Reviewed changes

Copilot reviewed 4 out of 4 changed files in this pull request and generated 8 comments.

File	Description
`src/linalg/stdlib_linalg_permutation.fypp`	New module implementing vector-backed permutation operations and apply routines.
`src/linalg/CMakeLists.txt`	Adds the new permutation module to the linalg build inputs.
`test/linalg/test_linalg_permutation.fypp`	New unit tests for the permutation module.
`test/linalg/CMakeLists.txt`	Registers the new permutation test for CMake/CTest execution.

💡 Add Copilot custom instructions for smarter, more guided reviews. Learn how to get started.

Copilot · 2026-03-25T18:11:34Z

src/linalg/stdlib_linalg_permutation.fypp

+!>   - 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


The header comment defines finalized permutations as “perm(i)=j means row i maps to row j”, but the implementation (apply_left/apply_vector and permutation_to_matrix) treats perm(i) as the source index for output position i (i.e., B(i,:)=A(perm(i),:) and b(i)=a(perm(i))). Please clarify/adjust the documentation so the meaning of perm(i) matches the actual semantics used throughout this module.

Suggested change

!> - 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

Copilot · 2026-03-25T18:11:35Z

src/linalg/stdlib_linalg_permutation.fypp

+  !> 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
+
+    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


permutation_is_valid currently always enforces “each integer 1..n appears exactly once”, which is only appropriate for finalized permutations. For a non-finalized/LAPACK-pivot representation (ipiv swap sequence), duplicates like [2,3,3] are valid but will be rejected. Consider branching on p%finalized (range-check only for ipiv), or rename/split the validator so callers can correctly validate both representations.

Copilot · 2026-03-25T18:11:35Z

src/linalg/stdlib_linalg_permutation.fypp

+    allocate(p%perm(p%n))
+
+    do i = 1, p%n
+      p%perm(i) = p1%perm(p2%perm(i))


permutation_compose appears to compose in the wrong order for the semantics used by permutation_apply_vector (b(i)=a(perm(i))). If p = p1 * p2 is meant to apply p2 then p1, the resulting permutation should satisfy p%perm(i)=p2%perm(p1%perm(i)), not p1%perm(p2%perm(i)). As written, composition will return incorrect results for general permutations.

Suggested change

p%perm(i) = p1%perm(p2%perm(i))

p%perm(i) = p2%perm(p1%perm(i))

Copilot · 2026-03-25T18:11:35Z

src/linalg/stdlib_linalg_permutation.fypp

+  !> 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
+
+    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
+
+    do j = 1, p%n
+      B(:, p%perm(j)) = A(:, j)
+    end do


The apply/to-matrix routines don’t account for p%finalized = .false. (swap-sequence representation). As a result, a permutation created with permutation_init(..., finalized=.false.) will be interpreted as a direct mapping and produce incorrect indexing (and potentially out-of-bounds). Either finalize internally (or expose a public finalize operation), or explicitly reject non-finalized permutations in these routines with a clear error.

Copilot · 2026-03-25T18:11:35Z

src/linalg/stdlib_linalg_permutation.fypp

+    ! 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


permutation_from_lapack silently ignores out-of-range pivot indices (it just skips the swap). This can hide corrupted ipiv inputs and yield a permutation that looks “valid” but is wrong. Consider validating that all used ipiv(i) are within 1..n (and raising an error/returning an invalid permutation) instead of skipping invalid entries.

Copilot · 2026-03-25T18:11:36Z

test/linalg/test_linalg_permutation.fypp

+    ! perm = [3,1,2] means: pos1->3, pos2->1, pos3->2
+    ! inverse: pos3->1, pos1->2, pos2->3 => inv = [2,3,1]


The comments describing the meaning of perm = [3,1,2] (“pos1->3, …”) don’t match the semantics exercised by permutation_apply_vector/permutation_apply_left (which use b(i)=a(perm(i)) and B(i,:)=A(perm(i),:)). Please update these comments to avoid implying the opposite mapping direction.

Suggested change

! 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: src1 -> dest2, src2 -> dest3, src3 -> dest1 => inv = [2,3,1]

Copilot · 2026-03-25T18:11:36Z

test/linalg/test_linalg_permutation.fypp

+    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)")


test_compose uses p2 that is the inverse of p1, so the test passes even if the composition order in permutation_compose is reversed. Add at least one composition test with two non-inverse permutations (and verify against explicit application) to catch order/definition errors.

Copilot · 2026-03-25T18:11:36Z

test/linalg/test_linalg_permutation.fypp

+  !> 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)
+
+    ! 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_type) :: 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_type) :: 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


The tests here exercise only real(dp) arrays. If the implementation is intended to support complex and integer kinds (as described in the PR), consider adding parameterized/fypp-generated cases (or additional tests) that call apply/to_matrix for representative complex and integer kinds as well.

…back

…er support

codecov · 2026-03-26T07:41:53Z

Codecov Report

✅ All modified and coverable lines are covered by tests.
✅ Project coverage is 68.01%. Comparing base (ac98d3b) to head (4780ae2).
⚠️ Report is 5 commits behind head on master.

Additional details and impacted files

@@           Coverage Diff           @@
##           master    #1155   +/-   ##
=======================================
  Coverage   68.00%   68.01%           
=======================================
  Files         404      404           
  Lines       12935    12948   +13     
  Branches     1392     1395    +3     
=======================================
+ Hits         8797     8806    +9     
- Misses       4138     4142    +4

☔ View full report in Codecov by Sentry.
📢 Have feedback on the report? Share it here.

🚀 New features to boost your workflow:

❄️ Test Analytics: Detect flaky tests, report on failures, and find test suite problems.

jalvesz · 2026-03-29T15:56:48Z

src/linalg/stdlib_linalg_permutation.fypp

+
+  ! Generic interfaces for multi-type support
+  interface permutation_apply_left
+#:for k1, t1 in RCI_KINDS_TYPES


Please indent the fypp preprocessing macros to the entity scope.

interface xyz #:for ... module procedure xyz_<> ... #:endfor

Please check across the PR, not only this line.

jalvesz · 2026-03-29T16:04:32Z

test/linalg/test_linalg_permutation.fypp

+    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)")


When possible, please check on whole arrays, not individual terms

call check(error, abs( all( b - [val1,val2,...]) < tol ), "some error comment" )

you can also store the reference values in a reference array and use it for the error check

loiseaujc

Here is my first pass.

I do have one major question though for @jalvesz, @jvdp1 and @perazz : at the moment, the permutation_apply_* are performing out-of-place computations. In many situations however, I guess permutations should be applied in-place for performance-sake. How about having a subroutine permutation_apply_*(A, p) performing in-place, and a function permutation_apply_*(A, p) result(B) for out-of-place (I know subs and functions can't have same name, just lazy here) ?

loiseaujc · 2026-03-30T06:00:15Z

src/linalg/stdlib_linalg_permutation.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


I would suggest making this contribution a submodule of linalg rather than a standalone module. It is just standard linear algebra after all. What's your take @jalvesz ?

Yes, agree!

loiseaujc · 2026-03-30T06:04:40Z

src/linalg/stdlib_linalg_permutation.fypp

+    allocate(p%perm(n))
+    do i = 1, n
+      p%perm(i) = i
+    end do


MInor detail, but you could use the arange function from stdlib_math to collapse these four lines into a single one.

loiseaujc · 2026-03-30T06:06:01Z

src/linalg/stdlib_linalg_permutation.fypp

+      p%finalized = finalized
+    else
+      p%finalized = .true.
+    end if


You can use optval from the stdlib_optval module.

p%finalize = optval(finalize, .true.)

src/linalg/stdlib_linalg_permutation.fypp

loiseaujc · 2026-03-30T06:09:32Z

src/linalg/stdlib_linalg_permutation.fypp

+
+    ! 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"


Please use linalg_state_type and linalg_handling_error for error handling.

src/linalg/stdlib_linalg_permutation.fypp

loiseaujc · 2026-03-30T06:14:33Z

src/linalg/stdlib_linalg_permutation.fypp

+    ${t1}$, intent(out) :: B(:,:)
+    integer(ilp) :: i
+
+    if (.not. p%finalized) error stop "Permutation must be finalized before applying."


Use linalg_state_type and linalg_error_handling.

src/linalg/stdlib_linalg_permutation.fypp

edish-github added 2 commits March 25, 2026 17:50

feat(linalg): add permutation matrix type with efficient vector stora…

9f6c2c1

…ge (fortran-lang#891)

refactor: remove recursive keyword from permutation_invert

e1fb158

jalvesz changed the title ~~Feat/linalg permutation matrices~~ feat (linalg): permutation matrices Mar 25, 2026

jalvesz requested changes Mar 25, 2026

View reviewed changes

refactor(linalg): adressing the PR review

6406811

jalvesz requested a review from Copilot March 25, 2026 18:05

Copilot started reviewing on behalf of jalvesz March 25, 2026 18:06 View session

Copilot AI reviewed Mar 25, 2026

View reviewed changes

edish-github added 2 commits March 25, 2026 23:51

fix(linalg): address Copilot AI semantic and mathematical review feed…

5dd11ea

…back

test(linalg): add non-inverse composition test and test complex/integ…

02d9f77

…er support

test(linalg): fix uninitialized array in compose test

894733c

jalvesz requested changes Mar 29, 2026

View reviewed changes

jalvesz requested review from jvdp1, loiseaujc and perazz March 29, 2026 16:05

fix fypp macro indents and use array checks in tests

83334b0

loiseaujc requested changes Mar 30, 2026

View reviewed changes

convert permutation to submodule of stdlib_linalg

4780ae2

		!> - 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 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

	p%perm(i) = p1%perm(p2%perm(i))
	p%perm(i) = p2%perm(p1%perm(i))

		! perm = [3,1,2] means: pos1->3, pos2->1, pos3->2
		! inverse: pos3->1, pos1->2, pos2->3 => inv = [2,3,1]

Conversation

edish-github commented Mar 25, 2026 • edited Loading Uh oh! There was an error while loading. Please reload this page.

Uh oh!

Summary

Changes

Benchmark

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

jalvesz Mar 25, 2026 • edited Loading Uh oh! There was an error while loading. Please reload this page.

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

edish-github commented Mar 25, 2026

Uh oh!

Copilot AI left a comment

Choose a reason for hiding this comment

Pull request overview

Reviewed changes

Uh oh!

Copilot AI Mar 25, 2026

Choose a reason for hiding this comment

Uh oh!

Copilot AI Mar 25, 2026

Choose a reason for hiding this comment

Uh oh!

Copilot AI Mar 25, 2026

Choose a reason for hiding this comment

Uh oh!

Copilot AI Mar 25, 2026

Choose a reason for hiding this comment

Uh oh!

Copilot AI Mar 25, 2026

Choose a reason for hiding this comment

Uh oh!

Copilot AI Mar 25, 2026

Choose a reason for hiding this comment

Uh oh!

Copilot AI Mar 25, 2026

Choose a reason for hiding this comment

Uh oh!

Copilot AI Mar 25, 2026

Choose a reason for hiding this comment

Uh oh!

codecov bot commented Mar 26, 2026 • edited Loading Uh oh! There was an error while loading. Please reload this page.

Uh oh!

Codecov Report

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

loiseaujc left a comment

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Uh oh!

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Uh oh!

Uh oh!

Reviewers

Assignees

edish-github commented Mar 25, 2026 •

edited

Loading

jalvesz Mar 25, 2026 •

edited

Loading

codecov bot commented Mar 26, 2026 •

edited

Loading