SUBROUTINE tred2(a, d, e, novectors) USE nrtype; USE nrutil, ONLY : assert_eq, outerprod IMPLICIT NONE REAL(sp), DIMENSION(:,:), INTENT(inout) :: a REAL(sp), DIMENSION(:), INTENT(out) :: d, e LOGICAL(lgt), OPTIONAL, INTENT(in) :: novectors INTEGER(i4b) :: i, j, l, n REAL(sp) :: f, g, h, hh, scale REAL(sp), DIMENSION(SIZE(a, 1)) :: g LOGICAL(lgt), SAVE :: yesvec = .TRUE. n = assert_eq(SIZE(a, 1), SIZE(a, 2), SIZE(d), SIZE(e), 'tred2') IF(PRESENT(novectors)) yesvec = .NOT. novectors DO i = n, 2, -1 l = i-1 h = 0.0 IF(l > 1) THEN scale = SUM(ABS(a(i, 1:l))) IF (scale == 0.0) THEN e(i) = a(i, l) ELSE a(i, 1:l) = a(i, 1:l) / scale h = SUM(a(i, 1:l) ** 2) f = a(i, l) g = -SIGN(SQRT(h), f) e(i) = scale*g h = h - f*g a(i, l) = f - g IF (yesvec) a(1:l, i) = a(i, 1:l)/h DO j = 1, l e(j) = (DOT_PRODUCT(a(j, 1:j), a(i, 1:j)) & + DOT_PRODUCT(a(j+1:l, j), a(i, j+1:l)))/h END DO f = DOT_PRODUCT(e(1:l), a(i,1:l)) hh=f/(h+h) e(1:l)=e(1:l) - hh*a(i,1:l) DO j=1, l a(j,1:j) = a(j,1:j)-a(i,j)*e(1:j)-e(j)*a(i,1:j) END DO END IF ELSE e(i)=a(i,l) END IF d(i) = h END DO IF (yesvec) d(1)=0.0 e(1) = 0.0 DO i=1, n IF(yesvec) THEN l=i-1 IF (d(i) /= 0.0) THEN gg(1:l) = MATMUL(a(i,1:l),a(1:l,1:l)) a(1:l,1:l) = a(1:l,1:l) - outerprod(a(1:l,i),gg(1:l)) END IF d(i) = a(i, i) a(i, i) = 1.0 a(i, 1:l) = 0.0 a(1:l, i) = 0.0 ELSE d(i) = a(i, i) END IF END DO END SUBROUTINE tred2

Subroutine `tred2`

performs a Housholder reduction to a triadiagonal
form of a real symmetric matrix `a`

, whose dimensions are
.
On output matrix `a`

is replaced with
matrix
,
the diagonal elements are written on array `d`

and
the *n*-1 off-diagonal elements are written on array `e`

with
`e(1)=0`

. The argument `novectors`

is optional:

LOGICAL(lgt), OPTIONAL, INTENT(in) :: novectorsIf it is set to

`.TRUE.`

then matrix
`a`

contains garbage.
You can test for the presence of optional parameters with the keyword
`PRESENT`

:

IF(PRESENT(novectors)) yesvec = .NOT. novectors

The action begins by checking dimensions of variables passed to the
subroutine as parameters, and if they are correct, i.e., if `a`

is ,
then *n* is extracted and written on `n`

.
Otherwise an error message which contains the word `'tred2'`

is
printed and the subroutine exits.

The main loop has the form

DO i = n, 2, -1 l = i-1 h = 0.0 IF(l > 1) THEN blah... blah... blah... ELSE e(i)=a(i,l) END IF d(i) = h END DOThe first thing to observe is that in this program we're moving not from the upper left corner towards the lower right corner, but the other way round, i.e., our index

`i`

starts at the last column
number `n`

, and then moves down towards the beginning of the matrix.
In other words our first evaluation will yield a new value for
e(i)=a(i,l)But note that by this stage

`a(i,l)`

is going to be something
quite different from what it was at the beginning. The first element
of array `e`

is set to `0`

after the main `DO`

loop
exits following the `IF`

statement:IF (yesvec) d(1)=0.0 e(1) = 0.0

Now let us have a look at the `blah... blah... blah...`

clause
in the main `DO`

loop. What we find inside there is another
`IF`

statement:

scale = SUM(ABS(a(i, 1:l))) IF (scale == 0.0) THEN e(i) = a(i, l) ELSE blah... blah... blah... END IFFirst we evaluate , i.e., we sum absolute values of all elements in row

`e(i)`

and go on to rotate the next row out of existence.
At this stage we have exhausted all possible avenues for avoidance
behaviour and have to do some real work, i.e., have to fill the
`blah... blah... blah...`

bit above.

The first thing we do is to scale :

a(i, 1:l) = a(i, 1:l) / scaleThis we do for the sake of improving numerical accuracy and stability. Is this kosher though? Remember that

(3.54) |

where is a symmetric tensor product, and

The next step evaluates what we have called , i.e., :

h = SUM(a(i, 1:l) ** 2)

Now we calculate
and the sign we choose is the
opposite to the sign of
*a*_{i, i-1}:

f = a(i, l) g = -SIGN(SQRT(h), f)The off-diagonal element that we're going to save on the array

`e`

is now going to be that `g`

scaled back up to the original
value. Remember
we have shown that
:e(i) = scale*g

Now recall that . This is implemented in the code as:

h = h - f*gConsequently, now

`h`

becomes
Vector
**u** is equal to the original vector
**x**everywhere with the exception of its head. So, if we just modify the
latter, we'll have
**u** stored in the
row of matrix `a`

:

a(i, l) = f - gwhere

`f`

is `g`

is
.
At the same time we store
**u**/*H* is the
column of matrix `a`

for future use in case the operator
**Q** is to be evaluated:

IF (yesvec) a(1:l, i) = a(i, 1:l)/h

The next step is to evaluate . This is done by the following code fragment:

DO j = 1, l e(j) = (DOT_PRODUCT(a(j, 1:j), a(i, 1:j)) & + DOT_PRODUCT(a(j+1:l, j), a(i, j+1:l)))/h END DOThe reason for such equilibristics is that by now we have already corrupted matrix

`e`

.
The next two lines evaluate :

f = DOT_PRODUCT(e(1:l), a(i,1:l)) hh=f/(h+h)And once we have

e(1:l)=e(1:l) - hh*a(i,1:l)And the result is once again stored on the unused portion of array

`e`

. We won't need vector
Finally we are ready to perform the final :

DO j=1, l a(j,1:j) = a(j,1:j)-a(i,j)*e(1:j)-e(j)*a(i,1:j) END DO

Remember that we have promised to return the diagonal elements on
array `d`

, but so far we haven't made any use of that array.
Luckily, for the time being the new diagonal elements are stored
safely on *a*_{i,i}. If matrix
**Q** is not wanted,
then we're done and we return with just that:

e(1) = 0.0 DO i=1, n IF(yesvec) THEN blah... blah... blah... ELSE d(i) = a(i, i) END IF END DO

If
**Q** is wanted though, then we make use of vectors
**u** and
**u**/*H*, which are now stored on
rows and on columns of matrix `a`

. The operation
**P**has the following form:
.
We continue to use matrix `a`

as a working
storage.

The logic of this part of the code is as follows:

d(1)=0.0 DO i=1, n l=i-1 IF (d(i) /= 0.0) THEN gg(1:l) = MATMUL(a(i,1:l),a(1:l,1:l)) a(1:l,1:l) = a(1:l,1:l) - outerprod(a(1:l,i),gg(1:l)) END IF blah... blah... blah... END DORemember that before we got to this place

`d(i)`

was set to
`h`

, which, in turn, was set to `d(1)=0`

and
then using the condition `IF(d(i) /= 0.0)`

protects us from
ever reaching for `a(0,...)`

. So the computation within the
`IF`

statement is performed only for .
The computation, in turn, reflects the following recursive relation:

Thus the first statement of the code in the

`IF`

brackets
simply calculates
and then
the second statement evaluates
.
Now we clean up for the next iteration:

d(i) = a(i, i) a(i, i) = 1.0 a(i, 1:l) = 0.0 a(1:l, i) = 0.0The diagonal element is finally transferred to array

`d`

, and
the corner portion of