svdmodule Module Reference

Functions/Subroutines
pure real(dp) function, dimension(:, :), allocatable	house_holder_matrix (x)
pure subroutine, public	bidiagonal_decomposition (A, P, Qt)
pure real(dp) function, dimension(2, 2)	givens_rotation (a, b)
pure real(dp) function	compute_shift (A)
pure subroutine, public	bidiagonal_qr_decomposition (A, U, VT)
pure subroutine	handle_zero_diagonal (A, U, VT)
pure real(dp) function	superdiagonal_norm (A)
pure subroutine	find_nonzero_superdiagonal (A, p, q)
pure logical(lgp) function	has_zero_diagonal (A)
pure subroutine	make_matrix_square (A, Qt)
pure subroutine	clean_superdiagonal (A)
pure subroutine, public	svd (A, U, S, VT)
	Singular Value Decomposition. More...

Function/Subroutine Documentation

bidiagonal_decomposition()

pure subroutine, public svdmodule::bidiagonal_decomposition	(	real(dp), dimension(:, :), intent(inout)	A,
		real(dp), dimension(:, :), intent(out), allocatable	P,
		real(dp), dimension(:, :), intent(out), allocatable	Qt
	)

Definition at line 64 of file SingularValueDecomposition.f90.

    ! -- dummy
    real(DP), intent(inout), dimension(:, :) :: A
    real(DP), intent(out), dimension(:, :), allocatable :: P, Qt
    ! -- locals
    integer(I4B) :: M, N
    integer(I4B) :: I
    real(DP), allocatable, dimension(:, :) :: Qi, Pi
    real(DP), allocatable, dimension(:) :: h
 
    m = size(a, dim=1) ! Number of rows
    n = size(a, dim=2) ! Number of columns
 
    qt = eye(n)
    p = eye(m)
 
    do i = 1, min(m, n)
      ! columns
      h = a(i:m, i)
      pi = house_holder_matrix(h)
      a(i:m, :) = matmul(pi, a(i:m, :)) ! Apply householder transformation from left
      p(:, i:m) = matmul(p(:, i:m), pi)
 
      ! rows
      if (i < n) then
        h = a(i, i + 1:n)
        qi = transpose(house_holder_matrix(h))
        a(:, i + 1:n) = matmul(a(:, i + 1:n), qi) ! Apply householder transformation from right
        qt(i + 1:n, :) = matmul(qi, qt(i + 1:n, :))
      end if
 
    end do
 

Here is the call graph for this function:

Here is the caller graph for this function:

bidiagonal_qr_decomposition()

pure subroutine, public svdmodule::bidiagonal_qr_decomposition	(	real(dp), dimension(:, :), intent(inout)	A,
		real(dp), dimension(:, :), intent(inout)	U,
		real(dp), dimension(:, :), intent(inout)	VT
	)

Definition at line 183 of file SingularValueDecomposition.f90.

    ! -- dummy
    real(DP), intent(inout), dimension(:, :) :: A
    real(DP), intent(inout), dimension(:, :) :: U, Vt
    ! -- locals
    integer(I4B) :: m, n, I, J
    real(DP), dimension(2, 2) :: G
    real(DP) :: t11, t12, mu
 
    m = size(a, dim=1) ! Number of rows
    n = size(a, dim=2) ! Number of columns
 
    DO i = 1, n - 1
      j = i + 1
      if (i == 1) then
        ! For the first iteration use the implicit shift
        mu = compute_shift(a)
        ! Apply the shift to the full tri-diagonal matrix.
        t11 = a(1, 1)**2 - mu
        t12 = a(1, 1) * a(1, 2)
 
        g = givens_rotation(t11, t12)
      else
        g = givens_rotation(a(i - 1, i), a(i - 1, j))
      end if
 
      a(:, i:j) = matmul(a(:, i:j), transpose(g))
      vt(i:j, :) = matmul(g, vt(i:j, :))
 
      if (j > m) cycle
      g = givens_rotation(a(i, i), a(j, i))
      a(i:j, :) = matmul(g, a(i:j, :))
      u(:, i:j) = matmul(u(:, i:j), transpose(g))
    END DO
 

Here is the call graph for this function:

Here is the caller graph for this function:

clean_superdiagonal()

pure subroutine svdmodule::clean_superdiagonal ( real(dp), dimension(:, :), intent(inout)  A )

private

Definition at line 385 of file SingularValueDecomposition.f90.

    ! -- dummy
    real(DP), intent(inout), dimension(:, :) :: A
    ! -- locals
    integer(I4B) :: n, m, j
 
    m = size(a, dim=1) ! Number of rows
    n = size(a, dim=2) ! Number of columns
 
    do j = 1, min(n, m) - 1
      if (abs(a(j, j + 1)) <= dprec * (abs(a(j, j)) + abs(a(j + 1, j + 1)))) then
        a(j, j + 1) = 0.0_dp
      end if
    end do
 

Here is the caller graph for this function:

compute_shift()

pure real(dp) function svdmodule::compute_shift ( real(dp), dimension(:, :), intent(in)  A )

private

Definition at line 135 of file SingularValueDecomposition.f90.

    ! -- dummy
    real(DP), intent(in), dimension(:, :) :: A
    real(DP) :: mu
    ! -- locals
    integer(I4B) :: m, n, min_mn
    real(DP) T11, T12, T21, T22
    real(DP) dm, fmmin, fm, dn
    real(DP) :: mean, product, mean_product, mu1, mu2
 
    m = size(a, dim=1) ! Number of rows
    n = size(a, dim=2) ! Number of columns
 
    min_mn = min(m, n)
 
    dn = a(min_mn, min_mn)
    dm = a(min_mn - 1, min_mn - 1)
    fm = a(min_mn - 1, min_mn)
    if (min_mn > 2) then
      fmmin = a(min_mn - 2, min_mn - 1)
    else
      fmmin = 0.0_dp
    end if
 
    t11 = dm**2 + fmmin**2
    t12 = dm * fm
    t21 = t12
    t22 = fm**2 + dn**2
 
    mean = (t11 + t22) / 2.0_dp
    product = t11 * t22 - t12 * t21
    mean_product = mean**2 - product
    if (abs(mean_product) < dem6) then
      mean_product = 0.0_dp
    end if
    mu1 = mean - sqrt(mean_product)
    mu2 = mean + sqrt(mean_product)
    if (abs(t22 - mu1) < abs(t22 - mu2)) then
      mu = mu1
    else
      mu = mu2
    end if
 

Here is the caller graph for this function:

find_nonzero_superdiagonal()

pure subroutine svdmodule::find_nonzero_superdiagonal ( real(dp), dimension(:, :), intent(in)  A, integer(i4b), intent(out)  p, integer(i4b), intent(out)  q  )

private

Definition at line 296 of file SingularValueDecomposition.f90.

    ! -- dummy
    real(DP), intent(in), dimension(:, :) :: A
    integer(I4B), intent(out) :: p ! Index of the first nonzero superdiagonal (start of active block)
    integer(I4B), intent(out) :: q ! Index of the last nonzero superdiagonal (end of active block)
    ! -- locals
    integer(I4B) :: m, n, j, min_mn
 
    m = size(a, dim=1) ! Number of rows
    n = size(a, dim=2) ! Number of columns
 
    min_mn = min(m, n)
    p = 1
    q = min_mn
 
    do j = min_mn, 2, -1
      if (abs(a(j - 1, j)) > dprec) then
        q = j
        exit
      end if
    end do
 
    do j = q - 1, 2, -1
      if (abs(a(j - 1, j)) < dprec) then
        p = j
        exit
      end if
    end do

Here is the caller graph for this function:

givens_rotation()

pure real(dp) function, dimension(2, 2) svdmodule::givens_rotation ( real(dp), intent(in)  a, real(dp), intent(in)  b  )

private

Definition at line 102 of file SingularValueDecomposition.f90.

    ! -- dummy
    real(DP), intent(in) :: a, b
    real(DP), dimension(2, 2) :: G
    ! -- locals
    real(DP) :: c, s, h, d
 
    if (abs(b) < dprec) then
      g = eye(2)
      return
    end if
 
    h = hypot(a, b)
    d = 1.0 / h
    c = abs(a) * d
    s = sign(d, a) * b
 
    g(1, 1) = c
    g(1, 2) = s
    g(2, 1) = -s
    g(2, 2) = c
 

Here is the call graph for this function:

Here is the caller graph for this function:

handle_zero_diagonal()

pure subroutine svdmodule::handle_zero_diagonal ( real(dp), dimension(:, :), intent(inout)  A, real(dp), dimension(:, :), intent(inout)  U, real(dp), dimension(:, :), intent(inout)  VT  )

private

Definition at line 230 of file SingularValueDecomposition.f90.

    ! -- dummy
    real(DP), intent(inout), dimension(:, :) :: A
    real(DP), intent(inout), dimension(:, :) :: U, Vt
    ! -- locals
    integer(I4B) :: m, n, I
    real(DP), dimension(2, 2) :: G
    integer(I4B) :: zero_index
 
    m = size(a, dim=1) ! Number of rows
    n = size(a, dim=2) ! Number of columns
 
    do i = 1, min(m, n)
      if (abs(a(i, i)) < dprec) then
        zero_index = i
        exit
      end if
    end do
 
    if (zero_index == min(n, m)) then
      ! If the zero index is the last element of the diagonal then zero out the column
      do i = zero_index - 1, 1, -1
        g = givens_rotation(a(i, i), a(i, zero_index))
        a(:, [i, zero_index]) = matmul(a(:, [i, zero_index]), transpose(g))
        vt([i, zero_index], :) = matmul(g, vt([i, zero_index], :))
      end do
    else
      ! Else zero out the row
      do i = zero_index + 1, n
        g = givens_rotation(a(i, i), a(zero_index, i))
        a([zero_index, i], :) = matmul(transpose(g), a([zero_index, i], :))
        u(:, [zero_index, i]) = matmul(u(:, [zero_index, i]), g)
      end do
    end if

Here is the call graph for this function:

Here is the caller graph for this function:

has_zero_diagonal()

pure logical(lgp) function svdmodule::has_zero_diagonal ( real(dp), dimension(:, :), intent(in)  A )

private

Definition at line 329 of file SingularValueDecomposition.f90.

    ! -- dummy
    real(DP), intent(in) :: A(:, :)
    logical(LGP) :: has_zero
    ! -- locals
    integer(I4B) :: m, n, I
 
    m = size(a, dim=1) ! Number of rows
    n = size(a, dim=2) ! Number of columns
 
    has_zero = .false.
    do i = 1, min(m, n)
      if (abs(a(i, i)) < dprec) then
        has_zero = .true.
        exit
      end if
    end do
 

Here is the caller graph for this function:

house_holder_matrix()

pure real(dp) function, dimension(:, :), allocatable svdmodule::house_holder_matrix ( real(dp), dimension(:), intent(in)  x )

private

Definition at line 18 of file SingularValueDecomposition.f90.

    ! -- dummy
    real(DP), intent(in) :: x(:)
    real(DP), allocatable, dimension(:, :) :: Q
    ! -- locals
    real(DP) :: x_norm, y_norm
    real(DP), allocatable :: v(:)
    real(DP), allocatable :: w(:)
 
    x_norm = norm2(x)
    y_norm = norm2(x(2:))
 
    q = eye(size(x))
 
    if (dabs(y_norm) < dsame) then
      return
    end if
 
    v = x
    if (x(1) > 0) then
      ! Parlett (1971) suggested this modification to avoid cancellation
      v(1) = -(y_norm**2) / (x(1) + x_norm)
    else
      v(1) = v(1) - x_norm
    end if
 
    w = v / norm2(v)
 
    q = q - 2.0_dp * outer_product(w, w)
 

Here is the call graph for this function:

Here is the caller graph for this function:

make_matrix_square()

pure subroutine svdmodule::make_matrix_square ( real(dp), dimension(:, :), intent(inout)  A, real(dp), dimension(:, :), intent(inout), allocatable  Qt  )

private

Definition at line 355 of file SingularValueDecomposition.f90.

    ! -- dummy
    real(DP), intent(inout), dimension(:, :) :: A
    real(DP), intent(inout), dimension(:, :), allocatable :: Qt
    ! -- locals
    real(DP), dimension(2, 2) :: G
    integer(I4B) :: m, n, I
 
    m = size(a, dim=1) ! Number of rows
    n = size(a, dim=2) ! Number of columns
 
    do i = m, 1, -1
      g = givens_rotation(a(i, i), a(i, m + 1))
      a(:, [i, m + 1]) = matmul(a(:, [i, m + 1]), transpose(g))
      qt([i, m + 1], :) = matmul(g, qt([i, m + 1], :))
    end do
 

Here is the call graph for this function:

Here is the caller graph for this function:

superdiagonal_norm()

pure real(dp) function svdmodule::superdiagonal_norm ( real(dp), dimension(:, :), intent(in)  A )

private

Definition at line 269 of file SingularValueDecomposition.f90.

    ! -- dummy
    real(DP), intent(in) :: A(:, :)
    real(DP) :: norm
    ! -- locals
    integer(I4B) :: m, n, I
 
    m = size(a, dim=1) ! Number of rows
    n = size(a, dim=2) ! Number of columns
 
    norm = 0.0_dp
    do i = 1, min(m, n) - 1
      norm = max(norm, abs(a(i, i + 1)))
    end do
 

Here is the caller graph for this function:

svd()

pure subroutine, public svdmodule::svd	(	real(dp), dimension(:, :), intent(in)	A,
		real(dp), dimension(:, :), intent(out), allocatable	U,
		real(dp), dimension(:, :), intent(out), allocatable	S,
		real(dp), dimension(:, :), intent(out), allocatable	VT
	)

This method decomposes the matrix A into U, S and VT. It follows the algorithm as described by Golub and Reinsch.

The first step is to decompose the matrix A into a bidiagonal matrix. This is done using Householder transformations. Then second step is to decompose the bidiagonal matrix into U, S and VT by repetitively applying the QR algorithm. If there is a zero on the diagonal or superdiagonal the matrix can be split into two smaller matrices and the QR algorithm can be applied to the smaller matrices.

The matrix U is the eigenvectors of A*A^T The matrix VT is the eigenvectors of A^T*A The matrix S is the square root of the eigenvalues of A*A^T or A^T*A

Definition at line 420 of file SingularValueDecomposition.f90.

    ! -- dummy
    real(DP), intent(in), dimension(:, :) :: A
    real(DP), intent(out), dimension(:, :), allocatable :: U
    real(DP), intent(out), dimension(:, :), allocatable :: S
    real(DP), intent(out), dimension(:, :), allocatable :: VT
    ! -- locals
    integer(I4B) :: i ! Loop counter for QR iterations
    integer(I4B) :: max_itr ! Maximum number of QR iterations allowed for convergence
    real(DP) :: error ! Convergence criterion (superdiagonal norm)
 
    integer(I4B) :: m ! Number of rows in input matrix A
    integer(I4B) :: n ! Number of columns in input matrix A
 
    integer(I4B) :: r ! Start index of active submatrix (for QR step)
    integer(I4B) :: q ! End index of active submatrix (for QR step)
 
    max_itr = 100
 
    m = size(a, dim=1) ! Number of rows
    n = size(a, dim=2) ! Number of columns
    s = a
 
    call bidiagonal_decomposition(s, u, vt)
 
    if (n > m) then
      call make_matrix_square(s, vt)
    end if
 
    do i = 1, max_itr
      call find_nonzero_superdiagonal(s, r, q)
 
      ! find zero diagonal
      if (has_zero_diagonal(s(r:q, r:q))) then
        ! write(*,*) 'Iteration: ', i, ' handle zero diagonal element'
        call handle_zero_diagonal(s(r:q, r:q), u(:, r:q), vt(r:q, :))
      else
        call bidiagonal_qr_decomposition(s(r:q, r:q), u(:, r:q), vt(r:q, :))
        call clean_superdiagonal(s(r:q, r:q))
 
        error = superdiagonal_norm(s)
        ! write(*,*) 'Iteration: ', i, ' Error: ', error
        if (error < dprec) exit
      end if
    end do
 

Here is the call graph for this function:

Here is the caller graph for this function: