subroutine div(ah,al,bh,bl,ch,cl)
        implicit none
        double precision, intent(in) :: ah, al, bh, bl
        double precision, intent(out) :: ch, cl
        double precision s1,s2,q1,q2,r1,r2,p1,p2,t1,t2,u1
        q1 = ah/bh
        ch = bh * q1
        u1 = 134217729.0*bh
        p1 = u1 - (u1-bh)
        p2 = bh - p1
        u1 = 134217729.0*q1
        r1 = u1 - (u1-q1)
        r2 = q1 - r1
        cl = ((p1*r1-ch) + p1*r2 + p2*r1) + p2*r2
        cl = cl + (bl * q1)
        s1 = ch + cl
        cl = cl - (s1 - ch)
        ch = s1
        s1 = ah - ch
        r1 = s1 - ah
        s2 = (ah - (s1 - r1)) - (ch + r1)
        s2 = s2 - cl
        s2 = s2 + al
        q2 = (s1 + s2)/bh
        ch = q1 + q2
        cl = q2 - (ch - q1)
      end subroutine div

      program hilbert5_sample_dd
      use ddcommunication_s, only : eigen_init, eigen_free
      implicit double precision (a-h,o-v,x-z)

      double precision, allocatable :: ah(:,:), al(:,:)
      double precision, pointer :: bh(:), bl(:)
      double precision, pointer :: zh(:), zl(:)
      double precision, pointer :: wh(:), wl(:)
      double precision temph, templ
      logical iexist
!
      include 'mpif.h'
      include 'trd.h'

      real(16)::whl16
!
      call mpi_init(ierr)
      call mpi_comm_rank(mpi_comm_world,i$inod,ierr)
      call mpi_comm_size(mpi_comm_world,i$nnod,ierr)

      n=5

      allocate(ah(n,n),al(n,n))

      call eigen_init(2)
      NPROW = size_of_col
      NPCOL = size_of_row
      nx = ((n-1)/NPROW+1)
      call CSTAB_get_optdim(nx, 2, 1, 2, nm)
!      call CSTAB_get_optdim(nx, 6, 16*4, 16*4*2, nm)
      call eigen_free(0)
      NB  = 32
!      NB  = 64+32
      nmz = ((n-1)/NPROW+1)
      nmz = ((nmz-1)/NB+1)*NB+1
      nmw = ((n-1)/NPCOL+1)
      nmw = ((nmw-1)/NB+1)*NB+1

      larray = MAX(nmz,nm)*nmw

      allocate( bh(larray), bl(larray),zh(larray),zl(larray), 
     &         wh(n),wl(n), stat=istat)
      if(istat.ne.0) then
         print*,"Memory exhausted"
         call flush(6)
         stop
      endif

      do i=1,n
         do j=1,n
            call div(1.0d0,0.0d0,dble(i+j-1),0.0d0,temph,templ)
            ah(i,j) = temph
            al(i,j) = templ
         end do
      end do

      call ddmatrix_adjust_s(n, ah, al, bh, bl, nm)

      call ddeigen_s(n, bh, bl, nm, wh, wl, zh, zl, nm, m, 1)

      write(*,*) 'EigenValue=', wh

      deallocate(ah,al)
      deallocate(bh,bl)
      deallocate(zh,zl)
      deallocate(wh,wl)
*
      call MPI_Finalize(ierr)
      end program