Update

CMakeLists.txt

@@ -3,4 +3,4 @@ project(GaussSeidel Fortran)
 
 
 enable_language(Fortran)
-add_executable(GaussSeidel main.f gauss_seidel.f)
\ No newline at end of file
+add_executable(GaussSeidel main.f gauss_seidel.f tril.f90 triu.f90)
\ No newline at end of file

Makefile

 FC = gfortran
 OPT = -g
 
-main: main.f gauss_seidel.o
+gss: main.f gauss_seidel.o
 	$(FC) $(OPT) main.f gauss_seidel.o -o main
 
 gauss_seidel.o: gauss_seidel.f

dgss

+2  1.e-6  100
+3 1 2 4
+3 3
+1 1
+n, prec,nmax
+A
+b
+x

gauss_seidel.f

-      module gauss_seidel
-        implicit none
-        private
-        public :: solve_system
-
-        interface solve_system
-          procedure solve_system_double, solve_system_real
-        endinterface
-      contains
-
-        subroutine solve_system_double(A, b, x)
-          implicit none
-          integer :: An(2), n(1), i, j
-          double precision :: A(:,:), b(:), x(:), tmp
-          n = shape(b)
-
-          do i=1,n(1)
-            tmp = 0
-            do j=1,n(1)
-              if (j /= i) then
-                tmp = tmp + A(i,j) * x(j)
-              endif
-            enddo
-            x(i) = (b(i) - tmp) / A(i, i)
-          enddo
-
-        end subroutine solve_system_double
-
-        subroutine solve_system_real(A, b, x)
+      subroutine gauss_seidel(A, n, nmax, b, x)
           implicit none
-          integer :: n(1), i, j
-          real :: A(:,:), b(:), x(:), tmp
-          n = shape(b)
+          integer :: n, nmax, i, j
+          double precision :: A(nmax,nmax), b(nmax), x(nmax), tmp
 
           do i=1,n(1)
             tmp = 0
@@ -41,6 +12,4 @@
             enddo
             x(i) = (b(i) - tmp) / A(i, i)
           enddo
-
-        end subroutine solve_system_real
-      end module
\ No newline at end of file
+      end
\ No newline at end of file

inv_trig_inf.f90

+function inv_trig_inf(L, nmax, n)
+    integer :: nmax, n, i, j
+    real :: L(nmax, nmax), inv_trig_inf(nmax, nmax), diag
+    inv_trig_inf = 0
+    ! inv(1,1) = 1/L(1,1)
+    ! inv(1,2) =
+    ! inv(2,1) = 0
+    diag = 1
+    do i = 1,n
+        inv_trig_inf(i,i) = 1/L(i,i)
+        diag = diag * L(i,i)
+        x(l) = b(l)
+        do j = 1, n
+            x(l) = x(l) - L(l, c) * x(c)
+        enddo
+        x(l) = x(l) / L(l, l)
+    enddo
+end
\ No newline at end of file

main.f

       program main
-      use gauss_seidel
       implicit none
-      integer i, n
-      real, allocatable :: A(:,:), b(:), x(:)
-	  real :: prec
-	  
-	  read (*,*) n
-	  read (*,*) prec
-	  allocate(A(n,n))
-	  allocate(b(n))
-	  allocate(x(n))
-	  read (*,*) A
-	  read (*,*) b
-	  read (*,*) x
+      integer, parameter :: nmax=5
+      integer i, j, n
+      real :: A(nmax,nmax), F(nmax,nmax), G(nmax,nmax), D1(nmax,nmax), L(nmax,nmax), LI(nmax,nmax), U(nmax,nmax), Id(nmax,nmax), b(nmax), x(nmax)
+      real :: prec,res
+  
+      read (*,*) n, prec
 
-      do i=1,3000000
+      read (*,*) A
+      read (*,*) b
+      read (*,*) x
+
+      D1 = 0
+      Id = 0
+
+      do i=1,n
+        D1(i,i) = 1.0/A(i,i)
+        Id(i,i) = 1
+      enddo
+      do i=1,n
+        do j=1,n
+          L(i, j) = A
+        enddo
+      enddo
+
+      L = tril(matmul(D1, A), nmax, n) - Id
+      U = triu(matmul(D1, A), nmax, n) - Id
+
+      LI = inv(L+Id)
+
+      F = -matmul(LI, U)
+      G = matmul(matmul(LI, D1), b)
+
+      do i=1,nmax
         call solve_system(A, b, x)
-        if (norm2(matmul(A, x) - b) < prec) then
+        res = norm2(matmul(A, x) - b)
+        write(*,*) i, res
+        if ( res < prec) then
           exit
         endif
       enddo
       
-	  write(*,*) x
-      write(*,*) matmul(A, x)
-	  end
+      end

src/Makefile

+FC = gfortran
+OPT = -g
+
+gss : main.f gauss_seidel.o
+	$(FC) $(OPT) main.f gauss_seidel.o -o gss 
+
+gauss_seidel.o: gauss_seidel.f
+	$(FC) $(OPT) gauss_seidel.f -c
+
+clean:
+	rm gss  *.o *.mod *.exe

src/dgss

+2  1.e-6  100
+3 1 2 4
+3 3
+1 1
+n, prec,nmax
+A
+b
+x

src/gauss_seidel.f

+      module gauss_seidel
+        implicit none
+        private
+        public :: solve_system
+
+        interface solve_system
+          procedure solve_system_double, solve_system_real
+        endinterface
+      contains
+
+        subroutine solve_system_double(A, b, x)
+          implicit none
+          integer :: An(2), n(1), i, j
+          double precision :: A(:,:), b(:), x(:), tmp
+          n = shape(b)
+
+          do i=1,n(1)
+            tmp = 0
+            do j=1,n(1)
+              if (j /= i) then
+                tmp = tmp + A(i,j) * x(j)
+              endif
+            enddo
+            x(i) = (b(i) - tmp) / A(i, i)
+          enddo
+
+        end subroutine solve_system_double
+
+        subroutine solve_system_real(A, b, x)
+          implicit none
+          integer :: n(1), i, j
+          real :: A(:,:), b(:), x(:), tmp
+          n = shape(b)
+
+          do i=1,n(1)
+            tmp = 0
+            do j=1,n(1)
+              if (j /= i) then
+                tmp = tmp + A(i,j) * x(j)
+              endif
+            enddo
+            x(i) = (b(i) - tmp) / A(i, i)
+          enddo
+
+        end subroutine solve_system_real
+      end module
\ No newline at end of file

src/gauss_seidel.m

+%A = [1 2 3 4 5 ; 5 1 2 3 4 ; 4 5 1 2 3 ; 3 4 5 1 2 ; 2 3 4 5 1];
+%Abis = [1 -1 -200 0 ; 0 1 100 -100 ; 0 1 101 -101 ; 0 0 0 100];
+%b = [1 ; 1 ; 3 ; 4 ; 2];
+%bbis = [1 ; 1 ; 2 ; 4];
+%x = [0 ; 0 ; 0 ; 0 ; 0];
+%xbis = [0 ; 0 ; 0 ; 0];
+
+
+tol = 1e-10; 
+cont = 0;
+delta = 10; 
+
+A = [3 1;2 4];
+b = [3;3];
+x = [1;1];
+
+n = size(b);
+
+while delta > tol  & cont < 200
+  ref = x;
+  for i = 1:n(1)
+    tmp = 0;
+    for j = 1:n(1)
+      if j != i
+        tmp = tmp + A(i,j)*x(j);
+      end
+      x(i) = (b(i) - tmp) / A(i,i);
+    end
+  end
+  delta = norm(ref - x);
+  cont = cont + 1;
+end
+
+fprintf('%f\n%f\n%f\n%f\n en %d boucles',x, cont);
\ No newline at end of file

src/main.f

+      program main
+      use gauss_seidel
+      implicit none
+      integer i, n,nmax
+      real, allocatable :: A(:,:), b(:), x(:)
+      real :: prec,res
+  
+      read (*,*) n, prec,nmax
+      allocate(A(n,n))
+      allocate(b(n))
+      allocate(x(n))
+
+      read (*,*) A
+      read (*,*) b
+      read (*,*) x
+
+      do i=1,nmax
+        call solve_system(A, b, x)
+        res = norm2(matmul(A, x) - b)
+        write(*,*) i, res
+        if ( res < prec) then
+          exit
+        endif
+      enddo
+      
+      end

tril.f90

+function tril(M, nmax, n)
+    integer :: nmax, n
+    real M(nmax,nmax), tril(nmax, nmax)
+    tril = 0
+    do i=1,n
+        do j=1,i
+            tril(i,j) = M(i,j)
+        end do
+    end do
+end function tril
\ No newline at end of file

triu.f90

+function triu(M, nmax, n)
+    integer :: nmax, n
+    real M(nmax,nmax), triu(nmax, nmax)
+    triu = 0
+    do i=1,n
+        do j=i,n
+            triu(i,j) = M(i,j)
+        end do
+    end do
+end function triu
\ No newline at end of file