MODULE Eig_Procedures¶

Module EIG_Procedures exports a large set of procedures for computing (selected) eigenvalues and/or (selected) eigenvectors of a symmetric (tridiagonal) matrix [Lawson_Hanson:1974] [Golub_VanLoan:1996] [Parlett:1998].

The standard real symmetric eigenvalue problem is to find eigenvalues $\lambda$ and eigenvectors $u$ such that

$MAT * u = \lambda * u$

where MAT is a n-by-n real symmetric matrix.

For an input n-by-n dense matrix MAT, this module provides routines for:

the transformation of MAT to tridiagonal form T,

$Q^T * MAT * Q = T$

where Q is a n-by-n orthogonal matrix;
the computation of the eigenvalues $\lambda_i$ and eigenvectors $p_i$ of T,

$P^T * T * P = S$

where S is a n-by-n diagonal matrix with $S(i,i)=\lambda_i$ and P is the n-by-n matrix of eigenvectors of T;
the back-transformation of the eigenvectors $p_i$ of T to eigenvectors $u_i$ of MAT,

$MAT = (Q * P) * S * (Q * P)^T = U * S * U^T$

where U is the n-by-n matrix of eigenvectors of MAT and the eigenvalues $S(i,i)=\lambda_i$ of T are also the eigenvalues of MAT.

The transformation of MAT to tridiagonal form T, the generation of the associated orthogonal factor Q, the computation of the eigenvectors of T and the back-transformation the eigenvectors of T to eigenvectors of MAT are done with multi-threaded and blocked algorithms [Dongarra_etal:1989] [Golub_VanLoan:1996] [Walker:1988]. All algorithms are parallelized with OpenMP [openmp]. Depending on the situation and the algorithm used, it is also possible to compute selected, or only the largest or smallest eigenvalues and the associated eigenvectors.

Currently, STATPACK includes four different algorithms for computing (selected) eigenvalues of a tridiagonal matrix T:

the implicit QR method [Lawson_Hanson:1974] [Golub_VanLoan:1996] [Parlett:1998] ;
the fast Pal-Walker-Kahan variant of the implicit QR method, which does not used square roots [Parlett:1998] ;
the bisection method, which is based on Sturm sequences and requires $O(nk)$ operations to compute k eigenvalues of T [Golub_VanLoan:1996] [Parlett:1998] ;
and a rational QR algorithm, where the shift is determined using Newton’s method and makes it possible to steer the iteration toward desired eigenvalues [Reinsch_Bauer:1968]. This method also allows computation of subset of eigenvalues.

These STATPACK eigenvalues routines generally compute eigenvalues of T (or MAT) to an absolute accuracy of $\epsilon ||T||_2$ (or $\epsilon ||MAT||_2$ ), where $\epsilon$ is the machine precision. If higher accuracy is wanted, subroutine symtrid_bisect() (or select_eigval_cmp3() in case of a dense matrix MAT) should be used with the optional argument ABSTOL set to sqrt(lamch("S")).

Currently, STATPACK includes three different methods for computing eigenvectors of a tridiagonal matrix T:

implicit QR iteration [Golub_VanLoan:1996] [Parlett:1998];
inverse iteration combined with Fernando’s method [Golub_VanLoan:1996] [Ipsen:1997] [Dhillon:1998] [Fernando:1997] [Bini_etal:2005];
and a deflation algorithm inspired from the work of Godunov and collaborators, which is also related to Fernando’s method for computing eigenvectors [Godunov_etal:1993] [Malyshev:2000] and [Fernando:1997].

The implicit QR algorithm applies a sequence of similarity transformations to the tridiagonal matrix T until its off-diagonal elements become negligible and the diagonal elements have converged to the eigenvalues of T [Golub_VanLoan:1996]. It consists of a bulge-chasing procedure that implicitly includes shifts and use plane rotations (e.g. Givens rotations) which preserve the tridiagonal form of T. High performance is obtained by restructuring the QR Algorithm with a wave-front algorithm to accumulate the Givens rotations for computing the eigenvectors [VanZee_etal:2011]. Subset computations are not possible with the QR algorithm, but it is possible to compute only all the eigenvalues or both all the eigenvalues and associated eigenvectors.

If the distance between the eigenvalues of T is sufficient relative to the norm of T, then computing eigenvectors by inverse iteration is a $O(nk)$ process, where k is the number of eigenvectors to compute [Ipsen:1997] [Dhillon:1998]. However, if the eigenvalues of T are too close, the eigenvectors must be orthogonalized by the modified Gram-Schmidt algorithm, which is much more expensive.

The deflation method combines Fernando’s method for the computation of eigenvectors [Fernando:1997] [Parlett_Dhillon:1997] with deflation procedures by Givens rotations, see [Godunov_etal:1993] [Parlett_Dhillon:1997] [Malyshev:2000] for more details. QR iterations with a perfect shift strategy are also used as a back-up procedure if the deflation technique fails [Mastronardi_etal:2006]. If the eigenvalues are well-separated, the deflation method is also a $O(nk)$ process, where k is the number of eigenvectors to compute. It is highly recommended to compute the eigenvalues to high accuracy (e.g. with symtrid_bisect() or select_eigval_cmp3()) for the success of the deflation technique.

Parallelism concerns only the computation of eigenvectors in the QR method, but both the computation of the eigenvalues and the eigenvectors in the bisection-inverse iteration and bisection-deflation methods.

Finally, as explained above, subset computations are not possible in the standard implicit QR algorithm, but is possible with the two other methods for computing eigenvectors and for the bisection method for computing eigenvalues. The m largest or smallest eigenvalues of a symmetric n-by-n tridiagonal matrix T can also be computed using the rational QR method [Reinsch_Bauer:1968]. This rational QR method is sequential.

Please note that driver and computational routines provided in this module apply only to real data of kind stnd. The real kind type stnd is defined in module Select_Parameters. Computation of eigenvalues and eigenvectors for a complex matrix are not provided in this release of STATPACK.

The driver and computational routines provided in this module are different from the corresponding routines provided by LAPACK [Anderson_etal:1999] and are (much) faster if OpenMP is used, but less accurate for the same precision.

In order to use one of these routines, you must include an appropriate use Eig_Procedures or use Statpack statement in your Fortran program, like:

use EIG_Procedures, only: eig_cmp

or :

use Statpack, only: eig_cmp

Here is the list of the public routines exported by module EIG_Procedures:

symtrid_cmp()¶

Purpose:

symtrid_cmp() reduces a real n-by-n symmetric matrix MAT (eventually stored in in packed form) to symmetric tridiagonal form T by an orthogonal similarity transformation:

$Q^T * MAT * Q = T$

The transformation to tridiagonal form T, which use Householder reflectors, is blocked and parallelized with OpenMP if the UPPER argument is not used [Dongarra_etal:1989].

On the other hand, if the UPPER argument is used, the standard sequential and non-blocked algorithm is used [Golub_VanLoan:1996] [Parlett:1998].

Synopsis:

call symtrid_cmp( mat(:n,:n) ,         d(:n) , e(:n) , store_q , upper )
call symtrid_cmp( mat(:n,:n) ,         d(:n) , e(:n) , store_q         )
call symtrid_cmp( matp(:(n*(n+1)/2)) , d(:n) , e(:n) , store_q , upper )
call symtrid_cmp( matp(:(n*(n+1)/2)) , d(:n) , e(:n) , store_q         )

Examples:

ex1_symtrid_cmp.F90

ex2_symtrid_cmp.F90

ex2_trid_deflate.F90

ortho_gen_symtrid()¶

Purpose:

ortho_gen_symtrid() generates a real orthogonal matrix Q which is defined as the product of n-1 elementary reflectors of order n, as returned by symtrid_cmp().

The generation of the real orthogonal matrix Q is blocked and parallelized with OpenMP in all cases [Walker:1988].

Synopsis:

call ortho_gen_symtrid( mat(:n,:n) , upper )
call ortho_gen_symtrid( mat(:n,:n)         )

Examples:

ex1_symtrid_cmp.F90

ex2_symtrid_cmp.F90

apply_q_symtrid()¶

Purpose:

apply_q_symtrid() overwrites the general real m-by-n matrix C with:

$Q * C$ if LEFT = true and TRANS = false ;
$Q^T * C$ if LEFT = true and TRANS = true ;
$C * Q$ if LEFT = false and TRANS = false ;
$C * Q^T$ if LEFT = false and TRANS = true .

where Q is a real orthogonal matrix of order nq and is defined as the product of nq-1 elementary reflectors

$Q = H(nq-1) * ... * H(2) * H(1)$ , if UPPER = true or is not used
$Q = H(1)*H(2) * ... * H(nq-1)$ , if UPPER = false

as returned by symtrid_cmp().

Q is of order m (=:data:nq) and is the product of m-1 reflectors if LEFT = true. Q is of order n (=:data:nq) and is the product of n-1 reflectors if LEFT = false.

The procedure is is blocked and parallelized with OpenMP in all cases [Walker:1988].

Synopsis:

call apply_q_symtrid( mat(:nq,:nq) ,         c(:,:) , left , trans , upper )
call apply_q_symtrid( mat(:nq,:nq) ,         c(:,:) , left , trans         )
call apply_q_symtrid( matp(:(nq*(nq+1)/2)) , c(:,:) , left , trans , upper )

The shape of C must verify:

size( C, 1 ) = nq if LEFT = true
size( C, 2 ) = nq if LEFT = false

symtrid_qri()¶

Purpose:

symtrid_qri() computes all eigenvalues and eigenvectors of a symmetric n-by-n tridiagonal matrix using the implicit QR method [Golub_VanLoan:1996] [Parlett:1998] [Greenbaum_Dongarra:1989] .

The eigenvalues and eigenvectors of a full symmetric matrix can also be found if symtrid_cmp() and ortho_gen_symtrid() have been used to reduce this matrix to tridiagonal form.

The computation of the eigenvectors are parallelized with OpenMP [Demmel_etal:1993].

If eigenvectors are not requested, symtrid_qri() computes all eigenvalues of a symmetric tridiagonal matrix using the fast Pal-Walker-Kahan variant of the QR algorithm, which does not use square roots [Parlett:1998].

Synopsis:

call symtrid_qri( d(:n) , e(:n) , failure , mat(:n,:n) , init_mat=init_mat , sort=sort , maxiter=maxiter )
call symtrid_qri( d(:n) , e(:n) , failure , sort=sort , maxiter=maxiter )

Examples:

ex1_symtrid_qri.F90

ex2_symtrid_qri.F90

symtrid_qri2()¶

Purpose:

symtrid_qri2() computes all eigenvalues and eigenvectors of a symmetric n-by-n tridiagonal matrix with a perfect shift strategy for the eigenvectors [Godunov_etal:1993] [Malyshev:2000] .

If the perfect shift strategy failed for some eigenvectors, these eigenvectors are computed with the implicit QR method [Golub_VanLoan:1996] [Parlett:1998].

The eigenvalues and eigenvectors of a full symmetric matrix can also be found if symtrid_cmp() and ortho_gen_symtrid() have been used to reduce this matrix to tridiagonal form.

Furthermore, the computation of the eigenvectors is blocked with a wave-front algorithm for applying Givens rotations [VanZee_etal:2011] and parellelized with OpenMP [Demmel_etal:1993].

If eigenvectors are not requested, symtrid_qri2() computes all eigenvalues of a symmetric tridiagonal matrix using the fast Pal-Walker-Kahan variant of the QR algorithm, which does not use square roots [Parlett:1998].

symtrid_qri2() is much faster than symtrid_qri() for computing eigenvectors of large matrices, but may be less efficient for small matrices.

Synopsis:

call symtrid_qri2( d(:n) , e(:n) , failure , mat(:n,:n) , init_mat=init_mat , sort=sort , maxiter=maxiter, max_francis_steps=max_francis_steps )
call symtrid_qri2( d(:n) , e(:n) , failure , sort=sort , maxiter=maxiter )

Examples:

ex1_symtrid_qri2.F90

ex2_symtrid_qri2.F90

symtrid_qri3()¶

Purpose:

symtrid_qri3() computes all eigenvalues and eigenvectors of a symmetric n-by-n tridiagonal matrix with the implicit QR method [Golub_VanLoan:1996] [Parlett:1998].

The eigenvalues and eigenvectors of a full symmetric matrix can also be found if symtrid_cmp() and ortho_gen_symtrid() have been used to reduce this matrix to tridiagonal form.

The computation of the eigenvectors is blocked with a wave-front algorithm for applying Givens rotations [VanZee_etal:2011] and parellelized with OpenMP [Demmel_etal:1993].

If eigenvectors are not requested, symtrid_qri3() computes all eigenvalues of a symmetric tridiagonal matrix using the implicit QR method [Golub_VanLoan:1996] [Parlett:1998].

symtrid_qri3() is (slightly) less faster than symtrid_qri2() for computing eigenvectors of large matrices, but may be more robust.

Synopsis:

call symtrid_qri3( d(:n) , e(:n) , failure , mat(:n,:n) , init_mat=init_mat , sort=sort , maxiter=maxiter, max_francis_steps=max_francis_steps )
call symtrid_qri3( d(:n) , e(:n) , failure , sort=sort , maxiter=maxiter )

Examples:

ex1_symtrid_qri3.F90

symtrid_ratqri()¶

Purpose:

symtrid_ratqri() computes the m largest or smallest eigenvalues of a real symmetric n-by-n tridiagonal matrix using a rational QR method [Reinsch_Bauer:1968].

The m largest or smallest eigenvalues of a full symmetric matrix can also be found if symtrid_cmp() and ortho_gen_symtrid() have been used to reduce this matrix to tridiagonal form.

This subroutine is not parallelized, but allows computation of a subset of the eigenvalues of a symmetric n-by-n tridiagonal matrix, a possibility, which is not offered by symtrid_qri(), symtrid_qri2() and symtrid_qri3().

Synopsis:

call symtrid_ratqri( d(:n) , e(:n) , m , failure , small=small , tol=tol )

Examples:

ex1_symtrid_ratqri.F90

symtrid_ratqri2()¶

Purpose:

symtrid_ratqri2() computes the largest or smallest eigenvalues of a symmetric n-by-n tridiagonal matrix in algebraic value whose sum (e.g. sum of the absolute values of the eigenvalues) exceeds a given value. A rational QR method is used [Reinsch_Bauer:1968].

The largest or smallest eigenvalues of a full symmetric matrix whose sum exceeds a given threshold in algebraic value can also be found, if symtrid_cmp() and ortho_gen_symtrid() have been used to reduce this matrix to tridiagonal form.

This subroutine is not parallelized, but allows computation of a subset of the eigenvalues of a symmetric n-by-n tridiagonal matrix, a possibility, which is not offered by symtrid_qri(), symtrid_qri2() and symtrid_qri3().

Synopsis:

call symtrid_ratqri2( d(:n) , e(:n) , val , failure , m , small=small , tol=tol )

Examples:

ex1_symtrid_ratqri2.F90

symtrid_bisect()¶

Purpose:

symtrid_bisect() computes all or some of the largest or smallest eigenvalues of a real n-by-n symmetric tridiagonal matrix T using a bisection method [Golub_VanLoan:1996].

The full set, largest or smallest eigenvalues of a full symmetric matrix can also be found if symtrid_cmp() has been used to reduce this matrix to tridiagonal form.

This subroutine is parallelized if OPENMP is used and allows the computation of eigenvalues to high (relative) accuracy if desired (e.g. with the optional argument ABSTOL).

Synopsis:

call symtrid_bisect( d(:n) , e(:n) , neig , w(:n) , failure , small=small , sort=sort , vector=vector , abstol=abstol , le=le , theta=theta , scaling=scaling , init=init )

Examples:

ex1_symtrid_bisect.F90

ex2_symtrid_bisect.F90

ex2_trid_deflate.F90

eig_cmp()¶

Purpose:

eig_cmp() computes all eigenvalues and eigenvectors of a n-by-n real symmetric matrix MAT.

The matrix MAT is first transformed to tridiagonal form T, then the eigenvalues and the eigenvectors are computed by the QR implicit algorithm [Golub_VanLoan:1996].

The transformation to tridiagonal form T is blocked and parallelized with OpenMP if the UPPER argument is not used [Dongarra_etal:1989]. Otherwise, a sequential and non-blocked algorithm is used for this transformation [Golub_VanLoan:1996] [Parlett:1998].

The computation of the eigenvectors by the QR implicit algorithm is parallelized with OpenMP in all cases [Demmel_etal:1993].

Synopsis:

call eig_cmp( mat(:n,:n) , eigval(:n) , failure , upper , sort=sort , maxiter=maxiter )
call eig_cmp( mat(:n,:n) , eigval(:n) , failure ,         sort=sort , maxiter=maxiter )

Examples:

ex1_eig_cmp.F90

ex2_eig_cmp.F90

eig_cmp2()¶

Purpose:

eig_cmp2() computes all eigenvalues and eigenvectors of a n-by-n real symmetric matrix MAT.

The matrix MAT is first transformed to tridiagonal form T, then the eigenvalues are computed by the Pal-Walker-Kahan variant of the QR algorithm [Parlett:1998] and the eigenvectors are computed with a perfect shift strategy [Godunov_etal:1993] [Malyshev:2000], or the implicit QR algorithm [Parlett:1998] if the perfect shift strategy fails. The Givens rotations, which are used in the perfect shift strategy or the implicit QR algorithm, are accumulated with a fast wavefront algorithm for computing the eigenvectors [VanZee_etal:2011].

The transformation to tridiagonal form T is blocked and parallelized with OpenMP if the UPPER argument is not used [Dongarra_etal:1989]. Otherwise, a sequential and non-blocked algorithm is used [Golub_VanLoan:1996] [Parlett:1998].

Furthermore, the computation of the eigenvectors is blocked with a wave-front algorithm for applying Givens rotations [VanZee_etal:2011] and parellelized with OpenMP in all cases [Demmel_etal:1993].

eig_cmp2() is much faster than eig_cmp() for large matrices, but may be less efficient for small matrices.

Synopsis:

call eig_cmp2( mat(:n,:n) , eigval(:n) , failure , upper , sort=sort , maxiter=maxiter , max_francis_steps=max_francis_steps )
call eig_cmp2( mat(:n,:n) , eigval(:n) , failure ,         sort=sort , maxiter=maxiter , max_francis_steps=max_francis_steps )

Examples:

ex1_eig_cmp2.F90

ex2_eig_cmp2.F90

eig_cmp3()¶

Purpose:

eig_cmp3() computes all eigenvalues and eigenvectors of a n-by-n real symmetric matrix MAT.

The matrix MAT is first transformed to tridiagonal form T, then the eigenvalues/eigenvectors are computed by the implicit QR algorithm [Parlett:1998]. The Givens rotations used in the implicit QR algorithm are accumulated with a fast wavefront algorithm for computing the eigenvectors [VanZee_etal:2011].

The transformation to tridiagonal form T is blocked and parallelized with OpenMP if the UPPER argument is not used [Dongarra_etal:1989]. Otherwise, a sequential and non-blocked algorithm is used [Golub_VanLoan:1996] [Parlett:1998].

Furthermore, the computation of the eigenvectors is blocked with a wave-front algorithm for applying Givens rotations [VanZee_etal:2011] and parellelized with OpenMP [Demmel_etal:1993].

eig_cmp3() is much faster than eig_cmp() for large matrices and less faster than eig_cmp2(), but may be more robust.

Synopsis:

call eig_cmp3( mat(:n,:n) , eigval(:n) , failure , upper , sort=sort , maxiter=maxiter , max_francis_steps=max_francis_steps )
call eig_cmp3( mat(:n,:n) , eigval(:n) , failure ,         sort=sort , maxiter=maxiter , max_francis_steps=max_francis_steps )

Examples:

ex1_eig_cmp3.F90

ex2_eig_cmp3.F90

laev2()¶

Purpose:

laev2() computes the eigendecomposition of a 2-by-2 real symmetric matrix

$\left( \begin{matrix} a & b \\ b & c \end{matrix} \right)$

On return, RT1 is the eigenvalue of larger absolute value, RT2 is the eigenvalue of smaller absolute value, and (CS1,SN1) is the unit right eigenvector for RT1, giving the decomposition

$\left( \begin{matrix} cs1 & sn1 \\ -sn1 & cs1 \end{matrix} \right) \left( \begin{matrix} a & b \\ b & c \end{matrix} \right) \left( \begin{matrix} cs1 & -sn1 \\ sn1 & cs1 \end{matrix} \right) = \left( \begin{matrix} rt1 & O \\ 0 & rt2 \end{matrix} \right)$

This subroutine is from LAPACK.

Synopsis:

call laev2( a , b , c , rt1 , rt2 , cs1 , sn1 )

eig_sort()¶

Purpose:

Given the eigenvalues D and, eventually, eigenvectors U as output from eig_cmp(), eig_cmp2(), eig_cmp3(), symtrid_qri(), symtrid_qri2(), symtrid_qri3() or other STATPACK routines which compute eigenvalues and eigenvectors, eig_sort() sorts the eigenvalues D into ascending or descending order, and, rearranges the columns of U correspondingly.

Synopsis:

call eig_sort( sort , d(:m), u(:n,:m) )
call eig_sort( sort , d(:m)           )

eigvalues()¶

Purpose:

eigvalues() computes all eigenvalues of a n-by-n real symmetric matrix MAT.

The matrix MAT is first transformed to symmetric tridiagonal form T by an orthogonal similarity transformation:

$Q^T * MAT * Q = T$

then the eigenvalues of T, which are also the eigenvalues of MAT, are computed by the Pal-Walker-Kahan variant of the QR algorithm [Parlett:1998].

If the QR algorithm fails to converge eigvalues() returns a n-vector filled with NANs.

Synopsis:

eigval(:n) = eigvalues( mat(:n,:n) , upper , sort=sort , maxiter=maxiter  )
eigval(:n) = eigvalues( mat(:n,:n) ,         sort=sort , maxiter=maxiter  )

Examples:

ex1_eigvalues.F90

eigval_cmp()¶

Purpose:

eigval_cmp() computes all eigenvalues of a n-by-n real symmetric matrix MAT, eventually stored in packed form.

The matrix MAT is first transformed to symmetric tridiagonal form T by an orthogonal similarity transformation:

$Q^T * MAT * Q = T$

then the eigenvalues of T, which are also the eigenvalues of MAT, are computed by the Pal-Walker-Kahan variant of the QR algorithm [Parlett:1998].

Synopsis:

call eigval_cmp( mat(:n,:n) ,         eigval(:n) , failure , upper , sort=sort , maxiter=maxiter , d_e=d_e(:n,:2) )
call eigval_cmp( mat(:n,:n) ,         eigval(:n) , failure ,         sort=sort , maxiter=maxiter , d_e=d_e(:n,:2) )
call eigval_cmp( matp(:(n*(n+1)/2)) , eigval(:n) , failure , upper , sort=sort , maxiter=maxiter , d_e=d_e(:n,:2) )
call eigval_cmp( matp(:(n*(n+1)/2)) , eigval(:n) , failure ,         sort=sort , maxiter=maxiter , d_e=d_e(:n,:2) )

Examples:

ex1_eigval_cmp.F90

ex2_eigval_cmp.F90

eigval_cmp2()¶

Purpose:

eigval_cmp2() computes all eigenvalues of a n-by-n real symmetric matrix MAT, eventually stored in packed form.

The matrix MAT is first transformed to symmetric tridiagonal form T by an orthogonal similarity transformation:

$Q^T * MAT * Q = T$

then the eigenvalues of T, which are also the eigenvalues of MAT, are computed by the Pal-Walker-Kahan variant of the QR algorithm [Parlett:1998].

A slightly different version of the Pal-Walker-Kahan algorithm is used compared to the eigval_cmp() generic subroutine.

Synopsis:

call eigval_cmp2( mat(:n,:n) ,         eigval(:n) , failure , upper , sort=sort , maxiter=maxiter , d_e=d_e(:n,:2) )
call eigval_cmp2( mat(:n,:n) ,         eigval(:n) , failure ,         sort=sort , maxiter=maxiter , d_e=d_e(:n,:2) )
call eigval_cmp2( matp(:(n*(n+1)/2)) , eigval(:n) , failure , upper , sort=sort , maxiter=maxiter , d_e=d_e(:n,:2) )
call eigval_cmp2( matp(:(n*(n+1)/2)) , eigval(:n) , failure ,         sort=sort , maxiter=maxiter , d_e=d_e(:n,:2) )

Examples:

ex1_eigval_cmp2.F90

ex2_eigval_cmp2.F90

eigval_cmp3()¶

Purpose:

eigval_cmp3() computes all eigenvalues of a n-by-n real symmetric matrix MAT, eventually stored in packed form.

The matrix MAT is first transformed to symmetric tridiagonal form T by an orthogonal similarity transformation:

$Q^T * MAT * Q = T$

then the eigenvalues of T, which are also the eigenvalues of MAT, are computed by the implicit QR algorithm [Parlett:1998] [Golub_VanLoan:1996] .

Synopsis:

call eigval_cmp3( mat(:n,:n) ,         eigval(:n) , failure , upper , sort=sort , maxiter=maxiter , d_e=d_e(:n,:2) )
call eigval_cmp3( mat(:n,:n) ,         eigval(:n) , failure ,         sort=sort , maxiter=maxiter , d_e=d_e(:n,:2) )
call eigval_cmp3( matp(:(n*(n+1)/2)) , eigval(:n) , failure , upper , sort=sort , maxiter=maxiter , d_e=d_e(:n,:2) )
call eigval_cmp3( matp(:(n*(n+1)/2)) , eigval(:n) , failure ,         sort=sort , maxiter=maxiter , d_e=d_e(:n,:2) )

Examples:

ex1_eigval_cmp3.F90

ex2_eigval_cmp3.F90

select_eigval_cmp()¶

Purpose:

select_eigval_cmp() computes the m = size(eigval) largest or smallest eigenvalues of a n-by-n real symmetric matrix MAT. MAT can be stored in packed form.

The matrix MAT is first transformed to symmetric tridiagonal form T by an orthogonal similarity transformation:

$Q^T * MAT * Q = T$

then the eigenvalues of T, which are also the eigenvalues of MAT, are computed by a rational QR method [Reinsch_Bauer:1968].

Synopsis:

call select_eigval_cmp( mat(:n,:n) ,         eigval(:m) , small , failure , upper , d_e=d_e(:n,:2) )
call select_eigval_cmp( mat(:n,:n) ,         eigval(:m) , small , failure ,         d_e=d_e(:n,:2) )
call select_eigval_cmp( matp(:(n*(n+1)/2)) , eigval(:m) , small , failure , upper , d_e=d_e(:n,:2) )
call select_eigval_cmp( matp(:(n*(n+1)/2)) , eigval(:m) , small , failure ,         d_e=d_e(:n,:2) )

Examples:

ex1_select_eigval_cmp.F90

ex2_select_eigval_cmp.F90

select_eigval_cmp2()¶

Purpose:

select_eigval_cmp2() computes the largest or smallest eigenvalues of a n-by-n real symmetric matrix MAT whose sum in algebraic value (e.g. the sum of the absolute values) exceeds a given value. MAT can be stored in packed form.

The matrix MAT is first transformed to symmetric tridiagonal form T by an orthogonal similarity transformation:

$Q^T * MAT * Q = T$

then the eigenvalues of T, which are also the eigenvalues of MAT, are computed by a rational QR method [Reinsch_Bauer:1968].

Synopsis:

call select_eigval_cmp2( mat(:n,:n) ,         eigval(:m) , small , val , failure , upper , d_e=d_e(:n,:2) )
call select_eigval_cmp2( mat(:n,:n) ,         eigval(:m) , small , val , failure ,         d_e=d_e(:n,:2) )
call select_eigval_cmp2( matp(:(n*(n+1)/2)) , eigval(:m) , small , val , failure , upper , d_e=d_e(:n,:2) )
call select_eigval_cmp2( matp(:(n*(n+1)/2)) , eigval(:m) , small , val , failure ,         d_e=d_e(:n,:2) )

Examples:

ex1_select_eigval_cmp2.F90

ex2_select_eigval_cmp2.F90

select_eigval_cmp3()¶

Purpose:

select_eigval_cmp3() computes the largest or smallest eigenvalues of a n-by-n real symmetric matrix MAT, eventually stored in packed form.

The matrix MAT is first transformed to symmetric tridiagonal form T by an orthogonal similarity transformation:

$Q^T * MAT * Q = T$

then the eigenvalues of T, which are also the eigenvalues of MAT, are computed by a bisection method [Golub_VanLoan:1996].

The eigenvalues of T can be computed to high accuracy with the optional argument ABSTOL.

Synopsis:

call select_eigval_cmp3( mat(:n,:n) ,         neig , eigval(:n) , small , failure , upper , sort=sort , vector=vector , scaling=scaling , init=init , abstol=abstol , le=le , theta=theta , d_e=d_e(:n,:2) )
call select_eigval_cmp3( mat(:n,:n) ,         neig , eigval(:n) , small , failure ,         sort=sort , vector=vector , scaling=scaling , init=init , abstol=abstol , le=le , theta=theta , d_e=d_e(:n,:2) )
call select_eigval_cmp3( matp(:(n*(n+1)/2)) , neig , eigval(:n) , small , failure , upper , sort=sort , vector=vector , scaling=scaling , init=init , abstol=abstol , le=le , theta=theta , d_e=d_e(:n,:2) )
call select_eigval_cmp3( matp(:(n*(n+1)/2)) , neig , eigval(:n) , small , failure ,         sort=sort , vector=vector , scaling=scaling , init=init , abstol=abstol , le=le , theta=theta , d_e=d_e(:n,:2) )

Examples:

ex1_select_eigval_cmp3.F90

ex2_select_eigval_cmp3.F90

lae2()¶

Purpose:

lae2() computes the eigenvalues of a 2-by-2 symmetric matrix

$\left( \begin{matrix} a & b \\ b & c \end{matrix} \right)$

On return, RT1 is the eigenvalue of larger absolute value, RT2 is the eigenvalue of smaller absolute value.

This subroutine is from LAPACK.

Synopsis:

call lae2( a , b , c , rt1 , rt2 )

eigval_sort()¶

Given the eigenvalues D as output from eigval_cmp(), eigval_cmp2(), eig_valcmp3(), symtrid_qri(), symtrid_qri2(), symtrid_qri3() or other STATPACK routines which compute eigenvalues, eigval_sort() sorts the eigenvalues D into ascending or descending order.

Synopsis:

call eigval_sort( sort , d(:m) )

dflgen()¶

Purpose:

dflgen() computes deflation parameters (e.g. a chain of Givens rotations) for a n-by-n symmetric unreduced tridiagonal matrix T and a given eigenvalue LAMBDA of T.

On output, the arguments CS and SN contain, respectively, the vectors of the cosines and sines coefficients of the chain of n-1 planar rotations that deflates the real n-by-n symmetric tridiagonal matrix T corresponding to the eigenvalue LAMBDA.

See [Malyshev:2000] for more details.

Synopsis:

call dflgen( d(:n) , e(:n-1) , lambda , cs(:n-1) , sn(:n-1) )

dflgen2()¶

Purpose:

dflgen2() computes and applies deflation parameters (e.g. a chain of Givens rotations) for a n-by-n symmetric unreduced tridiagonal matrix T and a given eigenvalue LAMBDA of T.

On output:

the arguments D and E contain, respectively, the new main diagonal and off-diagonal of the deflated symmetric tridiagonal matrix if DEFLATE is set to true.
the arguments CS and SN contain, respectively, the vectors of the cosines and sines coefficients of the chain of n-1 planar rotations that deflates the real n-by-n symmetric tridiagonal matrix T corresponding to the eigenvalue LAMBDA.

See [Malyshev:2000] for more details.

Synopsis:

call dflgen2( d(:n) , e(:n-1) , lambda , cs(:n-1) , sn(:n-1) , deflate )

dflapp()¶

Purpose:

dflapp() deflates a n-by-n real symmetric tridiagonal matrix T by a chain of planar rotations produced by dflgen().

On output, the arguments D and E contain, respectively, the new main diagonal and off-diagonal of the deflated symmetric tridiagonal matrix, if DEFLATE is set to true on output.

See [Malyshev:2000] for more details.

Synopsis:

call dflapp( d(:n) , e(:n-1) , cs(:n-1) , sn(:n-1) , deflate )

qrstep()¶

Purpose:

qrstep() performs one QR step with a given shift LAMBDA on a n-by-n real symmetric unreduced tridiagonal matrix T.

On output, the arguments D and E contain, respectively, the new main diagonal and off-diagonal of the deflated symmetric tridiagonal. The chain of n-1 planar rotations produced during the QR step are saved in the arguments CS and SN. See [Mastronardi_etal:2006] for more details.

Synopsis:

call qrstep( d(:n) , e(:n-1) , lambda , cs(:n-1) , sn(:n-1) , deflate )

prodgiv()¶

Purpose:

prodgiv() applies a chain of n-1 planar rotations produced by dflgen(), dflgen2() or qrstep() to the vector argument X of length n. See [Mastronardi_etal:2006] for more details.

Synopsis:

call prodgiv( cs(:n-1) , sn(:n-1) , x(:n) )

prodgiv_eigvec()¶

Purpose:

prodgiv_eigvec() computes an approximate eigenvector of a n-by-n symmetric tridiagonal matrix from a chain of n-1 planar rotations produced by dflgen(), dflgen2() or qrstep(). See [Mastronardi_etal:2006] for more details.

Synopsis:

eigvec(:n) = prodgiv_eigvec( cs(:n-1) , sn(:n-1) )

symtrid_deflate()¶

Purpose:

symtrid_deflate() computes eigenvectors of a real symmetric tridiagonal matrix T corresponding to specified eigenvalues, using sequential deflation techniques [Godunov_etal:1993] [Malyshev:2000] [Mastronardi_etal:2006].

symtrid_deflate() may fail if some the eigenvalues specified in parameter EIGVAL are nearly identical or for clusters of small eigenvalues or for some zero-diagonal tridiagonal matrices.

symtrid_deflate() is a low-level routine used by the trid_deflate() driver routine. Its direct use as a method for computing eigenvectors of a real symmetric tridiagonal matrix is not recommended.

Synopsis:

call symtrid_deflate( d(:n) , e(:n-1) , eigval     , eigvec(:n)    , failure     , max_qr_steps=max_qr_steps )
call symtrid_deflate( d(:n) , e(:n-1) , eigval(:p) , eigvec(:n,:p) , failure(:p) , max_qr_steps=max_qr_steps )

trid_deflate()¶

Purpose:

trid_deflate() computes all or selected eigenvectors of a n-by-n real symmetric tridiagonal T or full symmetric matrix MAT (eventually packed column-wise in a linear array MATP) corresponding to specified eigenvalues, using deflation techniques [Godunov_etal:1993] [Malyshev:2000] [Mastronardi_etal:2006].

If eigenvectors of a full symmetric matrix MAT are wanted, it is required that the original symmetric matrix MAT has been reduced to symmetric tridiagonal form T by an orthogonal similarity transformation:

$Q^T * MAT * Q = T$

with a call to symtrid_cmp() with parameter STORE_Q set to true, before calling trid_deflate().

The eigenvectors of T are computed using an efficient approach for the computation of (selected) eigenvectors of a tridiagonal matrix corresponding to (selected) eigenvalues by combining Fernando’s method for the computation of eigenvectors with deflation procedures by Givens rotations [Fernando:1997] [Parlett_Dhillon:1997] [Malyshev:2000]. QR iterations are also used as a back-up procedure if the deflation technique fails [Mastronardi_etal:2006].

If eigenvectors of a full symmetric matrix MAT are wanted, the computed eigenvectors of T are backed-transformed to the eigenvectors of MAT using the packed representation of Q as computed by symtrid_cmp().

Both the computation of the eigenvectors of T and the back-transformation phase are parallelized with OpenMP [Walker:1988].

It is essential that eigenvalues given on entry of trid_deflate() are computed to high (relative) accuracy. Subroutine symtrid_bisect() may be used for this purpose.

trid_deflate() may fail if some the eigenvalues specified in parameter EIGVAL are nearly identical for some pathological matrices or for clusters of small eigenvalues or for some zero-diagonal matrices.

The deflation algorithms used in trid_deflate() are competitive with the inverse iteration procedure implemented in trid_inviter() subroutine.

Synopsis:

call trid_deflate( d(:n) , e(:n) , eigval(:p) , eigvec(:n,:p) , failure ,                      max_qr_steps=max_qr_steps , ortho=ortho , inviter=inviter )
call trid_deflate( d(:n) , e(:n) , eigval(:p) , eigvec(:n,:p) , failure , mat(:n,:n)         , max_qr_steps=max_qr_steps , ortho=ortho , inviter=inviter )
call trid_deflate( d(:n) , e(:n) , eigval(:p) , eigvec(:n,:p) , failure , matp(:(n*(n+1)/2)) , max_qr_steps=max_qr_steps , ortho=ortho , inviter=inviter )

Examples:

ex1_trid_deflate.F90

ex2_trid_deflate.F90

ex3_trid_deflate.F90

maxdiag_tinv_qr()¶

Purpose:

maxdiag_tinv_qr() computes the index of the element of maximum absolute value in the diagonal entries of

$( T - lambda * I )^{-1}$

where T is a n-by-n symmetric tridiagonal matrix, I is the identity matrix and lambda is a scalar.

The diagonal entries of ( T - lambda * I )^-1 are computed by means of the QR factorization of T - lambda * I.

For more details, see [Bini_etal:2005].

It is assumed that T is unreduced, but no check is done in the subroutine to verify this assumption.

Synopsis:

maxdiag_tinv = maxdiag_tinv_qr( d(:n) , e(:n-1) , lambda )

maxdiag_tinv_ldu()¶

Purpose:

maxdiag_tinv_ldu() computes the index of the element of maximum absolute value in the diagonal entries of

$( T - lambda * I )^{-1}$

where T is a n-by-n symmetric tridiagonal matrix, I is the identity matrix and lambda is a scalar.

The diagonal entries of ( T - lambda * I )^-1 are computed by means of LDU and UDL factorization of T - lambda * I.

For more details, see [Fernando:1997].

It is assumed that T is unreduced, but no check is done in the subroutine to verify this assumption.

Synopsis:

maxdiag_tinv = maxdiag_tinv_ldu( d(:n) , e(:n-1) , lambda )

trid_qr_cmp()¶

Purpose:

trid_qr_cmp() factorizes the symmetric matrix T - lambda * I, where T is a n-by-n symmetric tridiagonal matrix, I is the identity matrix and lambda is a scalar. as

$T - lambda * I = Q * R$

where Q is an orthogonal matrix represented as the product of n-1 Givens rotations and R is an upper triangular matrix with at most two non-zero super-diagonal elements per column.

The parameter lambda is included in the routine so that trid_qr_cmp() may be used to obtain eigenvectors of T by inverse iteration.

The subroutine also computes the index of the entry of maximum absolute value in the diagonal of ( T - lambda * I )^-1, which provides a good initial approximation to start the inverse iteration process for computing the eigenvector associated with the eigenvalue lambda.

For further details, see [Bini_etal:2005] [Fernando:1997] [Parlett_Dhillon:1997].

Synopsis:

call trid_qr_cmp( d(:n) , e(:n-1) , lambda , cs(:n-1) , sn(:n-1) , diag(:n) , sup1(:n) , sup2(:n) , maxdiag_tinv )

trid_qr_solve()¶

Purpose:

trid_qr_solve() may be used to solve for the vector x the system of equations

$x(:) * (T - lambda * I) = scale * y(:)$

where T is a n-by-n symmetric tridiagonal matrix, I is the n-by-n identity matrix, lambda and scale are scalars, following the factorization of ( T - lambda * I ) by trid_qr_cmp() or gk_qr_cmp(), as

$T - lambda * I = Q * R$

where Q is an orthogonal matrix represented as the product of n-1 Givens rotations and R is an upper triangular matrix with at most two non-zero super-diagonal elements per column.

The matrix ( T - lambda * I ) is assumed to be ill-conditioned, and frequent rescalings are carried out in order to avoid overflow. However, no test for singularity or near-singularity is included in this routine. Such tests must be performed before calling this routine. The scalar scale is not output by this routine since this routine being intended for use in applications such as inverse iteration.

Synopsis:

call trid_qr_solve( cs(:n-1) , sn(:n-1) , diag(:n) , sup1(:n) , sup2(:n) , y(:n) )

trid_cmp()¶

Purpose:

trid_cmp() factorizes the symmetric matrix $( T - eigval * I )$ , where T is a n-by-n symmetric tridiagonal matrix, I is the n-by-n identity matrix, eigval is a scalar, as

$T - eigval * I = P * L * U$

where P is a permutation matrix, L is a unit lower tridiagonal matrix with at most one non-zero sub-diagonal elements per column and U is an upper triangular matrix with at most two non-zero super-diagonal elements per column.

The factorizations is obtained by Gaussian elimination with partial pivoting and implicit row scaling.

Several symmetric matrices $( T - eigval_i * I )$ can be handled in a single call.

The parameter EIGVAL is included in the routine so that trid_cmp() may be used to obtain eigenvectors of T by inverse iteration.

Synopsis:

call trid_cmp( d(:n) , e(:n) , eigval     , sub(:n)    , diag(:n)    , sup1(:n)    , sup2(:n)    , perm(:n)    , tol=tol )
call trid_cmp( d(:n) , e(:n) , eigval(:p) , sub(:p,:n) , diag(:p,:n) , sup1(:p,:n) , sup2(:p,:n) , perm(:p,:n) , tol=tol )

trid_cmp2()¶

Purpose:

trid_cmp2() factorizes the symmetric matrix $( T - eigval * I )$ , where T is a n-by-n symmetric tridiagonal matrix, I is the n-by-n identity matrix, eigval is a scalar, as

$T - eigval * I = P * L * U$

where P is a permutation matrix, L is a unit lower tridiagonal matrix with at most one non-zero sub-diagonal elements per column and U is an upper triangular matrix with at most two non-zero super-diagonal elements per column.

The factorization is obtained by Gaussian elimination with partial pivoting and row interchanges.

Several symmetric matrices $( T - eigval_i * I )$ can be handled in a single call.

The parameter EIGVAL is included in the routine so that trid_cmp2() may be used to obtain eigenvectors of T by inverse iteration.

trid_cmp2() is a simplified version of trid_cmp().

Synopsis:

call trid_cmp2( d(:n) , e(:n) , eigval     , sub(:n)    , diag(:n)    , sup1(:n)    , sup2(:n)    , perm(:n)    )
call trid_cmp2( d(:n) , e(:n) , eigval(:p) , sub(:p,:n) , diag(:p,:n) , sup1(:p,:n) , sup2(:p,:n) , perm(:p,:n) )

trid_solve()¶

Purpose:

trid_solve() may be used to solve for a vector x the system of equations

$x * (T - eigval * I) = scale*y$

where T is a n-by-n symmetric tridiagonal matrix, I is the n-by-n identity matrix, eigval and scale are scalars and y is a vector, following the factorization of $( T - eigval * I )$ by trid_cmp() or trid_cmp2() as

$T - eigval * I = P * L * U$

where P is a permutation matrix, L is a unit lower tridiagonal matrix with at most one non-zero sub-diagonal elements per column and U is an upper triangular matrix with at most two non-zero super-diagonal elements per column.

The matrix $( T - eigval * I )$ is assumed to be ill-conditioned, and frequent rescalings are carried out in order to avoid overflow. However, no test for singularity or near-singularity is included in this routine. Such tests must be performed before calling this routine.

Several systems $( T - eigval_i * I )$ can be handled in a single call.

The scalar scale is not output by this routine since this routine being intended for use in applications such as inverse iteration.

Synopsis:

call trid_solve( sub(:n)    , diag(:n)    , sup1(:n)    , sup2(:n)    , perm(:n)    , y(:n)    )
call trid_solve( sub(:p,:n) , diag(:p,:n) , sup1(:p,:n) , sup2(:p,:n) , perm(:p,:n) , y(:p,:n) )

trid_inviter()¶

Purpose:

trid_inviter() computes all or selected eigenvectors of a n-by-n real symmetric tridiagonal T or full symmetric matrix MAT (eventually packed column-wise in a linear array MATP) corresponding to specified eigenvalues, using inverse iteration and Fernando’s method [Golub_VanLoan:1996] [Ipsen:1997] [Fernando:1997] [Bini_etal:2005].

If eigenvectors of a full symmetric matrix MAT are wanted, it is required that the original symmetric matrix MAT has been reduced to symmetric tridiagonal form T by an orthogonal similarity transformation:

$Q^T * MAT * Q = T$

with a call to symtrid_cmp() with parameter STORE_Q set to true, before calling trid_inviter().

The eigenvectors of T are computed using inverse iteration, combined with Fernando’s method, for all the eigenvalues at one step. The eigenvectors are then orthogonalized by the Modified Gram-Schmidt algorithm if the eigenvalues are not well-separated [Golub_VanLoan:1996] [Ipsen:1997].

If eigenvectors of a full symmetric matrix MAT are wanted, the computed eigenvectors of T are backed-transformed to the eigenvectors of MAT using the packed representation of Q as computed by symtrid_cmp().

Both the computation of the eigenvectors of T and the back-transformation phase are parallelized with OpenMP [Walker:1988].

trid_inviter() may fail if some the eigenvalues specified in parameter EIGVAL are nearly identical for some pathological matrices.

Synopsis:

call trid_inviter( d(:n) , e(:n) , eigval     , eigvec(:n)    , failure ,                      maxiter=maxiter , scaling=scaling , initvec=initvec )
call trid_inviter( d(:n) , e(:n) , eigval(:p) , eigvec(:n,:p) , failure ,                      maxiter=maxiter , ortho=ortho ,  backward_sweep=backward_sweep , scaling=scaling , initvec=initvec )
call trid_inviter( d(:n) , e(:n) , eigval(:p) , eigvec(:n,:p) , failure , mat(:n,:n)         , maxiter=maxiter , ortho=ortho ,  backward_sweep=backward_sweep , scaling=scaling , initvec=initvec )
call trid_inviter( d(:n) , e(:n) , eigval(:p) , eigvec(:n,:p) , failure , matp(:(n*(n+1)/2)) , maxiter=maxiter , ortho=ortho ,  backward_sweep=backward_sweep , scaling=scaling , initvec=initvec )

Examples:

ex1_trid_inviter.F90

ex2_trid_inviter.F90

ex3_trid_inviter.F90

gen_symtrid_mat()¶

Purpose:

gen_symtrid_mat() generates different types of symmetric tridiagonal matrices with known eigenvalues or specific numerical properties such as clustered eigenvalues and so on for testing purposes of eigensolvers.

Optionally, the eigenvalues of the selected symmetric tridiagonal matrix can be computed analytically, if possible, or by a bisection algorithm with high absolute and relative accuracies.

Synopsis:

call gen_symtrid_mat( type , d(:n) , e(:n) , failure=failure , known_eigval=known_eigval , eigval=eigval , sort=sort , val1=val1 , val2=val2 , l0=l0 , glu0=glu0 )

Previous topic

Next topic

MODULE Eig_Procedures¶