#!/usr/bin/env python3
"``naginterfaces.library.sparseig.feast_poly_gen_solve`` Python Example."
# NAG Copyright 2020.
# pylint: disable=invalid-name,too-many-arguments,too-many-locals
import math
import numpy as np
from naginterfaces.library import machine, sparse, sparseig
from naginterfaces.base import utils as b_utils
def contour(handle):
""" Define the custom contour and create the nodes and weights """
zedge = np.array([-1.5 - 0.5j, 0.5 + 1.5j])
tedge = np.array([100, 0], dtype=b_utils.EngineIntNumPyType)
nedge = np.array([20, 10], dtype=b_utils.EngineIntNumPyType)
sparseig.feast_custom_contour(handle, nedge, tedge, zedge)
return handle
def input_data():
""" Get input data for the problem """
# The degree of the polynomial and the matrix size
deg = 2
n = 4
# The 3 columns of a, irowa and icola correspond to the matrices
# A0, A1 and A2, stored in sparse CS format. nnza contains the
# number of nonzero elements in A0, A1 and A2
nnza = np.array([5, 6, 4], dtype=b_utils.EngineIntNumPyType)
a = np.array([
[7.5 + 0.0j, 0.0 + 1.0j, 1.0 + 0.0j],
[-1.7 + 0.0j, 0.0 + 1.2j, 1.7 + 0.0j],
[0.0 + 0.4j, 0.0 + 2.9j, 0.2 + 0.0j],
[0.2 + 0.0j, 0.0 - 3.4j, -5.8 + 0.0j],
[-1.6 + 0.0j, 0.0 + 0.1j, 0.0 + 0.0j],
[0.0 + 0.0j, 0.0 + 1.2j, 0.0 + 0.0j]
])
irowa = np.array([
[1, 1, 1],
[2, 1, 2],
[3, 2, 3],
[3, 2, 4],
[4, 3, 0],
[0, 4, 0]
]).astype(int)
icola = np.array([
[2, 1, 1],
[1, 3, 2],
[3, 2, 3],
[4, 4, 4],
[2, 3, 0],
[0, 4, 0]
]).astype(int)
return deg, n, nnza, a, irowa, icola
def form_matrix_pencil(ze, deg, n, nnza, a, irowa, icola, conjg):
""" Form the sparse matrix sum ze^i A_i or (sum ze^i A_i)^H """
nnzsum = np.sum(nnza)
az_tmp = np.zeros(nnzsum, dtype=np.complex128)
irowz_tmp = np.zeros(nnzsum, dtype=b_utils.EngineIntNumPyType)
icolz_tmp = np.zeros(nnzsum, dtype=b_utils.EngineIntNumPyType)
ind = 0
for i in range(deg+1):
if conjg:
az_tmp[ind:ind+nnza[i]] = np.conj((ze**i)*a[0:nnza[i], i])
irowz_tmp[ind:ind+nnza[i]] = icola[0:nnza[i], i]
icolz_tmp[ind:ind+nnza[i]] = irowa[0:nnza[i], i]
else:
az_tmp[ind:ind+nnza[i]] = (ze**i)*a[0:nnza[i], i]
irowz_tmp[ind:ind+nnza[i]] = irowa[0:nnza[i], i]
icolz_tmp[ind:ind+nnza[i]] = icola[0:nnza[i], i]
ind = ind + nnza[i]
# Sort elements into correct coordinate storage format
az, irowz, icolz, _ = sparse.complex_gen_sort(
n, az_tmp, irowz_tmp, icolz_tmp, 'S', 'R',
)
return az, irowz, icolz
def incomplete_lu(az, irowz, icolz, n):
""" Form the incomplete LU factorization of a sparse matrix """
nnaz = np.size(az)
az = np.concatenate((az, np.zeros(2*nnaz)))
irowz = np.concatenate((
irowz,
np.zeros(2*nnaz, dtype=b_utils.EngineIntNumPyType)
))
icolz = np.concatenate((
icolz,
np.zeros(2*nnaz, dtype=b_utils.EngineIntNumPyType)
))
az, irowz, icolz, ipivp, ipivq, istr, idiag, _, _ \
= sparse.complex_gen_precon_ilu(
nnaz, az, irowz, icolz, n,
0, 0.0, 'P', 'N',
)
return az, irowz, icolz, ipivp, ipivq, istr, idiag
def solve_lin_sys(
nnz_az, az, irowz, icolz, ipivp, ipivq, istr,
idiag, y,
):
""" Solve a sparse linear system """
n = np.size(y, 0)
m0 = np.size(y, 1)
w = np.zeros(n, dtype=np.complex128)
for i in range(m0):
w[0:n] = y[0:n, i]
y[0:n, i] = 1.0
y[0:n, i], _, _ = sparse.complex_gen_solve_ilu(
'RGMRES', nnz_az, az, irowz, icolz, ipivp,
ipivq, istr, idiag, w, n,
math.sqrt(machine.precision()), 500, y[0:n, i])
def sparse_matmul(k, nnza, a, irowa, icola, x, z, transa):
""" Perform sparse matrix multiplcation with k^th matrix in the pencil """
n = np.size(x, 0)
m0 = int(np.size(x, 1)/2)
for i in range(m0):
x[0:n, i] = sparse.complex_gen_matvec(
transa, a[0:nnza[k], k], irowa[0:nnza[k], k],
icola[0:nnza[k], k], z[0:n, i], 'N',
)
[docs]def main():
"""
Example for :func:`naginterfaces.library.sparseig.feast_poly_gen_solve`.
Nonsymmetric Polynomial Eigenvalue Problem.
>>> main()
naginterfaces.library.sparseig.feast_poly_gen_solve Python Example Results.
Nonsymmetric Polynomial Eigenvalue Problem.
Eigenvalues:
[-1.01154017+0.74243953j -0.97330836+0.60663851j]
Right eigenvectors:
[[-0.27867645-0.55398883j 0.64961648+0.07275759j]
[ 0.15963493-0.05141476j -0.0367055 +0.19246955j]
[ 0.20169124+0.24677364j -0.02640785+0.10024232j]
[-0.022247 -0.02329904j 0.04061443-0.01888747j]]
Left eigenvectors:
[[-0.00184768-0.01019648j -0.03812357+0.02278326j]
[ 0.01537005-0.00110981j -0.02825777-0.04867545j]
[-0.7176863 -0.28399747j -0.5525527 +0.4415414j ]
[-0.01783012+0.00584953j -0.01808489+0.04505064j]]
"""
print(
'naginterfaces.library.sparseig.feast_poly_gen_solve '
'Python Example Results.'
)
print('Nonsymmetric Polynomial Eigenvalue Problem.')
# Get the input data for the problem
deg, n, nnza, a, irowa, icola = input_data()
# Initialize the FEAST data handle
handle = sparseig.feast_init()
# Set an option using the option setting routine
sparseig.feast_option(handle, 'Max Iter = 50')
# Get nodes and weights for the contour
handle = contour(handle)
# We need more arrays and variables before calling the solver
irevcm = 0
ze = 0.0 + 0.0j
m0 = n
eps = 0.0
itera = 0
x = np.zeros((n, 2*m0), dtype=np.complex128)
y = np.zeros((n, m0), dtype=np.complex128)
z = np.zeros((n, 2*m0), dtype=np.complex128)
d = np.zeros(m0, dtype=np.complex128)
resid = np.zeros(2*m0, dtype=np.float64)
# Call the reverse communication solver
while True:
irevcm, ze, k, m0, nconv, eps, itera = \
sparseig.feast_poly_gen_solve(
handle, irevcm, deg, ze, x,
y, m0, d, z, eps, itera, resid,
)
if irevcm == 0:
break
if irevcm == 1:
# Form the sparse matrix sum ze^i A_i
az, irowz, icolz = form_matrix_pencil(ze, deg, n, nnza, a,
irowa, icola, False)
nnz_az = np.size(az)
# Form incomplete LU factorization of sum ze^i A_i
az, irowz, icolz, ipivp, ipivq, istr, idiag = \
incomplete_lu(az, irowz, icolz, n)
if irevcm == 2:
# Solve (sum ze^i A_i) y = w, with m0 righthand sides
solve_lin_sys(nnz_az, az, irowz, icolz, ipivp, ipivq,
istr, idiag, y)
if irevcm == 3:
# Form the sparse matrix (sum ze^i A_i)^H
azh, irowzh, icolzh = form_matrix_pencil(ze, deg, n, nnza, a,
irowa, icola, True)
nnz_azh = np.size(azh)
# Form incomplete LU factorization of (sum ze^i A_i)^H
azh, irowzh, icolzh, ipivph, ipivqh, istrh, idiagh = \
incomplete_lu(azh, irowzh, icolzh, n)
if irevcm == 4:
# Solve (sum ze^i A_i)^H y = w, with m0 righthand sides
solve_lin_sys(nnz_azh, azh, irowzh, icolzh, ipivph, ipivqh,
istrh, idiagh, y)
if irevcm == 5:
# Compute x <- A_k z
sparse_matmul(k, nnza, a, irowa, icola, x, z, 'N')
if irevcm == 6:
# Compute x <- A_k^H z
sparse_matmul(k, nnza, a, irowa, icola, x, z, 'T')
print("Eigenvalues:")
print(d[0:nconv])
print("Right eigenvectors:")
print(z[0:n, 0:nconv])
print("Left eigenvectors:")
print(z[0:n, m0:m0+nconv])
if __name__ == '__main__':
import doctest
import sys
sys.exit(
doctest.testmod(
None, verbose=True, report=False,
optionflags=doctest.REPORT_NDIFF,
).failed
)
