Next: A Little Summary Up: The Householder Reduction Previous: The Householder Reduction

### Example Code

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) :: novectors

If it is set to .TRUE. then matrix Q is not computed and on exit 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 DO

The 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 an-1,n or an,n-1, depending on whether we're going to work on the lower or upper triangular part of the matrix. Remember that as we move towards the end, our operators (n-k)Pkbecome smaller and smaller, until the last one is simply the identity, which is why for i = 2 and l = 1 we have:
        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 IF

First we evaluate , i.e., we sum absolute values of all elements in row i up to, but excluding, the diagonal element ai,i. If that whole row adds up to zero, there is no point zeroing it further, so we simply transfer ai,i-1 to 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) / scale

This 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 q = p - K u. Now those vectors p and q are defined by dividing various combinations of u by H, which is . Consequently eventually divides by , which implies that it is OK to scale a column or a row before rotating it, because .

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 ai, 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*g

Consequently, now h becomes H, whereas previously it was .

Vector u is equal to the original vector xeverywhere 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 - g

where f is ai,i-1 and 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 DO

The reason for such equilibristics is that by now we have already corrupted matrix A and we must carefully pick up the right fragments of it in order to perform this multiplication. Remember that u is stored in row i. We store this result on an unused part of array e.

The next two lines evaluate :

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

And once we have K we can evaluate vector q = p - Ku:
           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 p any more so we can overwrite it with q.

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 ai,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 Phas 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 DO

Remember that before we got to this place d(i) was set to h, which, in turn, was set to H. Setting 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.0

The diagonal element is finally transferred to array d, and the corner portion of A that should be set to an identity matrix is set to an identity matrix: remember that the full matrix P comprised a smaller matrix (n-k)Pkand an identity matrix in the opposite corner.

Next: A Little Summary Up: The Householder Reduction Previous: The Householder Reduction
Zdzislaw Meglicki
2001-02-26