Bonjour à tous,

Je suis en train de réaliser un petit projet en Fortran 90 et pour l'instant j'en suis au calcul des valeurs propres.
Je calcule donc les valeurs propres d'une matrice de différentes façon.

Dans le code ci dessous, je cherche à obtenir la décomposition QR avec la méthode de Householder d'une matrice donné.
Malheureusement, mon programme ne fonctionne pas...

Pour info, j'ai testé l'algo utilisé à la main pour une matrice 3x3 et j'ai trouver que l'erreur venait pour k=2, c'est à dire lors du calcul de Q2 (l'algo trouve le bon résultat pour Q1).
(Sauf si je me suis tromper dans mes calculs ).

Voici le code : (pour une matrice 3x3, mais je le modifierais après pour une matrice n x n).
Code : Sélectionner tout - Visualiser dans une fenêtre à part
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
subroutine decomposition_QR(A,Q,R)
 
  real, dimension(:,:), intent(inout) :: A,Q,R 
  real, dimension(3,3) :: H
  real, dimension(3) :: v
  real :: alpha,b,c,d
  integer :: i,j,k,taille
 
 
  alpha=0
  b=0
  c=0
  d=0
 
  taille=3 !size(A)
 
 
!-----Initialisation de la matrice H, H=I------
  do i=1,taille
     do j=1,taille
        if( i .EQ. J) then
           H(i,j) = 1
        else
           H(i,j) = 0
        end if
     end do 
  end do
write(*,*)H
 
!-----------------------------------------------
 
!----Debut-----
 
do k=1,taille-1
   do i=k,taille
      alpha=alpha+(A(i,k)*A(i,k))
   end do
   alpha=sqrt(alpha)
   b=(alpha*alpha)-(alpha*A(k,k))
 
   !---Construction du vecteur v
   v(k)=A(k,k)-alpha
   do i=k+1,taille
      v(i)=A(i,k)
   end do
 
   do j=k,taille 
      !---Construction de la matrice A**(k+1)----
      do i=k,taille
         c=c+(v(i)*A(i,j))
      end do
      c=c*(1/b)
      do i=k,taille
         A(i,j)=A(i,j)-(c*v(i))
      end do
   end do
 
   do j=1,taille
      !---Construction de H = Hk x ... x H2 x H1
      do i=k,taille
         d=d+(v(i)*H(i,j))
      end do
      d=d*(1/b)
      do i=k,taille 
         H(i,j)=H(i,j)-(d*v(i))
      end do
   end do
end do
 
!---Fin----------
 
Q=transpose(H)
R=matmul(H,A)
 
end subroutine decomposition_QR
Si quelqu'un avait une solution ou une petite aide

Merci beaucoup