code de decomposition QR (Householder)

[invitee1cdc55e]

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).

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
-----

[invitee1cdc55e]

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é :

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.

[invitee1cdc55e]

C'est tout bon, le problème venait de l'initialisation des variables.

Merci.