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

DO 19 I=1,M
DO 19 J=NP1,M
VH(I,J)=VH(I,J)-2.D0*H1(I)*HU(J)
C(I,J)=C(I,J)-2.D0*H2(I)*HU(J)
19 CONTINUE
C Overwrite annihilated Nth row.
DO 20 J=1,M
C(N,J)=V(J)
20 CONTINUE
12 CONTINUE
C Convert upper diagonal matrix B to a real, positive matrix.
N=1
CONST=(1.D0,0.D0)
21 CONTINUE
NP1=N+1
C Find magnitudes of current row elements.
EMAG=DREAL(C(N,N)*DCONJG(C(N,N)))
EMAG=DSQRT(EMAG)
C Update unitary matrix U.
CONST=DCONJG(CONST)*DCONJG(C(N,N))/EMAG
C Scale current row.
DO 22 I=1,M
U(N,I)=CONST*U(N,I)
22 CONTINUE
C Replace diagonal element.
C(N,N)=DCMPLX(EMAG,0.D0)
IF(N.EQ.M) GO TO 23
EMAG=DREAL(C(N,NP1)*DCONJG(C(N,NP1)))
EMAG=DSQRT(EMAG)
C Update unitary matrix VH.
CONST=DCONJG(CONST)*DCONJG(C(N,NP1))/EMAG
C Scale next column.
DO 24 I=1,M
VH(I,NP1)=CONST*VH(I,NP1)
24 CONTINUE
C Replace above diagonal element.
C(N,NP1)=DCMPLX(EMAG,0.D0)
N=N+1
GO TO 21

23 CONTINUE

RETURN
END

5 complex matrix we are using as an
example into the following real, positive, upper bidiagonal matrix:

The subroutine BIDIAG converts the 5 ×

real part of matrix is,
0.63681D+03 0.12738D+04 0.00000D+00 0.00000D+00 0.00000D+00
0.00000D+00 0.45167D+00 0.20223D+00 0.00000D+00 0.00000D+00
0.00000D+00 0.00000D+00 0.10497D-01 0.48666D-02 0.00000D+00
0.00000D+00 0.00000D+00 0.00000D+00 0.51721D-02 0.72003D-02
0.00000D+00 0.00000D+00 0.00000D+00 0.00000D+00 0.31164D-02
imaginary part of matrix is,
0.00000D+00 0.00000D+00 0.00000D+00 0.00000D+00 0.00000D+00
0.00000D+00 0.00000D+00 0.00000D+00 0.00000D+00 0.00000D+00
0.00000D+00 0.00000D+00 0.00000D+00 0.00000D+00 0.00000D+00
0.00000D+00 0.00000D+00 0.00000D+00 0.00000D+00 0.00000D+00
0.00000D+00 0.00000D+00 0.00000D+00 0.00000D+00 0.00000D+00.