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

C
DOUBLE COMPLEX C(MMAX,MMAX),V(MMAX),HU(MMAX),
1 U(MMAX,MMAX),VH(MMAX,MMAX),H0(MMAX),H1(MMAX),H2(MMAX),CONST
DOUBLE PRECISION EMAG
C Initialize unitary matrices U and VH.
DO 10 I=1,M
DO 11 J=1,M
U(I,J)=(0.D0,0.D0)
VH(I,J)=(0.D0,0.D0)
11 CONTINUE
U(I,I)=(1.D0,0.D0)
VH(I,I)=(1.D0,0.D0)
10 CONTINUE
C Apply Householder annihilation.
MM1=M-1
DO 12 N=1,MM1
C Put Nth column in V.
DO 13 I=1,M
V(I)=C(I,N)
13 CONTINUE
C Annihilate Nth column.
CALL HHOLDER(V,HU,N,M,MMAX)
C Multiply Hermitian transpose of Householder vector HU
C into unitary matrix U and into transformed matrix B.
DO 14 I=1,M
H1(I)=(0.D0,0.D0)
H2(I)=(0.D0,0.D0)
DO 14 J=N,M
H0(J)=DCONJG(HU(J))
H1(I)=H1(I)+H0(J)*U(J,I)
H2(I)=H2(I)+H0(J)*C(J,I)
14 CONTINUE
C Update unitary matrix U and transformed matrix B.
DO 15 I=N,M
DO 15 J=1,M
U(I,J)=U(I,J)-2.D0*HU(I)*H1(J)
C(I,J)=C(I,J)-2.D0*HU(I)*H2(J)
15 CONTINUE
C Overwrite annihilated Nth column.
DO 16 I=1,M
C(I,N)=V(I)
16 CONTINUE
IF(N.EQ.MM1) GO TO 12
C Put Nth row in V.
DO 17 J=1,M
V(J)=C(N,J)
17 CONTINUE
NP1=N+1
C Annihilate Nth row.
CALL HHOLDER(V,HU,NP1,M,MMAX)
C Multiply unitary matrix VH and transformed matrix C
C into Householder vector HU.
DO 18 I=1,M
H1(I)=(0.D0,0.D0)
H2(I)=(0.D0,0.D0)
DO 18 J=NP1,M
H0(J)=DCONJG(HU(J))
H1(I)=H1(I)+VH(I,J)*H0(J)
H2(I)=H2(I)+C(I,J)*H0(J)
18 CONTINUE
C Update unitary matrix VH and transformed matrix B.