	PROGRAM TEST_LAPACK
!
!=================================================================================
!=================================================================================
! CRIHAN     11/2000
! Guy Moebs   Guy.Moebs@crihan.fr
!
! Utilisation de la bibliotheque scientifique LAPACK (www.netlib.org)
!
! Calcul de la matrice A <- ( U*(tU) + LAMBDA*V*(tV) + tridiag (-1,4,-1) )
!        du vecteur    V <- A  *  ( U + V )
! Resolution de A * U = V
!
! FUNCTION(S) APPELEE(S) :
!     DDOT    produit scalaire de 2 vecteurs                   en double precision 
!     DNRM2   norme L2 d'un vecteur                            en double precision
!
! SUBROUTINE(S) APPELEE(S) :
! BLAS
!     DGER    construction de la matrice A <- a * U * (tV) + A en double precision
!     DAXPY   construction du vecteur Y <- Y + a * X           en double precision
!     DGEMV   calcul du vecteur Y <- a * A * x + b * Y         en double precision
! LAPACK
!     DGESV   factorisation A=PLU et resolution de A * X = Y   en double precision
!
!=================================================================================
!=================================================================================
!
!     DECLARATIONS
!
!=================================================================================
!=================================================================================
!
! MODULE(S) NECESSAIRE(S)
!
      use modscal
!
      IMPLICIT NONE
!
! DECLARATION DES VARIABLES
!
      REAL(rp),  ALLOCATABLE,  DIMENSION(:,:)  ::  A
      REAL(rp),  ALLOCATABLE,  DIMENSION(:,:)  ::  B
      REAL(rp),  ALLOCATABLE,  DIMENSION(:)    ::  U
      REAL(rp),  ALLOCATABLE,  DIMENSION(:)    ::  V
      REAL(rp),  ALLOCATABLE,  DIMENSION(:)    ::  W
!
      REAL(rp)   ::  LAMBDA
      REAL(rp)   ::  ALPHA
      REAL(rp)   ::  BETA
      REAL(rp)   ::  scal
!
      INTEGER,  ALLOCATABLE,  DIMENSION(:)    :: IPIV
!
      INTEGER    ::  N
      INTEGER    ::  INCU
      INTEGER    ::  INCV
      INTEGER    ::  INCW
      INTEGER    ::  I
      INTEGER    ::  J
      INTEGER    ::  ierr = 0
!
      INTEGER    ::  LDA
      INTEGER    ::  LDB
      INTEGER    ::  INFO
      INTEGER    ::  NRHS
!
! DECLARATION DES FUNCTION(S) EXTERIEURE(S)
!
      REAL(rp)   ::  DDOT_
      EXTERNAL   ::  DDOT_
!
      REAL(rp)   ::  DNRM2_
      EXTERNAL   ::  DNRM2_
!
!=================================================================================
!=================================================================================
!
!     DEBUT DU PROGRAMME
!
!=================================================================================
!=================================================================================
!
!----Lecture des dimensions du probleme
!
      OPEN ( 19, FILE = 'lpk.dat', STATUS = 'OLD' )
      READ (19, *) N
      CLOSE (19)
!
!----Allocation dynamique
      ALLOCATE ( A(N,N), B(N,N), U(N), V(N), W(N), IPIV(N), STAT = ierr )
      IF (ierr /= 0) THEN
!
         WRITE (6,'(A)') ' Probleme allocation A, B, U, V, W, IPIV'
         WRITE (6,'(A)') ' ierr = ', ierr
         STOP
!
      END IF
!
!-----Initialisation des vecteurs U et V
!
      scal = 5.0_rp / 12.0_rp
      U (1:N) = scal
!
      DO I = 1, N
         V (I) = 0.2_rp + REAL (I, rp)
      END DO
!
!-----Calcul de LAMBDA = U'*V
!
      INCU = 1
      INCV = 1
!
      LAMBDA = DDOT_ ( N, U, INCU, V, INCV )
!
!-----Initialisation de la matrice A : tridiagonale (-1, 4, -1)
!
      A (1:N, 1:N) = zero
!
      DO I = 1, N
         A (I,I) = 4.0_rp
      ENDDO
!
      DO I = 1, N - 1
	 A (I, I+1) = - un
         A (I+1, I) = - un
      ENDDO
!
!-----Calcul de la matrice A = 1. U * (tU) + A
!
      ALPHA = un
      INCU = 1
!
      CALL DGER_ ( N, N, ALPHA, U, INCU, U, INCU, A, N )
!
!-----Calcul de la matrice A = LAMBDA * V * (tV) + A
!
      INCV = 1
      CALL DGER_ ( N, N, LAMBDA, V, INCV, V, INCV, A, N )
!
!-----Calcul du vecteur V = U + V
!
      ALPHA = un
!
      INCU = 1
      INCV = 1
!
      CALL DAXPY_ ( N, ALPHA, U, INCU, V, INCV )
!
!-----Calcul du U = ALPHA*A*V + BETA*U
!
      ALPHA = un
      BETA  = zero
!
      INCU = 1
      INCV = 1
!
      CALL DGEMV_ ('T', N, N, ALPHA, A, N, V, INCV, BETA, U, INCU)
!
!-----Sauvegarde de A
!
      B (1:N, 1:N) = A (1:N, 1:N)
!
!-----Sauvegarde de U
!
      W (1:N) = U (1:N)
!
!-----Appel a LAPACK pour retrouver la solution V tq A*V=U
!
      INFO = 0
      NRHS = 1
      LDA  = N
      LDB  = N
!
      CALL DGESV_ (N, NRHS, A, LDA, IPIV, U, LDB, INFO)
!
!-----Validation des resultats : residu sur la solution
!
      ALPHA = - un
!
      INCU = 1
      INCV = 1
!
      CALL DAXPY_ (N, ALPHA, U, INCU, V, INCV)
!
      LAMBDA = DNRM2_ ( N, V, INCV)
!
      write (6, 1010) ' NORME L2 DU RESIDU SUR LA SOLUTION : ', LAMBDA
 1010 FORMAT ( A, E12.4 )
!
!-----Validation des resultats : residu sur le second membre
!
      ALPHA = un
      BETA  = zero
!
      INCU = 1
      INCV = 1
!
      CALL DGEMV_ ('T', N, N, ALPHA, B, N, U, INCU, BETA, V, INCV)
!
      ALPHA = - un
!
      INCV = 1
      INCW = 1
!
      CALL DAXPY_ (N, ALPHA, V, INCV, W, INCW)
!
      LAMBDA = DNRM2_ ( N, W, INCW)
!
      WRITE (6, 1010) ' NORME L2 DU RESIDU SUR LE SECOND MEMBRE : ', LAMBDA
!
!-----Liberation de la memoire
!
      DEALLOCATE ( A, B, U, V, W, IPIV )
!
!=================================================================================
!=================================================================================
!
!     FIN DU PROGRAMME
!
!=================================================================================
!=================================================================================
!
      STOP
      END PROGRAM TEST_LAPACK

