Discussion :

[Glouglou69]

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

[FR119492]

Salut!
Pourquoi ne pas utiliser de routines toutes faites telles que celles qu'on trouve dans LAPack?
Jean-Marc Blanc

[Glouglou69]

Malheureusement, le projet que je réalise nécessite de détaillé différentes routines que utilisées.

Si ça peut aider, voilà le pseudo code de l'algorithme utilisé :







Merci bien

[Glouglou69]

Bonsoir à tous,

le problème commence à s'éclaircir
J'avais en effet posé dans la dernière ligne de mon code :

R=matmul(H,A)

Or, avec les transformation faite, il faut plutôt posé (comme l'indique le pseudo code):

R=A

En posant ceci, au lieu d'obtenir les résultat aberrant, j'obtiens une matrice plus ou moins cohérente.

Un nouveau problème ce pose .

En effet, je pose :

A=
12 -51 4
6 167 -68
-4 24 -41

je devrais donc obtenir :

R =
14 21 -14
0 175 -70
0 0 -35

Or avec mon programme j'obtiens :

R =

14.00000 21.07143 -11.42602
0.00000 175.11327 -62.83701
0.00000 0.09281 -41.05160


Il y a une bonne approximation pour les deux premières colonnes mais pour la troisième c'est pas encore ça...

Quelqu'un a-t-il une solution ?
(Je ne sais pas si c'est une erreur de calcul ou une imprécision)
Sinon existe-t-il un algorithme qui retourne les valeurs propre avec une très bonne approximation (même s'il me semblait que la méthode de décomposition QR en donnait déjà une très bonne).

Merci bien.

[Glouglou69]

Bonjour,

Le problème n'était en fait du qu'a la mauvaise initialisation des variables (il faut les réinitialiser à chaque boucle.