Matlab Function in Fortran - perms - permutation

Below is the Fortran implementation of the perms algorithm found in Matlab. It creates a matrix whose rows consist of all possible permutations of the n elements of vector vecIn.

There are many useful utility functions in Matlab and occassionally an important but rare (difficult to find source in the internet) algorithm such as this permutation algorithm. Hopefully there will be more Fortran version of Matlab utility functions appearing here in the future.

The one listed below is unpolished and unoptimized and not suitable for large inputs. It is a very basic implementation taken directly from the mathematical definition. It is also tested only for a small number of cases. Corrections, suggestions and comments are welcomed.




!/*
! * perms - All Possible permutations
! * 
! *     (C)  Copyright 2010 xTechNotes.blogspot.com
! * 
! *   This program is free software: you can redistribute it and/or modify
! *   it under the terms of the GNU General Public License as published by
! *   the Free Software Foundation, either version 3 of the License, or
! *   (at your option) any later version.
! *
! *    This program is distributed in the hope that it will be useful,
! *    but WITHOUT ANY WARRANTY; without even the implied warranty of
! *    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
! *    GNU General Public License for more details.
! *
! *    You should have received a copy of the GNU General Public License
! *    along with this program.  If not, see .
! */ 


    RECURSIVE subroutine perms(n, nfact, vecIn, vecOut)
        IMPLICIT NONE 
        INTEGER, INTENT(IN) :: n                ! number of elements of original vector
        INTEGER, INTENT(IN) :: nfact            ! Factorial(n)
        REAL, INTENT(IN) :: vecIn(n)            ! original vector to be permuted
        REAL, INTENT(OUT):: vecOut(nfact, n)    ! Permutations of original vector


        ! local variables         
        real :: vecInTmp(n)
        integer :: ii , mfact
        real :: dtmp
        
        vecOut = 0.0
        if (n .le. 0) then 
            return
        elseif ( n .eq. 1 ) then            
            vecOut(1,1) = vecIn(1)
        elseif ( n .eq. 2 ) then
            vecOut(1, :) = vecIn(:)
            vecOut(2, :) = (/vecIn(2), vecIn(1)/)
        else        
            ! ii = 1
            call factorial(n-1, mfact)
            vecOut(1:mfact, 1) = vecIn(1)
            call perms(n-1, mfact,  vecIn(2:n), vecOut(1:mfact, 2:n))
            
            do ii = 2, n
                vecInTmp = VecIn
                vecInTmp(1) = VecIn(ii)
                vecInTmp(ii) = VecIn(1)
                vecOut((ii-1)*mfact+1 : ii*mfact, 1) = vecInTmp(1)                
                call perms(n-1, mfact, vecInTmp(2:n), vecOut( (ii-1)*mfact+1 : ii*mfact, 2:n))
            enddo 
        
        endif    
        
    end subroutine perms