dgectrf.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       DGECTRF  -  Crout Factorization for A = L * U
C     ====================================================================
C       Description:
C       ------------------------------------------------------------------
C         Computes the LU decomposition of a general N x N matrix A
C         using the Crout algorithm. A is stored as a flat row-major
C         array. The lower triangular matrix L and the upper triangular
C         matrix U (with ones on the diagonal) are stored in L and U arrays.
C
C     ====================================================================
C       Arguments:
C       ------------------------------------------------------------------
C         N      : INTEGER            -> matrix size (N x N)
C         A(*)   : DOUBLE PRECISION   -> flat row-major matrix A
C         L(*)   : DOUBLE PRECISION   -> lower triangular matrix L
C         U(*)   : DOUBLE PRECISION   -> upper triangular matrix U
C         INFO   : INTEGER            -> return code:
C                                           0 = success
C                                          <0 = zero diagonal in L detected
C     ====================================================================
      SUBROUTINE dgectrf(N, A, L, U, 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(*), L(*), U(*)
 
C   L o c a l   V a r i a b l e s
C   ------------------------------------------------------------------
      INTEGER            :: I, J, K
      DOUBLE PRECISION   :: SUM
 
C   I n i t i a l i z e
C   ------------------------------------------------------------------
      info = 0
      DO i = 1, n * n
         l(i) = 0.0d0
         u(i) = 0.0d0
      END DO
 
C   M a i n   D e c o m p o s i t i o n   L o o p
C   ------------------------------------------------------------------
      DO j = 1, n
C        L o w e r   T r i a n g u l a r
C        --------------------------------------------------------------
         DO i = j, n
            sum = 0.0d0
            DO k = 1, j - 1
               sum = sum + l((i-1)*n + k) * u((k-1)*n + j)
            END DO
            l((i-1)*n + j) = a((i-1)*n + j) - sum
         END DO
 
C        U p p e r   T r i a n g u l a r
C        --------------------------------------------------------------
         DO i = j, n
            IF (i .EQ. j) THEN
               u((j-1)*n + i) = 1.0d0
            ELSE
               sum = 0.0d0
               DO k = 1, j - 1
                  sum = sum + l((j-1)*n + k) * u((k-1)*n + i)
               END DO
 
C              Check for zero diagonal to avoid division by zero
               IF (l((j-1)*n + j) .EQ. 0.0d0) THEN
                  info = -j
                  RETURN
               END IF
 
               u((j-1)*n + i) = (a((j-1)*n + i) - sum) / l((j-1)*n + j)
            END IF
         END DO
      END DO
 
C   S u c c e s s f u l   e x i t
C   ------------------------------------------------------------------
      RETURN
      END