c************************************************************************

	SUBROUTINE v_propres(a,n,v,x,ok)

c subroutine PUBLIC du module mod_numerique

c calcul des valeurs propres et vecteurs propres de matrices (REAL) 2X2 ou 3X3
c méthode calquée sur le calcul 'à la main'
c ne traite pas le cas de valeurs propres complexes ou multiples

c entrées:
c	A(n,n): matrice
c	1 < n <= 3: dimensions

c sorties
c	v(n): valeurs propres
c	x(n,n): vecteurs propres
c	ok=.TRUE. : C bon

c Auteur: P.Morel, Département Lagrange, O.C.A., CESAM2k

c----------------------------------------------------------------------

	USE mod_kind

	IMPLICIT NONE

	REAL(kind=dp), INTENT(in), DIMENSION(n,n) :: a
	INTEGER, INTENT(in) :: n

	REAL(kind=dp), INTENT(out), DIMENSION(n,n) :: x
	REAL(kind=dp), INTENT(out), DIMENSION(n) :: v
	LOGICAL, INTENT(out) :: ok

	REAL(kind=dp), DIMENSION(n+1,n) :: w
	REAL(kind=dp), DIMENSION(0:n) :: poly
        REAL(kind=dp), DIMENSION(1,n+1) :: b

        INTEGER, DIMENSION(n) :: indpc
	INTEGER :: i, j

	LOGICAL :: inversible

c---------------------------------------------------------------------

2000	FORMAT(8es10.3)

c matrices 2X2
	SELECT CASE(n)
	CASE(2)

c polynôme caractéristique
	 poly(0)=a(1,1)*a(2,2)-a(1,2)*a(2,1)
	 poly(1)=-(a(1,1)+a(2,2))
	 poly(2)=1.d0

c matrices 3X3
	CASE(3)

c polynôme caractéristique
         poly(0)=a(1,1)*(a(2,2)*a(3,3)-a(2,3)*a(3,2))-a(2,1)*(a(1,2)*a(3,3)-a(1,3)*a(3,2))
	1 +a(1,3)*(a(1,2)*a(2,3)-a(1,3)*a(2,1))
	 poly(1)=(a(1,1)*a(2,2)-a(2,1)*a(1,2))+(a(1,1)*a(3,3)-a(1,3)*a(3,1))
	1 +(a(2,2)*a(3,3)-a(2,3)*a(3,2))
	 poly(1)=-poly(1)
	 poly(2)=a(1,1)+a(2,2)+a(3,3)
	 poly(3)=-1.d0

c cas n /= é ou 3
	CASE DEFAULT
	 WRITE(*,1)n
1	 FORMAT('n=',i3,", v_propres n'est prévu que pour les  matrices 2X2 ou 3X3, désolé")
	 v=0.d0 ; x=0.d0 ; ok=.FALSE.
	 RETURN
	END SELECT

c valeurs propres
	CALL maehly(poly,n,n,v)
c	PRINT* ; PRINT*,'valeurs propres'
c	WRITE(*,2000)v

c x vecteur propre associé à v(k)
	DO i=1,n
	 w=0.d0
	 w(1:n,:)=a(1:n,:) ; b=0.d0 ; indpc=1
	 DO j=1,n
          w(j,j)=a(j,j)-v(i)
         ENDDO

c on impose 1 pour la 1-ière coordonnée des vecteurs propres
c le rang de la matrice aux vecteurs propres étant inférieur à n on ajoute une ligne fictive
         w(n+1,1)=1.d0 ; b(1,n+1)=1.d0
c        PRINT* ; PRINT*,'w,b pour calcul du vecteur propre n°',i
c	 DO j=1,n
c	  WRITE(*,2000)w(j,:),b(1,j)
c	 ENDDO

c gauss_band permet de résoudre des systèmes dont le rang est inférieur à n
	 CALL gauss_band(w,b,indpc,n+1,n,n,1,inversible)
c        PRINT* ; PRINT*,'vecteur propre n°',i
         x(i,1:n)=b(1,1:n)
        ENDDO

c C bon
        ok=.TRUE.

 	RETURN

	END SUBROUTINE v_propres