dot-product routines

c     Error checking routine
c     call as check_error (MPI_COMM_WORLD, ierror, 'MPI function name')
c     useless as MPICH handles errors internally

      subroutine check_error (comm, ierror, funct)

      implicit none

      include 'mpif.h'

      integer comm, ierror, lierror
      character*20 funct

      if (ierror .ne. 0) then
         write (6,*) 'Error number ', ierror, ' in ', funct
         write (6,*) 'Exiting'
         call MPI_ABORT (comm, ierror, lierror)
         stop
      endif
      return
      end


c**************************************************************************
c     Fortran loop based dot product


      double precision function dotprod (n, a, b, comm) 

      implicit none

      include 'mpif.h'
      
      double precision a(*), b(*)
      double precision temp1, temp2, lab, gab
      integer n, comm, i, ierror

      temp1 = 0.0d0
      temp2 = 0.0d0
      lab = 0.0d0
      gab = 0.0d0
c     local dot product 
c     main loop (doubly unrolled)
      do i = 1, n-1, 2
         temp1 = temp1 + a(i) * b(i)
         temp2 = temp2 + a(i+1) * b(i+1)
      enddo
c     cleanup loop
      if (mod (n, 2) .eq. 1) then
         lab = temp1 + temp2 + a(n) * b(n)
      else
         lab = temp1 + temp2
      endif
c     Sum local dot products to get global value on all processors
      call MPI_ALLREDUCE (lab, gab, 1, MPI_DOUBLE_PRECISION, MPI_SUM,
     $     comm, ierror)

      dotprod = gab
      return 
      end


c**************************************************************************
c     ESSL based dot product


      double precision function blas_dotprod (n, a, b, comm) 

      implicit none

      include 'mpif.h'

      double precision a(*), b(*)
      double precision ddot
      double precision lab, gab
      integer n, comm, ierror

      external ddot

      lab = 0.0d0
      gab = 0.0d0
c     local dot product
      lab = ddot (n, a, 1, b, 1)
c     Sum local dot products to get global value on all processors
      call MPI_ALLREDUCE (lab, gab, 1, MPI_DOUBLE_PRECISION, MPI_SUM,
     $     comm, ierror)

      blas_dotprod = gab
      return 
      end


c**************************************************************************
c     PESSL based dot product


c$$$      double precision function pblas_dotprod (nglobal, a, b, desc) 
c$$$
c$$$      integer desc(*), nglobal
c$$$      double precision lab
c$$$      double precision a(*), b(*)
c$$$
c$$$      pblas_dotprod = PDDOT (nglobal, lab, a, 1, 1, desc, 1, b, 1, 1,
c$$$     $     desc, 1)
c$$$      return
c$$$      end

c**************************************************************************
c     SCALAPACK support routine 

c$$$      subroutine descinit (desc, m, n, mb, nb, rsrc, csrc, ctxt, lda)
c$$$
c$$$      integer desc(8), m, n, mb, nb, rsrc, csrc, ctxt, lda
c$$$
c$$$c     minimal error checking
c$$$
c$$$      if (m .gt. 0) desc(1) = m
c$$$      if (n .gt. 0) desc(2) = n
c$$$      if (mb .gt. 0) desc(3) = mb
c$$$      if (nb .gt. 0) desc(4) = nb
c$$$      if (rsrc .ge. 0) desc(5) = rsrc
c$$$      if (csrc .ge. 0) desc(6) = csrc
c$$$
c$$$      desc(7) = ctxt
c$$$      desc(8) = lda
c$$$      return
c$$$      end