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 $n\times n$. On output matrix a is replaced with matrix $\boldsymbol{Q} = \boldsymbol{P}_1\cdot\boldsymbol{P}_2\cdot\ldots\boldsymbol{P}_{n-1}$, 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 $n\times n$, 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 a_n-1,n or a_n,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)P_kbecome 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 $\sum_{k=1}^{i-1} \vert a_{i,k}\vert$, i.e., we sum absolute values of all elements in row i up to, but excluding, the diagonal element a_i,i. If that whole row adds up to zero, there is no point zeroing it further, so we simply transfer a_i,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}, a_{i,2}, \ldots, a_{i, i-1}$:

          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

$$\boldsymbol{A}' = \boldsymbol{A} - \boldsymbol{q}\odot\boldsymbol{u}$$ (3.54)

where $\odot$ 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 $\boldsymbol{u}\cdot\boldsymbol{u}$. Consequently $\boldsymbol{q}\odot\boldsymbol{u}$ eventually divides $\boldsymbol{u}\otimes\boldsymbol{u}$ by $\boldsymbol{u}\cdot\boldsymbol{u}$, which implies that it is OK to scale a column or a row before rotating it, because $\boldsymbol{u} = \boldsymbol{x}\mp \vert\boldsymbol{x}\vert \boldsymbol{e}_1$.

The next step evaluates what we have called $\sigma$, i.e., $\sigma = a_{i, 1}^2 + a_{i,2}^2 + \ldots + a_{i, i-1}^2 = \boldsymbol{x} \cdot\boldsymbol{x}$:

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

Now we calculate $\mp \vert\boldsymbol{x}\vert$ 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 $\boldsymbol{P}\cdot\boldsymbol{x} = \pm\vert\boldsymbol{x}\vert \boldsymbol{e}_1$:

           e(i) = scale*g

Now recall that $\boldsymbol{u}\cdot\boldsymbol{u} = 2\vert\boldsymbol{x}\vert^2 \mp 2\vert\boldsymbol{x}\vert x_1$. This is implemented in the code as:

           h = h - f*g

Consequently, now h becomes H, whereas previously it was $\sigma$.

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 $i^{\hbox{\scriptsize {th}}}$row of matrix a:

           a(i, l) = f - g

where f is a_i,i-1 and g is $\mp \vert\boldsymbol{x}\vert x_1$.

At the same time we store u/H is the $i^{\hbox{\scriptsize {th}}}$ 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 $\boldsymbol{p} = \boldsymbol{A}\cdot\boldsymbol{u} / H$. 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 $K = \boldsymbol{u}\cdot\boldsymbol{p}/(2H)$:

           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 $\boldsymbol{A}' = \boldsymbol{A} - \boldsymbol{q}\otimes\boldsymbol{u} - \boldsymbol{u}\otimes\boldsymbol{q}$:

           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 Phas the following form: $\boldsymbol{P} = \boldsymbol{1} - \boldsymbol{u}\otimes\boldsymbol{u}/ H$. 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 $i \ge 2$.

The computation, in turn, reflects the following recursive relation:

$$\begin{eqnarray*}\boldsymbol{Q} &=& \boldsymbol{P}_1\cdot\boldsymbol{P}_2\cdot... ...cdot\boldsymbol{Q}_{j+1}\\ \boldsymbol{Q} &=& \boldsymbol{Q}_1 \end{eqnarray*}$$

Thus the first statement of the code in the IF brackets simply calculates $\frac{\boldsymbol{u}_j}{H}\cdot\boldsymbol{Q}_{j+1}$ and then the second statement evaluates $\boldsymbol{Q}_{j+1} - \boldsymbol{u}_j \otimes \hbox{(first statement)}$.

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)P_kand an identity matrix in the opposite corner.