dgegssv.f

Go to the documentation of this file.

C     ====================================================================
C                            N U L A P A C K
C                            U U L A P A C K
C                            L L L A P A C K
C                            A A A A P A C K
C                            P P P P P A C K
C                            A A A A A A C K
C                            C C C C C C C K
C                            K K K K K K K K
C
C     This file is part of NULAPACK - NUmerical Linear Algebra PACKage
C
C     Copyright (C) 2025  Saud Zahir
C
C     NULAPACK is free software: you can redistribute it and/or modify
C     it under the terms of the GNU General Public License as published by
C     the Free Software Foundation, either version 3 of the License, or
C     (at your option) any later version.
C
C     NULAPACK is distributed in the hope that it will be useful,
C     but WITHOUT ANY WARRANTY; without even the implied warranty of
C     MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
C     GNU General Public License for more details.
C
C     You should have received a copy of the GNU General Public License
C     along with NULAPACK.  If not, see <https://www.gnu.org/licenses/>.
C
C     ====================================================================
C       DGEGSSV  -   Gauss-Seidel Solver for A * X = B
C     ====================================================================
C       Description:
C       ------------------------------------------------------------------
C         Iterative Gauss-Seidel solver for solving linear systems of
C         equations A * X = B, where A is a square N x N matrix in
C         row-major flat array format. Double precision version.
C
C         On input:  X contains initial guess
C         On output: X contains solution
C
C         Convergence is based on maximum absolute difference per iteration.
C     ====================================================================
C       Arguments:
C       ------------------------------------------------------------------
C         N         : INTEGER            -> size of the matrix (N x N)
C         A(*)      : DOUBLE PRECISION   -> flat array, row-major matrix A
C         B(N)      : DOUBLE PRECISION   -> right-hand side vector
C         X(N)      : DOUBLE PRECISION   -> input: initial guess, output: solution
C         MAX_ITER  : INTEGER            -> max number of iterations
C         TOL       : DOUBLE PRECISION   -> convergence tolerance
C         OMEGA     : DOUBLE PRECISION   -> relaxation coefficient
C         STATUS    : INTEGER            -> return code:
C                                                 0 = success
C                                                >0 = did not converge
C                                                <0 = illegal or zero diagonal
C     ====================================================================
      SUBROUTINE dgegssv(N, A, B, X, MAX_ITER, TOL, OMEGA, STATUS)
 
C   I m p l i c i t   T y p e s
C   ------------------------------------------------------------------
      IMPLICIT NONE
 
C   D u m m y   A r g u m e n t s
C   ------------------------------------------------------------------
      INTEGER            :: N, MAX_ITER, STATUS
      DOUBLE PRECISION   :: A(*), B(N), X(N), TOL, OMEGA
 
C   L o c a l   V a r i a b l e s
C   ------------------------------------------------------------------
      INTEGER            :: I, J, K, INDEX
      DOUBLE PRECISION   :: X_NEW(N)
      DOUBLE PRECISION   :: S1, S2, DIFF, MAX_DIFF
 
C   I n i t i a l   S t a t u s
C   ------------------------------------------------------------------
      status = 1   ! Default: did not converge
 
C   M a i n   I t e r a t i o n   L o o p
C   ------------------------------------------------------------------
      DO k = 1, max_iter
         DO i = 1, n
            s1 = 0.0d0
            s2 = 0.0d0
 
C           Compute sum: S1 = sum_{j=1}^{i-1} A(i,j) * X_NEW(j)
            DO j = 1, i - 1
               index = (i - 1) * n + j
               s1 = s1 + a(index) * x_new(j)
            END DO
 
C           Compute sum: S2 = sum_{j=i+1}^{N} A(i,j) * X(j)
            DO j = i + 1, n
               index = (i - 1) * n + j
               s2 = s2 + a(index) * x(j)
            END DO
 
C           Check diagonal element A(i,i)
            index = (i - 1) * n + i
            IF (a(index) .EQ. 0.0d0) THEN
               status = -i
               RETURN
            END IF
 
C           Update X_NEW(i) with relaxation
            x_new(i) = (b(i) - s1 - s2) / a(index)
            x_new(i) = x(i) + omega * (x_new(i) - x(i))
         END DO
 
C        C o n v e r g e n c e   C h e c k
C        ----------------------------------------------------------------
         max_diff = 0.0d0
         DO i = 1, n
            diff = abs(x_new(i) - x(i))
            IF (diff .GT. max_diff) max_diff = diff
            x(i) = x_new(i)
         END DO
 
         IF (max_diff .LT. tol) THEN
            status = 0  ! Success
            RETURN
         END IF
      END DO
 
C   N o n - C o n v e r g e n c e   E x i t
C   ------------------------------------------------------------------
      RETURN
      END