SUBROUTINE MMUL (NA,NB,NC,L,M,N,A,B,C)
C
C     *****PARAMETERS:
      INTEGER NA,NB,NC,L,M,N
      DOUBLE PRECISION A(NA,N),B(NB,L),C(NC,L)
C
C     *****LOCAL VARIABLES:
      INTEGER I,J,K
C
C     *****SUBROUTINES CALLED:
C     NONE
C
C     ------------------------------------------------------------------
C
C     *****PURPOSE:
C     THIS SUBROUTINE COMPUTES THE MATRIX PRODUCT A*B AND STORES THE
C        RESULT IN THE ARRAY C.  A IS M X N, B IS N X L, AND C IS
C        M X L.  THE ARRAY C MUST BE DISTINCT FROM BOTH A AND B.
C
C     *****PARAMETER DESCRIPTION:
C     ON INPUT:
C        NA    ROW DIMENSION OF THE ARRAY CONTAINING A AS DECLARED
C              IN THE CALLING PROGRAM DIMENSION STATEMENT;
C
C        NB    ROW DIMENSION OF THE ARRAY CONTAINING B AS DECLARED
C              IN THE CALLING PROGRAM DIMENSION STATEMENT;
C
C        NC    ROW DIMENSION OF THE ARRAY CONTAINING C AS DECLARED
C              IN THE CALLING PROGRAM DIMENSION STATEMENT;
C
C        L     NUMBER OF COLUMNS OF THE MATRICES B AND C;
C
C        M     NUMBER OF ROWS OF THE MATRICES A AND C;
C
C        N     NUMBER OF COLUMNS OF THE MATRIX A AND NUMBER OF ROWS
C              OF THE MATRIX B;
C
C        A     AN M X N MATRIX;
C
C        B     AN N X L MATRIX.
C
C     ON OUTPUT:
C
C        C     AN M X L ARRAY CONTAINING A*B.
C
C     *****HISTORY:
C     WRITTEN BY ALAN J. LAUB (ELEC. SYS. LAB., M.I.T., RM. 35-331,
C     CAMBRIDGE, MA 02139,  PH.: (617)-253-2125), SEPTEMBER 1977.
C     MOST RECENT VERSION: SEP. 21, 1977.
C
C     ------------------------------------------------------------------
C
      DO 40 J=1,L
         DO 10 I=1,M
            C(I,J)=0.0D0
10       CONTINUE
         DO 30 K=1,N
            DO 20 I=1,M
               C(I,J)=C(I,J)+A(I,K)*B(K,J)
20          CONTINUE
30       CONTINUE
40    CONTINUE
      RETURN
C
C     LAST CARD OF MMUL
C
      END