Introduction

Our example diffusion code, based on Fourier's original paper and on the IBM example code, which we have discussed in P573, called two important subroutines from the PESSL library: the first one was the FFT subroutine, and we have dedicated the previous chapter, 2, to Fast Fourier Transform, the second one was the Eigenvalue problem routine.

At this stage you should review section ``The Eigenvalue Problem'' in P573 before proceeding.

The Eigenvalue problem, in short, is about finding a non-zero solution to the following equation:

$$\boldsymbol{A}\cdot\boldsymbol{x} = \lambda\boldsymbol{x}$$ (3.1)

As I have already mentioned before, if A is a normal matrix, i.e., if

$$\boldsymbol{A}\cdot\boldsymbol{A}^\dagger = \boldsymbol{A}^\dagger\cdot\boldsymbol{A}$$ (3.2)

then there exists a rotation $\boldsymbol{\Lambda}\in\hbox{SO}(n)$such that

$$\boldsymbol{\Lambda}^{-1}\cdot\boldsymbol{A}\cdot\boldsymbol{\Lambda} = \hbox{diag}\left(\lambda_1, \ldots,\lambda_n\right)$$ (3.3)

If A is defective, i.e., if it is not normal, then there is still a linear transformation X such that

$$\boldsymbol{X}^{-1}\cdot\boldsymbol{A}\cdot\boldsymbol{X} = \hbox{diag}\left(\lambda_1, \ldots,\lambda_n\right)$$ (3.4)

but this transformation is no longer of $\hbox{SO}(n)$.

But whichever is the case the columns of either $\boldsymbol{\Lambda}$ or X are simply eigenvectors of A that correspond to appropriate $\lambda$s.

There are actually two types of eigenvectors, which are referred to as left and right eigenvectors. If you took to heart what I have told you about vectors and forms in the past, you will recognise now that right eigenvectors are vectors, and left eigenvectors are forms, and that an orthogonality relationship exists between the two:

$$\begin{eqnarray*}\boldsymbol{A}\cdot\boldsymbol{x}_R &=& \lambda\boldsymbol{x}_R... ...ldsymbol{\chi}^i,\boldsymbol{x}_k\right\rangle &=& \delta^i{}_k \end{eqnarray*}$$

This is the case even for the defective matrices, although for defective matrices it may happen that neither $\{\boldsymbol{x}_k\}$nor $\{\boldsymbol{\chi}^i\}$ form an (orthogonal) basis.

The way most eigenproblem routines work is that matrix Ais nudged, usually iteratively, towards a diagonal form by a sequence of linear transformations:

$$\begin{eqnarray*}&&\boldsymbol{A} \\ &&\boldsymbol{P}^{-1}_1\cdot \boldsymbol{... ...P}_1\cdot \boldsymbol{P}_2 \cdot \ldots\cdot \boldsymbol{P}_k \end{eqnarray*}$$

Once the resulting matrix is diagonal, up to the required accuracy, eigenvectors are returned in

$$\boldsymbol{P}_1\cdot \boldsymbol{P}_2 \cdot \ldots\cdot \boldsymbol{P}_k$$

The solution of an eigenproblem may get quite complicated and costly. Public domain eigenproblem routines are available within LAPACK, e.g.,:

CGEEV

compute for an $N\times N$ complex nonsymmetric matrix A the eigenvalues and, optionally, the left and/or right eigenvectors

CHBEV

compute all the eigenvalues and, optionally, eigenvectors of a complex Hermitian band matrix A

CPTEQR

compute all eigenvalues and, optionally, eigenvectors of a symmetric positive definite tridiagonal matrix

DGEEV

compute for an $N\times N$ real nonsymmetric matrix A, the eigenvalues, and, optionally, the lef and/or right eigenvectors

DSBEV

compute all the eigenvalues and, optionally, eigenvectors of a real symmetric band matrix A

DSPEV

compute all the eigenvalues and, optionally, eigenvectors of a real symmetric matrix A in packed storage

DSYEV

compute all eigenvalues and, optionally, eigenvectors of a real symmetric tridiagonal matrix A

There are many, many more in there specialised for various combinations of symmetries and requirements.

Normally you would want such specialised routines to cover at least the following:

• real, symmetric, tridiagonal
• real, symmetric, banded
• real, symmetric
• real, nonsymmetric
• complex, Hermitian
• complex, non-Hermitian

The purpose of those specialisations is to save time and storage, an important considerations especially when the matrices grow very large.