Time sequence and spectral analysis - Earth Dynamics: Deformations and Oscillations of the Rotating Earth - page 151

The subroutine SVD begins with a call to the subroutine BIDIAG to construct the
real, positive, upper bidiagonal matrix B . It then performs an implicit QR iteration,
through a series of Givens rotations, to reduce B to diagonal form containing the
singular values of C . At each step of the iteration, the upper diagonal element of
the leading row of the working matrix is tested for negligibility. If it is negligibly
small compared to the smallest of its nearest diagonal neighbours, the matrix is
deflated and a new leading row of the working matrix is established. Similarly, a
new trailing row of the working matrix is established as the first row below the
leading row with a negligible upper diagonal element compared to its smallest
diagonal neighbour. The QR iteration is then repeated on the smaller split matrix.
The final iteration is on a 2

2 matrix. At each stage the unitary matrices U and
V H are updated, and at the end of the iteration their Hermitian transposes are taken
to give the unitary matrices appearing in the SVD factored form (2.199).

×

SUBROUTINE SVD(C,U,S,VH,M,MMAX)

C
C The subroutine SVD calls the subroutine BIDIAG to convert the complex
C matrix C to real, positive, upper bidiagonal form. The diagonal elements
C are stored in the vector S and its upper diagonal is stored in
C the vector UD. An implicit QR iteration with Wilkinson shift is then
C performed on the bidiagonal matrix through a series of Givens rotations
C to chase a rogue element, RE, from the matrix. After each QR iteration,
C the upper diagonal element in the leading row is tested for negligibility.
C The matrix is deflated when this element is found negligible, to
C computational accuracy, compared to its smallest diagonal neighbour,
C establishing a new leading row K. Similarly, a new trailing row, N, is
C established as the first row below the leading row with an upper diagonal
C element that is negligible compared to its smallest diagonal neighbour.
C The final iteration is on a 2 x 2 matrix. The singular values are on
C the diagonal of the final matrix. The singular values and their
C associated eigenvectors are then sorted so that they are in descending
C order down the diagonal of the final matrix. The Hermitian transposes of
C the accumulated matrices U and VH are taken to convert them to
C the matrices appearing in the SVD factorization.
C

IMPLICIT DOUBLE PRECISION(A-H,O-Z)
DOUBLE COMPLEX C(MMAX,MMAX),U(MMAX,MMAX),VH(MMAX,MMAX),HOLD
DIMENSION S(MMAX),UD(MMAX)
C Set maximum number of iterations.
MAXI=100
CALL BIDIAG(U,C,VH,M,MMAX)
C Put diagonal of real, positive bidiagonal matrix in S and its upper
C diagonal in UD.
DO 10 I=1,M
IP1=I+1
S(I)=DREAL(C(I,I))
IF(IP1.GT.M) GO TO 10
UD(I)=DREAL(C(I,IP1))
10 CONTINUE
C Begin diagonalization of matrix B.
C Set initial values of iteration number, ITER, leading row, K, and
C trailing row, N, of working matrix.
ITER=0
K=1
N=M