Geology Reference
In-Depth Information
The solution of the full equations for
m
=
n
+
1 is then derived from the solution
for
m
=
n
through the recurrence relations,
γ
n
=
f
n
,0
r
n
+
1
+
f
n
,1
r
n
+···+
f
n
,
n
r
1
,
g
n
+
1
−
γ
n
α
∗
n
+
1
,
q
n
=
q
n
a
∗
n
+
1,
n
+
1
,
f
n
+
1,0
=
f
n
,0
+
q
n
a
∗
n
+
1,
n
,
f
n
+
1,1
=
f
n
,1
+
.
.
.
(2.133)
q
n
a
∗
n
+
1,1
,
f
n
+
1,
n
=
f
n
,
n
+
q
n
a
∗
n
+
1,0
.
By making use of the Toeplitz form of the coe
f
n
+
1,
n
+
1
=
cient matrix, the Levinson
ort from the usual
N
3
, for linear algebraic
algorithm reduces the computational e
ff
N
, to only
N
2
.
The subroutine LEVSOL implements the Levinson algorithm through the suc-
cessive recursions (2.123) and (2.133).
systems of size
N
×
SUBROUTINE LEVSOL(LR,R,G,F,A)
C
C LEVSOL solves a linear system with a Toeplitz coefficient matrix by
C the Levinson recursive method. LR is the dimension of the system and R is
C the LR-length vector of autocorrelations for zero and positive lags,
C or in the general case, the first column vector of the coefficient matrix.
C G is the right-hand side vector, F is the solution vector and A is
C a vector representing the solution of an auxiliary prediction error
C system used in the recursion. Since A is overwritten in the recursion,
C its previous values are stored in the vector AO.
C
IMPLICIT DOUBLE COMPLEX(A-H,O-Z)
DIMENSION R(LR),G(LR),F(LR),A(LR),AO(LR)
C Solve 1x1 system.
F(1)=G(1)/R(1)
IF(LR.EQ.1) RETURN
C Set first element of the auxiliary vector to unity.
A(1)=(1.D0,0.D0)
C Initialize alpha, beta and gamma.
ALPHA=R(1)
BETA=R(2)
GAMMA=F(1)*R(2)
C Begin loop over remaining equation systems.
DO 10 L=2,LR
C Find k for L-1.
AK=-BETA/DCONJG(ALPHA)
C Find last element of new auxiliary vector.
LM1=L-1
A(L)=AK
C 2x2 system is a special case.
IF(L.EQ.2) GO TO 20
C Store old auxiliary vector elements.
DO 40 I=2,LM1
Search WWH ::
Custom Search