Geology Reference
In-Depth Information
DO 15 J=1,M
HOLD=VH(J,IP1)
VH(J,IP1)=CO*HOLD-SI*VH(J,IP2)
VH(J,IP2)=SI*HOLD+CO*VH(J,IP2)
15 CONTINUE
C Post-multiply matrix B by rotation.
H1=CMAG
H2=CO*S(IP2)
H3=CO*S(IP1)-SI*UD(IP1)
H4=CO*UD(IP1)+SI*S(IP1)
UD(I)=H1
S(IP1)=H3
UD(IP1)=H4
RE=-SI*S(IP2)
S(IP2)=H2
13 CONTINUE
C Find leading row of working matrix.
J=K
MM1=M-1
DO 16 I=K,MM1
IP1=I+1
C Find smallest neighbouring diagonal element.
SMIN=DMIN1(DABS(S(I)),DABS(S(IP1)))
C Compare to see if upper diagonal element is non-negligible.
SIZE=SMIN+DABS(UD(I))
IF(SIZE.NE.SMIN) GO TO 17
J=J+1
16 CONTINUE
17 CONTINUE
C Set new leading row number.
K=J
C Test to see if diagonalization is complete.
IF(K.EQ.M) GO TO 20
C Find trailing row of working matrix.
J=K+1
KKP1=K+1
DO 18 I=KKP1,MM1
IP1=I+1
C Find smallest neighbouring diagonal element.
SMIN=DMIN1(DABS(S(I)),DABS(S(IP1)))
C Compare to see if upper diagonal element is negligible.
SIZE=SMIN+DABS(UD(I))
IF(SIZE.EQ.SMIN) GO TO 19
J=J+1
18 CONTINUE
19 CONTINUE
C Set new trailing row number.
N=J
C Test for failure to converge.
NNM1=N-1
IF(KKP1.NE.KP1.OR.NNM1.NE.NM1) ITER=0
IF(ITER.EQ.MAXI) PAUSE 'Maximum iterations reached.'
GO TO 11
20 CONTINUE
C Take Hermitian transpose of U and VH to get their final forms.
DO 21 I=1,M
DO 21 J=1,I
HOLD=U(I,J)
U(I,J)=DCONJG(U(J,I))
U(J,I)=DCONJG(HOLD)
21
CONTINUE
Search WWH ::




Custom Search