Discussion :

[cedrix57]

Bonjour,

Dans le cadre d'un projet scolaire de calcul numérique, je suis en train d'implémenter l'algorithme GMRES.
J'ai un problème de mémoire que je n'arrive pas à résoudre. Quelqu'un aurait une idée de ce qui ne va pas dans mon programme :
Voicl la version de gfortran que j'utilise :

Code :

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

1
2
3
4
5
6
7
8
~/projet/fortran/src$ gfortran --version
GNU Fortran (Debian 4.3.2-1) 4.3.2
Copyright (C) 2008 Free Software Foundation, Inc.
 
GNU Fortran comes with NO WARRANTY, to the extent permitted by law.
You may redistribute copies of GNU Fortran
under the terms of the GNU General Public License.
For more information about these matters, see the file named COPYING

Voicl ce qu'affiche mon programme et l'erreur que j'obtiens :

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
~/projet/fortran/src$ ./arnoldi 
 INFO : variables intialisee
 A :
  1.00000000000000000        2.0000000000000000        3.0000000000000000     
   4.0000000000000000        5.0000000000000000        6.0000000000000000     
   7.0000000000000000        8.0000000000000000        9.0000000000000000     
 A :
  1.00000000000000000        2.0000000000000000        3.0000000000000000     
   4.0000000000000000        5.0000000000000000        6.0000000000000000     
   7.0000000000000000        8.0000000000000000        9.0000000000000000     
 Q :
  0.19611613513818404       0.11623161232356750       0.97366764033479380     
  0.78446454055273618      -0.61436709371028553      -8.46667513334609795E-002
  0.58834840541455213       0.78041225417252458      -0.21166687833365103     
 H :
   14.769230769230770        3.5885676911335862        4.2507561078333183     
   5.4512362204697169        1.1483319547835669        1.9710023043995237     
   0.0000000000000000       0.40207322329405470      -0.91756272401433692     
   0.0000000000000000        0.0000000000000000       9.63326575041553881E-016
 *****************************************************
 do while :            1
 end do 1
 do while :            2
*** glibc detected *** ./arnoldi: free(): invalid next size (fast): 0x080665b8 ***
======= Backtrace: =========
/lib/libc.so.6[0xb6f37905]
/lib/libc.so.6[0xb6f391a3]
/lib/libc.so.6(cfree+0x6d)[0xb6f3c1fd]
./arnoldi[0x804a915]
./arnoldi[0x8048df6]
./arnoldi[0x804b099]
/lib/libc.so.6(__libc_start_main+0xe5)[0xb6ee2b55]
./arnoldi[0x80488f1]
======= Memory map: ========
08048000-0804c000 r-xp 00000000 00:12 88728223   /home/gis3/cbetting/projet/fortran/src/arnoldi
0804c000-0804d000 rw-p 00003000 00:12 88728223   /home/gis3/cbetting/projet/fortran/src/arnoldi
0804d000-08081000 rw-p 00000000 00:00 0          [heap]
b6d00000-b6d21000 rw-p 00000000 00:00 0 
b6d21000-b6e00000 ---p 00000000 00:00 0 
b6e43000-b6e45000 rw-p 00000000 00:00 0 
b6e45000-b6ecb000 r-xp 00000000 08:01 427491     /usr/lib/libblas.so.3gf.0
b6ecb000-b6ecc000 rw-p 00086000 08:01 427491     /usr/lib/libblas.so.3gf.0
b6ecc000-b700d000 r-xp 00000000 08:01 474448     /lib/libc-2.10.2.so
b700d000-b700f000 r--p 00141000 08:01 474448     /lib/libc-2.10.2.so
b700f000-b7010000 rw-p 00143000 08:01 474448     /lib/libc-2.10.2.so
b7010000-b7013000 rw-p 00000000 00:00 0 
b7013000-b701f000 r-xp 00000000 08:01 407475     /lib/libgcc_s.so.1
b701f000-b7020000 rw-p 0000b000 08:01 407475     /lib/libgcc_s.so.1
b7020000-b7044000 r-xp 00000000 08:01 474444     /lib/libm-2.10.2.so
b7044000-b7045000 r--p 00023000 08:01 474444     /lib/libm-2.10.2.so
b7045000-b7046000 rw-p 00024000 08:01 474444     /lib/libm-2.10.2.so
b7046000-b7047000 rw-p 00000000 00:00 0 
b7047000-b70f7000 r-xp 00000000 08:01 407519     /usr/lib/libgfortran.so.3.0.0
b70f7000-b70f8000 rw-p 000af000 08:01 407519     /usr/lib/libgfortran.so.3.0.0
b70f8000-b70f9000 rw-p 00000000 00:00 0 
b70f9000-b772a000 r-xp 00000000 08:01 388109     /usr/lib/liblapack.so.3gf.0
b772a000-b772d000 rw-p 00631000 08:01 388109     /usr/lib/liblapack.so.3gf.0
b772d000-b783a000 rw-p 00000000 00:00 0 
b7854000-b7856000 rw-p 00000000 00:00 0 
b7856000-b7872000 r-xp 00000000 08:01 474453     /lib/ld-2.10.2.so
b7872000-b7873000 r--p 0001b000 08:01 474453     /lib/ld-2.10.2.so
b7873000-b7874000 rw-p 0001c000 08:01 474453     /lib/ld-2.10.2.so
bfaa4000-bfab9000 rw-p 00000000 00:00 0          [stack]
ffffe000-fffff000 r-xp 00000000 00:00 0          [vdso]
Aborted

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
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
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
      program main
      REAL*8 Q(3,3), H(4,3)
      REAL*8 A(3,3), b(3), x(3), res(3)
      integer n_max, n_end, m
      real*8 eps_max, eps_end
      n_max = 3
      m = 3
      eps_max = 0.001
      A(1,:)=(/ 1, 2, 3 /)
      A(2,:)=(/ 4, 5, 6 /)
      A(3,:)=(/ 7, 8, 9 /)
      b = (/ 1, 4, 3 /)
      x = (/ 0, 0, 0 /)
      write(*,*) "INFO : variables intialisee"
      write(*,*) "A :"
      call printMatrix(A, 3, 3)
      call arnoldi(A, Q, H, b, 3, 3)
      write(*,*) "A :"
      call printMatrix(A, 3, 3)
      write(*,*) "Q :"
      call printMatrix(Q, 3, 3)
      write(*,*) "H :"
      call printMatrix(H, 4, 3)
      write(*,*) "*****************************************************"
      call gmres(A, b, x, eps_max, n_max, eps_end, n_end, m)
      write(*,*) "INFO : gmres execute"
      res = matmul(A,x)
      write(*,*) "A vaut :"
      call printMatrix(A, 3, 3)
      write(*,*) "res vaut :"
      write(*,*) (res(i), i=1,m)
 
      write(*,*) "INFO : end program"
      stop
      end
 
c     -------------------------------------------------------------
c     ----------------------------arnoldi--------------------------
c     -------------------------------------------------------------
      subroutine arnoldi(A, Q, H, b, m, n)
      integer i, j, k, m, n
      REAL*8 A(m,m), Q(m,n), H(n+1,n), b(m), omega(m)
      Q(:,1) = b/sqrt(dot_product(b,b))
      H = 0
      do k = 1,n
      	omega = MATMUL(A, Q(:,k))
      	do j = 1,k
      		H(j, k) = DOT_PRODUCT(Q(:,j), omega)
      		omega = omega - H(j, k)*Q(:,j)
      	end do
      	H(k+1, k) = SQRT(DOT_PRODUCT(omega, omega))
      	if (abs(H(k+1, k)) .GT. 0.000000000000000000000000001) then
      		omega = omega/H(k+1,k)
      	else
      		write(*,*) "WARNING : H(k+1, k) :", H(k+1, k), "k :", k
      		H(k+1, k) = 0
      	end if
      	if (k.LT.n) then
      		Q(:,k+1)=omega
      	end if
      end do
      return
      end
 
c     -------------------------------------------------------------
c     ----------------------------GMRES----------------------------
c     -------------------------------------------------------------
      subroutine gmres(A, b, x, eps_max, n_max, eps_end, n_end, m)
      real*8 A(m,m), b(m), x(m), r(m), y(n_max+1)
      real*8 A_x(m), ecart(m), be(n_max+1)
      real*8 Q(m,n_max), H(n_max+1, n_max)
      real eps_max, eps_end, eps, norme_b
      integer n_max, n_end, m, i
c     Variables pour la resolution du systeme lineaire
      integer info, lda, ldb
      real*8 work(10000)
c     LWORK >= 3*min(M,N)*min(M,N) + max(max(M,N),4*min(M,N)*min(M,N)+4*min(M,N)).
      lwork=10000
      lda=n_max+1
      ldb=n_max+1
      i = 1
c     on normalise b si necessaire
      norme_b = sqrt(dot_product(b, b))
      if (abs(norme_b-1).GT.0.000001) then
      	b=b/norme_b
      end if
      A_x = matmul(A, x)
      ecart = A_x-b
      eps = sqrt(dot_product(ecart, ecart))
      be = 0
      be(1) = eps
      do while(.NOT.((eps.LE.eps_max).OR.(i.EQ.n_max+1)))
      	write(*,*) "do while : ", i
      	call arnoldi(A, Q(:,:i), H(:i+1,:i), ecart/eps, m, i)
      	y=be
      	call dgels('No transpose',i+1,i,1,H(:i+1,:i),lda,y(:i+1),ldb,
     s              work, lwork,info)
      	if(info.LT.0) then
      		write(*,*) "ERROR : sgels info<0, argument pb :", -info
      	end if
      	x=matmul(Q(:,:i), y(:i))
      	A_x=matmul(A, x)
      	ecart = A_x-b
      	eps = sqrt(dot_product(ecart, ecart))
      	be(1) = eps
      	i=i+1
      	write(*,*) "end do 1"
      end do
      eps_end=eps
      n_end=i
c     on ajuste la valeur de x si b avait ete normalise
      if (abs(norme_b-1).GT.0.000001) then
      	x=x*norme_b
      end if
      write(*,*) "INFO : fin de la methode"
      return
      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

edit : j'ai oublié de dire que j'utilise la méthode dgels de LAPACK