Arnoldi iteration: Eigenvalues

Study performance of Arnoldi approximation of eigenvalues on matrices from the Harwell-Boeing collection. Demonstrate use of mixed Fortran-Python coding.

Problem setup

Define Python environment

In [1]:
import numpy as np
import pylab as plt
import scipy as sci
import scipy.io as sio
import scipy.sparse as sp
import os

---------------------------------------------------------------------------
ModuleNotFoundError                       Traceback (most recent call last)
<ipython-input-1-93a97396dbed> in <module>
----> 1 import numpy as np
      2 import pylab as plt
      3 import scipy as sci
      4 import scipy.io as sio
      5 import scipy.sparse as sp

ModuleNotFoundError: No module named 'numpy'
In [2]:
os.getcwd()
Out[2]:
'/home/student/courses/MATH662'

Load a matrix

In [3]:
Asp=sio.mmread("./examples/Arnoldi/tols4000.mtx");
plt.spy(Asp)
Out[3]:
<matplotlib.lines.Line2D at 0x7fcc3fddf940>
In [4]:
sp.issparse(Asp)
Out[4]:
True
In [5]:
A = np.array(Asp.todense('F'))
A.dtype
Out[5]:
dtype('float64')
In [6]:
[m,mcol]=np.shape(A)
[m,mcol]
Out[6]:
[4000, 4000]

Find eigenvalues through np.eigvals for comparison with Arnoldi-Ritz estimates

In [7]:
Lambda=np.linalg.eigvals(A)
In [8]:
Lambda[0:10]
Out[8]:
array([-12.098+0.j, -12.098+0.j, -12.098+0.j, -12.098+0.j, -12.098+0.j,
       -12.098+0.j, -12.098+0.j, -12.098+0.j, -12.098+0.j, -12.098+0.j])
In [9]:
Lambda[3990:4000]
Out[9]:
array([-483.66159+2784.3037669j , -483.66159-2784.3037669j ,
       -517.96745+2875.72490352j, -517.96745-2875.72490352j,
       -504.9651 +2841.47958778j, -504.9651 -2841.47958778j,
       -513.615  +2864.31486254j, -513.615  -2864.31486254j,
       -490.00447+2801.4695464j , -490.00447-2801.4695464j ])

Define Fortran computational routine

If the instruction below leads to errors, the fortran_magic module is probably not yet install. Do so in a terminal by:

sudo pip install -U fortran-magic

In [10]:
%load_ext fortranmagic

In the background, FortranMagic invokes the f2py utility to compile Fortran code and transparently provide an interface for calls from Python

In [11]:
%%fortran -v
subroutine arnoldi(A,b,H,Q,v,nK,m,n)
  integer, intent(in) :: m,n
  integer, intent(out) :: nK
  real*8, intent(in) :: A(m,m), b(m)
  real*8, intent(out) :: H(n+1,n), Q(m,n+1), v(m)
  H = 0.d0
  Q(:,1) = b/NORM2(b)
  do k=1,n
    v = MATMUL(A,Q(:,k)); nK=k
    do j=1,k
      H(j,k) = DOT_PRODUCT(Q(:,j),v)      
      v = v - H(j,k)*Q(:,j)
    end do
    H(k+1,k) = NORM2(v)
    if (H(k+1,k)<1.0d-14) exit      
    Q(:,k+1) = v/H(k+1,k)
  end do
end subroutine arnoldi

Ok. The following fortran objects are ready to use: arnoldi

Invoke Fortran arnoldi routine

Krylov space of dimension $n=10$

Invoke the Arnoldi algorithm to return

In [12]:
n=10
b=np.array(np.ones(m)       , dtype='float64', order='F')
H=np.array(np.zeros([n+1,n]), dtype='float64', order='F')
Q=np.array(np.zeros([m,n+1]), dtype='float64', order='F')
v=np.array(np.zeros([m,1])  , dtype='float64', order='F')
[H,Q,v,nK]=arnoldi(A,b,n)
nK
Out[12]:
10
In [13]:
Ritz = np.linalg.eigvals(H[0:nK,0:nK])
In [14]:
plt.plot(Lambda.real,Lambda.imag,'.',Ritz.real,Ritz.imag,'or')
Out[14]:
[<matplotlib.lines.Line2D at 0x7fcc3ca89250>,
 <matplotlib.lines.Line2D at 0x7fcc3ca89340>]
In [ ]:

Krylov space of dimension $n=20$

Invoke the Arnoldi algorithm to return

In [23]:
n=20
b=np.array(np.random.randn(m)       , dtype='float64', order='F')
H=np.array(np.zeros([n+1,n]), dtype='float64', order='F')
Q=np.array(np.zeros([m,n+1]), dtype='float64', order='F')
v=np.array(np.zeros([m,1])  , dtype='float64', order='F')
[H,Q,v,nK]=arnoldi(A,b,n)
nK
Out[23]:
20
In [24]:
Ritz = np.linalg.eigvals(H[0:nK,0:nK])
In [25]:
plt.plot(Lambda.real,Lambda.imag,'.',Ritz.real,Ritz.imag,'or')
Out[25]:
[<matplotlib.lines.Line2D at 0x7fcc3c922c40>,
 <matplotlib.lines.Line2D at 0x7fcc3c922c70>]
In [ ]:

Krylov space of dimension $n=200$

Invoke the Arnoldi algorithm to return

In [15]:
n=200
b=np.array(np.ones(m)       , dtype='float64', order='F')
H=np.array(np.zeros([n+1,n]), dtype='float64', order='F')
Q=np.array(np.zeros([m,n+1]), dtype='float64', order='F')
v=np.array(np.zeros([m,1])  , dtype='float64', order='F')
[H,Q,v,nK]=arnoldi(A,b,n)
nK
Out[15]:
200
In [16]:
Ritz = np.linalg.eigvals(H[0:nK,0:nK])
In [17]:
plt.plot(Lambda.real,Lambda.imag,'.',Ritz.real,Ritz.imag,'or')
Out[17]:
[<matplotlib.lines.Line2D at 0x7fcc3ca68a00>,
 <matplotlib.lines.Line2D at 0x7fcc3ca68a30>]
In [ ]:

Krylov space of dimension $n=400$

Invoke the Arnoldi algorithm to return

In [18]:
n=400
b=np.array(np.ones(m)       , dtype='float64', order='F')
H=np.array(np.zeros([n+1,n]), dtype='float64', order='F')
Q=np.array(np.zeros([m,n+1]), dtype='float64', order='F')
v=np.array(np.zeros([m,1])  , dtype='float64', order='F')
[H,Q,v,nK]=arnoldi(A,b,n)
nK
Out[18]:
400
In [19]:
Ritz = np.linalg.eigvals(H[0:nK,0:nK])
In [20]:
plt.plot(Lambda.real,Lambda.imag,'.',Ritz.real,Ritz.imag,'or')
Out[20]:
[<matplotlib.lines.Line2D at 0x7fcc3ca3fac0>,
 <matplotlib.lines.Line2D at 0x7fcc3ca3faf0>]
In [ ]: