!***************************************************************************
!*               Test program for Bauhuber's method                        *
!* ----------------------------------------------------------------------- *
!*  This program uses Bauhuber's Method to find all real or complex        *
!*  roots of a polynomial of degree n:                                     *
!*                                               n-1           n           *
!*      P(x) = a[0] + a[1] * x + ... + a[n-1] * x    + a[n] * x ,          *
!*                                                                         *
!*  where a[i], i=0, ..., n, can be complex.                               *
!* ----------------------------------------------------------------------- *
!* SAMPLE RUN:                                                             *
!* (Find all roots of polynomial of degree 4:                              *
!*  P(x) = 0.000089248 - 0.04368 X + 2.9948 X2 - 6.798 X3 + X4)            *
!*                                                                         *
!*   ----------------------------------------------------------            *
!*   Complex and Real Roots of a Polynomial (Bauhuber's method)            *
!*   ----------------------------------------------------------            *
!*   Polynomial coefficients:                                              *
!*   a( 0) =  0.892480E-04   0.000000E+00                                  *
!*   a( 1) = -0.436800E-01   0.000000E+00                                  *
!*   a( 2) =  0.299480E+01   0.000000E+00                                  *
!*   a( 3) = -0.679800E+01   0.000000E+00                                  *
!*   a( 4) =  0.100000E+01   0.000000E+00                                  *
!*   Roots (without scaling)                                               *
!*     No   Real part    Imaginary part   Function value                   *
!*      0   0.245372E-02   0.000000E+00   0.205566E-14                     *
!*      1   0.125733E-01   0.000000E+00   0.205559E-14                     *
!*      2   0.457319E+00   0.000000E+00   0.206118E-14                     *
!*      3   0.632565E+01   0.000000E+00   0.514305E-13                     *
!*   Roots (with scaling)                                                  *
!*     No   Real part    Imaginary part   Function value                   *
!*      0   0.245372E-02   0.000000E+00   0.205566E-14                     *
!*      1   0.125733E-01   0.000000E+00   0.205565E-14                     *
!*      2   0.457319E+00   0.000000E+00   0.224209E-14                     *
!*      3   0.632565E+01   0.000000E+00   0.514305E-13                     *
!*   ----------------------------------------------------------            *
!*                                                                         *
!* ----------------------------------------------------------------------- *
!* Ref.: "Numerical algorithms with C, By Gisela Engeln-Muellges and       *
!*        Frank Uhlig, Springer-Verlag, 1996" [BIBLI 11].                  *
!*                                                                         *
!*                                       F90 Release By J-P Moreau, Paris. *
!*                                              (www.jpmoreau.fr)          *
!***************************************************************************
Program TBauhube

Use FBauhube

IMPLICIT REAL*8 A-H,O-Z

Dimension ar(0:NMAX), ai(0:NMAX), rootr(0:NMAX), rooti(0:NMAX), val(0:NMAX)
Integer rc, skala 


  print *,' '
  print *,'  ----------------------------------------------------------'
  print *,'  Complex and Real Roots of a Polynomial (Bauhuber''s method)'
  print *,'  ----------------------------------------------------------'

  n = 4      !order of polynomial

  !define ar vector
  ar(0) =  0.000089248d0
  ar(1) = -0.04368d0
  ar(2) =  2.9948d0
  ar(3) = -6.798d0
  ar(4) =  1.d0

  !ai vector is null
  do i=0, n 
    ai(i) = 0.d0
  end do

  print *,'  Polynomial coefficients:'
  do i = 0, n
    write(*,10)  i, ar(i), ai(i)
  end do

  do j = 0, 1
    skala = j
    rc = bauhub (0, skala, n, ar, ai, rootr, rooti, val)
    if (rc == 0) then
      if (skala == 0) then
        print *,'  Roots (without scaling)'
      else
        print *,'  Roots (with scaling)'
      end if

      print *,'    No   Real part    Imaginary part   Function value'

      do i = 0, n-1
        write(*,20) i, rootr(i), rooti(i), val(i)
      end do
    else
      print *,'  *** Error in bauhube, rc=', rc
    end if
  end do

  print *,'  ----------------------------------------------------------'
  print *,' '
  stop

10 Format('   a(',I2,') = ', E13.6,'  ',E13.6)
20 Format('   ',I4,'  ',E13.6,'  ',E13.6,'  ',E13.6)

END

!end of file tbauhube.f90