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
| program test_gc !!gradient conjugue
implicit none
integer, parameter :: pr=selected_real_kind(10,70)
real(pr),dimension(:,:), allocatable :: AA
real(pr),dimension(:), allocatable :: f,x0,sol, gradconj
integer :: i,j
Allocate( AA(4,4),f(4),x0(4),sol(4),gc(4) )
! Affectation des AA,f,x0
do i=1,4
f(i) = 1.0
x0(i) = 4.0
do j=1,4
AA(i,j) = 4*(i-1)+j
end do
end do
! Printing for check
print*, "AA est :"
do i=1,4
print*, AA(i,:)
end do
print*, "f est : ", f
print*, "x0 est :", x0
! Solving by Gradient Conjugue ---> Problems here !!!
sol = gradconj(AA,f,x0)
print*, "la solution est :"
print*, "Sol = ", sol
Deallocate( AA,f,x0,sol,gc )
end program test_gc
!contains
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
function gradconj(A,b,x0) result (gc)
implicit none
integer, parameter :: pr=selected_real_kind(10,70)
real(pr), dimension(:,:), intent(in) :: A
real(pr), dimension(:), intent(in) :: b, x0
real(pr), dimension(size(x0)) :: gc
integer :: i,j,k
real(pr) :: ff
real(pr), dimension(:,:), allocatable :: r,p,q, x
real(pr), dimension(:), allocatable :: alfa,beta
!! rO = b - (A|x0) ET p0 = r0
do i=1, size(b)
do j= 1, size(A,1)
ff = A(i,j)*x0(j)
end do
r(i,0) = b(i) - ff
p(i,0) = r(i,0)
end do
k=0 !! Nombre d'iterations
do while ( sum(r(:,k)) > 0.001 )
k = k+1
if (k==1) then
p(:,k) = r(:,0)
else
beta(k-1) = dot_product ( r(:,k-1),r(:,k-1) ) / dot_product ( r(:,k-2),r(:,k-2) )
p(:,k) = r(:,k-1) + beta(k-1)*p(:,k-1)
end if
q(i,k) = dot_product( A(i,:),p(:,k) )
alfa(k) = dot_product( r(:,k-1),r(:,k-1) ) / dot_product ( p(:,k),q(:,k) )
x(:,k) = x(:,k-1) + alfa(k)*p(:,k)
r(:,k) = r(:,k-1) - alfa(k)*q(:,k)
end do
gc = x(:,k)
print*, "Solution trouvee! Le nombre d'iterations = ",k
end function gradconj
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! |
