      subroutine daxa(a,b,c,m,n,l)
      implicit none
      double precision a(m,n), b(n,l), c(m,l)
      integer m,n,l
      integer i,j,k
      
      do i=1,m
         do k=1,l
            c(i,k)=0
            do j=1,n
               c(i,k)=c(i,k)+a(i,j)*b(j,k)
            end do
         end do
      end do
      end ! subroutine mmult

c     ..................................................................
      double precision function d_dot(a, inca, b, incb, n)
      implicit none
      integer inca, incb, n
      double precision a(n), b(n)
      integer i, aend, bend
      
      d_dot = 0d0
      if (n .le. 0) return
      if (inca .lt. 0) then
         aend=n
      else
         aend=1
      end if

      if (incb .lt. 0) then
         bend=n
      else
         bend=1
      end if

      do i=1,n
         d_dot = d_dot + a( (i-aend)*inca+1 ) * b( (i-bend)*incb+1 )
      end do

      end ! double precision function d_dot
