Discussion :

[cedrix57]

Bonjour,

J'ai un problème pour utiliser la fonction dgels, qui est censé résoudre un système des moindres carrées. Voici mon code source :

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
 
      program main
      REAL*8 A(4,3), b(4), res(4)
      integer n_max, n_end, m
 
      integer info, lda, ldb 
c     LWORK >= 3*min(M,N)*min(M,N) + max(max(M,N),4*min(M,N)*min(M,N)+4*min(M,N)).
      real*8 work(1000)
      lwork=1000
      m = 3 
      lda=m+1
      ldb=m+1
      m = 2 
      A(1,:)=(/ 1, 2, 3 /)
      A(2,:)=(/ 4, 5, 6 /)
      A(3,:)=(/ 7, 8, 9 /)
      A(4,:)=(/ 9, 9, 9 /)
      b = (/ 1, 4, 3, 2 /)
      write(*,*) "A :"
      call printMatrix(A(:m+1,:m), m+1, m)
      write(*,*) "b vaut :"
      write(*,*) (b(i), i=1,m+1)
      call dgels('No transpose',m+1,m,1,A,lda,b,ldb,
     s              work, lwork,info)
      write(*,*) "INFO : dgels end"
      write(*,*) "res vaut :"
      write(*,*) (b(i), i=1,m)
 
      write(*,*) "INFO : end program"
      stop
      end 
 
      subroutine printMatrix(mat, m, n)
      integer i, j, m, n
      real*8 mat(m,n)
      do i =1,m
      >-WRITE(*,*) (mat(i,j), j=1,n)
      end do
      return
      end

J'essai de trouver un ymin tel que ||A*ymin - b|| soit minimiser. Si je prend A en entier, en mettant m = 3, au lieu de m = 2, ça marche, j'obtiens ce résultat qui semble être cohérent :

Code :

Sélectionner tout - Visualiser dans une fenêtre à part

1
2
3
4
5
6
7
8
9
10
11
12
~/projet/fortran/test$ gfortran -o miny miny.f -llapack;./miny
 A :
  1.00000000000000000        2.0000000000000000        3.0000000000000000     
   4.0000000000000000        5.0000000000000000        6.0000000000000000     
   7.0000000000000000        8.0000000000000000        9.0000000000000000     
   9.0000000000000000        9.0000000000000000        9.0000000000000000     
 b vaut :
  1.00000000000000000        4.0000000000000000        3.0000000000000000        2.0000000000000000     
 INFO : dgels end
 res vaut :
   306553334728091.56       -613106669456184.75        306553334728093.44     
 INFO : end program

Par contre, quand je met m = 2, le résultat n'est plus cohérent :

Code :

Sélectionner tout - Visualiser dans une fenêtre à part

1
2
3
4
5
6
7
8
9
10
11
~/projet/fortran/test$ gfortran -o miny miny.f -llapack;./miny
 A :
  1.00000000000000000        2.0000000000000000     
   4.0000000000000000        5.0000000000000000     
   7.0000000000000000        8.0000000000000000     
 b vaut :
  1.00000000000000000        4.0000000000000000        3.0000000000000000     
 INFO : dgels end
 res vaut :
 -0.99999999999999567        1.3333333333333295     
 INFO : end program

La solution est loin d'être optimale, par exemple le vecteur (0, 0.5) est une bien meilleur solution.

J'ai trouvé un exemple d'utilisation de cette fonction ici :
http://www.nag.co.uk/lapack-ex/examples/source/dgels-ex.f
Il y a aussi quelque infos ici:
http://www.netlib.org/lapack/lug/node27.html
http://download.oracle.com/docs/cd/B...a.htm#CIADBIAE


Auriez vous une idée de ce qui ne va pas ?

Merci d'avance,

Cédric

[cedrix57]

Désolé c'est moi que m'étais trompé lors de la vérification des calculs.