SUBROUTINE tqli(d,e,z) USE nrtype; USE nrutil, ONLY : assert_eq, nrerror USE nr, ONLY : pythag IMPLICIT NONE REAL(sp), DIMENSION(:), INTENT(inout) :: d, e REAL(sp), DIMENSION(:,:), OPTIONAL, INTENT(inout) :: z INTEGER(i4b) :: i, iter, l, m, n, ndum REAL(sp) :: b, c, dd, f, g, p, r, s REAL(sp), DIMENSION(SIZE(e)) :: ff n = assert_eq(SIZE(d), SIZE(e), 'tqli: n') IF(PRESENT(z)) ndum=assert_eq(n, SIZE(z, 1), SIZE(z, 2), 'tqli: ndum') e(:)=EOSHIFT(e(:),1) DO l=1, n iter=0 iterate: DO DO m=l, n-1 dd=ABS(d(m)) + ABS(d(m+1)) IF(ABS(e(m))+dd == dd) EXIT END DO IF (m==l) EXIT iterate IF(iter == 30) CALL nrerror('too many iterations in tqli') iter = iter+1 g=(d(l+1)-d(l))/(2.0_sp*e(l)) r = pythag(g, 1.0_sp) g = d(m) - d(l) + e(l) / (g + SIGN(r, g)) s = 1.0 c = 1.0 p = 0.0 DO i = m-1, l, -1 f = s*e(i) b = c*e(i) r = pythag(f, g) e(i+1) = r IF(r == 0.0) THEN d(i+1) = d(i+1) - p e(m) = 0.0 CYCLE iterate END IF s = f/r c = g/r g = d(i + 1) - p r = (d(i) - g) * s + 2.0_sp * c * b p = s * r d(i + 1) = g + p g = c * r - b IF(PRESENT(z)) THEN ff(1:n) = z(1:n,i+1) z(1:n, i+1) = s*z(1:n,i) + c*ff(1:n) z(1:n,i) = c*z(1:n, i) - s*ff(1:n) END IF END DO d(l) = d(l) - p e(l) = g e(m) = 0.0 END DO iterate END DO END SUBROUTINE tqli

This subroutine takes two arrays of length *n* on input: `d`

,
which contains diagonal terms of a tridiagonal symmetric
matrix
**A**, and `e`

, which contains
sub- (and super-) diagonal terms of
**A**. The third parameter, `z`

, is
optional. If it is provided the eigenvectors that correspond to
the diagonalising transformation
**Q** (and thus the
transformation itself) are written on it. Its
dimensions must therefore be .
On output the subroutine
returns diagonal terms on `d`

, whereas vector `e`

is
destroyed. Because there are only *n*-1 sub-diagonal terms,
`e(1)`

is ignored and the valid terms are assumed to
live on `e(2)`

through `e(n)`

.

The first two statements check if arrays `d`

and `e`

are
properly dimensioned. If they are the dimension is written on
`n`

, if not, an error message is written on standard output,
which contains a string `'tqli: n'`

and the subroutine aborts.
Then, if the third argument `z`

is present, we check if that
matrix is correctly dimensioned, i.e.,
and write the
dimension on `ndum`

, which is a dummy argument that's not used
by the subroutine any more, because the real *n* already lives in
`n`

:

n = assert_eq(SIZE(d), SIZE(e), 'tqli: n') IF(PRESENT(z)) ndum=assert_eq(n, SIZE(z, 1), SIZE(z, 2), 'tqli: ndum')

The next step is to shift array `e`

to the left by one place
for convenience:

e(:)=EOSHIFT(e(:),1)so that the first valid element of it goes into

`e(1)`

and the last goes into `e(n-1)`

, i.e., the matrix
Now the real fun begins. The main body of the subroutine is a large
`DO`

loop, with another `DO`

loop inside and a number of
conditionals:

DO l=1, n iter=0 iterate: DO DO m=l, n-1 dd=ABS(d(m)) + ABS(d(m+1)) IF(ABS(e(m))+dd == dd) EXIT END DO IF (m==l) EXIT iterate IF(iter == 30) CALL nrerror('too many iterations in tqli') iter = iter+1 blah... blah... blah... DO i = m-1, l, -1 blah... blah... blah... IF(r == 0.0) THEN d(i+1) = d(i+1) - p e(m) = 0.0 CYCLE iterate END IF blah... blah... blah... END DO blah... blah... blah... d(l) = d(l) - p e(l) = g e(m) = 0.0 END DO iterate END DOThere are several ideas here. First we move downwards, i.e., we start with

As you remember the
**QL** algorithm generates
eigenvalues in the order of diminishing absolute value. So once the first
such value pops up in the upper left corner, that part of the job
is finished, and the matrix can be deflated.

Within each such deflated submatrix, as *l* moves towards *n*, we try to
*split* the submatrix of
**A** on
a very small subdiagonal element *e*_{m}, so small that

|*d*_{m}| + |*d*_{m+1}| + |*e*_{m}| = |*d*_{m}| + |*d*_{m+1}|

within the accuracy of floating point operations. Assuming that we have found such an element, we are eventually going to set it to zero formally (at the end of the

`iterate`

loop),
while performing some mathematics on the elements between
Now, recall the shifting strategy:

OK, so we have this little submatrix:Choosek_{s}equal to the eigenvalue of the leadingsubmatrix that is closer tod_{l}.

and how do we find its eigenvalue? Well, we simply hit it with the Jacobi rotation:

and this

In the code this is implemented as:

g=(d(l+1)-d(l))/(2.0_sp*e(l)) r = pythag(g, 1.0_sp) g = d(m) - d(l) + e(l) / (g + SIGN(r, g))where

`g`

is initially ,
`r`

is
,
the shift
`d(l) - e(l) / (g + SIGN(r, g))`

,
and the second instance of `g`

becomes
Now let us have a look at the lower right corner of our submatrix
of a deflated original
**A** that results from a successful
splitting:

We choose our
and
to be

It is easy to see that their squares add up to 1, so they are good and .

But what if
? This
can happen only if
*e*_{m-1} = 0 and
*d*_{m} - *k*_{s} = 0, or so close
to zero that the machine can't make a difference. A condition like that
should not really
happen on the first sweep, because we have already located
such *m* that
,
and that *m* wasn't *m*-1!

This part of the computation is implemented thusly:

f = s*e(i) b = c*e(i) r = pythag(f, g) e(i+1) = r IF(r == 0.0) THEN d(i+1) = d(i+1) - p e(m) = 0.0 CYCLE iterate END IF s = f/r c = g/rwhere

`s`

and `c`

stand for
and
accordingly.
What are
and
chosen thusly going
to do to *e*_{m-1}?

Remember that we are now looking for a Jacobi rotation that kills
element
*A*_{m-1, m} when applied to matrix
**A** from
the left, i.e.,
:

You can now see that choosing and as given by equations (3.59) and (3.60) will indeed do the job.

Now, let us have a look at the next part of the code:

g = d(i + 1) - p r = (d(i) - g) * s + 2.0_sp * c * b p = s * r d(i + 1) = g + p g = c * r - bTo understand what happens here we

d(i + 1) = g + pwhere

`p = s * r`

, and Looking 3 lines up we see that

In turn

so

Now, if you compare this with equation (3.15) on page in section about Jacobi rotations (section 3.2), you will see that this is simply:

i.e., a plain Jacobi rotation. This is then followed with the corresponding generation of eigenvectors:

IF(PRESENT(z)) THEN ff(1:n) = z(1:n,i+1) z(1:n, i+1) = s*z(1:n,i) + c*ff(1:n) z(1:n,i) = c*z(1:n, i) - s*ff(1:n) END IFwhich mimics the following Jacobi formulae:

But even though for *i* = *m* - 1 the operations correspond indeed
to a Jacobi rotation, on the following iterations for
the operations are somewhat different.
They are the Givens rotations, whose purpose is to restore the tridiagonal
form.

At the end of this process, as we go all the way back to *l*, *d*_{l}and *e*_{l} themselves will become affected and *e*_{m} will become annihilated,
whereas the matrix will remain tridiagonal. This part of the code closes
the `iterate`

loop:

d(l) = d(l) - p e(l) = g e(m) = 0.0

The innermost `DO`

loop corresponds to subjecting an ever shrinking
submatrix of
**A** to a series of *QL* iterations. The iterations
stop once the element *e*_{l} has been annihilated.