Geology Reference
In-Depth Information
C Continue reduction of working matrix.
11 CONTINUE
C Advance iteration counter.
ITER=ITER+1
KP1=K+1
NM1=N-1
C Find Wilkinson shift.
P=S(N)
Q=S(NM1)
R=UD(NM1)
H1=DSQRT((P+Q)*(P+Q)+R*R)
H2=DSQRT((P-Q)*(P-Q)+R*R)
SIG1=(H1-H2)/2.D0
C Find first Givens rotation.
T11=(S(K)-SIG1)*(S(K)+SIG1)
T12=S(K)*UD(K)
CMAG=DSQRT(T11*T11+T12*T12)
CO=T11/CMAG
SI=-T12/CMAG
C Update unitary matrix VH.
DO 12 I=1,M
HOLD=VH(I,K)
VH(I,K)=CO*HOLD-SI*VH(I,KP1)
VH(I,KP1)=SI*HOLD+CO*VH(I,KP1)
12 CONTINUE
C Post-multiply bidiagonal matrix B by rotation.
H1=CO*S(K)-SI*UD(K)
H2=CO*UD(K)+SI*S(K)
S(K)=H1
UD(K)=H2
C Find rogue element RE.
RE=-SI*S(KP1)
S(KP1)=CO*S(KP1)
C Find remaining Givens rotations.
DO 13 I=K,NM1
IP1=I+1
IP2=I+2
CMAG=DSQRT(S(I)*S(I)+RE*RE)
CO=S(I)/CMAG
SI=-RE/CMAG
C Update unitary matrix U.
DO 14 J=1,M
HOLD=U(I,J)
U(I,J)=CO*HOLD-SI*U(IP1,J)
U(IP1,J)=SI*HOLD+CO*U(IP1,J)
14 CONTINUE
C Pre-multiply matrix B by rotation.
H1=CMAG
H2=CO*UD(I)-SI*S(IP1)
H4=CO*S(IP1)+SI*UD(I)
S(I)=H1
UD(I)=H2
S(IP1)=H4
IF(IP2.GT.N) GO TO 13
RE=-SI*UD(IP1)
UD(IP1)=CO*UD(IP1)
C Find next Givens rotation.
CMAG=DSQRT(UD(I)*UD(I)+RE*RE)
CO=UD(I)/CMAG
SI=-RE/CMAG
C Update unitary matrix VH.
Search WWH ::




Custom Search