1. Solve the shift linear system (A - sigma B) x = b

2. FFT for matrix-vector multiplication with matrix T

3. Solve FFT-based preconditioning linear system

4. Solve the generalized eigenvalue problem A x = lambda B x by Krylov method

5. Solve the linear system K x = b in the null-space free eigenvalue problem (NFSEP)

6. Solve NFSEP K u = lambda u by Lanczos method

7. Solve NFSEP K u = lambda u by SIRA

8. MATLAB codes

License: Copyright (c) 2022 by Tsung-Ming Huang and Wen-Wei Lin

Reference: Xing-Long Lyu, Tiexiang Li, Tsung-Ming Huang, Jia-Wei Lin, Wen-Wei Lin, and Sheng Wang, FAME: fast algorithms for Maxwell's equations for three-dimensional photonic crystals, ACM Trans. Math. Softw. 47, 3, Article 26 (June), 2021.

This program computes the periodic Schur decomposition (PSD) of the matrix product

Bp^{-1} * Ap * ... * B1^{-1} * A1, where p is the period.

In addition, we can also select to accumulate the orthogonal matices Q_i's or Z_i's. That is

Q_j * B_j * Z_j = Bj, Q_j * A_j * Z_{j-1} = Aj, Z_0 = Z_p, for j = 1: p,

where A1 is qusi-upper triangular and A2, ..., Ap, B1, ..., Bp are upper triangular matrices.

• Download (PQZ.zip)

• License: Copyright (c) 2016 by Hung-Yuan Fan, Tsung-Ming Huang and Wen-Wei Lin

This program is a structure-preserving solver for finding the eigenpairs with the largest magnitude eigenvalues of the eigenproblem

B x = lambda x,

where B is T-skew-Hamiltonian. It is applied T-skew-Hamiltonian Arnoldi algorithm to compute the target eigenpairs with the largest magnitude eigenvalues.

Syntax:

[ev, ew, inv_u, hess_h, outer_it] = SHIRA_Gen(mtx_B, nn, eigswanted, opts)

Input:

• mtx_B : coefficient matrix or the function handle instead of the explicitly matrix form with dim = 2*nn.

• nn : 2*nn is the matrix dimension of the coefficient matrix.

• eigswanted : the number of the wanted largest magnitude eigenvalues.

• opts : specify options

• opts.tol : tolerance of the convergence (default is eps)

• opts.p : number of Arnoldi vectors to restart, eigswanted+1<p<=N (default is 3*eigswanted)

• opts.maxit : maximum number of iterations (default is 300)

• opts.v0 : initial starting vector (default is random vector)

Output:

• ev : convergence right eigenvectors

• ew : convergence eigenvalues

• inv_u : T-symplectic matrix spanning the invariant subspace

• hess_h : upper Hessenberg matrix for mtx_B * inv_u = inv_u * hess_h

• outer_it : number of the iterations of Arnoldi method

II. Second type of T-skew-Hamiltonian B:

B = N1^{-1} K N2^{-1}

• K = - tau L J L^T J^T,

• N1 = ( M - tau L ),

• N2 = J ( M^T - tau L^T ) J^T,

• M = [ A1, 0; -A0, -I ],

• L = [ 0, I; A1^T, 0],

• J = [0, I; -I, 0].

Note that (M, L) is a T-sympletic pair.

Syntax:

[ev, ew, inv_u, hess_h, outer_it] = SHIRA_TSymplP(mtx_A0, mtx_A1, mtx_A1T, linear_solver_A1A0A1T, linear_solver_A1TA0A1, nn, tau, eigswanted, opts)

Input:

• mtx_A0, mtx_A1 : coefficient matrices or the function handles of matrix-vector multiplications.

• mtx_A1T : the transpose matrix of mtx_A1, i.e., mtx_A1T = mtx_A1.' or the function handles of matrix-vector with matrix mtx_A1.'

• linear_solver_A1A0A1T : the function handle for solving the linear system

(tau^2 mtx_A1 + tau mtx_A0 + mtx_A1.') x = b.

Syntax: x = linear_solver_A1A0A1T( tau, b)

• linear_solver_A1TA0A1 : the function handle for solving the linear system

(tau^2 mtx_A1.' + tau mtx_A0 + mtx_A1) x = b.

Syntax: x = linear_solver_A1TA0A1( tau, b)

• nn : the matrix dimension of the coefficient matrices.

• tau : shift value

• eigswanted : the number of the wanted largest magnitude eigenvalues.

• opts : specify options

• opts.tol : tolerance of the convergence (default is eps)

• opts.p : number of Arnoldi vectors to restart, eigswanted+1<p<=N (default is 3*eigswanted)

• opts.maxit : maximum number of iterations (default is 300)

• opts.v0 : initial starting vector (default is random vector)

Output:

• ev : convergence right eigenvectors

• ew : convergence eigenvalues

• inv_u : T-symplectic matrix spanning the invariant subspace

• hess_h : upper Hessenberg matrix for mtx_B * inv_u = inv_u * hess_h

• outer_it : number of the iterations of Arnoldi method

• Download (SHIRA.zip)

• License: Copyright (c) 2018 by Tsung-Ming Huang and Wen-Wei Lin

• Reference:

1. Tsung-Ming Huang, Wen-Wei Lin and Jiang Qian, Structure-preserving algorithms for palindromic quadratic eigenvalue problems arising from vibration on fast trains, SIAM J. Matrix Anal. Appl., Vol. 30, pp. 1566-1592, 2009.

This program is a structure-preserving solver for finding the eigenpairs with the eigenvalues closest to a given shift of the eigenvalue problem

M * z = lambda * L * z,

where

M = [ mtx_A, 0; mtx_Q, -I ], L = [ 0, I; mtx_A.', 0 ], z = [ x; y ].

Note that (M, L) is a T-sympletic pair and (lambda, x) is an eigenpair of the T-Palindromic Quadratic Eigenvalue Problem

P(lambda) x = (lambda^2 * mtx_A.' - lambda * mtx_Q + mtx_A) x = 0.

The program is applied Generalized T-skew-Hamiltonian Arnoldi algorithm to compute the target eigenpairs with the eigenvalues closest to a given shift.

Syntax:

[ ew_new, ev_new, outer_it ] = TSymplecticPair_GTSHIRA( n, mtx_A, mtx_Q, solve_LS_PQEP, shift, eigenwanted, restart_info, opts, R, H, Y, Z)

Input:

• n : the matrix dimension of the coefficient matrices mtx_A and mtx_Q.

• mtx_A : the function handle of matrix-vector multiplications mtx_A*u and mtx_A.'*u.

Syntax: x = mtx_A( u, transp_flag)

if transp_flag = 'transp', then the function handle computes

x = mtx_A.' * u

else (transp_flag = 'notransp') it computes

x = mtx_A * u

end

• mtx_Q : coefficient matrix or the function handle of matrix-vector multiplication mtx_Q*v.

• solve_LS_PQEP : the function handle for solving the linear system

Syntax: x = solve_LS_PQEP( b, transp_flag)

if transp_flag = 'transp', then the function handle solves

(shift^2 * mtx_A.' - shift * mtx_Q + mtx_A).' x = b

else (transp_flag = 'notransp') it solves

(shift^2 * mtx_A.' - shift * mtx_Q + mtx_A) x = b

end

• shift : shift value

• eigswanted : the number of the wanted largest magnitude eigenvalues

• restart_info : 'no' or 'yes'

'yes' : need to input Y, Z, H, R which satisfy the initial generalized T-isotropic Arnoldi factorization, i.e.,

mtx_K * Z = Y * H, mtx_N * Z = Y * R,

where (K,N) is T-skew-Hamiltonian.

'no' : using random initial vector as an initial vector of GTSHIRA

• Y, Z : orthonormal matrices, Y and Z are T-bi-isotropic (option)

• H : upper Hessenberg matrix (option)

• R : upper triangular matrix (option)

• opts : specify options

• opts.tol : tolerance of the convergence (default is eps)

• opts.p : number of Arnoldi vectors to restart, eigswanted+1<p<=N (default is 3*eigswanted)

• opts.maxit : maximum number of iterations (default is 300)

• opts.v0 : initial starting vector (default is random vector)

Output:

• ew_new : one dimensional array (convergence eigenvalues)

ew_new(1:2:end) save the eigenvalues lambda with abs(lambda) <= 1

ew_new(2:2:end) save the eigenvalues lambda with abs(lambda) >= 1

ew_new(2*i-1) * ew_new(2*i) = 1

• ev_new : two-dimensional array (convergence right eigenvectors)

ev_new(:,1:2:end) save the eigenvectors corresponding to ew_new(1:2:end)

ev_new(:,2:2:end) save the eigenvectors corresponding to ew_new(2:2:end)

• outer_it : number of the iterations of Arnoldi method

Remark:

( ew_new, ev_new(1:n,:) ) are the eigenpairs of T-Palindromic Quadratic Eigenvalue Problem.

• Download (GTSHIRA.zip)

• License: Copyright (c) 2019 by Tsung-Ming Huang and Wen-Wei Lin

• Reference:

1. Tsung-Ming Huang, Wen-Wei Lin and Jiang Qian, Structure-preserving algorithms for palindromic quadratic eigenvalue problems arising from vibration on fast trains, SIAM J. Matrix Anal. Appl., Vol. 30, pp. 1566-1592, 2009.

2. Tsung-Ming Huang, Wen-Wei Lin, Heng Tian and Guan-Hua Chen, Computing the Full Spectrum of Large Sparse Palindromic Quadratic Eigenvalue Problems arising from Surface Green's Function calculations, Journal of Computational Physics, Vol. 356, 340-355, 2018.

We provide a generator of the coefficient matrices of the generalized eigenvalue problems that arise from modeling leaky surface wave propagation in an acoustic resonator with an infinite amount of periodically arranged interdigital transducers. The constitutive equations are discretized by finite element methods with mesh refinements along the electrode interfaces and corners.

Syntax:

[mtx_M1, mtx_M2, mtx_F, mtx_G] = mtx_SAW_filter(test_case, omega, mesh_size, refine_level)

Input variables:

• test_case : benchmark problem ('LiNbO_3', or, 'LiTaO_3', or 'quartz')

• omega : frequency

• mesh_size : length of the mesh

• refine_level : level of the refining (0: uniform mesh, 2: non-uniform mesh)

• Download (mtx_SAW_filter.zip)

• License: Copyright (c) 2018 by Tsung-Ming Huang, Wen-Wei Lin and Chin-Tien Wu

• Reference:

1. Tsung-Ming Huang, Wen-Wei Lin and Chin-Tien Wu, Structure-preserving Arnoldi-type algorithm for solving eigenvalue problems in Leaky surface wave propagation, Applied Mathematics and Computation, Vol. 219, 9947-9958, 2013.

2. Tsung-Ming Huang, Tiexiang Li, Wen-Wei Lin and Chin-Tien Wu, Numerical studies on structure-preserving algorithms for surface acoustic wave simulations, Journal of Computational and Applied Mathematics, Vol. 244, pp. 140-154, 2013.

SAW_filter_TSHIRA and SAW_filter_GTSHIRA are matlab functions which are used T-skew-Hamiltonian Arnoldi algorithm and Generalized T-skew-Hamiltonian Arnoldi algorithm, respectively, to compute the target eigenpairs of the generalized eigenvalue problem arising leaky surface wave propagation in an acoustic resonator with an infinite amount of periodically arranged interdigital transducers.

Syntax:

[EW, EV, res_norm, NRes, err_recip, outer_it] = SAW_filter_TSHIRA(mtx_G, mtx_F, mtx_M1, mtx_M2, tau, eigswanted, opts)

Input variables:

• mtx_G . : complex matrix

• mtx_F . : complex matrix

• mtx_M1 . : complex symmetric matrix with dimension n

• mtx_M2 . : complex symmetric matrix with dimension m where m < n

• tau . : the shift value

• eigswanted : number of wanted eigenvalues which are closest to tau

• opts : specify options

• opts.tol : tolerance of the convergence (default is eps)

• opts.p : number of Arnoldi vectors to restart, eigswanted+1<p<=N (default is 3*eigswanted)

• opts.maxit : maximum number of iterations (default is 300)

• opts.v0 : initial starting vector (default is random vector)

Output variables:

• EW : Matrix with two columns. The entries in the first and second column of this matrix are the target eigenvalues which are located inside and outside, respectively, the unit circle.

• EV : three dimensional array. The columns of EV(:,:,1) and EV(:,:,2) are the eigenvectors corresponding to eigenvalues EW(:,1) and EW(:,2), respectively.

• res_nrom : two dimensional array. The first and second columns are represented the residuals of (EW(:,1),EV(:,:,1)) and (EW(:,2),EV(:,:,2)), respectively.

• NRes . : two dimensional array. The first and second columns are represented the normalized residuals of (EW(:,1),EV(:,:,1)) and (EW(:,2),EV(:,:,2)), respectively.

• err_recip : the reciprocal error of (EW(:,1), EW(:,2)).

• outer_it . : number of outer iterations of TSHIRA.

• Download (SAW_filter_TSHIRA.zip)

• License: Copyright (c) 2018 by Tsung-Ming Huang and Wen-Wei Lin