File:Arnoldi.f: Difference between revisions

$ \renewcommand{\vec}[1]{ {\bf #1} } \newcommand{\bnabla}{ \vec{\nabla} } \newcommand{\Rey}{Re} \def\vechat#1{ \hat{ \vec{#1} } } \def\mat#1{#1} $

The Arnoldi Method

The Arnoldi method is a method for calculating the eigenvalues and eigenvectors of a matrix, i.e. for calculating the scalar $\sigma$ and $n$-vectors $\vec{x}$ that satisfy \[

  \sigma\,A = A\,\vec{x}

\] for a given $n\times n$ matrix $A$.

The main advantage of the method is that it only requires calculations of multiplies by $A$ for a given $\vec{x}$ -- it does not need to know $A$ itself. This means that $A$ need not even be stored, and could correspond to a very complex linear 'action' on $\vec{x}$, e.g. a time integral with initial condition $\vec{x}$. The method seeks eigenvectors in $\mathrm{span}\{\vec{x},\,A\vec{x},\,A^2\vec{x},...\}$, but uses Gram-Schmidt orthogonalisation to improve the suitability of the basis. An extra vector must be stored after each iteration. It is possible to restart without losing information. This has not been implemented.

The method finds the largest eigenvalues and those most separated in the complex plane first. If $A$ is expected to have many negative eigenvalues of little interest, it may be better to work with $\tilde{A}=\mathrm{e}^A=1+A+\frac{1}{2}A^2+...$, which has no effect on the eigenvectors but exponentiates the eigenvalues, $\tilde{\sigma}=\mathrm{e}^\sigma$. The negative eigenvalues then correspond to bunched eigenvalues close to the origin. The Arnoldi method then favours the largest $|\tilde{\sigma}|$, being the $\sigma$ with largest real parts.

Note that time integration corresponds to exponentiation: An eigenvector of a differential equation with a strongly negative eigenvalue decays. Upon time integration, the decay corresponds to multiplication of the initial condition by a value close to zero, i.e. the eigenvalue of the action of the time integral is close to zero. Time integration of a marginal eigenvector corresponds to multiplication by an eigenvalue on the unit circle.

How to use the code

To download, click the link above.

The subroutine arnold(...) needs to be passed a subroutine that calculates the dot product of two eigenvectors. It should look like, for example,

double precision function dotprod(n,a,b)
   implicit none
   integer :: n
   double precision :: a(n), b(n)
   dotprod = sum(a*b)
end function dotprod

arnold(...) needs to be called repeatedly. It communicates the status of the computation via the flag ifail, which tells the user how many eigenvalues are converged up to a given tolerance, to multiply a vector by $A$ again, or tells the user if the method has failed, e.g. reached maximum number of vector that can be stored.

An example of use of the code:

  ! declare workspace vectors, h, q, b... - see header of arnoldi.f
  sv = ... ! random initial vector x
  k = 0    ! initialise iteration counter
  do while(.true.)
     call arnold(n,k,kmax,ncgd,dotprod,tol,sv,h,q,b,wr,wi,ifail)
     if(ifail==-1) then
        print*, ' arnoldi converged!'
        exit
     else if(ifail==0) then
        call multA(sv, sv)      ! possibly complicated routine that multiplies sv by A
     else if(ifail==1) then
        print*, 'WARNING: arnoldi reached max its'
     else if(ifail>=2) then
        print*, 'WARNING: arnoldi error:', ifail
        exit
     end if   
  end do

On exit, the eigenvectors are stored in columns of b, in order corresponding to eigenvalues in wr and wi. If the first eigenvalue is real (wi(1)==0.), the eigenvector occupies the first column of b only. If the next eigenvalue is complex, the real and imaginary parts will occupy the next two columns of b.