spoctrf.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       SPOCTRF  -  Cholesky Factorization for A = L * L^T
C     ====================================================================
C       Description:
C       ------------------------------------------------------------------
C         Computes the Cholesky factorization of a real symmetric
C         positive-definite matrix A stored in a flat row-major array.
C         The lower-triangular matrix L overwrites the lower triangle of A.
C
C         On output: A contains L in the lower triangle, upper triangle is
C         not referenced.
C
C     ====================================================================
C       Arguments:
C       ------------------------------------------------------------------
C         N    : INTEGER       -> order of the matrix
C         A(*) : REAL          -> flat row-major matrix, size (LDA*N)
C         L(*) : REAL          -> flat row-major matrix, size (LDA*N)
C         LDA  : INTEGER       -> leading dimension of A (usually N)
C         INFO : INTEGER       -> return code:
C                                0 = success
C                               <0 = illegal argument
C                               >0 = matrix not positive definite
C     ====================================================================
      SUBROUTINE spoctrf(N, A, L, LDA, 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, LDA, INFO
      REAL          :: A(*), L(*)
 
C   L o c a l   v a r i a b l e s
C   ------------------------------------------------------------------
      INTEGER       :: I, J, K, INDEX
      REAL          :: SUM
 
C   I n i t i a l i z e
C   ------------------------------------------------------------------
      info = 0
 
C   M a i n   l o o p   o v e r   r o w s
C   ------------------------------------------------------------------
      DO i = 1, n
 
C       Compute diagonal element L(I,I)
         sum = 0.0
         DO k = 1, i-1
            index = (i-1)*lda + k
            sum = sum + l(index)*l(index)
         END DO
 
         index = (i-1)*lda + i
         IF (a(index) - sum .LE. 0.0) THEN
            info = i
            RETURN
         END IF
         l(index) = sqrt(a(index) - sum)
 
C       Compute off-diagonal elements L(J,I), J = I+1:N
         DO j = i+1, n
            sum = 0.0
            DO k = 1, i-1
               sum = sum + l((j-1)*lda + k) * l((i-1)*lda + k)
            END DO
            index = (j-1)*lda + i
            l(index) = (a(index) - sum)/l((i-1)*lda + i)
         END DO
 
      END DO
 
C   S u c c e s s f u l   e x i t
C   ------------------------------------------------------------------
      RETURN
      END