schinx.f90

Go to the documentation of this file.

! Copyright (C) 2008, 2009  VZLU Prague, a.s., Czech Republic, Jaroslav Hajek <highegg@gmail.com>
! Copyright (C) 2026 Martin Köhler <koehlerm(AT)mpi-magdeburg.mpg.de>
!
! This file is part of qrupdate-ng.
!
! qrupdate is free software; you can redistribute it and/or modify
! it under the terms of the GNU General Public License as published by
! the Free Software Foundation; either version 3 of the License, or
! (at your option) any later version.
!
! This program is distributed in the hope that it will be useful,
! but WITHOUT ANY WARRANTY; without even the implied warranty of
! MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
! GNU General Public License for more details.
!
! You should have received a copy of the GNU General Public License
! along with this software; see the file COPYING.  If not, see
! <http://www.gnu.org/licenses/>.
!
!> \brief Updates a Cholesky factorization after inserting a row and column.
!>
!> \par Definition:
! =============
!> \verbatim
!>       subroutine schinx(n,R,ldr,j,u,w,info)
!>
!>       .. Scalar Arguments ..
!>       integer            n, j, ldr, info
!>       ..
!>       .. Array Arguments ..
!>       real               R(ldr,*), u(*), w(*)
!>       ..
!> \endverbatim
!>
!> \par Purpose:
! =============
!> \verbatim
!>
!> SCHINX updates the Cholesky factorization of a symmetric
!> positive definite matrix A after inserting a row and column.
!> Given an upper triangular matrix R that is a Cholesky factor of
!> A, i.e., A = R.'*R, where R.' denotes the transpose of R, this
!> SCHINX updates R -> R1 so that R1.'*R1 = A1, where
!> A1(jj,jj) = A, A1(j,:) = u.', A1(:,j) = u, and
!> jj = [1:j-1, j+1:n+1].
!>
!> On exit, u is destroyed and R is extended by one row and column.
!> The insertion is performed by first solving R.'*u = v, checking
!> positive definiteness, and then retriangularizing.
!> \endverbatim
!>
!> \param[in] n
!> \verbatim
!>          n is INTEGER
!>          The order of matrix R.  n >= 0.
!> \endverbatim
!>
!> \param[in,out] R
!> \verbatim
!>          R is REAL array, dimension (ldr,n+1)
!>          On entry, the upper triangular matrix R, the Cholesky
!>          factor of A.  On exit, the updated upper triangular
!>          matrix R1, the Cholesky factor of A1.
!> \endverbatim
!>
!> \param[in] ldr
!> \verbatim
!>          ldr is INTEGER
!>          The leading dimension of the array R.  ldr >= n+1.
!> \endverbatim
!>
!> \param[in] j
!> \verbatim
!>          j is INTEGER
!>          The position of the inserted row and column.
!>          1 <= j <= n+1.
!> \endverbatim
!>
!> \param[in,out] u
!> \verbatim
!>          u is REAL array, dimension (n+1)
!>          On entry, the vector defining the inserted row/column.
!>          On exit, u is destroyed.
!> \endverbatim
!>
!> \param[out] w
!> \verbatim
!>          w is REAL array, dimension (n+1)
!>          Workspace vector used during the retriangularization.
!> \endverbatim
!>
!> \param[out] info
!> \verbatim
!>          info is INTEGER
!>          = 0:  successful exit
!>          = 1:  the update would violate positive-definiteness
!>          = 2:  R is singular
!> \endverbatim
!>
!> \ingroup choldecomp
subroutine schinx(n,R,ldr,j,u,w,info)
    use iso_fortran_env
    use qrupdate_error
    integer, intent(in) :: n, j, ldr
    real(real32), intent(inout) :: R(ldr,*)
    real(real32), intent(inout) :: u(*)
    real(real32), intent(out) :: w(*)
    integer, intent(out) :: info
    external scopy,snrm2,strsv,sqrtv1,sqrqh
    real(real32) snrm2,rho,t
    integer i
    ! check arguments
    info = 0
    if (n < 0) then
        info = -1
    else if (j < 1 .or. j > n+1) then
        info = -4
    end if
    if (info /= 0) then
        call qrupdate_xerror('SCHINX',info)
        return
    end if
    ! shift vector.
    t = u(j)
    do i = j,n
        u(i) = u(i+1)
    end do
    ! check for singularity of R.
    do i = 1,n
        if (r(i,i) == 0e0) goto 20
    end do
    ! form R' \ u
    call strsv('U','T','N',n,r,ldr,u,1)
    rho = snrm2(n,u,1)
    ! check positive definiteness.
    rho = t - rho**2
    if (rho <= 0e0) goto 10
    ! shift columns
    do i = n,j,-1
        call scopy(i,r(1,i),1,r(1,i+1),1)
        r(i+1,i+1) = 0e0
    end do
    call scopy(n,u,1,r(1,j),1)
    r(n+1,j) = sqrt(rho)
    ! retriangularize
    if (j < n+1) then
        ! eliminate the introduced spike.
        call sqrtv1(n+2-j,r(j,j),w)
        ! apply rotations to R
        call sqrqh(n+2-j,n+1-j,r(j,j+1),ldr,w,r(j+1,j))
        ! zero spike.
        do i = j+1,n+1
            r(i,j) = 0e0
        end do
    end if
    ! normal return.
    return
    ! error returns.
    10 info = 1
    return
    20 info = 2
    return
end subroutine