C.hr dgmprd
C@
*....*...1.........2.........3.........4.........5.........6.........7.*
*     dgmprd  6/15/83, 12/27/94
*
*     purpose
*     multiply two matrices: R = AB
*
*     usage
*     call dgmprd(A,B,R,n,m,l)
*
*     arguments
*     A - input n by m matrix stored columwise, i.e. A must either be
*         dimensioned as A(n,m) in the calling program or be a vector
*         that is indexed as A(-n+n*j+i) where i is the row index and
*         j is the column index.
*         input, real*8
*     B - input m by l matrix stored columnwise.
*         input, real*8
*     R - output n by l matrix stored columnwise.
*         output, real*8
*     n - number of rows in A and R.
*         input, integer*4
*     m - number of columns in A and number of rows in B.
*         input, integer*4
*     l - number of columns in B and R.
*         input, integer*4

      subroutine dgmprd(A,B,R,n,m,l)
      implicit none

*     arguments

      real*8 A,B,R
      integer*4 n,m,l

      dimension A(n,m),B(m,l),R(n,l)

*     local variables

      integer*4 i,j,k

*     clear

      do j=1,l
      do i=1,n
        R(i,j)=0.d0
      end do
      end do

*     accumulate

      do j=1,l
      do k=1,m
      do i=1,n
        R(i,j)=R(i,j)+A(i,k)*B(k,j)
      end do
      end do
      end do

      return
      end