3
\$\begingroup\$

I have written a quick sort routine to inplace sort an array in descending order:

module sort
  contains
    recursive subroutine qsort_self(array)
      real(kind=8)    :: array(:), temp
      integer :: length, pivot, i, j

      i = 0
      j = 0
      length = 0
      pivot = 1
      temp = 0
      length = size(array)
      if (length > 1) then
        i = 2
        j = size(array)
        do while (i <= j)
           do while (array(j) < array(pivot))
            j = j - 1
          enddo
          do while ((array(i) > array(pivot)) .and. (i<=length))
            i = i + 1
          enddo
         if (i > j) exit
          temp = array(j)
          array(j) = array(i)
          array(i) = temp
          i = i + 1
          j = j - 1
        enddo
          temp = array(pivot)
          array(pivot) = array(j)
          array(j) = temp
          pivot = j
          call qsort_self(array(1:pivot-1))
          call qsort_self(array(pivot + 1 : length))
      endif
    end subroutine qsort_self
  end module sort


program p
  use sort
  implicit none
  integer :: size_,i
  real(kind=8), allocatable    :: b6(:)

  size_ = 100
  allocate(b6(size_))
  do i = 1, size_
    b6(i) = real(rand(), 8)
  enddo 
  print *, "unsorted: ", b6
  call qsort_self(b6)
  print *,"sorted a= ", b6
end program p

It works fine for now, but without .and. (i<=length) part of code, i run into errors where if array is in ascending order, i goes out of bound. I checked several online implementations, such as this:

https://stackoverflow.com/questions/11196571/quick-sort-sorts-descending-not-ascending

and I believe they suffer from same flaw. Is here any better way of doing it? I Always want to take first element as pivot because of the array structure.

\$\endgroup\$
  • 1
    \$\begingroup\$ Welcome to CodeReview@SE. Could you mention the Fortan version/dialect? (from Fortran 90, inline comments should be possible) (punched cards, anyone?) \$\endgroup\$ – greybeard Dec 24 '19 at 15:03
  • \$\begingroup\$ @greybeard at least it doesn't have line numbers and gotos. \$\endgroup\$ – pacmaninbw Dec 24 '19 at 15:11
  • \$\begingroup\$ F90. Comments on <40 lines, one of the most famous algorithms? \$\endgroup\$ – ipcamit Dec 25 '19 at 3:51
  • \$\begingroup\$ Documenting comments with every single piece of code. \$\endgroup\$ – greybeard Dec 25 '19 at 9:34
  • \$\begingroup\$ I really wont argue further unless their is some actual advice included. But, I am sorry, this is somewhat stupid. Comments are not some holy relic that makes any code more readable. Uselessly commenting a 40 line code with three variables and 3 while loops can break the flow of reader too. So i disagree on commenting every single piece of code. Its useless and leads to situation like: i.redd.it/iuy9fxt300811.png. \$\endgroup\$ – ipcamit Dec 25 '19 at 11:37
3
\$\begingroup\$

Errors in the code

This piece

array(i) > array(pivot) .and. i <= length

is not valid code. If i > length you will get an out of bounds access with array(i). If you rewrite it as

i <= length .and. array(i) > array(pivot)

it is still wrong, because there is no short circuiting in Fortran. Or to be more precise: Fortran operators are neither short-circuit nor eager. The language specification allows the compiler to select the method for optimization. In your case the compiler reorders and short circuits apparently. The only way to write it reliably correctly is unfortunately:

if (i <= length) then
    if (array(i) > array(pivot)) then
        ...

A declaration with a number literal like real(kind=8) is not portable. You might think, that it is 8 bytes and/or double precision, but the actual integer value for different kinds is implementation dependent. It just happens to be the case, that gfortran and ifort use the kind value 8 to denote double precision, but on another compiler this might be wrong. You should always use

integer, parameter :: dp = selected_real_kind(15, 307)

if you want to have a certain range and precision or

use, intrinsic :: iso_fortran_env, only: real64

if you want to have a float with 8 bytes. (Of course these both will be the same on most architectures today.)


There is no implicit none in the module. Strictly speaking this is not an error and you declared all variables correctly, but code without implicit none should be considered faulty.

Major Code improvements

You should always add one of intent(in), intent(out), intent(inout) to the dummy arguments of your procedures.


Fortran arrays can be declared with arbitrary numerical indices. You implicitly assumed in your code, that they start at 1 (which is the default). If you want your routine to work with differently indexed arrays, it is necessary to query the upper and lower bounds.

do i = lbound(array, dim=1), ubound(array, dim=1)

This makes it also explicit when you really mean the size, and when the bounds.

What I initially wrote is wrong. If you use assumed-shape arrays, the dummy argument becomes again 1-indexed by default. Had to learn that myself the hard way. Sorry. :-/


I agree you with you that these <40 lines don't require comments for the most parts, but you can write it more readable without comments, by defining a swap subroutine. This would also reduce the clutter of local variables.

If you define a get_pivot function, it will be easy to change to different pivoting strategies.


You have superfluous assignments to length in your code. IMHO you don't need this variable anyway, because you can insert size(array) at every appearance. Any decent compiler will do common-subexpression-elimination on these calls.


Where possible I would use pure and elemental on your procedures. They are amazing. ;-)

Architecture

I would always make all names in a module private and export explicitly with public. On the other hand I would only import used names.


In the long run, you want to have a second argument where one can pass a pure logical function that compares two elements. Doing so, the user can sort according to his/her constraints.

Formatting

Thank you for writing Fortran and not FORTRAN code. If you don't scream at people it is much easier to read your code. ;-)


Of course this is opinion based so, whatever I say, stick to the style of your environment. But for example the whitespace around your binary operators was not consistent. An indentation of two spaces is IMHO too short.


I assume that you are from science and you and your peers use python. In this case it is quite easy to stick to the python formatting guidelines in Fortran.

Minor Code improvements

I agree that the code mostly does not require comments, but at least public functions of a module should get a docstring on how to use them.


size in the program can be made constant.


The upper and lower bound in array slices can be ommited A(lbound(A, 1) : 5) == A(: 5)


I would not use p as program name. If you do so, you can't use it as variable name anymore.


There is random_number to generate an array of real numbers.

Possible improved version

module sort
    implicit none
    private
    public :: dp, qsort
    integer, parameter :: dp = selected_real_kind(15, 307)

    abstract interface
        logical pure function compare_t(a, b)
            import :: dp
            real(dp), intent(in) :: a, b
        end function
    end interface

contains

    !> @brief
    !>      Sort array inplace according to comparison function comp.
    !>
    !> @details
    !> The comparison function has to be a pure, logical function
    !> that compares two elements of array in a strict partial order.
    !> That means equality has to evaluate to false.
    pure recursive subroutine qsort(array, comp)
        real(kind=dp), intent(inout) :: array(:)
        procedure(compare_t) :: comp

        integer :: pivot, i, j

        i = 0
        j = 0
        pivot = get_pivot(array)
        if (size(array) > 1) then
            i = lbound(array, dim=1)
            j = ubound(array, dim=1)
            do while (i <= j)
                do while (comp(array(i), array(pivot)))
                    i = i + 1
                end do
                do while (comp(array(pivot), array(j)))
                    j = j - 1
                end do
                if (i >= j) exit
                call swap(array(i), array(j))
                i = i + 1
                j = j - 1
            end do
            call swap(array(pivot), array(j))
            call qsort(array(: j - 1), comp)
            call qsort(array(j + 1 : ), comp)
        end if
    end subroutine qsort

    pure subroutine swap(a, b)
        real(dp), intent(inout) :: a, b
        real(dp) :: temp
        temp = a
        a = b
        b = temp
    end subroutine

    pure function get_pivot(array) result(res)
        real(dp), intent(in) :: array(:)
        integer :: res
        res = lbound(array, dim=1)
    end function
  end module sort


program test_qsort_procedures
    use sort, only: dp, qsort
    implicit none
    integer, parameter :: size_ = 100
    real(kind=dp), allocatable :: b6(:)

    allocate(b6(size_))
    call random_number(b6)

    print *, "unsorted: ", b6
    call qsort(b6, increasing)
    print *,"increasing a= ", b6
    call qsort(b6, decreasing)
    print *,"decreasing a= ", b6

    contains

    logical pure function increasing(a, b)
        real(dp), intent(in) :: a, b
        increasing = a < b
    end function

    logical pure function decreasing(a, b)
        real(dp), intent(in) :: a, b
        decreasing = a > b
    end function
end program test_qsort_procedures
\$\endgroup\$
  • \$\begingroup\$ Thank you for such comprehensive answer. I was not aware of lbounds and definitely let the short circuiting error slip. Regarding 2-space indent, program name and docstring, this function is supposed to be a part of bigger submodule, hence the short program name; and already existing indent was 2 :'(. However thank you again for such comprehensive review, I would definitely include the changes \$\endgroup\$ – ipcamit Jan 15 at 11:00
  • 1
    \$\begingroup\$ Welcome. ;-) I just wrote down the stuff, that I would have liked to hear in the beginning. Especially the short circuiting can become really nasty, as you can imagine. \$\endgroup\$ – mcocdawc Jan 15 at 12:18

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.