Earth’s rotation: observations and theory - Earth Dynamics: Deformations and Oscillations of the Rotating Earth - page 283

C Perform sums.
A(J,1)=A(J,1)+DCMPLX(DCOS(ARG1),DSIN(ARG1))/SIGS(I)
B(J)=B(J)+G(I)*DCMPLX(DCOS(ARG2),DSIN(ARG2))/SIGS(I)
20 CONTINUE
B(J)=T*B(J)
19 CONTINUE
C Complete construction of conditional equations matrix A.
DO 21 J=2,M
DO 23 I=1,M
K=I-J+1
L=2-K
IF(K.LE.0)GO TO 22
A(I,J)=A(K,1)
GO TO 23
22 A(I,J)=DCONJG(A(L,1))
23 CONTINUE
21 CONTINUE
C Do singular value decomposition of conditional equations matrix A=U*S*VH.
CALL SVD(A,U,S,VH,M,4000)
C Multiply the Hermitian transpose of U into right-hand side vector B.
DO 24 I=1,M
C Accumulate product vector in SV.
SV(I)=(0.D0,0.D0)
DO 25 J=1,M
SV(I)=SV(I)+DCONJG(U(J,I))*B(J)
25 CONTINUE
24 CONTINUE
C Replace right-hand side vector B by product vector SV.
DO 26 I=1,M
B(I)=SV(I)
26 CONTINUE
C Write out singular values.
WRITE(4,27)
27 FORMAT(1X,'Singular values are:')
WRITE(4,28)(S(I),I=1,M)
28 FORMAT(1X,D23.15)
C Begin solution of reduced conditional equations.
C Set number of singular values to be eliminated.
WRITE(6,29)
29 FORMAT(1X,'Enter number of singular values to be eliminated.')
READ(5,*)NSVE
C Find number of remaining singular values.
30 NMAX=M-NSVE
C Multiply V into the inverse of truncated S.
DO 31 I=1,M
DO 32 J=1,M
IF(J.GT.NMAX)GO TO 33
US(I,J)=DCONJG(VH(J,I))/S(J)
GO TO 32
33 US(I,J)=(0.D0,0.D0)
32 CONTINUE
31 CONTINUE
C Multiply result into modified right-hand side to find solution vector SV
C and solution vector length SVL.
SVL=0.D0
DO 34 I=1,M
SV(I)=(0.D0,0.D0)
DO 35 J=1,M
SV(I)=SV(I)+US(I,J)*B(J)

35 CONTINUE

SVL=SVL+SV(I)*DCONJG(SV(I))