C MULMAT    SOURCE    PV        18/10/17    11:46:14     9965           
      SUBROUTINE MULMAT(A,B,C,I,J,K)
C=======================================================================
C
C  S.P. PRODUIT DE 2 MATRICES A(I,J)=B(I,K)*C(K,J)
C
C    AOUT 84
C
C  Version optimise pour la bande passante memoire
C  Unroll sur 4 et 6 car ce sont les dimensions des contraintes
C=======================================================================
      IMPLICIT INTEGER(I-N)
      IMPLICIT REAL*8(A-H,O-Z)
      
      DIMENSION A(I,*),B(I,*),C(K,*)
      
      PARAMETER (XZERO=0.D0)

**    print *, ' mulmat i j k ',i,j,k
      DO JJ=1,J
        DO II=1,I
          A(II,JJ)=XZERO
        ENDDO
        KK=1
        DO KK=1,K-5,6
          XFLOT1=C(KK  ,JJ)
          XFLOT2=C(KK+1,JJ)
          XFLOT3=C(KK+2,JJ)
          XFLOT4=C(KK+3,JJ)
          XFLOT5=C(KK+4,JJ)
          XFLOT6=C(KK+5,JJ)

          DO II=1,I
            A(II,JJ)=A(II,JJ)+B(II,KK  )*XFLOT1
     >                       +B(II,KK+1)*XFLOT2
     >                       +B(II,KK+2)*XFLOT3
     >                       +B(II,KK+3)*XFLOT4
     >                       +B(II,KK+4)*XFLOT5
     >                       +B(II,KK+5)*XFLOT6
          ENDDO          
        ENDDO

         KKK=KK                
         DO KK=KKK,K-3,4
          XFLOT1=C(KK  ,JJ)
          XFLOT2=C(KK+1,JJ)
          XFLOT3=C(KK+2,JJ)
          XFLOT4=C(KK+3,JJ)

          DO II=1,I
            A(II,JJ)=A(II,JJ)+B(II,KK  )*XFLOT1
     >                       +B(II,KK+1)*XFLOT2
     >                       +B(II,KK+2)*XFLOT3
     >                       +B(II,KK+3)*XFLOT4
          ENDDO          
         ENDDO

         KKK=KK                
         DO KK=KKK,K-2,3
          XFLOT1=C(KK  ,JJ)
          XFLOT2=C(KK+1,JJ)
          XFLOT3=C(KK+2,JJ)

          DO II=1,I
            A(II,JJ)=A(II,JJ)+B(II,KK  )*XFLOT1
     >                       +B(II,KK+1)*XFLOT2
     >                       +B(II,KK+2)*XFLOT3
          ENDDO          
         ENDDO

         KKK=KK                
         DO KK=KKK,K-1,2
          XFLOT1=C(KK  ,JJ)
          XFLOT2=C(KK+1,JJ)

          DO II=1,I
            A(II,JJ)=A(II,JJ)+B(II,KK  )*XFLOT1
     >                       +B(II,KK+1)*XFLOT2
          ENDDO          
         ENDDO

        KKK=kk
        DO KK=kkK,K
          XFLOT1=C(KK,JJ)
          DO II=1,I
            A(II,JJ)=A(II,JJ)+B(II,KK)*XFLOT1
          ENDDO          
        ENDDO
      ENDDO
         
      END