dgttsv.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       DGTTSV  -  Thomas Algorithm for (tridiagonal system)
C     ====================================================================
C       Description:
C       ------------------------------------------------------------------
C         Direct solver for tridiagonal linear systems A * X = B using
C         the Thomas algorithm.  A is supplied as a FULL N x N matrix in
C         flat row-major storage and is overwritten during computation.
C
C         On output: X contains the solution vector.
C
C         No pivoting is performed. Zero diagonal entries produce failure.
C     ====================================================================
C       Arguments:
C       ------------------------------------------------------------------
C         N      : INTEGER            -> matrix size (N x N)
C         A(*)   : DOUBLE PRECISION   -> flat row-major matrix A (modified)
C         B(N)   : DOUBLE PRECISION   -> right-hand side (modified)
C         X(N)   : DOUBLE PRECISION   -> solution vector (output)
C         INFO   : INTEGER            -> return code:
C                                           0 = success
C                                          <0 = zero diagonal detected
C     ====================================================================
      SUBROUTINE dgttsv(N, A, B, X, INFO)
 
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, INFO
      DOUBLE PRECISION   :: A(*), B(N), X(N)
 
C   L o c a l   V a r i a b l e s
C   ------------------------------------------------------------------
      INTEGER            :: I, INDEX_D, INDEX_L, INDEX_U
      DOUBLE PRECISION   :: FACT
 
C   I n i t i a l   S t a t u s
C   ------------------------------------------------------------------
      info = 0
 
C   F o r w a r d   E l i m i n a t i o n
C   ------------------------------------------------------------------
      DO i = 2, n
 
C        Compute indices in flat row-major layout
         index_l = (i-1) * n + (i-1)
         index_d = (i-2) * n + (i-1)
         index_u = (i-2) * n + i
 
C        Check previous diagonal
         IF (a(index_d) .EQ. 0.0d0) THEN
            info = -(i-1)
            RETURN
         END IF
 
         fact = a(index_l) / a(index_d)
 
C        Update diagonal A(I,I)
         index_d = (i-1) * n + i
         a(index_d) = a(index_d) - fact * a(index_u)
 
C        Update RHS B(I)
         b(i) = b(i) - fact * b(i-1)
 
C        Check new diagonal
         IF (a(index_d) .EQ. 0.0d0) THEN
            info = -i
            RETURN
         END IF
      END DO
 
C   B a c k w a r d   S u b s t i t u t i o n
C   ------------------------------------------------------------------
      index_d = (n-1) * n + n
      x(n) = b(n) / a(index_d)
 
      DO i = n-1, 1, -1
         index_d = (i-1) * n + i
         index_u = (i-1) * n + (i+1)
         x(i) = (b(i) - a(index_u) * x(i+1)) / a(index_d)
      END DO
 
C   S u c c e s s
C   ------------------------------------------------------------------
      RETURN
      END