FUNCTION PDDOT(nmax,X,xn,Y,yn,myid,my_veclen,numprocs)
      implicit none
      include 'mpif.h'

      integer :: xn, yn
      integer :: i,nmax,my_veclen,nwork,mystart,myend
      integer :: ierr, myid, numprocs
      integer :: requests(2)
      integer :: statuses(MPI_STATUS_SIZE,2)

      real(kind=8) :: answer, total, totend, DD
      real(kind=8) :: x(nmax), y(nmax)
      DOUBLE PRECISION PDDOT

      DOUBLE PRECISION DDOT
      EXTERNAL DDOT
      total = 0.0
      totend = 0.0

      nwork = nmax/(my_veclen*numprocs)
      write(*,*) "nwork = ", nwork
      
      do i=1,nwork
        !
        mystart = 1 + (i-1)*(my_veclen*numprocs) + myid*my_veclen
        myend = mystart + my_veclen - 1
!        write(*,*)"myid mystart myend = ", myid, mystart, myend

        answer = DDOT( my_veclen, X(mystart:myend), xn, Y(mystart:myend), yn )
!        write(*,*) "myid answer = ", myid, answer

        call MPI_WAITALL(2,requests,statuses,ierr)
        call MPI_REDUCE(answer,totend,1,MPI_DOUBLE_PRECISION,MPI_SUM,0,MPI_COMM_WORLD,ierr)
        total = total + totend
      enddo


      DD = (my_veclen*numprocs)*nwork
!     Finish the rest
      answer = 0.0
      mystart = 1 + DD + myid*my_veclen
      myend = mystart + my_veclen - 1
      if(myend .gt. nmax) myend=nmax
      if (mystart.le.nmax) then
        answer = DDOT( my_veclen, X(mystart:myend), 1, Y(mystart:myend), 1 )
!        write(*,*) "myid answer = ", myid, answer
      endif
      call MPI_WAITALL(2,requests,statuses,ierr)
      call MPI_REDUCE(answer,totend,1,MPI_DOUBLE_PRECISION,MPI_SUM,0,MPI_COMM_WORLD,ierr)

      PDDOT = total + totend

      return
      end function PDDOT