      PROGRAM MATR_MULT
C
C.....Paradeigma xrnsns upoprogrammatwv tou BLAS 
C
      PARAMETER (NDIM = 4)
      IMPLICIT REAL*8 (A-H,O-Z)
      DIMENSION A(NDIM,NDIM),B(NDIM,NDIM),C(NDIM,NDIM),X(NDIM),Y(NDIM)
      DO 10 I = 1,NDIM
        X(I) = DFLOAT(I)
        DO 20 J = 1,NDIM
          A(I,J) = DFLOAT(I) * DFLOAT(J) - 8.0D0
          B(I,J) = DFLOAT(I) + DFLOAT(J)
 20     CONTINUE
 10   CONTINUE
      PRINT*,' O pivakas A eivai:'
      DO 30 I = 1,NDIM
        WRITE(*,*) (A(I,J),J = 1,NDIM)
 30   CONTINUE
      PRINT*
      PRINT*,' O pivakas B eivai:'
      DO 40 I = 1,NDIM
        WRITE(*,*) (B(I,J),J = 1,NDIM)
 40   CONTINUE
      PRINT*
      PRINT*,' To diavusma X eivai:'        
      WRITE(*,*) (X(I),I = 1,NDIM)
      PRINT*
      CALL DGEMM('N','N',NDIM,NDIM,NDIM,1.0d0,A,NDIM,B,NDIM,
     &           0.0d0,C,NDIM) 
      PRINT*,' O pivakas C = A*B eivai:'     
      DO 50 I = 1,NDIM
        WRITE(*,*) (C(I,J),J = 1,NDIM)
 50   CONTINUE
      PRINT*
      CALL DGEMV('N',NDIM,NDIM,1.0d0,C,NDIM,X,1,0.0d0,Y,1)
      PRINT*,' To diavusma Y = A*B eivai:'        
      WRITE(*,*) (Y(I),I = 1,NDIM)
      STOP
      END
