Jacobi Method Implementation in Fortran

This section details the implementation of the Jacobi method in Fortran for solving a system of linear equations. The subroutine solucion_Jacobi takes a matrix A, a vector b, and an initial guess X as input, and returns the solution vector X.

Subroutine: solucion_Jacobi(A, b, X)

The subroutine begins with the following declarations and initializations:

• use normas
• integer, parameter :: Kx = 100
• real(8), parameter :: tolerancia = 1E-6
• real(8), allocatable :: X0(:), ERROR(:)
• real(8), intent(in) :: A(:,:), b(:)
• real(8), allocatable, intent(out) :: X(:)
• integer :: i, j, k, n
• real(8) :: sum, zero, ek

Where:

• Kx is the maximum number of iterations.
• tolerancia is the error tolerance.
• X0 is the previous iteration's solution vector.
• ERROR is the error vector.
• A is the coefficient matrix.
• b is the right-hand side vector.
• X is the solution vector.

The size of the system n is determined, memory is allocated for X, X0, and ERROR, and X0 is initialized to zero:

n = size(A, dim = 1)
zero = 0d0
allocate(X(1:N), X0(1:N), ERROR(1:N))
X0 = 0d0

Condition Check

The code checks if the matrix A satisfies the condition for convergence. For each row i, it verifies that the absolute value of the diagonal element A(i,i) is greater than the sum of the absolute values of the other elements in that row:

do i = 1, N
 sum = 0d0
 do j = 1, N
  if (i == j) then
   cycle
  else
   sum = sum + (A(i,j))
  end if
 end do
 if ((A(i,i)) <= sum) then
  write(*,*) "No cumple la condicion"
  write(*,*) A(i,i), "<", sum
  !stop
 else
  write(*,*) "Cumple la condicion"
 end if
end do

Diagonal Invertibility Check

The code checks if the diagonal elements of matrix A are non-zero to ensure invertibility:

do i = 1, N
 if (abs(A(i,i)) <= spacing(zero)) then
  write(*,*) "Diagonal no inversible"
  stop
 else
  write(*,*)"Diagonal inversible"
 end if
end do

Jacobi Iteration

The Jacobi iteration is performed within a loop that runs up to Kx times:

do k = 1, kx
 !write(*,*)k
 do i = 1, N
  sum = 0d0
  do j = 1, N
   if (i == j) then
    !write(*,*)"cycle"
    cycle
   else
    sum = sum + A(i,j) * X0(j)
   end if
  end do
 end do

 ERROR = X0 - X
 call normavector2(ERROR, ek)
 if (ek < tolerancia) then
  write(*,*) "alcanzada tolerancia en k =", k
  exit
 end if

 X0 = X

end do

Inside the outer loop, the new value of X(i) is computed using the values from the previous iteration (X0). The error is calculated, and if it's less than the specified tolerance, the loop exits.