PROGRAM MAIN
      IMPLICIT NONE
      INTEGER I,J,K,NTIME,TYP,M
      DOUBLE PRECISION X,ALP
      PARAMETER MM = 500000
C      COMMON /XXVEC/ DOUBLE PRECISION V1(MM),V2(MM),V3(MM),V4(MM)
      DOUBLE PRECISION V1(MM),V2(MM),V3(MM),V4(MM)
 
    
      NTIME = 1000
C  TYP = 2  boucle normale - copy/addition 
C  TYP = 3  utilisation de blas 
      TYP = 2
      M = 50000
      ALP = 2.726
      PRINT *,' Remplissage V1,V2 M=', M, ' NTimes=', NTIME,' Typ=',TYP
      IF (TYP.EQ.2) PRINT *,' TYP=2 --- SANS BLAS'
      IF (TYP.EQ.3) PRINT *,' TYP=3 --- Appel BLAS'
      DO I=1,M
        X = FLOAT(I)*0.00314
        V1(I) = SIN(X)
        V2(I) = COS(X)
      ENDDO

      PRINT *, 'Appel Addition V3 = V1*V2'

      DO K=1,NTIME
        CALL VECADD(TYP,V1,V2,V3,M,ALP)
      ENDDO

      DO J=1,3
        PRINT *, 'J, V1,V2 , V3 = ', J,V1(J),V2(J),V3(J)
      ENDDO 
      DO J=991,993
        PRINT *, 'J, V1,V2 , V3 = ', J,V1(J),V2(J),V3(J)
      ENDDO 
      PRINT *, 'Fin programme Test-Addition'
      STOP
      END

      SUBROUTINE VECADD(TYP,V1,V2,V3,N,ALP)
      INTEGER V1(*), V2(*), V3(*)
      DOUBLE PRECISION ALP
      INTEGER TYP,N,TYPO
      TYPO = 0
      IF (TYP.EQ.1)  THEN
        CALL VECADDI(V1,V2,V3,N,TYPO,ALP)
      ELSE
        IF (TYP.EQ.2)  THEN
          CALL VECADDD2(V1,V2,V3,N,TYPO,ALP)
        ELSE 
          CALL VECADDD(V1,V2,V3,N,TYPO,ALP)
        ENDIF
      ENDIF
      IF ((TYPO.NE.TYP).OR.(TYPO.EQ.0))  THEN
        PRINT *, 'ERREUR TYP-OPeration ' , TYP, TYPO
        STOP
      ENDIF
      RETURN
      END

      SUBROUTINE VECADDD(V1,V2,V3,N,TYP,ALP)
      DOUBLE PRECISION V1(*), V2(*), V3(*),ALP
      INTEGER  N,I,TYP
      INTEGER INCX,INCY
      INCX = 1
      INCY = 1
      CALL DCOPY(N, V2, INCX, V3, INCY)
      CALL DAXPY(N, ALP, V1, INCX, V3, INCY)  
      TYP = 3
      RETURN
      END

      SUBROUTINE VECADDD2(V1,V2,V3,N,TYP,ALP)
      DOUBLE PRECISION V1(*), V2(*), V3(*)
      INTEGER  N,I,TYP
      DO I=1,N
C       V3(I) = V1(I)+V2(I)
       V3(I) = V2(I)
      ENDDO
      DO I=1,N
      V3(I) = V1(I)*ALP+V3(I)
      ENDDO
      TYP = 2
      RETURN
      END

      SUBROUTINE VECADDI(V1,V2,V3,N,TYP,ALP)
      INTEGER V1(*), V2(*), V3(*)
      INTEGER  N,I,TYP
      DOUBLE PRECISION ALP
      DO I=1,N
       V3(I) = V1(I)+V2(I)
      ENDDO
      TYP = 1
      RETURN
      END

      SUBROUTINE VECADDF(V1,V2,V3,N,TYP)
      REAL V1(*), V2(*), V3(*)
      INTEGER  N,I,TYP
      CALL SCOPY(N, V2, 1, V3, 1)
      CALL SAXPY(N, 1., V1, 1, V3, 1)  
      TYP = 2
      RETURN
      END