! Metodos Numericos y Simulaciones en Astrofisica
! FCEFN - UNSJ, 2010
!
! ejemplo 3b
!
! Ordenamiento de arrays - Metodo Quicksort
!
!
program eje3c
implicit none
  integer, parameter ::d=selected_real_kind(16,200)
  integer :: ii
  real(d), dimension(50) :: ori,dat
  real(d) :: val

! inicializa generador de numeros aleatorios
!
  call random_seed()

! arma un array de numeros reales al azar
!
  do ii=1,50
    call random_number(val)
    dat(ii)=float(int(10000._d*val))/10._d
  enddo

! guarda el array original
!
  ori=dat

! lo ordena usando insercion
!
  call quick(50,dat)

! salida a pantalla
!
  write(*,*)'Array original: '
  write(*,*)ori
  write(*,*)' '
  write(*,*)'Array ordenado: '
  write(*,*)dat

end program eje3c

! subrutina para efectuar un ordenamiento
! utilizando el metodo Quicksort
!
! se ingresa con:
!
! N   numero de elementos del array
! ARR array a ordenar
!
! en la salida se devuelve el array ordenado en ARR.
!
! El parametro MM indica el tamaño maximo del subarray que
! se ordena por insercion y NSTACK es el tamaño del array
! auxiliar.
!
subroutine quick(n,arr)
implicit none
  integer, parameter ::d=selected_real_kind(16,200)
  integer, parameter :: mm=7,nstack=50
  integer :: n,ii,jj,kk,ll,ir,jstack,flg
  integer,dimension(nstack) :: istack
  real(d),dimension(n) :: arr  
  real(d) :: aa,temp

  jstack=0
  ll=1
  ir=n

  do

! los subarrays son pequeños y usa insercion
!
    if((ir-ll) < mm)then
      jj=ll+1
      do while(jj <= ir)
        aa=arr(jj)
        ii=jj-1
        flg=0
        do while(ii >= ll)
          if(arr(ii) <= aa)flg=1
          if(flg == 1)exit
          arr(ii+1)=arr(ii)
          ii=ii-1
        enddo 
        if(flg == 0)ii=ll-1
        arr(ii+1)=aa
        jj=jj+1
      enddo

! punto de salida
!
      if(jstack == 0)return

! no sale e inicia una nuevo ordenamiento
!
      ir=istack(jstack)
      ll=istack(jstack-1)
      jstack=jstack-2
    else

! elije el pivot calculando la media del primero, el 
! elemento del centro y el ultimo. Tambien los ordena
!
      kk=(ll+ir)/2
      temp=arr(kk)
      arr(kk)=arr(ll+1)
      arr(ll+1)=temp
      if(arr(ll) > arr(ir))then
        temp=arr(ll)
        arr(ll)=arr(ir)
        arr(ir)=temp
      endif
      if(arr(ll+1) > arr(ir))then
        temp=arr(ll+1)
        arr(ll+1)=arr(ir)
        arr(ir)=temp
      endif
      if(arr(ll) > arr(ll+1))then
        temp=arr(ll)
        arr(ll)=arr(ll+1)
        arr(ll+1)=temp
      endif

! inicializa punteros para cada particion
!
      ii=ll+1
      jj=ir

! guarda pivot
!
      aa=arr(ll+1)

! inicio del bucle interno
!
      do

! busca un elemento mayor que el pivot
!
        do while(arr(ii) <= aa)
          ii=ii+1
        enddo

! busca un elemento menor que el pivot
!
        do while(arr(jj) >= aa)
          jj=jj-1
        enddo

! si los punteros no estan cruzados los elementos
! encontrados estan en la zona correcta
!
        if(jj < ii)exit

! no esta bien asi que los intercambia
!
        temp=arr(ii)
        arr(ii)=arr(jj)
        arr(jj)=temp
      enddo

! inserta el pivot
!
      arr(ll+1)=arr(jj)
      arr(jj)=aa

      jstack=jstack+2
      if(jstack > nstack)then
        write(*,*)'Array auxiliar muy chico'
        pause
      endif

      if(ir-ii+1 >= jj-1)then
        istack(jstack)=ir
        istack(jstack-1)=ii
        ir=jj-1
      else
        istack(jstack)=jj-1
        istack(jstack-1)=ll
        ll=ii
      endif
    endif
  enddo

end subroutine
