Module_Eig_Procedures¶

This file is part of statpack.

statpack is free software: you can redistribute it and/or modify it under the terms of the GNU Lesser General Public License as published by the Free Software Foundation, either version 3 of the License, or (at your option) any later version.

statpack is distributed in the hope that it will be useful, but WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public License for more details.

You can find a copy of the GNU Lesser General Public License in the statpack/doc directory.

Authors: Pascal Terray (LOCEAN/IPSL, Paris, France)

MODULE EXPORTING PROCEDURES FOR COMPUTING (SELECTED) EIGENVALUES AND/OR (SELECTED) EIGENVECTORS OF A SYMMETRIC (TRIDIAGONAL) MATRIX.

SUBROUTINES FOR COMPUTING A PARTIAL EIGENVALUE DECOMPOSITION OF SYMMETRIC MATRICES BASED ON RANDOMIZED ALGORITHMS ARE ALSO PROVIDED.

LATEST REVISION : 11/03/2024

`subroutine symtrid_cmp ( mat, d, e, store_q, upper )`¶

Purpose¶

SYMTRID_CMP reduces a real n-by-n symmetric matrix MAT to symmetric tridiagonal form T by an orthogonal similarity transformation:

Q’ * MAT * Q = T

Arguments¶

MAT (INPUT/OUTPUT) real(stnd), dimension(:,:)

On entry:

If UPPER = true : The leading n-by-n upper triangular part of MAT contains the upper triangular part of the symmetric matrix MAT, and the strictly lower triangular part of MAT is not referenced.

If UPPER = false : The leading n-by-n lower triangular part of MAT contains the lower triangular part of the symmetric matrix MAT, and the strictly upper triangular part of MAT is not referenced.

On exit:

If UPPER = true and STORE_Q = true : The leading n-by-n upper triangular part of MAT is overwritten by the matrix Q as a product of elementary reflectors.

If UPPER = false and STORE_Q = true : The leading n-by-n lower triangular part of MAT is overwritten by the matrix Q as a product of elementary reflectors.

If UPPER = true and STORE_Q = false : The leading n-by-n upper triangular part of MAT is destroyed.

If UPPER = false and STORE_Q = false : The leading n-by-n lower triangular part of MAT is destroyed.

D (OUTPUT) real(stnd), dimension(:)

The diagonal elements of the tridiagonal matrix T:

D(i) = T(i,i).

The size of D must verify: size( D ) = size( MAT, 1 ) = size( MAT, 2 ) = n .

E (OUTPUT) real(stnd), dimension(:)

The off-diagonal elements of the tridiagonal matrix T:

E(i) = T(i,i+1) = T(i+1,i)

E(n) is arbitrary.

The size of E must verify: size( E ) = size( MAT, 1 ) = size( MAT, 2 ) = n .

STORE_Q (INPUT) logical(lgl)

On exit:

If UPPER = true and STORE_Q = true : The leading n-by-n upper triangular part of MAT is overwritten by the matrix Q as a product of elementary reflectors.

If UPPER = false and STORE_Q = true : The leading n-by-n lower triangular part of MAT is overwritten by the matrix Q as a product of elementary reflectors.

If UPPER = true and STORE_Q = false : The leading n-by-n upper triangular part of MAT is destroyed.

If UPPER = false and STORE_Q = false : The leading n-by-n lower triangular part of MAT is destroyed.

UPPER (INPUT) logical(lgl)

Specifies whether the upper or lower triangular part of the symmetric matrix MAT is stored. If:

UPPER= true : Upper triangular is stored ;

UPPER= false: Lower triangular is stored .

Further Details¶

If UPPER = true and STORE_Q = true, the matrix Q is represented as a product of elementary reflectors

Q = H(n-1) * … * H(2) * H(1).

Each H(i) has the form

H(i) = I + tau * v * v’

where tau is a real scalar, and v is a real vector with v(i+1:n) = 0 ; v(1:i) is stored on exit in MAT(1:i,i+1), and tau in MAT(i+1,i+1).

If UPPER = false and STORE_Q = true, the matrix Q is represented as a product of elementary reflectors

Q = H(1) * H(2) * … * H(n-1).

Each H(i) has the form

H(i) = I + tau * v * v’

where tau is a real scalar, and v is a real vector with v(1:i) = 0 ; v(i+1:n) is stored on exit in MAT(i+1:n,i), and tau in MAT(i,i).

The contents of MAT on exit are illustrated by the following examples with n = 5:

if UPPER = true and STORE_Q = true :

( xx v1 v2 v3 v4 )

( yy t1 v2 v3 v4 )

( yy yy t2 v3 v4 )

( yy yy yy t3 v4 )

( yy yy yy yy t4 )

if UPPER = false and STORE_Q = true :

( t1 yy yy yy yy )

( v1 t2 yy yy yy )

( v1 v2 t3 yy yy )

( v1 v2 v3 t4 yy )

( v1 v2 v3 v4 xx )

where vi and ti denote an element of the vector v and the scalar tau defining H(i), respectively. xx = machhuge and is used by the subroutine ORTHO_GEN_SYMTRID in order to verify that SYMTRID_CMP has been called before ORTHO_GEN_SYMTRID. Elements yy are not modified by the subroutine.

This subroutine is adapted from the routine DSYTD2 in LAPACK. Note that this subroutine is not blocked and not parallelized.

For more details on the reduction algorithm used in SYMTRID_CMP, see:

Golub, G.H., and Van Loan, C.F., 1996:

Matrix Computations. 3rd ed. The Johns Hopkins University Press, Baltimore.

`subroutine symtrid_cmp ( mat, d, e, store_q )`¶

Purpose¶

SYMTRID_CMP reduces a real n-by-n symmetric matrix MAT to symmetric tridiagonal form T by an orthogonal similarity transformation:

Q’ * MAT * Q = T

Arguments¶

MAT (INPUT/OUTPUT) real(stnd), dimension(:,:)

On entry:

The leading n-by-n upper triangular part of MAT contains the upper triangular part of the symmetric matrix MAT, and the strictly lower triangular part of MAT is not referenced.

On exit:

If STORE_Q = true : The leading n-by-n upper triangular part of MAT is overwritten by the matrix Q as a product of elementary reflectors.

If STORE_Q = false : The leading n-by-n upper triangular part of MAT is destroyed.

D (OUTPUT) real(stnd), dimension(:)

The diagonal elements of the tridiagonal matrix T:

D(i) = T(i,i).

The size of D must verify: size( D ) = size( MAT, 1 ) = size( MAT, 2 ) = n.

E (OUTPUT) real(stnd), dimension(:)

The off-diagonal elements of the tridiagonal matrix T:

E(i) = T(i,i+1) = T(i+1,i)

E(n) is arbitrary.

The size of E must verify: size( E ) = size( MAT, 1 ) = size( MAT, 2 ) = n.

STORE_Q (INPUT) logical(lgl)

On exit:

If STORE_Q = true : The leading n-by-n upper triangular part of MAT is overwritten by the matrix Q as a product of elementary reflectors and the lower triangular part of MAT is destroyed. See Further Details.

If STORE_Q = false : The symmetric matrix MAT is destroyed.

Further Details¶

If STORE_Q = true, the matrix Q is represented as a product of elementary reflectors

Q = H(n-1) * … * H(2) * H(1).

Each H(i) has the form

H(i) = I + tau * v * v’

where tau is a real scalar, and v is a real vector with v(i+1:n) = 0 ; v(1:i) is stored on exit in MAT(1:i,i+1), and tau in MAT(i+1,i+1).

The contents of MAT on exit are illustrated by the following example with n = 5:

( xx v1 v2 v3 v4 )

( yy t1 v2 v3 v4 )

( yy yy t2 v3 v4 )

( yy yy yy t3 v4 )

( yy yy yy yy t4 )

where vi and ti denote an element of the vector v and the scalar tau defining H(i), respectively. xx = machhuge and is used by the subroutine ORTHO_GEN_SYMTRID in order to verify that SYMTRID_CMP has been called before ORTHO_GEN_SYMTRID. Elements yy are not modified by the subroutine.

This subroutine is adapted from the routine DSYTD2 in LAPACK. An efficient blocked algorithm is used to reduced the n-by-n symmetric matrix MAT to tridiagonal form T. Furthermore, the computations are parallelized if OPENMP is used.

In other words, SYMTRID_CMP is much more efficient then SYMTRID_CMP with argument UPPER, which is not blocked and not parallelized.

For further details, see:

Golub, G.H., and Van Loan, C.F., 1996:

Matrix Computations. 3rd ed. The Johns Hopkins University Press, Baltimore.

Dongarra, J.J., Sorensen, D.C., and Hammarling, S.J., 1989:

Block reduction of matrices to condensed form for eigenvalue computations. J. of Computational and Applied Mathematics, Vol. 27, pp. 215-227.

`subroutine symtrid_cmp ( matp, d, e, store_q, upper )`¶

Purpose¶

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

Q’ * MAT * Q = T

Arguments¶

MATP (INPUT/OUTPUT) real(stnd), dimension(:)

On entry, the upper or lower triangle of the symmetric matrix MAT, packed column-wise in a linear array. The j-th column of MAT is stored in the array MATP as follows:

if UPPER = true, MATP(i + (j-1) * j/2) = MAT(i,j) for 1<=i<=j;

if UPPER = false, MATP(i + (j-1) * (2*n-j)/2) = MAT(i,j) for j<=i<=n.

On exit:

If UPPER = true and STORE_Q = true : The leading n-by-n upper triangular part of MAT is overwritten by the matrix Q as a product of elementary reflectors.

If UPPER = false and STORE_Q = true : The leading n-by-n lower triangular part of MAT is overwritten by the matrix Q as a product of elementary reflectors.

If UPPER = true and STORE_Q = false : The leading n-by-n upper triangular part of MAT is destroyed.

If UPPER = false and STORE_Q = false : The leading n-by-n lower triangular part of MAT is destroyed.

The size of MATP must verify: size( MATP ) = (n * (n+1)/2)

D (OUTPUT) real(stnd), dimension(:)

The diagonal elements of the tridiagonal matrix T:

D(i) = T(i,i).

The size of D must verify: size( D ) = n .

E (OUTPUT) real(stnd), dimension(:)

The off-diagonal elements of the tridiagonal matrix T:

E(i) = T(i,i+1) = T(i+1,i)

E(n) is arbitrary.

The size of E must verify: size( E ) = n .

STORE_Q (INPUT) logical(lgl)

On exit:

If UPPER = true and STORE_Q = true : The leading n-by-n upper triangular part of MAT is overwritten by the matrix Q as a product of elementary reflectors.

If UPPER = false and STORE_Q = true : The leading n-by-n lower triangular part of MAT is overwritten by the matrix Q as a product of elementary reflectors.

If UPPER = true and STORE_Q = false : The leading n-by-n upper triangular part of MAT is destroyed.

If UPPER = false and STORE_Q = false : The leading n-by-n lower triangular part of MAT is destroyed.

UPPER (INPUT) logical(lgl)

Specifies whether the upper or lower triangular part of the symmetric matrix MAT is stored. If:

UPPER = true : Upper triangle of MAT is stored;

UPPER = false: Lower triangle of MAT is stored.

Further Details¶

If UPPER = true and STORE_Q = true, the matrix Q is represented as a product of elementary reflectors

Q = H(n-1) * … * H(2) * H(1).

Each H(i) has the form

H(i) = I + tau * v * v’

where tau is a real scalar, and v is a real vector with v(i+1:n) = 0 ; v(1:i) and tau are stored on exit in MATP, overwriting MAT(1:i,i+1) and MAT(i+1,i+1), respectively.

If UPPER = false and STORE_Q = true, the matrix Q is represented as a product of elementary reflectors

Q = H(1) * H(2) * … * H(n-1).

Each H(i) has the form

H(i) = I + tau * v * v’

where tau is a real scalar, and v is a real vector with v(1:i) = 0 ; v(i+1:n) and tau are stored on exit in MATP, overwriting MAT(i+1:n,i) and MAT(i,i), respectively.

The contents of MATP on exit are illustrated by the following examples with n = 5:

if UPPER = true and STORE_Q = true, MAT is equal to :

( xx v1 v2 v3 v4 )

( yy t1 v2 v3 v4 )

( yy yy t2 v3 v4 )

( yy yy yy t3 v4 )

( yy yy yy yy t4 )

if UPPER = false and STORE_Q = true, MAT is equal to :

( t1 yy yy yy yy )

( v1 t2 yy yy yy )

( v1 v2 t3 yy yy )

( v1 v2 v3 t4 yy )

( v1 v2 v3 v4 xx )

where vi and ti denote an element of the vector v and the scalar tau defining H(i), respectively. Elements yy are not used and not stored in MATP. xx = machhuge and is used by other subroutines in order to verify that SYMTRID_CMP has been called.

This subroutine is adapted from the routine DSPTRD in LAPACK. Note that this subroutine is not blocked and not parallelized.

For more details on the reduction algorithm used in SYMTRID_CMP, see:

Golub, G.H., and Van Loan, C.F., 1996:

Matrix Computations. 3rd ed. The Johns Hopkins University Press, Baltimore.

`subroutine symtrid_cmp ( matp, d, e, store_q )`¶

Purpose¶

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

Q’ * MAT * Q = T

Arguments¶

MATP (INPUT/OUTPUT) real(stnd), dimension(:)

On entry, the upper triangle of the symmetric matrix MAT, packed column-wise in a linear array. The j-th column of MAT is stored in the array MATP as follows:

MATP(i + (j-1) * j/2) = MAT(i,j) for 1<=i<=j;

On exit:

If STORE_Q = true : The leading n-by-n upper triangular part of MAT is overwritten by the matrix Q as a product of elementary reflectors.

If STORE_Q = false : The leading n-by-n upper triangular part of MAT is destroyed.

The size of MATP must verify: size( MATP ) = (n * (n+1)/2)

D (OUTPUT) real(stnd), dimension(:)

The diagonal elements of the tridiagonal matrix T:

D(i) = T(i,i).

The size of D must verify: size( D ) = n .

E (OUTPUT) real(stnd), dimension(:)

The off-diagonal elements of the tridiagonal matrix T:

E(i) = T(i,i+1) = T(i+1,i)

E(n) is arbitrary.

The size of E must verify: size( E ) = n .

STORE_Q (INPUT) logical(lgl)

On exit:

If STORE_Q = true : The leading n-by-n upper triangular part of MAT is overwritten by the matrix Q as a product of elementary reflectors.

If STORE_Q = false : The leading n-by-n upper triangular part of MAT is destroyed.

Further Details¶

The matrix Q is represented as a product of elementary reflectors

Q = H(n-1) * … * H(2) * H(1).

Each H(i) has the form

H(i) = I + tau * v * v’

where tau is a real scalar, and v is a real vector with v(i+1:n) = 0 ; v(1:i) and tau are stored on exit in MATP, overwriting MAT(1:i,i+1) and MAT(i+1,i+1), respectively, if STORE_Q = true.

The contents of MATP (if STORE_Q = true) on exit are illustrated by the following example with n = 5 (giving the contents of MAT):

( xx v1 v2 v3 v4 )

( yy t1 v2 v3 v4 )

( yy yy t2 v3 v4 )

( yy yy yy t3 v4 )

( yy yy yy yy t4 )

where vi and ti denote an element of the vector v and the scalar tau defining H(i), respectively. Elements yy are not used and not stored in MATP. xx = machhuge and is used by other subroutines in order to verify that SYMTRID_CMP has been called.

This subroutine is adapted from the routine DSPTRD in LAPACK. An efficient blocked algorithm is used to reduced the n-by-n symmetric matrix MAT to tridiagonal form T. Furthermore, the computations are parallelized if OPENMP is used.

In other words, SYMTRID_CMP is much more efficient then SYMTRID_CMP with argument UPPER which is not blocked and not parallelized.

For further details, see:

Golub, G.H., and Van Loan, C.F., 1996:

Matrix Computations. 3rd ed. The Johns Hopkins University Press, Baltimore.

Dongarra, J.J., Sorensen, D.C., and Hammarling, S.J., 1989:

Block reduction of matrices to condensed form for eigenvalue computations. J. of Computational and Applied Mathematics, Vol. 27, pp. 215-227.

`subroutine symtrid_cmp2 ( mat, d, e, store_q )`¶

Purpose¶

SYMTRID_CMP2 reduces a real n-by-n symmetric matrix product

MAT’ * MAT

, where MAT is a m-by-n matrix with m >= n, to symmetric tridiagonal form T by an orthogonal similarity transformation:

Q’ * MAT’ * MAT * Q = T

where Q is an n-by-n orthogonal matrix.

SYMTRID_CMP2 computes T and Q, using the one-sided Ralha tridiagonal reduction algorithm without explicitly forming the matrix product MAT’ * MAT .

Arguments¶

MAT (INPUT/OUTPUT) real(stnd), dimension(:,:)

On entry, the general m-by-n matrix to be reduced.

On exit:

If STORE_Q = true : The leading n-by-n lower triangular part of MAT is overwritten by the matrix Q as a product of elementary Householder reflectors. The other part of MAT is also destroyed.

If STORE_Q = false : The m-by-n matrix MAT is destroyed.

The shape of MAT must verify: size( MAT, 1 ) >= size( MAT, 2 ) = n .

D (OUTPUT) real(stnd), dimension(:)

The diagonal elements of the tridiagonal matrix T:

D(i) = T(i,i).

The size of D must verify: size( D ) = size( MAT, 1 ) = size( MAT, 2 ) = n .

E (OUTPUT) real(stnd), dimension(:)

The off-diagonal elements of the tridiagonal matrix T:

E(i) = T(i,i+1) = T(i+1,i)

E(n) is arbitrary.

The size of E must verify: size( E ) = size( MAT, 1 ) = size( MAT, 2 ) = n .

STORE_Q (INPUT) logical(lgl)

On exit:

If STORE_Q = true : The leading n-by-n lower triangular part of MAT is overwritten by the matrix Q as a product of elementary Householder reflectors.

If STORE_Q = false : The m-by-n matrix MAT is destroyed.

See Further details.

Further Details¶

This subroutine is an implementation of the Ralha one-sided method to reduce implicitly a matrix product MAT’ * MAT to tridiagonal form T. Q is computed as a product of n-1 elementary reflectors (e.g., Householder transformations):

Q = H(1) * H(2) * … * H(n-1)

Each H(i) has the form:

H(i) = I + tau * v * v’

where tau is a real scalar, and v is a real vector. IF STORE_Q is set to true, the n-1 H(i) elementary reflectors are stored in the leading lower triangle of the array MAT. For the H(i) reflector, tau is stored in MAT(i,i) and v is stored in MATi+1:n,i). Note that if n is equal to 1, no elementary reflectors are needed.

This is the blocked version of the algorithm. See the references (1), (2) and (3) for further details. Furthermore the algorithm is parallelized if OPENMP is used.

The contents of MAT on exit are illustrated by the following example with n = 5 and m = 6:

if STORE_Q = true :

( t1 yy yy yy yy )

( v1 t2 yy yy yy )

( v1 v2 t3 yy yy )

( v1 v2 v3 t4 yy )

( v1 v2 v3 v4 xx )

( yy yy yy yy yy )

where vi and ti denote an element of the vector v and the scalar tau defining H(i), respectively. xx = machhuge and is used by the subroutine ORTHO_GEN_SYMTRID in order to verify that SYMTRID_CMP2 has been called before ORTHO_GEN_SYMTRID. Elements yy are destroyed by the subroutine.

Subroutines ORTHO_GEN_SYMTRID or APPLY_Q_SYMTRID, with logical argument UPPER set to false, can be used to generate the orthogonal matrix Q or to apply it to another matrix. See descriptions of these subroutines for further details.

For further details, see:

Hegland, M., Kahn, M., and Osborn, M., 1999:

A parallel algorithm for the reduction to tridiagonal form for eigendecomposition. SIAM Journal on Scientific Computing, 21:3, pp. 987-1005.

Ralha, R.M.S., 2003:

One-sided reduction to bidiagonal form. Linear Algebra Appl., No 358, pp. 219-238.

Barlow, J.L., Bosner, N., and Drmac, Z., 2005:

A new stable bidiagonal reduction algorithm. Linear Algebra Appl., No 397, pp. 35-84.

Bosner, N., and Barlow, J.L., 2007:

Block and Parallel versions of one-sided bidiagonalization. SIAM J. Matrix Anal. Appl., Volume 29, No 3, pp. 927-953.

`subroutine ortho_gen_symtrid ( mat, upper )`¶

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 with STORE_Q = true:

if UPPER = true, Q = H(n-1) * … * H(2) * H(1),

if UPPER = false, Q = H(1) * H(2) * … * H(n-1).

Arguments¶

MAT (INPUT/OUTPUT) real(stnd), dimension(:,:)

On entry, the vectors and the scalars, which define the elementary reflectors, as returned by SYMTRID_CMP with STORE_Q = true.

On exit, the n-by-n orthogonal matrix Q.

UPPER (INPUT) logical(lgl)

If:

UPPER= true : Upper triangle of MAT contains elementary reflectors from SYMTRID_CMP;

UPPER = false: Lower triangle of MAT contains elementary reflectors from SYMTRID_CMP.

Further Details¶

This subroutine is adapted from the routine SORGTR in LAPACK. A blocked algorithm is used for agregating the Householder transformations (e.g. the elementary reflectors) stored in the upper or lower triangle of MAT and generating the orthogonal matrix Q. Furthermore, the computations are parallelized if OPENMP is used.

For further details, see:

Golub, G.H., and Van Loan, C.F., 1996:

Matrix Computations. 3rd ed. The Johns Hopkins University Press, Baltimore.

Dongarra, J.J., Sorensen, D.C., and Hammarling, S.J., 1989:

Block reduction of matrices to condensed form for eigenvalue computations. J. of Computational and Applied Mathematics, Vol. 27, pp. 215-227.

Walker, H.F., 1988:

Implementation of the GMRES method using Householder transformations. Siam J. Sci. Stat. Comput., Vol. 9, No 1, pp. 152-163.

`subroutine ortho_gen_symtrid ( mat )`¶

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 with STORE_Q = true:

Q = H(n-1) * … * H(2) * H(1)

Arguments¶

MAT (INPUT/OUTPUT) real(stnd), dimension(:,:)

On entry, the vectors and the scalars, which define the elementary reflectors, as returned by SYMTRID_CMP with STORE_Q = true.

On exit, the n-by-n orthogonal matrix Q.

Further Details¶

This subroutine is adapted from the routine SORGTR in LAPACK. A blocked algorithm is used for agregating the Householder transformations (e.g. the elementary reflectors) stored in the upper triangle of MAT and generating the orthogonal matrix Q. Furthermore, the computations are parallelized if OPENMP is used.

For further details, see:

Golub, G.H., and Van Loan, C.F., 1996:

Matrix Computations. 3rd ed. The Johns Hopkins University Press, Baltimore.

Dongarra, J.J., Sorensen, D.C., and Hammarling, S.J., 1989:

Block reduction of matrices to condensed form for eigenvalue computations. J. of Computational and Applied Mathematics, Vol. 27, pp. 215-227.

Walker, H.F., 1988:

Implementation of the GMRES method using Householder transformations. Siam J. Sci. Stat. Comput., Vol. 9, No 1, pp. 152-163.

`subroutine apply_q_symtrid ( mat, c, left, trans, upper )`¶

Purpose¶

APPLY_Q_SYMTRID overwrites the general real m-by-n matrix C with:

Q * C if LEFT = true and TRANS = false, or

Q’* C if LEFT = true and TRANS = true, or

C * Q if LEFT = false and TRANS = false, or

C * Q’ 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

Q = H(1) * H(2) * … * H(nq-1), if UPPER = false

as returned by SYMTRID_CMP with STORE_Q = true.

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

Arguments¶

MAT (INPUT) real(stnd), dimension(:,:)

On entry, the vectors and the scalars, which define the elementary reflectors, as returned by SYMTRID_CMP with STORE_Q = true. MAT is not modified by the routine.

The shape of MAT must verify: size( MAT, 1 ) = size( MAT, 2 ) = nq.

C (INPUT/OUTPUT) real(stnd), dimension(:,:)

On entry, the m by n matrix C.

On exit, C is overwritten by Q * C or Q’ * C or C * Q’ or C * Q.

The shape of C must verify:

if LEFT = true, size( C, 1 ) = nq ;

if LEFT = false, size( C, 2 ) = nq .

LEFT (INPUT) logical(lgl)

If:

LEFT = true : apply Q or Q’ from the left ;

LEFT = false : apply Q or Q’ from the right .

TRANS (INPUT) logical(lgl)

If:

TRANS = false : apply Q (no transpose) ;

TRANS = true : apply Q’ (transpose) .

UPPER (INPUT) logical(lgl)

If:

UPPER = true : The upper triangle of MAT contains elementary reflectors generated by SYMTRID_CMP;

UPPER = false: The lower triangle of MAT contains elementary reflectors generated by SYMTRID_CMP.

Further Details¶

This subroutine is adapted from the routine SORMTR in LAPACK. A blocked algorithm is used to apply the Householder transformations (e.g. the elementary reflectors) stored in the upper or lower triangle of MAT (see the references (2) and (3) below).

Furthermore, the subroutine is parallelized if OPENMP is used.

For further details, see:

Golub, G.H., and Van Loan, C.F., 1996:

Matrix Computations. 3rd ed. The Johns Hopkins University Press, Baltimore.

Dongarra, J.J., Sorensen, D.C., and Hammarling, S.J., 1989:

Block reduction of matrices to condensed form for eigenvalue computations. J. of Computational and Applied Mathematics, Vol. 27, pp. 215-227.

Walker, H.F., 1988:

Implementation of the GMRES method using Householder transformations. Siam J. Sci. Stat. Comput., Vol. 9, No 1, pp. 152-163.

`subroutine apply_q_symtrid ( mat, c, left, trans )`¶

Purpose¶

APPLY_Q_SYMTRID overwrites the general real m-by-n matrix C with

Q * C if LEFT = true and TRANS = false, or

Q’* C if LEFT = true and TRANS = true, or

C * Q if LEFT = false and TRANS = false, or

C * Q’ 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)

as returned by SYMTRID_CMP (with UPPER = true or without this argument) and with STORE_Q = true.

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

Arguments¶

MAT (INPUT) real(stnd), dimension(:,:)

On entry, the vectors and the scalars, which define the elementary reflectors, as returned by SYMTRID_CMP (with UPPER = true or SYMTRID_CMP without argument UPPER) and with STORE_Q = true. MAT is not modified by the routine.

The shape of MAT must verify: size( MAT, 1 ) = size( MAT, 2 ) = nq.

C (INPUT/OUTPUT) real(stnd), dimension(:,:)

On entry, the m by n matrix C.

On exit, C is overwritten by Q * C or Q’ * C or C * Q’ or C * Q.

The shape of C must verify:

if LEFT = true, size( C, 1 ) = nq ;

if LEFT = false, size( C, 2 ) = nq .

LEFT (INPUT) logical(lgl)

If:

LEFT = true : apply Q or Q’ from the left ;

LEFT = false : apply Q or Q’ from the right .

TRANS (INPUT) logical(lgl)

If:

TRANS = false : apply Q (no transpose) ;

TRANS = true : apply Q’ (transpose) .

Further Details¶

This subroutine is adapted from the routine SORMTR in LAPACK. A blocked algorithm is used to apply the Householder transformations (e.g. the elementary reflectors) stored in the upper triangle of MAT (see the references (2) and (3) below).

Furthermore, the subroutine is parallelized if OPENMP is used.

For further details, see:

Golub, G.H., and Van Loan, C.F., 1996:

Matrix Computations. 3rd ed. The Johns Hopkins University Press, Baltimore.

Dongarra, J.J., Sorensen, D.C., and Hammarling, S.J., 1989:

Block reduction of matrices to condensed form for eigenvalue computations. J. of Computational and Applied Mathematics, Vol. 27, pp. 215-227.

Walker, H.F., 1988:

Implementation of the GMRES method using Householder transformations. Siam J. Sci. Stat. Comput., Vol. 9, No 1, pp. 152-163.

`subroutine apply_q_symtrid ( matp, c, left, trans, upper )`¶

Purpose¶

APPLY_Q_SYMTRID overwrites the general real m-by-n matrix C with

Q * C if LEFT = true and TRANS = false, or

Q’* C if LEFT = true and TRANS = true, or

C * Q if LEFT = false and TRANS = false, or

C * Q’ 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

Q = H(1) * H(2) * … * H(nq-1), if UPPER = false

as returned by SYMTRID_CMP with STORE_Q = true.

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

Arguments¶

MATP (INPUT) real(stnd), dimension(:)

On entry, the vectors and the scalars, which define the elementary reflectors, as returned by SYMTRID_CMP in argument MATP. MATP is not modified by the routine.

The size of MATP must verify size( MATP ) = (nq * (nq+1)/2)

C (INPUT/OUTPUT) real(stnd), dimension(:,:)

On entry, the m by n matrix C.

On exit, C is overwritten by Q * C or Q’ * C or C * Q’ or C * Q.

The shape of C must verify:

if LEFT = true, size( C, 1 ) = nq ;

if LEFT = false, size( C, 2 ) = nq .

LEFT (INPUT) logical(lgl)

If:

LEFT = true : apply Q or Q’ from the left ;

LEFT = false : apply Q or Q’ from the right .

TRANS (INPUT) logical(lgl)

If:

TRANS = false : apply Q (no transpose) ;

TRANS = true : apply Q’ (transpose) .

UPPER (INPUT) logical(lgl)

Specifies whether the upper or lower triangular part of the original symmetric matrix MAT was stored in packed form in MATP before the reduction by SYMTRID_CMP. If:

UPPER = true : Upper triangle of MAT was stored;

UPPER = false: Lower triangle of MAT was stored.

Further Details¶

This subroutine is adapted from the routine DORMTR in LAPACK and uses a blocked algorithm to apply the Householder transformations (e.g. the elementary reflectors) stored in packed form in the vector MATP.

Furthermore, the subroutine is parallelized if OPENMP is used.

For further details, see:

Golub, G.H., and Van Loan, C.F., 1996:

Matrix Computations. 3rd ed. The Johns Hopkins University Press, Baltimore.

Dongarra, J.J., Sorensen, D.C., and Hammarling, S.J., 1989:

Block reduction of matrices to condensed form for eigenvalue computations. J. of Computational and Applied Mathematics, Vol. 27, pp. 215-227.

Walker, H.F., 1988:

Implementation of the GMRES method using Householder transformations. Siam J. Sci. Stat. Comput., Vol. 9, No 1, pp. 152-163.

`subroutine apply_q_symtrid ( matp, c, left, trans )`¶

Purpose¶

APPLY_Q_SYMTRID overwrites the general real m-by-n matrix C with

Q * C if LEFT = true and TRANS = false, or

Q’* C if LEFT = true and TRANS = true, or

C * Q if LEFT = false and TRANS = false, or

C * Q’ 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)

as returned by SYMTRID_CMP (with UPPER = true or without this argument) and with STORE_Q = true.

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

Arguments¶

MATP (INPUT) real(stnd), dimension(:)

On entry, the vectors and the scalars, which define the elementary reflectors, as returned by SYMTRID_CMP in argument MATP. MATP is not modified by the routine.

The size of MATP must verify size( MATP ) = (nq * (nq+1)/2)

C (INPUT/OUTPUT) real(stnd), dimension(:,:)

On entry, the m by n matrix C.

On exit, C is overwritten by Q * C or Q’ * C or C * Q’ or C * Q.

The shape of C must verify:

if LEFT = true, size( C, 1 ) = nq ;

if LEFT = false, size( C, 2 ) = nq .

LEFT (INPUT) logical(lgl)

If:

LEFT = true : apply Q or Q’ from the left ;

LEFT = false : apply Q or Q’ from the right .

TRANS (INPUT) logical(lgl)

If:

TRANS = false : apply Q (no transpose) ;

TRANS = true : apply Q’ (transpose) .

Further Details¶

This subroutine is adapted from the routine DORMTR in LAPACK and uses a blocked algorithm to apply the Householder transformations (e.g. the elementary reflectors) stored in packed form in the vector MATP.

Furthermore, the subroutine is parallelized if OPENMP is used.

For further details, see:

Golub, G.H., and Van Loan, C.F., 1996:

Matrix Computations. 3rd ed. The Johns Hopkins University Press, Baltimore.

Dongarra, J.J., Sorensen, D.C., and Hammarling, S.J., 1989:

Block reduction of matrices to condensed form for eigenvalue computations. J. of Computational and Applied Mathematics, Vol. 27, pp. 215-227.

Walker, H.F., 1988:

Implementation of the GMRES method using Householder transformations. Siam J. Sci. Stat. Comput., Vol. 9, No 1, pp. 152-163.

`subroutine eig_sort ( sort, d, u )`¶

Purpose¶

Given the eigenvalues D and eigenvectors U as output from EIG_CMP, EIG_CMP2 EIG_CMP3 or SYMTRID_QRI, SYMTRID_QRI2 and SYMTRID_QRI3, this routine sorts the eigenvalues into ascending or descending order, and, rearranges the columns of U correspondingly.

Arguments¶

SORT (INPUT) character

Sort the eigenvalues into ascending order if SORT = ‘A’ or ‘a’, or in descending order if SORT = ‘D’ or ‘d’. The eigenvectors are reordered accordingly.

D (INPUT/OUTPUT) real(stnd), dimension(:)

On entry, the eigenvalues.

On exit, the eigenvalues in ascending or decreasing order.

U (INPUT/OUTPUT) real(stnd), dimension(:,:)

On entry, the columns of U are the eigenvectors.

On exit, U contains the reordered eigenvectors.

The shape of U must verify: size( U, 2 ) = size( D ) = m .

Further Details¶

The method is straight insertion.

`subroutine eig_abs_sort ( sort, d, u )`¶

Purpose¶

Given the eigenvalues D and eigenvectors U as output from EIG_CMP, EIG_CMP2 EIG_CMP3 or SYMTRID_QRI, SYMTRID_QRI2 and SYMTRID_QRI3, this routine sorts the eigenvalues into ascending or descending order of absolute magnitude, and, rearranges the columns of U correspondingly.

Arguments¶

SORT (INPUT) character

Sort the eigenvalues into ascending order if SORT = ‘A’ or ‘a’, or in descending order if SORT = ‘D’ or ‘d’ of absolute magnitude. The eigenvectors are reordered accordingly.

D (INPUT/OUTPUT) real(stnd), dimension(:)

On entry, the eigenvalues.

On exit, the eigenvalues in ascending or decreasing order of absolute magnitude.

U (INPUT/OUTPUT) real(stnd), dimension(:,:)

On entry, the columns of U are the eigenvectors.

On exit, U contains the reordered eigenvectors.

The shape of U must verify: size( U, 2 ) = size( D ) = m .

Further Details¶

The method is straight insertion.

`subroutine eigval_sort ( sort, d )`¶

Purpose¶

Given the eigenvalues D as output from EIGVAL_CMP, EIGVAL_CMP2, EIGVAL_CMP3 or SYMTRID_QRI, SYMTRID_QRI2 and SYMTRID_QRI3, this routine sorts the eigenvalues into ascending or descending order.

Arguments¶

SORT (INPUT) character

Sort the eigenvalues into ascending order if SORT = ‘A’ or ‘a’, or in descending order if SORT = ‘D’ or ‘d’.

D (INPUT/OUTPUT) real(stnd), dimension(:)

On entry, the eigenvalues.

On exit, the eigenvalues in ascending or decreasing order.

Further Details¶

The method is quick sort.

`subroutine eigval_abs_sort ( sort, d )`¶

Purpose¶

Given the eigenvalues D as output from EIGVAL_CMP, EIGVAL_CMP2, EIGVAL_CMP3 or SYMTRID_QRI, SYMTRID_QRI2 and SYMTRID_QRI3, this routine sorts the eigenvalues into ascending or descending order of absolute magnitude.

Arguments¶

SORT (INPUT) character

Sort the eigenvalues into ascending order if SORT = ‘A’ or ‘a’, or in descending order if SORT = ‘D’ or ‘d’ of absolute magnitude.

D (INPUT/OUTPUT) real(stnd), dimension(:)

On entry, the eigenvalues.

On exit, the eigenvalues in ascending or decreasing order of absolute magnitude.

Further Details¶

The method is straight insertion.

`subroutine symtrid_qri ( d, e, failure, mat, init_mat, sort, maxiter )`¶

Purpose¶

SYMTRID_QRI computes all eigenvalues and eigenvectors of a symmetric n-by-n tridiagonal matrix using the implicit QR method.

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 before calling SYMTRID_QRI.

Arguments¶

D (INPUT/OUTPUT) real(stnd), dimension(:)

On entry, the diagonal elements of the tridiagonal matrix.

On exit, the eigenvalues.

E (INPUT/OUTPUT) real(stnd), dimension(:)

On entry, the n-1 subdiagonal elements of the tridiagonal matrix. E(n) is arbitrary and is used as workspace.

On exit, E has been destroyed.

The size of E must verify: size( E ) = size( D ) = n .

FAILURE (OUTPUT) logical(lgl)

On exit:

FAILURE = false : indicates successful exit.

FAILURE = true : indicates that the algorithm did not converge and that full accuracy was not attained in the Schur decomposition of the tridiagonal matrix.

MAT (INPUT/OUTPUT) real(stnd), dimension(:,:)

On entry, if:

INIT_MAT is absent or if INIT_MAT = false, then MAT contains the orthogonal matrix used in the reduction to tridiagonal form as returned by ORTHO_GEN_SYMTRID.

INIT_MAT is present and INIT_MAT = true, MAT is set to the identity matrix of order n.

On exit, if FAILURE = false:

MAT contains the orthonormal eigenvectors of the original symmetric matrix if INIT_MAT is absent or if INIT_MAT = false,

MAT contains the orthonormal eigenvectors of the symmetric tridiagonal matrix if INIT_MAT is present and INIT_MAT = true.

The shape of MAT must verify: size( MAT, 1 ) = size( MAT, 2 ) = size( D ) = n .

INIT_MAT (INPUT, OPTIONAL) logical(lgl)

If:

INIT_MAT = false: The subroutine computes eigenvalues and eigenvectors of the original symmetric matrix. On entry, MAT must contain the orthogonal matrix used to reduce the original matrix to tridiagonal form.

INIT_MAT = true: The subroutine computes eigenvalues and eigenvectors of the tridiagonal matrix. MAT is initialized to the identity matrix of order n.

The default is false.

SORT (INPUT, OPTIONAL) character

Sort the eigenvalues into ascending order if SORT = ‘A’ or ‘a’, or in descending order if SORT = ‘D’ or ‘d’. The eigenvectors are reordered accordingly.

By default, the eigenvalues are not sorted.

MAXITER (INPUT, OPTIONAL) integer(i4b)

MAXITER controls the maximum number of QR sweeps in the Schur decomposition of the tridiagonal matrix. The algorithm fails to converge if the number of QR sweeps exceeds MAXITER * n. Convergence usually occurs in about 2 * n QR sweeps.

The default is 30.

Further Details¶

The eigenvalues and eigenvectors are computed by the implicit tridiagonal QR algorithm described in the reference (1) with modifications suggested in the reference (2).

This subroutine is adapted from the routine DSTEQR in LAPACK. Note that this subroutine is parallelized with OPENMP using the method described in the reference (3).

For further details, see:

Golub, G.H., and Van Loan, C.F., 1996:

Matrix Computations. 3rd ed. The Johns Hopkins University Press, Baltimore.

Greenbaum, A., and J. Dongarra, J., 1989:

Experiments with QR/QL Methods for the Symmetric Tridiagonal Eigenproblem. LAPACK Working Note No 17, November 1989.

Demmel, J., Heath, M.T., and Van Der Vorst, H., 1993:

Parallel numerical linear algebra. Acta Numerica 2, 111-197.

`subroutine symtrid_qri ( d, e, failure, sort, maxiter )`¶

Purpose¶

SYMTRID_QRI computes all eigenvalues of a symmetric n-by-n tridiagonal matrix using the Pal-Walker-Kahan variant of the QR algorithm.

The eigenvalues of a full symmetric matrix can also be found if SYMTRID_CMP has been used to reduce this matrix to tridiagonal form before calling SYMTRID_QRI.

Arguments¶

D (INPUT/OUTPUT) real(stnd), dimension(:)

On entry, the diagonal elements of the tridiagonal matrix.

On exit, the eigenvalues.

E (INPUT/OUTPUT) real(stnd), dimension(:)

On entry, the n-1 subdiagonal elements of the tridiagonal matrix. E(n) is arbitrary and is used as workspace.

On exit, E has been destroyed.

The size of E must verify: size( E ) = size( D ) = n .

FAILURE (OUTPUT) logical(lgl)

On exit:

FAILURE = false : this indicates successful exit.

FAILURE = true : this indicates that the algorithm did not converge and that full accuracy was not attained in the Schur decomposition of the tridiagonal matrix.

SORT (INPUT, OPTIONAL) character

Sort the eigenvalues into ascending order if SORT = ‘A’ or ‘a’, or in descending order if SORT = ‘D’ or ‘d’.

By default, the eigenvalues are not sorted.

MAXITER (INPUT, OPTIONAL) integer(i4b)

MAXITER controls the maximum number of QR sweeps in the Schur decomposition of the tridiagonal matrix. The algorithm fails to converge if the number of QR sweeps exceeds MAXITER * n. Convergence usually occurs in about 2 * n QR sweeps.

The default is 30.

Further Details¶

The eigenvalues are computed by the Pal-Walker-Kahan variant of the implicit tridiagonal QR algorithm described in the reference (1) with modifications suggested in the reference (2).

This subroutine is adapted from the routine DSTERF in LAPACK. This subroutine is not parallelized.

For further details, see:

Parlett, B.N., 1998:

The Symmetric Eigenvalue Problem. Revised edition, SIAM, Philadelphia.

Greenbaum, A., and J. Dongarra, J., 1989:

Experiments with QR/QL Methods for the Symmetric Tridiagonal Eigenproblem. LAPACK Working Note No 17, November 1989.

`subroutine symtrid_qri2 ( d, e, failure, mat, init_mat, sort, maxiter, max_francis_steps )`¶

Purpose¶

SYMTRID_QRI2 computes all eigenvalues and eigenvectors of a symmetric n-by-n tridiagonal matrix with a perfect shift strategy for the eigenvectors.

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 before calling SYMTRID_QRI2.

Arguments¶

D (INPUT/OUTPUT) real(stnd), dimension(:)

On entry, the diagonal elements of the tridiagonal matrix.

On exit, the eigenvalues.

E (INPUT/OUTPUT) real(stnd), dimension(:)

On entry, the n-1 subdiagonal elements of the tridiagonal matrix. E(n) is arbitrary and is used as workspace.

On exit, E has been destroyed.

The size of E must verify: size( E ) = size( D ) = n .

FAILURE (OUTPUT) logical(lgl)

On exit:

FAILURE = false : indicates successful exit.

FAILURE = true : indicates that the algorithm did not converge and that full accuracy was not attained in the Schur decomposition of the tridiagonal matrix.

MAT (INPUT/OUTPUT) real(stnd), dimension(:,:)

On entry, if:

INIT_MAT is absent or if INIT_MAT = false, then MAT contains the orthogonal matrix used in the reduction to tridiagonal form as returned by ORTHO_GEN_SYMTRID.

INIT_MAT is present and INIT_MAT = true, MAT is set to the identity matrix of order n.

On exit, if FAILURE = false:

MAT contains the orthonormal eigenvectors of the original symmetric matrix if INIT_MAT is absent or if INIT_MAT = false,

MAT contains the orthonormal eigenvectors of the symmetric tridiagonal matrix if INIT_MAT is present and INIT_MAT = true.

The shape of MAT must verify: size( MAT, 1 ) = size( MAT, 2 ) = size( D ) = n .

INIT_MAT (INPUT, OPTIONAL) logical(lgl)

If:

INIT_MAT = false: The subroutine computes eigenvalues and eigenvectors of the original symmetric matrix. On entry, MAT must contain the orthogonal matrix used to reduce the original matrix to tridiagonal form.

INIT_MAT = true: The subroutine computes eigenvalues and eigenvectors of the tridiagonal matrix. MAT is initialized to the identity matrix of order n.

The default is false.

SORT (INPUT, OPTIONAL) character

Sort the eigenvalues into ascending order if SORT = ‘A’ or ‘a’, or in descending order if SORT = ‘D’ or ‘d’. The eigenvectors are reordered accordingly.

By default, the eigenvalues are not sorted.

MAXITER (INPUT, OPTIONAL) integer(i4b)

MAXITER controls the maximum number of QR sweeps in the Schur decomposition of the tridiagonal matrix. The algorithm fails to converge if the number of QR sweeps exceeds MAXITER * n. Convergence usually occurs in about 2 * n QR sweeps.

The default is 30.

MAX_FRANCIS_STEPS (INPUT, OPTIONAL) integer(i4b)

MAX_FRANCIS_STEPS controls the maximum number of Francis sets (e.g. QR sweeps) of Givens rotations which must be saved before applying them with a wavefront algorithm to accumulate the eigenvectors in the QR algorithm. MAX_FRANCIS_STEPS is a strictly positive integer, otherwise the default value is used.

The default is the minimum of n and the integer parameter MAX_FRANCIS_STEPS_EIG specified in the module Select_Parameters.

Further Details¶

The eigenvalues are computed by the Pal-Walker-Kahan variant of the implicit tridiagonal QR algorithm described in the reference (1).

The eigenvectors are computed with a perfect shift strategy (see the references (1) and (2)) with modifications suggested in the references (3) and (4) for deflating a given eigenvalue from the tridiagonal matrix, and modifications described in references (5) and (6) for applying a set of Givens rotations.

Furthermore, the computation of the eigenvectors is parallelized if OPENMP is used.

With all these changes, SYMTRID_QRI2 is much more efficient than SYMTRID_QRI for computing the full set of eigenvectors of a real n-by-n symmetric tridiagonal matrix for large matrices.

For further details, see:

Parlett, B.N., 1998:

The Symmetric Eigenvalue Problem. Revised edition, SIAM, Philadelphia.

Greenbaum, A., and J. Dongarra, J., 1989:

Experiments with QR/QL Methods for the Symmetric Tridiagonal Eigenproblem. LAPACK Working Note No 17, November 1989.

Fernando, K.V., 1997:

On computing an eigenvector of a tridiagonal matrix. Part I: Basic results. Siam J. Matrix Anal. Appl., Vol. 18, 1013-1034.

Malyshev, A.N., 2000:

On deflation for symmetric tridiagonal matrices. Report 182 of the Department of Informatics, University of Bergen, Norway.

Lang, B., 1998:

Using level 3 BLAS in rotation-based algorithms. Siam J. Sci. Comput., Vol. 19, 626-634.

Van Zee, F.G., Van de Geijn, R., and Quintana-Orti, G., 2011:

Restructuring the QR Algorithm for High-Performance Application of Givens Rotations. FLAME Working Note 60. The University of Texas at Austin, Department of Computer Sciences. Technical Report TR-11-36.

`subroutine symtrid_qri2 ( d, e, failure, sort, maxiter )`¶

Purpose¶

SYMTRID_QRI2 computes all eigenvalues of a symmetric n-by-n tridiagonal matrix using the Pal-Walker-Kahan variant of the QR algorithm.

The eigenvalues of a full symmetric matrix can also be found if SYMTRID_CMP has been used to reduce this matrix to tridiagonal form before calling SYMTRID_QRI2.

Arguments¶

D (INPUT/OUTPUT) real(stnd), dimension(:)

On entry, the diagonal elements of the tridiagonal matrix.

On exit, the eigenvalues.

E (INPUT/OUTPUT) real(stnd), dimension(:)

On entry, the n-1 subdiagonal elements of the tridiagonal matrix. E(n) is arbitrary and is used as workspace.

On exit, E has been destroyed

The size of E must verify: size( E ) = size( D ) = n .

FAILURE (OUTPUT) logical(lgl)

On exit:

FAILURE = false : indicates successful exit.

FAILURE = true : indicates that the algorithm did not converge and that full accuracy was not attained in the Schur decomposition of the tridiagonal matrix.

SORT (INPUT, OPTIONAL) character

Sort the eigenvalues into ascending order if SORT = ‘A’ or ‘a’, or in descending order if SORT = ‘D’ or ‘d’.

By default, the eigenvalues are not sorted.

MAXITER (INPUT, OPTIONAL) integer(i4b)

MAXITER controls the maximum number of QR sweeps in the Schur decomposition of the tridiagonal matrix. The algorithm fails to converge if the number of QR sweeps exceeds MAXITER * n. Convergence usually occurs in about 2 * n QR sweeps.

The default is 30.

Further Details¶

The eigenvalues are computed by the Pal-Walker-Kahan variant of the implicit tridiagonal QR algorithm described in the reference (1).

This subroutine is adapted from the routine DSTERF in LAPACK. This subroutine is not parallelized.

For further details, see:

Parlett, B.N., 1998:

The Symmetric Eigenvalue Problem. revised edition, SIAM, Philadelphia.

`subroutine symtrid_qri3 ( d, e, failure, mat, init_mat, sort, maxiter, max_francis_steps )`¶

Purpose¶

SYMTRID_QRI3 computes all eigenvalues and eigenvectors of a symmetric n-by-n tridiagonal matrix using the implicit QR method.

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 before calling SYMTRID_QRI3.

Arguments¶

D (INPUT/OUTPUT) real(stnd), dimension(:)

On entry, the diagonal elements of the tridiagonal matrix.

On exit, the eigenvalues.

E (INPUT/OUTPUT) real(stnd), dimension(:)

On entry, the n-1 subdiagonal elements of the tridiagonal matrix. E(n) is arbitrary and is used as workspace.

On exit, E has been destroyed.

The size of E must verify: size( E ) = size( D ) = n .

FAILURE (OUTPUT) logical(lgl)

On exit:

FAILURE = false : indicates successful exit.

FAILURE = true : indicates that the algorithm did not converge and that full accuracy was not attained in the Schur decomposition of the tridiagonal matrix.

MAT (INPUT/OUTPUT) real(stnd), dimension(:,:)

On entry, if:

INIT_MAT is absent or if INIT_MAT = false, then MAT contains the orthogonal matrix used in the reduction to tridiagonal form as returned by ORTHO_GEN_SYMTRID.

INIT_MAT is present and INIT_MAT = true, MAT is set to the identity matrix of order n.

On exit, if FAILURE = false:

MAT contains the orthonormal eigenvectors of the original symmetric matrix if INIT_MAT is absent or if INIT_MAT = false,

MAT contains the orthonormal eigenvectors of the symmetric tridiagonal matrix if INIT_MAT is present and INIT_MAT = true.

The shape of MAT must verify: size( MAT, 1 ) = size( MAT, 2 ) = size( D ) = n .

INIT_MAT (INPUT, OPTIONAL) logical(lgl)

If:

INIT_MAT = false: The subroutine computes eigenvalues and eigenvectors of the original symmetric matrix. On entry, MAT must contain the orthogonal matrix used to reduce the original matrix to tridiagonal form.

INIT_MAT = true: The subroutine computes eigenvalues and eigenvectors of the tridiagonal matrix. MAT is initialized to the identity matrix of order n.

The default is false.

SORT (INPUT, OPTIONAL) character

Sort the eigenvalues into ascending order if SORT = ‘A’ or ‘a’, or in descending order if SORT = ‘D’ or ‘d’. The eigenvectors are reordered accordingly.

By default, the eigenvalues are not sorted.

MAXITER (INPUT, OPTIONAL) integer(i4b)

MAXITER controls the maximum number of QR sweeps in the Schur decomposition of the tridiagonal matrix. The algorithm fails to converge if the number of QR sweeps exceeds MAXITER * n. Convergence usually occurs in about 2 * n QR sweeps.

The default is 30.

MAX_FRANCIS_STEPS (INPUT, OPTIONAL) integer(i4b)

MAX_FRANCIS_STEPS controls the maximum number of Francis sets (e.g. QR sweeps) of Givens rotations which must be saved before applying them with a wavefront algorithm to accumulate the eigenvectors in the QR algorithm. MAX_FRANCIS_STEPS is a strictly positive integer, otherwise the default value is used.

The default is the minimum of n and the integer parameter MAX_FRANCIS_STEPS_EIG specified in the module Select_Parameters.

Further Details¶

The eigenvalues and eigenvectors are computed by the implicit tridiagonal QR algorithm described in the reference (1).

This subroutine is adapted from the routine DSTEQR in LAPACK with modifications suggested in the references (2) and (3) for the application of a set of Givens rotations.

Furthermore, the computations of the eigenvectors are parallelized if OPENMP is used.

SYMTRID_QRI3 is much more efficient than SYMTRID_QRI and only slightly less efficient than SYMTRID_QRI2 (but more robust) for computing the full set of eigenvectors of a real n-by-n symmetric tridiagonal matrix.

For further details, see:

Golub, G.H., and Van Loan, C.F., 1996:

Matrix Computations. 3rd ed. The Johns Hopkins University Press, Baltimore.

Lang, B., 1998:

Using level 3 BLAS in rotation-based algorithms. Siam J. Sci. Comput., Vol. 19, 626-634.

Van Zee, F.G., Van de Geijn, R., and Quintana-Orti, G., 2011:

Restructuring the QR Algorithm for High-Performance Application of Givens Rotations. FLAME Working Note 60. The University of Texas at Austin, Department of Computer Sciences. Technical Report TR-11-36.

`subroutine symtrid_qri3 ( d, e, failure, sort, maxiter )`¶

Purpose¶

SYMTRID_QRI3 computes all eigenvalues of a symmetric n-by-n tridiagonal matrix using the implicit QR method. The eigenvalues of a full symmetric matrix can also be found if SYMTRID_CMP has been used to reduce this matrix to tridiagonal form before calling SYMTRID_QRI3.

Arguments¶

D (INPUT/OUTPUT) real(stnd), dimension(:)

On entry, the diagonal elements of the tridiagonal matrix.

On exit, the eigenvalues.

E (INPUT/OUTPUT) real(stnd), dimension(:)

On entry, the n-1 subdiagonal elements of the tridiagonal matrix. E(n) is arbitrary and is used as workspace.

On exit, E has been destroyed.

The size of E must verify: size( E ) = size( D ) = n .

FAILURE (OUTPUT) logical(lgl)

On exit:

FAILURE = false : indicates successful exit.

FAILURE = true : indicates that the algorithm did not converge and that full accuracy was not attained in the Schur decomposition of the tridiagonal matrix.

SORT (INPUT, OPTIONAL) character

Sort the eigenvalues into ascending order if SORT = ‘A’ or ‘a’, or in descending order if SORT = ‘D’ or ‘d’.

By default, the eigenvalues are not sorted.

MAXITER (INPUT, OPTIONAL) integer(i4b)

MAXITER controls the maximum number of QR sweeps in the Schur decomposition of the tridiagonal matrix. The algorithm fails to converge if the number of QR sweeps exceeds MAXITER * n. Convergence usually occurs in about 2 * n QR sweeps.

The default is 30.

Further Details¶

The eigenvalues are computed by the implicit tridiagonal QR algorithm described in the reference (1).

This subroutine is adapted from the routine DSTEQR in LAPACK. This subroutine is not parallelized.

For further details, see:

Golub, G.H., and Van Loan, C.F., 1996:

Matrix Computations. 3rd ed. The Johns Hopkins University Press, Baltimore.

`subroutine lae2 ( a, b, c, rt1, rt2 )`¶

Purpose¶

LAE2 computes the eigenvalues of a 2-by-2 symmetric matrix

[ A B ]

[ B C ]

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

Arguments¶

A (INPUT) real(stnd)

The (1,1) element of the 2-by-2 matrix.

B (INPUT) real(stnd)

The (1,2) and (2,1) elements of the 2-by-2 matrix.

C (INPUT) real(stnd)

The (2,2) element of the 2-by-2 matrix.

RT1 (OUTPUT) real(stnd)

The eigenvalue of larger absolute value.

RT2 (OUTPUT) real(stnd)

The eigenvalue of smaller absolute value.

Further Details¶

RT1 is accurate to a few ulps barring over/underflow.

RT2 may be inaccurate if there is massive cancellation in the determinant A * C - B * B; higher precision or correctly rounded or correctly truncated arithmetic would be needed to compute RT2 accurately in all cases.

Overflow is possible only if RT1 is within a factor of 5 of overflow. Underflow is harmless if the input data is 0 or exceeds underflow_threshold / macheps.

This subroutine is translated from the routine DLAE2 in LAPACK.

`subroutine laev2 ( a, b, c, rt1, rt2, cs1, sn1 )`¶

Purpose¶

LAEV2 computes the eigendecomposition of a 2-by-2 symmetric matrix

[ A B ]

[ B C ]

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

[+CS1 +SN1 ] [ A B ] [ +CS1 -SN1 ] = [ RT1 000 ]

[-SN1 +CS1 ] [ B C ] [ +SN1 +CS1 ] _ [ 000 RT2 ].

Arguments¶

A (INPUT) real(stnd)

The (1,1) element of the 2-by-2 matrix.

B (INPUT) real(stnd)

The (1,2) and (2,1) elements of the 2-by-2 matrix.

C (INPUT) real(stnd)

The (2,2) element of the 2-by-2 matrix.

RT1 (OUTPUT) real(stnd)

The eigenvalue of larger absolute value.

RT2 (OUTPUT) real(stnd)

The eigenvalue of smaller absolute value.

CS1 (OUTPUT) real(stnd)

SN1 (OUTPUT) real(stnd)

The vector (CS1, SN1) is a unit right eigenvector for RT1.

Further Details¶

RT1 is accurate to a few ulps barring over/underflow.

RT2 may be inaccurate if there is massive cancellation in the determinant A * C - B * B; higher precision or correctly rounded or correctly truncated arithmetic would be needed to compute RT2 accurately in all cases.

CS1 and SN1 are accurate to a few ulps barring over/underflow.

Overflow is possible only if RT1 is within a factor of 5 of overflow. Underflow is harmless if the input data is 0 or exceeds underflow_threshold / macheps.

This subroutine is translated from the routine DLAE2 in LAPACK.

`subroutine symtrid_ratqri ( d, e, m, failure, small, tol )`¶

Purpose¶

SYMTRID_RATQRI computes the m largest or smallest eigenvalues of a symmetric n-by-n tridiagonal matrix using a rational QR method.

The m 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.

Arguments¶

D (INPUT/OUTPUT) real(stnd), dimension(:)

On entry, the diagonal elements of the tridiagonal matrix.

On exit, the computed eigenvalues replace the first m elements of D in decreasing sequence. The other elements are lost.

E (INPUT/OUTPUT) real(stnd), dimension(:)

On entry, the n-1 subdiagonal elements of the tridiagonal matrix. E(n) is arbitrary and is used as workspace.

On exit, E has been destroyed.

The size of E must verify: size( E ) = size( D ) = n .

M (INPUT) integer(i4b)

On entry, the number of smallest or largest eigenvalues wanted. M must be less than or equal to size( E ) = size( D ) = n.

FAILURE (OUTPUT) logical(lgl)

On exit:

FAILURE = false : indicates successful exit.

FAILURE = true : indicates that the algorithm did not converge and that full accuracy was not attained in the rational QR iterations for some eigenvalues.

SMALL (INPUT, OPTIONAL) logical(lgl)

On entry:

SMALL = false : indicates that the M largest eigenvalues are desired.

SMALL = true : indicates that the M smallest eigenvalues are desired.

The default is false.

TOL (INPUT, OPTIONAL) real(stnd)

On entry, TOL specifies a tolerance for the theoretical error of the computed eigenvalues. The theoretical error of the k-th eigenvalue is usually not greater than k * TOL.

The default is zero.

Further Details¶

This subroutine is not parallelized.

For further details, see:

Reinsch, C., and Bauer, F.L., 1968:

Rational QR transformation with Newton shift for symmetric tridiagonal matrices. Numerische Mathematik 11, 264-272.

`subroutine symtrid_ratqri2 ( d, e, val, failure, m, small, tol )`¶

Purpose¶

SYMTRID_RATQRI2 computes the largest or smallest eigenvalues of a symmetric n-by-n tridiagonal matrix whose sum of absolute values exceeds a given (positive) value. A rational QR method is used.

The largest or smallest eigenvalues of a full symmetric matrix whose sum of absolute values exceeds a given threshold can also be found, if SYMTRID_CMP has been used to reduce this matrix to tridiagonal form.

Arguments¶

D (INPUT/OUTPUT) real(stnd), dimension(:)

On entry, the diagonal elements of the tridiagonal matrix.

On exit, the computed eigenvalues replace the first M elements of D in decreasing sequence. The other elements are lost.

E (INPUT/OUTPUT) real(stnd), dimension(:)

On entry, the n-1 subdiagonal elements of the tridiagonal matrix. E(n) is arbitrary and is used as workspace.

On exit, E has been destroyed.

The size of E must verify: size( E ) = size( D ) = n .

VAL (INPUT) real(stnd)

On entry, the sum of the M eigenvalues (in absolute value) found will exceed abs(VAL) or M is equal to n.

FAILURE (OUTPUT) logical(lgl)

On exit:

FAILURE = false : indicates successful exit.

FAILURE = true : indicates that the algorithm did not converge and that full accuracy was not attained in the rational QR iterations for some eigenvalues.

M (OUTPUT) integer(i4b)

On exit, the number of eigenvalues found.

SMALL (INPUT, OPTIONAL) logical(lgl)

On entry:

SMALL = false : indicates that the M largest eigenvalues are desired.

SMALL = true : indicates that the M smallest eigenvalues are desired.

The default is false.

TOL (INPUT, OPTIONAL) real(stnd)

On entry, TOL specifies a tolerance for the theoretical error of the computed eigenvalues. The theoretical error of the k-th eigenvalue is usually not greater than k * TOL.

The default is zero.

Further Details¶

This subroutine is not parallelized.

For further details, see:

Reinsch, C., and Bauer, F.L., 1968:

Rational QR transformation with Newton shift for symmetric tridiagonal matrices. Numerische Mathematik 11, 264-272.

`subroutine symtrid_bisect ( d, e, neig, w, failure, small, sort, vector, abstol, le, theta, scaling, init )`¶

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.

The 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.

Arguments¶

D (INPUT) real(stnd), dimension(:)

On entry, D contains the diagonal elements of the tridiagonal matrix T.

E (INPUT) real(stnd), dimension(:)

On entry, E contains the n-1 off-diagonal elements of the tridiagonal matrix T whose eigenvalues are desired. E(n) is arbitrary.

The size of E must verify: size( E ) = size( D ) = n .

NEIG (OUTPUT) integer(i4b)

On output, NEIG specifies the number of eigenvalues which have been computed. Note that NEIG may be greater than the optional argument LE, if multiple eigenvalues at index LE make unique selection impossible.

If none of the optional arguments LE and THETA are used, NEIG is set to size(D) and all the eigenvalues of T are computed.

W (OUTPUT) real(stnd), dimension(:)

On exit, W(1:NEIG) contains the first NEIG largest (or smallest) eigenvalues of T. The other values in W (e.g. W(NEIG+1:size(D)) ) are flagged by a quiet NAN.

The size of W must verify: size( W ) = size( D ) = n .

FAILURE (OUTPUT) logical(lgl)

On exit:

FAILURE = false : indicates successful exit and the bisection algorithm converged for all the computed eigenvalues to the desired accuracy;

FAILURE = true : indicates that some or all of the eigenvalues failed to converge or were not computed. This is generally caused by unexpectedly inaccurate arithmetic.

SMALL (INPUT, OPTIONAL) logical(lgl)

On entry:

SMALL = false : indicates that the largest eigenvalues are desired;

SMALL = true : indicates that the smallest eigenvalues are desired.

The default is false.

SORT (INPUT, OPTIONAL) character

Sort the eigenvalues into ascending order if SORT = ‘A’ or ‘a’, or in descending order if SORT = ‘D’ or ‘d’. For other values of SORT nothing is done and W(:NEIG) may not be sorted.

VECTOR (INPUT, OPTIONAL) logical(lgl)

On entry, if VECTOR is set to true, a vectorized version of the bisection algorithm is used.

The default is VECTOR=false.

ABSTOL (INPUT, OPTIONAL) real(stnd)

On entry, the absolute tolerance for the eigenvalues. An eigenvalue (or cluster) is considered to be located if it has been determined to lie in an interval whose width is ABSTOL or less.

If ABSTOL is less than or equal to zero, or is not specified, then ULP * | T | will be used, where | T | means the 1-norm of T and ULP is the machine precision (distance from 1 to the next larger floating point number).

Eigenvalues will be computed most accurately when ABSTOL is set to the square root of the underflow threshold, sqrt(LAMCH(‘S’)), not zero.

LE (INPUT, OPTIONAL) integer(i4b)

On entry, LE specifies the number of eigenvalues which must be computed by the subroutine. On output, NEIG may be different than LE if multiple eigenvalues at index LE make unique selection impossible.

If:

SMALL=false, the subroutine computes the LE largest eigenvalues of T ;

SMALL=true, the subroutine computes the LE smallest eigenvalues of T .

Only one of the optional arguments LE and THETA must be specified, otherwise the subroutine will stop with an error message.

LE must be greater than 0 and less or equal to size( D ) .

The default is LE = size( D ).

THETA (INPUT, OPTIONAL) real(stnd)

On entry:

if SMALL=false, THETA specifies that the eigenvalues which are greater or equal to THETA must be computed. If none of the eigenvalues are greater or equal to THETA, NEIG is set to zero and W(:) to a quiet NAN;

if SMALL=true, THETA specifies that the eigenvalues which are less or equal to THETA must be computed. If none of the eigenvalues are smaller or equal to THETA, NEIG is set to zero and W(:) to a quiet NAN.

Only one of the optional arguments LE and THETA must be specified, otherwise the subroutine will stop with an error message.

SCALING (INPUT, OPTIONAL) logical(lgl)

On entry, if SCALING=true the tridiagonal matrix T is scaled before computing the eigenvalues.

The default is to scale the tridiagonal matrix.

INIT (INPUT, OPTIONAL) logical(lgl)

On entry, if INIT=true the initial intervals for the bisection steps are computed from estimates of the eigenvalues of the tridiagonal matrix obtained from the Pal-Walker-Kahan algorithm.

The default is not to use the Pal-Walker-Kahan algorithm.

Further Details¶

Let S(i), i=1,…,N=size( D ), be the N eigenvalues of the symmetric tridiagonal matrix T in decreasing order of magnitude. SYMTRID_BISECT then computes the LE largest or smallest eigenvalues ( or the eigenvalues which are greater/smaller or equal to THETA) of T by a bisection method (see the reference (1) below, Sec.8.5 ).

This subroutine is parallelized if OPENMP is used.

For further details, see:

Golub, G.H., and Van Loan, C.F., 1996:

Matrix Computations. 3rd ed. The Johns Hopkins University Press, Baltimore.

`subroutine dflgen ( d, e, lambda, cs, sn )`¶

Purpose¶

DFLGEN computes deflation parameters (e.g. a chain of Givens rotations) for a symmetric unreduced n-by-n tridiagonal matrix T and a given eigenvalue 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 an eigenvalue LAMBDA.

Arguments¶

D (INPUT) real(stnd), dimension(:)

On entry, the diagonal elements of the tridiagonal matrix.

E (INPUT) real(stnd), dimension(:)

On entry, the n-1 subdiagonal elements of the tridiagonal matrix.

The size of E must verify: size( E ) = size( D ) - 1.

LAMBDA (INPUT) real(stnd)

On entry, an eigenvalue of the tridiagonal matrix.

CS (OUTPUT) real(stnd), dimension(:)

On exit, the vector of the cosines coefficients of the chain of n-1 Givens rotations that deflates the symmetric tridiagonal matrix.

The size of CS must verify: size( CS ) = size( E ) = size( D ) - 1.

SN (OUTPUT) real(stnd), dimension(:)

On exit, the vector of the sines coefficients of the chain of n-1 Givens rotations that deflates the symmetric tridiagonal matrix.

The size of SN must verify: size( SN ) = size( E ) = size( D ) - 1.

Further Details¶

This subroutine is adapted from the matlab routine DFLGEN given in the reference (1). No check is done in the subroutine to verify that the input tridiagonal matrix is unreduced.

For further details, see:

Malyshev, A.N., 2000:

On deflation for symmetric tridiagonal matrices. Report 182 of the Department of Informatics, University of Bergen, Norway.

`subroutine dflgen2 ( d, e, lambda, cs, sn, deflate )`¶

Purpose¶

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

On output:

the arguments D and E contain, respectively, the new main diagonal and subdiagonal 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.

Arguments¶

D (INPUT/OUTPUT) real(stnd), dimension(:)

On entry, the diagonal elements of the tridiagonal matrix.

On exit, the new main diagonal of the symmetric tridiagonal matrix if DEFLATE=true. Otherwise, D is not changed.

E (INPUT/OUTPUT) real(stnd), dimension(:)

On entry, the n-1 subdiagonal elements of the tridiagonal matrix.

On exit, the new subdiagonal of the symmetric tridiagonal matrix if DEFLATE=true. Otherwise, E is not changed.

The size of E must verify: size( E ) = size( D ) - 1.

LAMBDA (INPUT) real(stnd)

On entry, an eigenvalue of the tridiagonal matrix.

CS (OUTPUT) real(stnd), dimension(:)

On exit, the vector of the cosines coefficients of the chain of n-1 Givens rotations that deflates the symmetric tridiagonal matrix.

The size of CS must verify: size( CS ) = size( E ) = size( D ) - 1.

SN (OUTPUT) real(stnd), dimension(:)

On exit, the vector of the sines coefficients of the chain of n-1 Givens rotations that deflates the symmetric tridiagonal matrix.

The size of SN must verify: size( SN ) = size( E ) = size( D ) - 1.

DEFLATE (OUTPUT) logical(lgl)

On exit:

DEFLATE = true : indicates successful exit.

DEFLATE = false: indicates that full accuracy was not attained in the deflation of the tridiagonal matrix.

Further Details¶

This subroutine is adapted from the matlab routine DFLGEN given in the reference (1). No check is done in the subroutine to verify that the input tridiagonal matrix is unreduced.

For further details, see:

Malyshev, A.N., 2000:

On deflation for symmetric tridiagonal matrices. Report 182 of the Department of Informatics, University of Bergen, Norway.

`subroutine dflapp ( d, e, cs, sn, deflate )`¶

Purpose¶

DFLAPP deflates a real symmetric n-by-n 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 subdiagonal of the deflated symmetric tridiagonal matrix if DEFLATE is set to true.

Arguments¶

D (INPUT/OUTPUT) real(stnd), dimension(:)

On entry, the diagonal elements of the tridiagonal matrix.

On exit, the new main diagonal of the symmetric tridiagonal matrix if DEFLATE=true. Otherwise, D is not changed.

E (INPUT/OUTPUT) real(stnd), dimension(:)

On entry, the n-1 subdiagonal elements of the tridiagonal matrix.

On exit, the new subdiagonal of the symmetric tridiagonal matrix if DEFLATE=true. Otherwise, E is not changed.

The size of E must verify: size( E ) = size( D ) - 1.

CS (INPUT) real(stnd), dimension(:)

On entry, the vector of the cosines coefficients of the chain of n-1 Givens rotations that deflates the symmetric tridiagonal matrix as computed by DFLGEN subroutine.

The size of CS must verify: size( CS ) = size( E ) = size( D ) - 1.

SN (INPUT) real(stnd), dimension(:)

On entry, the vector of the sines coefficients of the chain of n-1 Givens rotations that deflates the symmetric tridiagonal matrix as computed by DFLGEN subroutine.

The size of SN must verify: size( SN ) = size( E ) = size( D ) - 1.

DEFLATE (OUTPUT) logical(lgl)

On exit:

DEFLATE = true : indicates successful exit.

DEFLATE = false: indicates that full accuracy was not attained in the deflation of the tridiagonal matrix.

Further Details¶

This subroutine is adapted from the matlab routine DFLAPP given in the reference (1).

For further details, see:

Malyshev, A.N., 2000:

On deflation for symmetric tridiagonal matrices. Report 182 of the Department of Informatics, University of Bergen, Norway.

`subroutine qrstep ( d, e, lambda, cs, sn, deflate )`¶

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 subdiagonal 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.

Arguments¶

D (INPUT/OUTPUT) real(stnd), dimension(:)

On entry, the diagonal elements of the tridiagonal matrix.

On exit, the new main diagonal of the symmetric tridiagonal matrix after the QR step.

E (INPUT/OUTPUT) real(stnd), dimension(:)

On entry, the n-1 subdiagonal elements of the tridiagonal matrix.

On exit, the new subdiagonal of the symmetric tridiagonal matrix after the QR step.

The size of E must verify: size( E ) = size( D ) - 1.

LAMBDA (INPUT) real(stnd)

On entry, the shift used in the current QR step.

CS (OUTPUT) real(stnd), dimension(:)

On exit, the vector of the cosines coefficients of the chain of n-1 Givens rotations for the current QR step.

The size of CS must verify: size( CS ) = size( E ) = size( D ) - 1.

SN (OUTPUT) real(stnd), dimension(:)

On exit, the vector of the sines coefficients of the chain of n-1 Givens rotations for the current QR step.

The size of SN must verify: size( SN ) = size( E ) = size( D ) - 1.

DEFLATE (OUTPUT) logical(lgl)

On exit:

DEFLATE = true : indicates that deflation occured at the end of the step.

DEFLATE = false: indicates that the last subdiagonal element of the tridiagonal matrix is not small.

Further Details¶

This subroutine is adapted from the matlab routine QRSTEP given in the reference (1). No check is done in the subroutine to verify that the input tridiagonal matrix is unreduced.

For further details, see:

Mastronardi, M., Van Barel, M., Van Camp, E., and Vandebril, R., 2006:

On computing the eigenvectors of a class of structured matrices. Journal of Computational and Applied Mathematics, 189, 580-591.

`subroutine prodgiv ( cs, sn, x )`¶

Purpose¶

PRODGIV applies a chain of n-1 planar rotations produced by DFLGEN, DFLGEN2 or QRSTEP to a vector of length n.

Arguments¶

X (INPUT/OUTPUT) real(stnd), dimension(:)

On entry, the input vector of length n.

On exit, the product of the chain of the n-1 planar rotations by the input vector.

CS (INPUT) real(stnd), dimension(:)

On entry, the vector of the cosines coefficients of the chain of n-1 Givens rotations as computed by DFLGEN or DFLGEN2 subroutines.

The size of CS must verify: size( CS ) = size( X ) - 1.

SN (INPUT) real(stnd), dimension(:)

On entry, the vector of the sines coefficients of the chain of n-1 Givens rotations as computed by DFLGEN or DFLGEN2 subroutines.

The size of SN must verify: size( SN ) = size( CS ) = size( X ) - 1.

Further Details¶

This subroutine is adapted from the matlab routine PRODGIV given in the reference (1).

For further details, see:

Mastronardi, M., Van Barel, M., Van Camp, E., and Vandebril, R., 2006:

On computing the eigenvectors of a class of structured matrices. Journal of Computational and Applied Mathematics, 189, 580-591.

`function prodgiv_eigvec ( cs, sn )`¶

Purpose¶

PRODGIV_EIGVEC computes an eigenvector of a n-by-n symmetric tridiagonal matrix from a chain of n-1 planar rotations produced by DFLGEN, DFLGEN2 or QRSTEP.

Arguments¶

CS (INPUT) real(stnd), dimension(:)

On entry, the vector of the cosines coefficients of the chain of n-1 Givens rotations as computed by DFLGEN or DFLGEN2 subroutines.

SN (INPUT) real(stnd), dimension(:)

On entry, the vector of the sines coefficients of the chain of n-1 Givens rotations as computed by DFLGEN or DFLGEN2 subroutines.

The size of SN must verify: size( SN ) = size( CS ) = n - 1.

Further Details¶

This subroutine is adapted from the matlab routine PRODGIV given in the reference (1).

For further details, see:

Mastronardi, M., Van Barel, M., Van Camp, E., and Vandebril, R., 2006:

On computing the eigenvectors of a class of structured matrices. Journal of Computational and Applied Mathematics, 189, 580-591.

`subroutine symtrid_deflate ( d, e, eigval, eigvec, failure, max_qr_steps )`¶

Purpose¶

SYMTRID_DEFLATE computes an eigenvector of a real symmetric tridiagonal matrix T corresponding to a specified eigenvalue, using a deflation technique.

Arguments¶

D (INPUT) real(stnd), dimension(:)

On entry, the diagonal elements of the symmetric tridiagonal matrix T.

E (INPUT) real(stnd), dimension(:)

On entry, the n-1 subdiagonal elements of the symmetric tridiagonal matrix T.

The size of E must verify: size( E ) = size( D ) - 1 = n - 1.

EIGVAL (INPUT) real(stnd)

On entry, an eigenvalue of the symmetric tridiagonal matrix.

EIGVEC (OUTPUT) real(stnd), dimension(:)

On exit, the computed eigenvector associated with the eigenvalue EIGVAL.

The sise of EIGVEC must verify: size( EIGVEC ) = size( D ) = n .

FAILURE (OUTPUT) logical(lgl)

On exit:

FAILURE = false : indicates successful exit.

FAILURE = true : indicates that the algorithm did not converge and that full accuracy was not attained in the deflation procedure of the tridiagonal matrix.

MAX_QR_STEPS (INPUT, OPTIONAL) integer(i4b)

MAX_QR_STEPS controls the maximum number of QR sweeps for deflating the tridiagonal matrix for a given eigenvalue. The algorithm fails to converge if the total number of QR sweeps exceeds MAX_QR_STEPS.

The default is 4.

Further Details¶

SYMTRID_DEFLATE may fail for some zero-diagonal tridiagonal matrices.

For further details, see:

Malyshev, A.N., 2000:

On deflation for symmetric tridiagonal matrices. Report 182 of the Department of Informatics, University of Bergen, Norway.

Mastronardi, M., Van Barel, M., Van Camp, E., and Vandebril, R., 2006:

On computing the eigenvectors of a class of structured matrices. Journal of Computational and Applied Mathematics, 189, 580-591.

`subroutine symtrid_deflate ( d, e, eigval, eigvec, failure, max_qr_steps )`¶

Purpose¶

SYMTRID_DEFLATE computes eigenvectors of a real symmetric tridiagonal matrix T corresponding to specified eigenvalues, using a deflation technique.

Arguments¶

D (INPUT) real(stnd), dimension(:)

On entry, the diagonal elements of the symmetric tridiagonal matrix T.

E (INPUT) real(stnd), dimension(:)

On entry, the n-1 subdiagonal elements of the symmetric tridiagonal matrix T.

The size of E must verify: size( E ) = size( D ) - 1 = n - 1.

EIGVAL (INPUT) real(stnd), dimension(:)

On entry, selected eigenvalues of the symmetric tridiagonal matrix. The eigenvalues can be given in any order.

The size of EIGVAL must verify: size( EIGVAL ) <= size( D ) = n .

EIGVEC (OUTPUT) real(stnd), dimension(:,:)

On exit, the computed eigenvectors. The eigenvector associated with the eigenvalue EIGVAL(j) is stored in the j-th column of EIGVEC.

The shape of EIGVEC must verify:

size( EIGVEC, 1 ) = size( D ) = n ;

size( EIGVEC, 2 ) = size( EIGVAL ).

FAILURE (OUTPUT) logical(lgl), dimension(:)

On exit:

FAILURE(j) = false : indicates successful exit for the jth eigenvector.

FAILURE(j) = true : indicates that the algorithm did not converge and full accuracy was not attained in the deflation procedure of the tridiagonal matrix for the jth eigenvector.

The size of FAILURE must verify: size( FAILURE ) = size( EIGVAL ) .

MAX_QR_STEPS (INPUT, OPTIONAL) integer(i4b)

MAX_QR_STEPS controls the maximum number of QR sweeps for deflating the tridiagonal matrix for a given eigenvalue. The algorithm fails to converge if the total number of QR sweeps for all eigenvalues exceeds MAX_QR_STEPS * size(EIGVAL).

The default is 4.

Further Details¶

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.

For further details, see:

Malyshev, A.N., 2000:

On deflation for symmetric tridiagonal matrices. Report 182 of the Department of Informatics, University of Bergen, Norway.

Mastronardi, M., Van Barel, M., Van Camp, E., and Vandebril, R., 2006:

On computing the eigenvectors of a class of structured matrices. Journal of Computational and Applied Mathematics, 189, 580-591.

`subroutine trid_deflate ( d, e, eigval, eigvec, failure, max_qr_steps, ortho, inviter )`¶

Purpose¶

TRID_DEFLATE computes the eigenvectors of a real symmetric tridiagonal matrix T corresponding to specified eigenvalues, using deflation techniques.

Arguments¶

D (INPUT) real(stnd), dimension(:)

On entry, the diagonal elements of the symmetric tridiagonal matrix T.

E (INPUT) real(stnd), dimension(:)

On entry, the n-1 subdiagonal elements of the symmetric tridiagonal matrix T, E(n) is arbitrary and is not used .

The size of E must verify: size( E ) = size( D ) = n.

EIGVAL (INPUT) real(stnd), dimension(:)

On entry, selected eigenvalues of the symmetric tridiagonal matrix. The eigenvalues must be given in decreasing order.

The size of EIGVAL must verify: size( EIGVAL ) <= size( D ) = n .

EIGVEC (OUTPUT) real(stnd), dimension(:,:)

On exit, the computed eigenvectors. The eigenvector associated with the eigenvalue EIGVAL(j) is stored in the j-th column of EIGVEC.

The shape of EIGVEC must verify:

size( EIGVEC, 1 ) = size( D ) = n ;

size( EIGVEC, 2 ) = size( EIGVAL ) .

FAILURE (OUTPUT) logical(lgl)

On exit:

FAILURE = false : indicates successful exit.

FAILURE = true : indicates that the algorithm did not converge and that full accuracy was not attained in the deflation procedure of the tridiagonal matrix.

MAX_QR_STEPS (INPUT, OPTIONAL) integer(i4b)

MAX_QR_STEPS controls the maximum number of QR sweeps for deflating the tridiagonal matrix for a given eigenvalue. The algorithm fails to converge if the total number of QR sweeps for all eigenvalues exceeds MAX_QR_STEPS * size(EIGVAL).

The default is 4.

ORTHO (INPUT, OPTIONAL) logical(lgl)

On entry, if:

ORTHO=true, the tridiagonal matrix is deflated sequentially for all the specified eigenvalues; this implies that the eigenvectors will be automatically orthogonal on exit.

ORTHO=false, the tridiagonal matrix is deflated in parallel for the different clusters of eigenvalues or isolated eigenvalues; this implies that orthogonality is preserved inside each cluster, but not automatically between clusters.

The default is ORTHO=false.

INVITER (INPUT, OPTIONAL) logical(lgl)

On entry, if INVITER=true eigenvectors corresponding to isolated eigenvalues are computed by inverse iteration instead of deflation.

The default is INVITER=true.

Further Details¶

TRID_DEFLATE uses an efficient and robust 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 (see the references (1), (2) and (3) below). QR iterations are also used as a back-up procedure if the deflation technique fails (see the reference (4)).

Optionally, eigenvectors corresponding to isolated eigenvalues may be also computed by inverse iteration on the tridiagonal matrix T. This is the default for eigenvectors associated with isolated eigenvalues since in this case inverse iteration is safe and faster than the deflation algorithms.

The computation of the eigenvectors is parallelized if OPENMP is used.

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 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.

For further details, see:

Fernando, K.V., 1997:

On computing an eigenvector of a tridiagonal matrix. Part I: Basic results. Siam J. Matrix Anal. Appl., Vol. 18, pp. 1013-1034.

Parlett, B.N., and Dhillon, I.S., 1997:

Fernando’s solution to Wilkinsin’s problem: An application of double factorization. Linear Algebra and its Appl., 267, pp.247-279.

Malyshev, A.N., 2000:

On deflation for symmetric tridiagonal matrices. Report 182 of the Department of Informatics, University of Bergen, Norway.

Mastronardi, M., Van Barel, M., Van Camp, E., and Vandebril, R., 2006:

On computing the eigenvectors of a class of structured matrices. Journal of Computational and Applied Mathematics, 189, 580-591.

`subroutine trid_deflate ( d, e, eigval, eigvec, failure, mat, max_qr_steps, ortho, inviter )`¶

Purpose¶

TRID_DEFLATE computes the eigenvectors of a full real n-by-n symmetric matrix MAT corresponding to specified eigenvalues, using deflation techniques applied to a symmetric tridiagonal matrix T followed by a back-transformation procedure.

It is required that the original symmetric matrix MAT has been reduced to symmetric tridiagonal form T by an orthogonal similarity transformation:

Q’ * MAT * Q = T

with a call to SYMTRID_CMP with parameter STORE_Q set to true, before calling TRID_DEFLATE.

Arguments¶

D (INPUT) real(stnd), dimension(:)

On entry, the diagonal elements of the symmetric tridiagonal form T of MAT.

E (INPUT) real(stnd), dimension(:)

On entry, the n-1 subdiagonal elements of the symmetric tridiagonal form T of MAT. E(n) is arbitrary .

The size of E must verify: size( E ) = size( D ) = n .

EIGVAL (INPUT) real(stnd), dimension(:)

On entry, selected eigenvalues of the symmetric matrix MAT. The eigenvalues must be given in decreasing order.

The size of EIGVAL must verify: size( EIGVAL ) <= size( D ) = n .

EIGVEC (OUTPUT) real(stnd), dimension(:,:)

On exit, the computed eigenvectors. The eigenvector associated with the eigenvalue EIGVAL(j) is stored in the j-th column of EIGVEC.

The shape of EIGVEC must verify:

size( EIGVEC, 1 ) = size( D ) = n ;

size( EIGVEC, 2 ) = size( EIGVAL ).

FAILURE (OUTPUT) logical(lgl)

On exit:

FAILURE = false : indicates successful exit.

FAILURE = true : indicates that the algorithm did not converge and that full accuracy was not attained in the deflation procedure of the tridiagonal matrix.

MAT (INPUT) real(stnd), dimension(:,:)

On entry, the vectors and the scalars which define the elementary reflectors used to reduce the full real n-by-n symmetric matrix MAT to symmetric tridiagonal form T, as returned by SYMTRID_CMP with STORE_Q=true, in its argument MAT. MAT is not modified by the routine.

Back-transformation is used to find the selected eigenvectors of the original matrix MAT and these eigenvectors are stored in argument EIGVEC.

The shape of MAT must verify:

size( MAT, 1 ) = size( D ) = n ;

size( MAT, 2 ) = size( D ) = n .

MAX_QR_STEPS (INPUT, OPTIONAL) integer(i4b)

MAX_QR_STEPS controls the maximum number of QR sweeps for deflating the tridiagonal matrix for a given eigenvalue. The algorithm fails to converge if the total number of QR sweeps for all eigenvalues exceeds MAX_QR_STEPS * size(EIGVAL).

The default is 4.

ORTHO (INPUT, OPTIONAL) logical(lgl)

On entry, if:

ORTHO=true, the tridiagonal matrix is deflated sequentially for all the specified eigenvalues; this implies that the eigenvectors will be automatically orthogonal on exit.

ORTHO=false, the tridiagonal matrix is deflated in parallel for the different clusters of eigenvalues or isolated eigenvalues; this implies that orthogonality is preserved inside each cluster, but not automatically between clusters.

The default is ORTHO=false.

INVITER (INPUT, OPTIONAL) logical(lgl)

On entry, if INVITER=true eigenvectors corresponding to isolated eigenvalues are computed by inverse iteration instead of deflation.

The default is INVITER=true.

Further Details¶

TRID_DEFLATE uses an efficient and robust 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 (see the references (1), (2) and (3) below). QR iterations are also used as a back-up procedure if the deflation technique fails (see the reference (4)).

Optionally, eigenvectors corresponding to isolated eigenvalues may be also computed by inverse iteration on the tridiagonal matrix T. This is the default for eigenvectors associated with isolated eigenvalues since in this case inverse iteration is safe and faster than the deflation algorithms.

In a second step, the corresponding (selected) eigenvectors of the full real n-by-n symmetric matrix MAT are computed by a blocked back-transformation algorithm using the Householder transformations used to reduce the full real n-by-n symmetric matrix MAT to symmetric tridiagonal form T (see the references (5) and (6)).

The computation of the eigenvectors is parallelized if OPENMP is used.

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 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.

For further details on the deflation algorithm and the blocked backed-transformation algorithm, see:

Fernando, K.V., 1997:

On computing an eigenvector of a tridiagonal matrix. Part I: Basic results. Siam J. Matrix Anal. Appl., Vol. 18, 1013-1034.

Parlett, B.N., and Dhillon, I.S., 1997:

Fernando’s solution to Wilkinsin’s problem: An application of double factorization. Linear Algebra and its Appl., 267, pp.247-279.

Malyshev, A.N., 2000:

On deflation for symmetric tridiagonal matrices. Report 182 of the Department of Informatics, University of Bergen, Norway.

Mastronardi, M., Van Barel, M., Van Camp, E., and Vandebril, R., 2006:

On computing the eigenvectors of a class of structured matrices. Journal of Computational and Applied Mathematics, 189, 580-591.

Dongarra, J.J., Sorensen, D.C., and Hammarling, S.J., 1989:

Block reduction of matrices to condensed form for eigenvalue computations. J. of Computational and Applied Mathematics, Vol. 27, pp. 215-227.

Walker, H.F., 1988:

Implementation of the GMRES method using Householder transformations. Siam J. Sci. Stat. Comput., Vol. 9, No 1, pp. 152-163.

`subroutine trid_deflate ( d, e, eigval, eigvec, failure, matp, max_qr_steps, ortho, inviter )`¶

Purpose¶

TRID_DEFLATE computes the eigenvectors of a full real n-by-n symmetric matrix MAT, packed column-wise in a linear array MATP, corresponding to specified eigenvalues, using deflation techniques applied to a symmetric tridiagonal matrix T followed by a back-transformation procedure.

It is required that the original packed symmetric matrix mat has been reduced to symmetric tridiagonal form T by an orthogonal similarity transformation:

Q’ * MAT * Q = T

with a call to SYMTRID_CMP with parameter STORE_Q set to true, before calling TRID_DEFLATE.

Arguments¶

D (INPUT) real(stnd), dimension(:)

On entry, the diagonal elements of the symmetric tridiagonal form T of MAT.

E (INPUT) real(stnd), dimension(:)

On entry, the n-1 subdiagonal elements of the symmetric tridiagonal form T of MAT. E(n) is arbitrary .

The size of E must verify: size( E ) = size( D ) = n .

EIGVAL (INPUT) real(stnd), dimension(:)

On entry, selected eigenvalues of the symmetric matrix MAT. The eigenvalues must be given in decreasing order.

The size of EIGVAL must verify: size( EIGVAL ) <= size( D ) = n .

EIGVEC (OUTPUT) real(stnd), dimension(:,:)

On exit, the computed eigenvectors. The eigenvector associated with the eigenvalue EIGVAL(j) is stored in the j-th column of EIGVEC.

The shape of EIGVEC must verify:

size( EIGVEC, 1 ) = size( D ) = n ;

size( EIGVEC, 2 ) = size( EIGVAL ) .

FAILURE (OUTPUT) logical(lgl)

On exit:

FAILURE = false : indicates successful exit.

FAILURE = true : indicates that the algorithm did not converge and that full accuracy was not attained in the deflation procedure of the tridiagonal matrix.

MATP (INPUT) real(stnd), dimension(:)

On entry, the vectors and the scalars which define the elementary reflectors used to reduce the packed real n-by-n symmetric matrix MAT to symmetric tridiagonal form T, as returned by SYMTRID_CMP with STORE_Q=true, in its argument MATP. MATP is not modified by the routine.

Back-transformation is used to find the selected eigenvectors of the original matrix MAT and these eigenvectors are stored in argument EIGVEC.

The size of MATP must verify: size( MATP ) = (n * (n+1)/2)

MAX_QR_STEPS (INPUT, OPTIONAL) integer(i4b)

MAX_QR_STEPS controls the maximum number of QR sweeps for deflating the tridiagonal matrix for a given eigenvalue. The algorithm fails to converge if the total number of QR sweeps for all eigenvalues exceeds MAX_QR_STEPS * size(EIGVAL).

The default is 4.

ORTHO (INPUT, OPTIONAL) logical(lgl)

On entry, if:

ORTHO=true, the tridiagonal matrix is deflated sequentially for all the specified eigenvalues; this implies that the eigenvectors will be automatically orthogonal on exit.

ORTHO=false, the tridiagonal matrix is deflated in parallel for the different clusters of eigenvalues or isolated eigenvalues; this implies that orthogonality is preserved inside each cluster, but not automatically between clusters.

The default is ORTHO=false.

INVITER (INPUT, OPTIONAL) logical(lgl)

On entry, if INVITER=true eigenvectors corresponding to isolated eigenvalues are computed by inverse iteration instead of deflation.

The default is INVITER=true.

Further Details¶

TRID_DEFLATE uses an efficient and robust 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 (see the references (1), (2) and (3) below). QR iterations are also used as a back-up procedure if the deflation technique fails (see the reference (4)).

Optionally, eigenvectors corresponding to isolated eigenvalues may be also computed by inverse iteration on the tridiagonal matrix T. This is the default for eigenvectors associated with isolated eigenvalues since in this case inverse iteration is safe and faster than the deflation algorithms.

In a second step, the corresponding (selected) eigenvectors of the full real n-by-n symmetric matrix MAT are computed by a blocked back-transformation algorithm with the Householder transformations used to reduce the full real n-by-n symmetric matrix MAT to symmetric tridiagonal form T (see the references (5) and (6)). These Householder transformations must be packed in the linear array MATP (as returned by SYMTRID_CMP) on entry of TRID_DEFLATE.

The computation of the eigenvectors is parallelized if OPENMP is used.

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 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.

For further details, on the deflation algorithm and the blocked backed-transformation algorithm, see:

Fernando, K.V., 1997:

On computing an eigenvector of a tridiagonal matrix. Part I: Basic results. Siam J. Matrix Anal. Appl., Vol. 18, 1013-1034.

Parlett, B.N., and Dhillon, I.S., 1997:

Fernando’s solution to Wilkinsin’s problem: An application of double factorization. Linear Algebra and its Appl., 267, pp.247-279.

Malyshev, A.N., 2000:

On deflation for symmetric tridiagonal matrices. Report 182 of the Department of Informatics, University of Bergen, Norway.

Mastronardi, M., Van Barel, M., Van Camp, E., and Vandebril, R., 2006:

On computing the eigenvectors of a class of structured matrices. Journal of Computational and Applied Mathematics, 189, 580-591.

Dongarra, J.J., Sorensen, D.C., and Hammarling, S.J., 1989:

Block reduction of matrices to condensed form for eigenvalue computations. J. of Computational and Applied Mathematics, Vol. 27, pp. 215-227.

Walker, H.F., 1988:

Implementation of the GMRES method using Householder transformations. Siam J. Sci. Stat. Comput., Vol. 9, No 1, pp. 152-163.

`subroutine eig_cmp ( mat, eigval, failure, upper, sort, maxiter )`¶

Purpose¶

EIG_CMP computes all eigenvalues and eigenvectors of a n-by-n real symmetric matrix MAT.

Arguments¶

MAT (INPUT/OUTPUT) real(stnd), dimension(:,:)

On entry, the symmetric matrix MAT.

If:

UPPER = true: the leading n-by-n upper triangular part of MAT contains the upper triangular part of the matrix MAT.

UPPER = false: the leading n-by-n lower triangular part of MAT contains the lower triangular part of the matrix MAT.

On exit, MAT contains the orthonormal eigenvectors of the matrix MAT.

The shape of MAT must verify: size( MAT, 1 ) = size( MAT, 2 ) = n .

EIGVAL (OUTPUT) real(stnd), dimension(:)

On exit, the eigenvalues.

The size of EIGVAL must verify: size( EIGVAL ) = n .

FAILURE (OUTPUT) logical(lgl)

On exit:

FAILURE = false : indicates successful exit.

FAILURE = true : indicates that the algorithm did not converge and that full accuracy was not attained in the Schur decomposition of an intermediate tridiagonal form T of MAT .

UPPER (INPUT) logical(lgl)

Specifies whether the upper or lower triangular part of the symmetric matrix MAT is stored. If:

UPPER= true : Upper triangular is stored ;

UPPER= false: Lower triangular is stored .

SORT (INPUT, OPTIONAL) character

Sort the eigenvalues into ascending order if SORT = ‘A’ or ‘a’, or in descending order if SORT = ‘D’ or ‘d’. The eigenvectors are reordered accordingly.

MAXITER (INPUT, OPTIONAL) integer(i4b)

MAXITER controls the maximum number of QR sweeps in the Schur decomposition of an intermediate tridiagonal form T of MAT.

The algorithm fails to converge if the number of QR sweeps exceeds MAXITER * size(EIGVAL). Convergence usually occurs in about 2 * size(EIGVAL) QR sweeps.

The default is 30.

Further Details¶

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

For further details, see:

Parlett, B.N., 1998:

The Symmetric Eigenvalue Problem. Revised edition, SIAM, Philadelphia.

`subroutine eig_cmp ( mat, eigval, failure, sort, maxiter )`¶

Purpose¶

EIG_CMP computes all eigenvalues and eigenvectors of a n-by-n real symmetric matrix MAT.

Arguments¶

MAT (INPUT/OUTPUT) real(stnd), dimension(:,:)

On entry, the symmetric matrix MAT. The leading n-by-n upper triangular part of MAT contains the upper triangular part of the matrix MAT. The strictly n-by-n lower triangular part of MAT is not referenced.

On exit, MAT contains the orthonormal eigenvectors of the matrix MAT.

The shape of MAT must verify: size( MAT, 1 ) = size( MAT, 2 ) = n .

EIGVAL (OUTPUT) real(stnd), dimension(:)

On exit, the eigenvalues.

The size of EIGVAL must verify: size( EIGVAL ) = n .

FAILURE (OUTPUT) logical(lgl)

On exit:

FAILURE = false : indicates successful exit.

FAILURE = true : indicates that the algorithm did not converge and that full accuracy was not attained in the Schur decomposition of an intermediate tridiagonal form T of MAT .

SORT (INPUT, OPTIONAL) character

Sort the eigenvalues into ascending order if SORT = ‘A’ or ‘a’, or in descending order if SORT = ‘D’ or ‘d’. The eigenvectors are reordered accordingly.

MAXITER (INPUT, OPTIONAL) integer(i4b)

MAXITER controls the maximum number of QR sweeps in the Schur decomposition of an intermediate tridiagonal form T of MAT.

The algorithm fails to converge if the number of QR sweeps exceeds MAXITER * size(EIGVAL). Convergence usually occurs in about 2 * size(EIGVAL) QR sweeps.

The default is 30.

Further Details¶

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

For further details, see

Parlett, B.N., 1998:

The Symmetric Eigenvalue Problem. Revised edition, SIAM, Philadelphia.

`subroutine eig_cmp2 ( mat, eigval, failure, upper, sort, maxiter, max_francis_steps )`¶

Purpose¶

EIG_CMP2 computes all eigenvalues and eigenvectors of a n-by-n real symmetric matrix MAT.

Arguments¶

MAT (INPUT/OUTPUT) real(stnd), dimension(:,:)

On entry, the symmetric matrix MAT.

If:

UPPER = true: the leading n-by-n upper triangular part of MAT contains the upper triangular part of the matrix MAT.

UPPER = false: the leading n-by-n lower triangular part of MAT contains the lower triangular part of the matrix MAT.

On exit, MAT contains the orthonormal eigenvectors of the matrix MAT.

The shape of MAT must verify: size( MAT, 1 ) = size( MAT, 2 ) = n .

EIGVAL (OUTPUT) real(stnd), dimension(:)

On exit, the eigenvalues.

The size of EIGVAL must verify: size( EIGVAL ) = n .

FAILURE (OUTPUT) logical(lgl)

On exit:

FAILURE = false : indicates successful exit.

FAILURE = true : indicates that the algorithm did not converge and that full accuracy was not attained in the Schur decomposition of an intermediate tridiagonal form T of MAT .

UPPER (INPUT) logical(lgl)

Specifies whether the upper or lower triangular part of the symmetric matrix MAT is stored. If:

UPPER = true : Upper triangular is stored ;

UPPER = false: Lower triangular is stored .

SORT (INPUT, OPTIONAL) character

Sort the eigenvalues into ascending order if SORT = ‘A’ or ‘a’, or in descending order if SORT = ‘D’ or ‘d’. The eigenvectors are reordered accordingly.

MAXITER (INPUT, OPTIONAL) integer(i4b)

MAXITER controls the maximum number of QR sweeps in the Schur decomposition of an intermediate tridiagonal form T of MAT .

The algorithm fails to converge if the number of QR sweeps exceeds MAXITER * size(EIGVAL). Convergence usually occurs in about 2 * size(EIGVAL) QR sweeps.

The default is 30.

MAX_FRANCIS_STEPS (INPUT, OPTIONAL) integer(i4b)

MAX_FRANCIS_STEPS controls the maximum number of Francis sets (e.g. QR sweeps) of Givens rotations which must be saved before applying them with a wavefront algorithm to accumulate the eigenvectors in the QR algorithm. MAX_FRANCIS_STEPS is a strictly positive integer, otherwise the default value is used.

The default is the minimum of n and the integer parameter MAX_FRANCIS_STEPS_EIG specified in the module Select_Parameters.

Further Details¶

This driver subroutine is adapted from the routine DSYEV in LAPACK.

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 and the eigenvectors are computed with a perfect shift strategy.

For further details, see:

Parlett, B.N., 1998:

The Symmetric Eigenvalue Problem. Revised edition, SIAM, Philadelphia.

Greenbaum, A., and J. Dongarra, J., 1989:

Experiments with QR/QL Methods for the Symmetric Tridiagonal Eigenproblem. LAPACK Working Note No 17, November 1989.

Van Zee, F.G., van de Geijn, R., and Quintana-Orti, G., 2011:

Restructuring the QR Algorithm for High-Performance Application of Givens Rotations. FLAME Working Note 60. The University of Texas at Austin, Department of Computer Sciences. Technical Report TR-11-36.

`subroutine eig_cmp2 ( mat, eigval, failure, sort, maxiter, max_francis_steps )`¶

Purpose¶

EIG_CMP2 computes all eigenvalues and eigenvectors of a n-by-n real symmetric matrix MAT.

Arguments¶

MAT (INPUT/OUTPUT) real(stnd), dimension(:,:)

On entry, the symmetric matrix MAT. The leading n-by-n upper triangular part of MAT contains the upper triangular part of the matrix MAT. The strictly n-by-n lower triangular part of MAT is not referenced.

On exit, MAT contains the orthonormal eigenvectors of the matrix MAT.

The shape of MAT must verify: size( MAT, 1 ) = size( MAT, 2 ) = n .

EIGVAL (OUTPUT) real(stnd), dimension(:)

On exit, the eigenvalues.

The size of EIGVAL must verify: size( EIGVAL ) = n .

FAILURE (OUTPUT) logical(lgl)

On exit:

FAILURE = false : indicates successful exit.

FAILURE = true : indicates that the algorithm did not converge and that full accuracy was not attained in the Schur decomposition of an intermediate tridiagonal form T of MAT .

SORT (INPUT, OPTIONAL) character

Sort the eigenvalues into ascending order if SORT = ‘A’ or ‘a’, or in descending order if SORT = ‘D’ or ‘d’. The eigenvectors are reordered accordingly.

MAXITER (INPUT, OPTIONAL) integer(i4b)

MAXITER controls the maximum number of QR sweeps in the Schur decomposition of an intermediate tridiagonal form T of MAT .

The algorithm fails to converge if the number of QR sweeps exceeds MAXITER * size(EIGVAL). Convergence usually occurs in about 2 * size(EIGVAL) QR sweeps.

The default is 30.

MAX_FRANCIS_STEPS (INPUT, OPTIONAL) integer(i4b)

MAX_FRANCIS_STEPS controls the maximum number of Francis sets (e.g. QR sweeps) of Givens rotations which must be saved before applying them with a wavefront algorithm to accumulate the eigenvectors in the QR algorithm. MAX_FRANCIS_STEPS is a strictly positive integer, otherwise the default value is used.

The default is the minimum of n and the integer parameter MAX_FRANCIS_STEPS_EIG specified in the module Select_Parameters.

Further Details¶

This driver subroutine is adapted from the routine DSYEV in LAPACK.

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 and the eigenvectors are computed with a perfect shift strategy.

For further details, see:

Parlett, B.N., 1998:

The Symmetric Eigenvalue Problem. Revised edition, SIAM, Philadelphia.

Greenbaum, A., and J. Dongarra, J., 1989:

Experiments with QR/QL Methods for the Symmetric Tridiagonal Eigenproblem. LAPACK Working Note No 17, November 1989.

Van Zee, F.G., van de Geijn, R., and Quintana-Orti, G., 2011:

Restructuring the QR Algorithm for High-Performance Application of Givens Rotations. FLAME Working Note 60. The University of Texas at Austin, Department of Computer Sciences. Technical Report TR-11-36.

`subroutine eig_cmp3 ( mat, eigval, failure, upper, sort, maxiter, max_francis_steps )`¶

Purpose¶

EIG_CMP3 computes all eigenvalues and eigenvectors of a n-by-n real symmetric matrix MAT.

Arguments¶

MAT (INPUT/OUTPUT) real(stnd), dimension(:,:)

On entry, the symmetric matrix MAT.

If:

UPPER = true: the leading n-by-n upper triangular part of MAT contains the upper triangular part of the matrix MAT.

UPPER = false: the leading n-by-n lower triangular part of MAT contains the lower triangular part of the matrix MAT.

On exit, MAT contains the orthonormal eigenvectors of the matrix MAT.

The shape of MAT must verify: size( MAT, 1 ) = size( MAT, 2 ) = n .

EIGVAL (OUTPUT) real(stnd), dimension(:)

On exit, the eigenvalues.

The size of EIGVAL must verify: size( EIGVAL ) = n .

FAILURE (OUTPUT) logical(lgl)

On exit:

FAILURE = false : indicates successful exit.

FAILURE = true : indicates that the algorithm did not converge and that full accuracy was not attained in the Schur decomposition of an intermediate tridiagonal form T of MAT .

UPPER (INPUT) logical(lgl)

Specifies whether the upper or lower triangular part of the symmetric matrix MAT is stored. If:

UPPER = true : Upper triangular is stored ;

UPPER = false: Lower triangular is stored .

SORT (INPUT, OPTIONAL) character

Sort the eigenvalues into ascending order if SORT = ‘A’ or ‘a’, or in descending order if SORT = ‘D’ or ‘d’. The eigenvectors are reordered accordingly.

MAXITER (INPUT, OPTIONAL) integer(i4b)

MAXITER controls the maximum number of QR sweeps in the Schur decomposition of an intermediate tridiagonal form T of MAT .

The algorithm fails to converge if the number of QR sweeps exceeds MAXITER * size(EIGVAL). Convergence usually occurs in about 2 * size(EIGVAL) QR sweeps.

The default is 30.

MAX_FRANCIS_STEPS (INPUT, OPTIONAL) integer(i4b)

MAX_FRANCIS_STEPS controls the maximum number of Francis sets (e.g. QR sweeps) of Givens rotations which must be saved before applying them with a wavefront algorithm to accumulate the eigenvectors in the QR algorithm. MAX_FRANCIS_STEPS is a strictly positive integer, otherwise the default value is used.

The default is the minimum of n and the integer parameter MAX_FRANCIS_STEPS_EIG specified in the module Select_Parameters.

Further Details¶

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

For further details, see:

Parlett, B.N., 1998:

The Symmetric Eigenvalue Problem. Revised edition, SIAM, Philadelphia.

`subroutine eig_cmp3 ( mat, eigval, failure, sort, maxiter, max_francis_steps )`¶

Purpose¶

EIG_CMP3 computes all eigenvalues and eigenvectors of a n-by-n real symmetric matrix MAT.

Arguments¶

MAT (INPUT/OUTPUT) real(stnd), dimension(:,:)

On entry, the symmetric matrix MAT. The leading n-by-n upper triangular part of MAT contains the upper triangular part of the matrix MAT. The strictly n-by-n lower triangular part of MAT is not referenced.

On exit, MAT contains the orthonormal eigenvectors of the matrix MAT.

The shape of MAT must verify: size( MAT, 1 ) = size( MAT, 2 ) = n .

EIGVAL (OUTPUT) real(stnd), dimension(:)

On exit, the eigenvalues.

The size of EIGVAL must verify: size( EIGVAL ) = n .

FAILURE (OUTPUT) logical(lgl)

On exit:

FAILURE = false : indicates successful exit.

FAILURE = true : indicates that the algorithm did not converge and that full accuracy was not attained in the Schur decomposition of an intermediate tridiagonal form T of MAT .

SORT (INPUT, OPTIONAL) character

Sort the eigenvalues into ascending order if SORT = ‘A’ or ‘a’, or in descending order if SORT = ‘D’ or ‘d’. The eigenvectors are reordered accordingly.

MAXITER (INPUT, OPTIONAL) integer(i4b)

MAXITER controls the maximum number of QR sweeps in the Schur decomposition of an intermediate tridiagonal form T of MAT .

The algorithm fails to converge if the number of QR sweeps exceeds MAXITER * size(EIGVAL). Convergence usually occurs in about 2 * size(EIGVAL) QR sweeps.

The default is 30.

MAX_FRANCIS_STEPS (INPUT, OPTIONAL) integer(i4b)

MAX_FRANCIS_STEPS controls the maximum number of Francis sets (e.g. QR sweeps) of Givens rotations which must be saved before applying them with a wavefront algorithm to accumulate the eigenvectors in the QR algorithm. MAX_FRANCIS_STEPS is a strictly positive integer, otherwise the default value is used.

The default is the minimum of n and the integer parameter MAX_FRANCIS_STEPS_EIG specified in the module Select_Parameters.

Further Details¶

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

For further details, see:

Parlett, B.N., 1998:

The Symmetric Eigenvalue Problem. Revised edition, SIAM, Philadelphia.

`function eigvalues ( mat, upper, sort, maxiter )`¶

Purpose¶

Function EIGVALUES computes all eigenvalues of a n-by-n real symmetric matrix MAT.

Arguments¶

MAT (INPUT) real(stnd), dimension(:,:)

On entry, the symmetric matrix MAT.

If:

UPPER = true: The leading n-by-n upper triangular part of MAT contains the upper triangular part of the matrix MAT.

UPPER = false: The leading n-by-n lower triangular part of MAT contains the lower triangular part of the matrix MAT.

The shape of MAT must verify: size( MAT, 1 ) = size( MAT, 2 ) = n .

UPPER (INPUT) logical(lgl)

Specifies whether the upper or lower triangular part of the symmetric matrix MAT is stored. If:

UPPER = true : Upper triangular is stored ;

UPPER = false: Lower triangular is stored .

SORT (INPUT, OPTIONAL) character

Sort the eigenvalues into ascending order if SORT = ‘A’ or ‘a’, or in descending order if SORT = ‘D’ or ‘d’.

MAXITER (INPUT, OPTIONAL) integer(i4b)

MAXITER controls the maximum number of QR sweeps in the Schur decomposition of an intermediate tridiagonal form T of MAT .

The algorithm fails to converge if the number of QR sweeps exceeds MAXITER * size(MAT,1). Convergence usually occurs in about 2 * size(MAT,1) QR sweeps.

The default is 30.

Further Details¶

This function is adapted from the routine DSYEV in LAPACK77 (version 3).

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. If the QR algorithm fails to converge EIGVALUES returns a n-vector filled with NAN() function.

For further details, see:

Parlett, B.N., 1998:

The Symmetric Eigenvalue Problem. Revised edition, SIAM, Philadelphia.

`function eigvalues ( mat, sort, maxiter )`¶

Purpose¶

Function EIGVALUES computes all eigenvalues of a n-by-n real symmetric matrix MAT.

Arguments¶

MAT (INPUT) real(stnd), dimension(:,:)

On entry, the symmetric matrix MAT. The leading n-by-n upper triangular part of MAT contains the upper triangular part of the matrix MAT. The strictly n-by-n lower triangular part of MAT is not used.

The shape of MAT must verify: size( MAT, 1 ) = size( MAT, 2 ) = n .

SORT (INPUT, OPTIONAL) character

Sort the eigenvalues into ascending order if SORT = ‘A’ or ‘a’, or in descending order if SORT = ‘D’ or ‘d’.

MAXITER (INPUT, OPTIONAL) integer(i4b)

MAXITER controls the maximum number of QR sweeps in the Schur decomposition of an intermediate tridiagonal form T of MAT .

The algorithm fails to converge if the number of QR sweeps exceeds MAXITER * size(MAT,1). Convergence usually occurs in about 2 * size(MAT,1) QR sweeps.

The default is 30.

Further Details¶

This function is adapted from the routine DSYEV in LAPACK77 (version 3).

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. If the QR algorithm fails to converge EIGVALUES returns a n-vector filled with NAN() function.

For further details, see:

Parlett, B.N., 1998:

The Symmetric Eigenvalue Problem. Revised edition, SIAM, Philadelphia.

`subroutine eigval_cmp ( mat, eigval, failure, upper, sort, maxiter, d_e )`¶

Purpose¶

EIGVAL_CMP 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’ * MAT * Q = T

,then the eigenvalues are computed by the Pal-Walker-Kahan variant of the QR algorithm.

Arguments¶

MAT (INPUT/OUTPUT) real(stnd), dimension(:,:)

On entry, the symmetric matrix MAT.

If:

UPPER = true: the leading n-by-n upper triangular part of MAT contains the upper triangular part of the matrix MAT.

UPPER = false: the leading n-by-n lower triangular part of MAT contains the lower triangular part of the matrix MAT.

On exit:

If UPPER = true and D_E is present : The leading n-by-n upper triangular part of MAT is overwritten by the matrix Q as a product of elementary reflectors.

If UPPER = false and D_E is present : The leading n-by-n lower triangular part of MAT is overwritten by the matrix Q as a product of elementary reflectors.

If UPPER = true and D_E is absent : The leading n-by-n upper triangular part of MAT is destroyed.

If UPPER = false and D_E is absent : The leading n-by-n lower triangular part of MAT is destroyed.

The shape of MAT must verify: size( MAT, 1 ) = size( MAT, 2 ) = n .

EIGVAL (OUTPUT) real(stnd), dimension(:)

On exit, the eigenvalues.

The size of EIGVAL must verify: size( EIGVAL ) = n .

FAILURE (OUTPUT) logical(lgl)

On exit:

FAILURE = false : indicates successful exit.

FAILURE = true : indicates that the algorithm did not converge and that full accuracy was not attained in the Schur decomposition of the intermediate tridiagonal form T of MAT .

UPPER (INPUT) logical(lgl)

Specifies whether the upper or lower triangular part of the symmetric matrix MAT is stored. If:

UPPER = true : Upper triangular is stored ;

UPPER = false: Lower triangular is stored .

SORT (INPUT, OPTIONAL) character

Sort the eigenvalues into ascending order if SORT = ‘A’ or ‘a’, or in descending order if SORT = ‘D’ or ‘d’.

MAXITER (INPUT, OPTIONAL) integer(i4b)

MAXITER controls the maximum number of QR sweeps in the Schur decomposition of the intermediate tridiagonal form T of MAT .

The algorithm fails to converge if the number of QR sweeps exceeds MAXITER * size(EIGVAL). Convergence usually occurs in about 2 * size(EIGVAL) QR sweeps.

The default is 30.

D_E (OUTPUT, OPTIONAL) real(stnd), dimension(:,2)

On exit, the first column of D_E contains the n diagonal elements of the intermediate tridiagonal form T of MAT. The n-1 first elements of the second column of D_E contains the n-1 subdiagonal elements of the intermediate tridiagonal form T of MAT. D_E(n,2) is arbitrary.

The shape of D_E must verify: size( D_E, 1 ) = n and size( D_E, 2 ) = 2

Further Details¶

This driver subroutine is adapted from the routine DSYEV in LAPACK.

`subroutine eigval_cmp ( mat, eigval, failure, sort, maxiter, d_e )`¶

Purpose¶

EIGVAL_CMP 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’ * MAT * Q = T

,then the eigenvalues are computed by the Pal-Walker-Kahan variant of the QR algorithm.

Arguments¶

MAT (INPUT/OUTPUT) real(stnd), dimension(:,:)

On entry, the symmetric matrix MAT. The leading n-by-n upper triangular part of MAT contains the upper triangular part of the matrix MAT. The strictly n-by-n lower triangular part of MAT is not referenced.

On exit:

If D_E is present: The leading n-by-n upper triangular part of MAT is overwritten by the matrix Q as a product of elementary reflectors.

If D_E is absent : The leading n-by-n upper triangular part of MAT is destroyed.

The shape of MAT must verify: size( MAT, 1 ) = size( MAT, 2 ) = n .

EIGVAL (OUTPUT) real(stnd), dimension(:)

On exit, the eigenvalues.

The size of EIGVAL must verify: size( EIGVAL ) = n .

FAILURE (OUTPUT) logical(lgl)

On exit:

FAILURE = false : indicates successful exit.

FAILURE = true : indicates that the algorithm did not converge and that full accuracy was not attained in the Schur decomposition of an intermediate tridiagonal form T of MAT .

SORT (INPUT, OPTIONAL) character

Sort the eigenvalues into ascending order if SORT = ‘A’ or ‘a’, or in descending order if SORT = ‘D’ or ‘d’.

MAXITER (INPUT, OPTIONAL) integer(i4b)

MAXITER controls the maximum number of QR sweeps in the Schur decomposition of an intermediate tridiagonal form T of MAT .

The algorithm fails to converge if the number of QR sweeps exceeds MAXITER * size(EIGVAL). Convergence usually occurs in about 2 * size(EIGVAL) QR sweeps.

The default is 30.

D_E (OUTPUT, OPTIONAL) real(stnd), dimension(:,2)

On exit, the first column of D_E contains the n diagonal elements of the intermediate tridiagonal form T of MAT. The n-1 first elements of the second column of D_E contains the n-1 subdiagonal elements of the intermediate tridiagonal form T of MAT. D_E(n,2) is arbitrary.

The shape of D_E must verify: size( D_E, 1 ) = n and size( D_E, 2 ) = 2

Further Details¶

This driver subroutine is adapted from the routine DSYEV in LAPACK.

`subroutine eigval_cmp ( matp, eigval, failure, upper, sort, maxiter, d_e )`¶

Purpose¶

EIGVAL_CMP computes all eigenvalues of a n-by-n real symmetric matrix mat stored in packed form in a linear array MATP.

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

Q’ * mat * Q = T

,then the eigenvalues are computed by the Pal-Walker-Kahan variant of the QR algorithm.

Arguments¶

MATP (INPUT/OUTPUT) real(stnd), dimension(:)

On entry, the upper or lower triangle of the symmetric matrix mat, packed column-wise in a linear array. The j-th column of mat is stored in the array MATP as follows:

if UPPER = true, MATP(i + (j-1) * j/2) = mat(i,j) for 1<=i<=j;

if UPPER = false, MATP(i + (j-1) * (2 * n-j)/2) = mat(i,j) for j<=i<=n.

On exit:

If D_E is present : MATP is overwritten by the matrix Q as a product of elementary reflectors.

If D_E is absent : MATP is destroyed.

The size of MATP must verify: size( MATP ) = (n * (n+1)/2)

EIGVAL (OUTPUT) real(stnd), dimension(:)

On exit, the eigenvalues.

The size of EIGVAL must verify: size( EIGVAL ) = n .

FAILURE (OUTPUT) logical(lgl)

On exit:

FAILURE = false : indicates successful exit.

FAILURE = true : indicates that the algorithm did not converge and that full accuracy was not attained in the Schur decomposition of an intermediate tridiagonal form T of mat .

UPPER (INPUT) logical(lgl)

Specifies whether the upper or lower triangular part of the symmetric matrix mat is stored in the linear array MATP. If:

UPPER = true : Upper triangle of mat is stored;

UPPER = false: Lower triangle of mat is stored.

SORT (INPUT, OPTIONAL) character

Sort the eigenvalues into ascending order if SORT = ‘A’ or ‘a’, or in descending order if SORT = ‘D’ or ‘d’.

MAXITER (INPUT, OPTIONAL) integer(i4b)

MAXITER controls the maximum number of QR sweeps in the Schur decomposition of an intermediate tridiagonal form T of mat .

The algorithm fails to converge if the number of QR sweeps exceeds MAXITER * size(EIGVAL). Convergence usually occurs in about 2 * size(EIGVAL) QR sweeps.

The default is 30.

D_E (OUTPUT, OPTIONAL) real(stnd), dimension(:,2)

On exit, the first column of D_E contains the n diagonal elements of the intermediate tridiagonal form T of mat. The n-1 first elements of the second column of D_E contains the n-1 subdiagonal elements of the intermediate tridiagonal form T of mat. D_E(n,2) is arbitrary.

The shape of D_E must verify: size( D_E, 1 ) = n and size( D_E, 2 ) = 2

Further Details¶

This driver subroutine is adapted from the routine DSYEV in LAPACK.

`subroutine eigval_cmp ( matp, eigval, failure, sort, maxiter, d_e )`¶

Purpose¶

EIGVAL_CMP computes all eigenvalues of a n-by-n real symmetric matrix mat stored in packed form in a linear array MATP.

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

Q’ * mat * Q = T

,then the eigenvalues are computed by the Pal-Walker-Kahan variant of the QR algorithm.

Arguments¶

MATP (INPUT/OUTPUT) real(stnd), dimension(:)

On entry, the upper triangle of the symmetric matrix mat, packed column-wise in a linear array. The j-th column of mat is stored in the array MATP as follows:

MATP(i + (j-1) * j/2) = mat(i,j) for 1<=i<=j;

On exit:

If D_E is present : MATP is overwritten by the matrix Q as a product of elementary reflectors.

If D_E is absent : MATP is destroyed.

The size of MATP must verify: size( MATP ) = (n * (n+1)/2)

EIGVAL (OUTPUT) real(stnd), dimension(:)

On exit, the eigenvalues.

The size of EIGVAL must verify: size( EIGVAL ) = n .

FAILURE (OUTPUT) logical(lgl)

On exit:

FAILURE = false : indicates successful exit.

FAILURE = true : indicates that the algorithm did not converge and that full accuracy was not attained in the Schur decomposition of an intermediate tridiagonal form T of mat .

SORT (INPUT, OPTIONAL) character

Sort the eigenvalues into ascending order if SORT = ‘A’ or ‘a’, or in descending order if SORT = ‘D’ or ‘d’.

MAXITER (INPUT, OPTIONAL) integer(i4b)

MAXITER controls the maximum number of QR sweeps in the Schur decomposition of an intermediate tridiagonal form T of mat .

The algorithm fails to converge if the number of QR sweeps exceeds MAXITER * size(EIGVAL). Convergence usually occurs in about 2 * size(EIGVAL) QR sweeps.

The default is 30.

D_E (OUTPUT, OPTIONAL) real(stnd), dimension(:,2)

On exit, the first column of D_E contains the n diagonal elements of the intermediate tridiagonal form T of mat. The n-1 first elements of the second column of D_E contains the n-1 subdiagonal elements of the intermediate tridiagonal form T of mat. D_E(n,2) is arbitrary.

The shape of D_E must verify: size( D_E, 1 ) = n and size( D_E, 2 ) = 2

Further Details¶

This driver subroutine is adapted from the routine DSYEV in LAPACK.

`subroutine eigval_cmp2 ( mat, eigval, failure, upper, sort, maxiter, d_e )`¶

Purpose¶

EIGVAL_CMP2 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’ * MAT * Q = T

,then the eigenvalues are computed by the Pal-Walker-Kahan variant of the QR algorithm.

Arguments¶

MAT (INPUT/OUTPUT) real(stnd), dimension(:,:)

On entry, the symmetric matrix MAT.

If:

UPPER = true: the leading n-by-n upper triangular part of MAT contains the upper triangular part of the matrix MAT.

UPPER = false: the leading n-by-n lower triangular part of MAT contains the lower triangular part of the matrix MAT.

On exit:

If UPPER = true and D_E is present : The leading n-by-n upper triangular part of MAT is overwritten by the matrix Q as a product of elementary reflectors.

If UPPER = false and D_E is present : The leading n-by-n lower triangular part of MAT is overwritten by the matrix Q as a product of elementary reflectors.

If UPPER = true and D_E is absent : The leading n-by-n upper triangular part of MAT is destroyed.

If UPPER = false and D_E is absent : The leading n-by-n lower triangular part of MAT is destroyed.

The shape of MAT must verify: size( MAT, 1 ) = size( MAT, 2 ) = n .

EIGVAL (OUTPUT) real(stnd), dimension(:)

On exit, the eigenvalues.

The size of EIGVAL must verify: size( EIGVAL ) = n .

FAILURE (OUTPUT) logical(lgl)

On exit:

FAILURE = false : indicates successful exit.

FAILURE = true : indicates that the algorithm did not converge and that full accuracy was not attained in the Schur decomposition of an intermediate tridiagonal form T of MAT .

UPPER (INPUT) logical(lgl)

Specifies whether the upper or lower triangular part of the symmetric matrix MAT is stored. If:

UPPER = true : Upper triangular is stored ;

UPPER = false: Lower triangular is stored .

SORT (INPUT, OPTIONAL) character

Sort the eigenvalues into ascending order if SORT = ‘A’ or ‘a’, or in descending order if SORT = ‘D’ or ‘d’.

MAXITER (INPUT, OPTIONAL) integer(i4b)

MAXITER controls the maximum number of QR sweeps in the Schur decomposition of an intermediate tridiagonal form T of MAT .

The algorithm fails to converge if the number of QR sweeps exceeds MAXITER * size(EIGVAL). Convergence usually occurs in about 2 * size(EIGVAL) QR sweeps.

The default is 30.

D_E (OUTPUT, OPTIONAL) real(stnd), dimension(:,2)

On exit, the first column of D_E contains the n diagonal elements of the intermediate tridiagonal form T of MAT. The n-1 first elements of the second column of D_E contains the n-1 subdiagonal elements of the intermediate tridiagonal form T of MAT. D_E(n,2) is arbitrary.

The shape of D_E must verify: size( D_E, 1 ) = n and size( D_E, 2 ) = 2

Further Details¶

This driver subroutine is adapted from the routine DSYEV in LAPACK.

`subroutine eigval_cmp2 ( mat, eigval, failure, sort, maxiter, d_e )`¶

Purpose¶

EIGVAL_CMP2 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’ * MAT * Q = T

,then the eigenvalues are computed by the Pal-Walker-Kahan variant of the QR algorithm.

Arguments¶

MAT (INPUT/OUTPUT) real(stnd), dimension(:,:)

On entry, the symmetric matrix MAT. The leading n-by-n upper triangular part of MAT contains the upper triangular part of the matrix MAT. The strictly n-by-n lower triangular part of MAT is not referenced.

On exit:

If D_E is present: The leading n-by-n upper triangular part of MAT is overwritten by the matrix Q as a product of elementary reflectors.

If D_E is absent : The leading n-by-n upper triangular part of MAT is destroyed.

The shape of MAT must verify: size( MAT, 1 ) = size( MAT, 2 ) = n .

EIGVAL (OUTPUT) real(stnd), dimension(:)

On exit, the eigenvalues.

The size of EIGVAL must verify: size( EIGVAL ) = n .

FAILURE (OUTPUT) logical(lgl)

On exit:

FAILURE = false : indicates successful exit.

FAILURE = true : indicates that the algorithm did not converge and that full accuracy was not attained in the Schur decomposition of an intermediate tridiagonal form T of MAT .

SORT (INPUT, OPTIONAL) character

Sort the eigenvalues into ascending order if SORT = ‘A’ or ‘a’, or in descending order if SORT = ‘D’ or ‘d’.

MAXITER (INPUT, OPTIONAL) integer(i4b)

MAXITER controls the maximum number of QR sweeps in the Schur decomposition of an intermediate tridiagonal form T of MAT .

The algorithm fails to converge if the number of QR sweeps exceeds MAXITER * size(EIGVAL). Convergence usually occurs in about 2 * size(EIGVAL) QR sweeps.

The default is 30.

D_E (OUTPUT, OPTIONAL) real(stnd), dimension(:,2)

On exit, the first column of D_E contains the n diagonal elements of the intermediate tridiagonal form T of MAT. The n-1 first elements of the second column of D_E contains the n-1 subdiagonal elements of the intermediate tridiagonal form T of MAT. D_E(n,2) is arbitrary.

The shape of D_E must verify: size( D_E, 1 ) = n and size( D_E, 2 ) = 2

Further Details¶

This driver subroutine is adapted from the routine DSYEV in LAPACK.

`subroutine eigval_cmp2 ( matp, eigval, failure, upper, sort, maxiter, d_e )`¶

Purpose¶

EIGVAL_CMP2 computes all eigenvalues of a n-by-n real symmetric matrix mat stored in packed form in a linear array MATP.

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

Q’ * mat * Q = T

,then the eigenvalues are computed by the Pal-Walker-Kahan variant of the QR algorithm.

Arguments¶

MATP (INPUT/OUTPUT) real(stnd), dimension(:)

On entry, the upper or lower triangle of the symmetric matrix mat , packed column-wise in a linear array. The j-th column of mat is stored in the array MATP as follows:

if UPPER = true, MATP(i + (j-1) * j/2) = mat(i,j) for 1<=i<=j;

if UPPER = false, MATP(i + (j-1) * (2 * n-j)/2) = mat(i,j) for j<=i<=n.

On exit:

If D_E is present : MATP is overwritten by the matrix Q as a product of elementary reflectors.

If D_E is absent : MATP is destroyed.

The size of MATP must verify: size( MATP ) = (n * (n+1)/2)

EIGVAL (OUTPUT) real(stnd), dimension(:)

On exit, the eigenvalues.

The size of EIGVAL must verify: size( EIGVAL ) = n .

FAILURE (OUTPUT) logical(lgl)

On exit:

FAILURE = false : indicates successful exit.

FAILURE = true : indicates that the algorithm did not converge and that full accuracy was not attained in the Schur decomposition of an intermediate tridiagonal form T of mat .

UPPER (INPUT) logical(lgl)

Specifies whether the upper or lower triangular part of the symmetric matrix mat is stored in the linear array MATP. If:

UPPER = true : Upper triangle of mat is stored;

UPPER = false: Lower triangle of mat is stored.

SORT (INPUT, OPTIONAL) character

Sort the eigenvalues into ascending order if SORT = ‘A’ or ‘a’, or in descending order if SORT = ‘D’ or ‘d’.

MAXITER (INPUT, OPTIONAL) integer(i4b)

MAXITER controls the maximum number of QR sweeps in the Schur decomposition of an intermediate tridiagonal form T of mat .

The algorithm fails to converge if the number of QR sweeps exceeds MAXITER * size(EIGVAL). Convergence usually occurs in about 2 * size(EIGVAL) QR sweeps.

The default is 30.

D_E (OUTPUT, OPTIONAL) real(stnd), dimension(:,2)

On exit, the first column of D_E contains the n diagonal elements of the intermediate tridiagonal form T of mat. The n-1 first elements of the second column of D_E contains the n-1 subdiagonal elements of the intermediate tridiagonal form T of mat. D_E(n,2) is arbitrary.

The shape of D_E must verify: size( D_E, 1 ) = n and size( D_E, 2 ) = 2

Further Details¶

This driver subroutine is adapted from the routine DSYEV in LAPACK.

`subroutine eigval_cmp2 ( matp, eigval, failure, sort, maxiter, d_e )`¶

Purpose¶

EIGVAL_CMP2 computes all eigenvalues of a n-by-n real symmetric matrix mat stored in packed form in a linear array MATP.

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

Q’ * mat * Q = T

,then the eigenvalues are computed by the Pal-Walker-Kahan variant of the QR algorithm.

Arguments¶

MATP (INPUT/OUTPUT) real(stnd), dimension(:)

On entry, the upper triangle of the symmetric matrix mat, packed column-wise in a linear array. The j-th column of mat is stored in the array MATP as follows:

MATP(i + (j-1) * j/2) = mat(i,j) for 1<=i<=j;

On exit:

If D_E is present : MATP is overwritten by the matrix Q as a product of elementary reflectors.

If D_E is absent : MATP is destroyed.

The size of MATP must verify: size( MATP ) = (n * (n+1)/2)

EIGVAL (OUTPUT) real(stnd), dimension(:)

On exit, the eigenvalues.

The size of EIGVAL must verify: size( EIGVAL ) = n .

FAILURE (OUTPUT) logical(lgl)

On exit:

FAILURE = false : indicates successful exit.

FAILURE = true : indicates that the algorithm did not converge and that full accuracy was not attained in the Schur decomposition of an intermediate tridiagonal form T of mat .

SORT (INPUT, OPTIONAL) character

Sort the eigenvalues into ascending order if SORT = ‘A’ or ‘a’, or in descending order if SORT = ‘D’ or ‘d’.

MAXITER (INPUT, OPTIONAL) integer(i4b)

MAXITER controls the maximum number of QR sweeps in the Schur decomposition of an intermediate tridiagonal form T of mat .

The algorithm fails to converge if the number of QR sweeps exceeds MAXITER * size(EIGVAL). Convergence usually occurs in about 2 * size(EIGVAL) QR sweeps.

The default is 30.

D_E (OUTPUT, OPTIONAL) real(stnd), dimension(:,2)

On exit, the first column of D_E contains the n diagonal elements of the intermediate tridiagonal form T of mat. The n-1 first elements of the second column of D_E contains the n-1 subdiagonal elements of the intermediate tridiagonal form T of mat. D_E(n,2) is arbitrary.

The shape of D_E must verify: size( D_E, 1 ) = n and size( D_E, 2 ) = 2

Further Details¶

This driver subroutine is adapted from the routine DSYEV in LAPACK.

`subroutine eigval_cmp3 ( mat, eigval, failure, upper, sort, maxiter, d_e )`¶

Purpose¶

EIGVAL_CMP3 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’ * MAT * Q = T

,then the eigenvalues are computed by the implicit QR algorithm.

Arguments¶

MAT (INPUT/OUTPUT) real(stnd), dimension(:,:)

On entry, the symmetric matrix MAT.

If:

UPPER = true: the leading n-by-n upper triangular part of MAT contains the upper triangular part of the matrix MAT.

UPPER = false: the leading n-by-n lower triangular part of MAT contains the lower triangular part of the matrix MAT.

On exit:

If UPPER = true and D_E is present : The leading n-by-n upper triangular part of MAT is overwritten by the matrix Q as a product of elementary reflectors.

If UPPER = false and D_E is present : The leading n-by-n lower triangular part of MAT is overwritten by the matrix Q as a product of elementary reflectors.

If UPPER = true and D_E is absent : The leading n-by-n upper triangular part of MAT is destroyed.

If UPPER = false and D_E is absent : The leading n-by-n lower triangular part of MAT is destroyed.

The shape of MAT must verify: size( MAT, 1 ) = size( MAT, 2 ) = n .

EIGVAL (OUTPUT) real(stnd), dimension(:)

On exit, the eigenvalues.

The size of EIGVAL must verify: size( EIGVAL ) = n .

FAILURE (OUTPUT) logical(lgl)

On exit:

FAILURE = false : indicates successful exit.

FAILURE = true : indicates that the algorithm did not converge and that full accuracy was not attained in the Schur decomposition of an intermediate tridiagonal form T of MAT .

UPPER (INPUT) logical(lgl)

Specifies whether the upper or lower triangular part of the symmetric matrix MAT is stored. If:

UPPER = true : Upper triangular is stored ;

UPPER = false: Lower triangular is stored .

SORT (INPUT, OPTIONAL) character

Sort the eigenvalues into ascending order if SORT = ‘A’ or ‘a’, or in descending order if SORT = ‘D’ or ‘d’.

MAXITER (INPUT, OPTIONAL) integer(i4b)

MAXITER controls the maximum number of QR sweeps in the Schur decomposition of an intermediate tridiagonal form T of MAT .

The algorithm fails to converge if the number of QR sweeps exceeds MAXITER * size(EIGVAL). Convergence usually occurs in about 2 * size(EIGVAL) QR sweeps.

The default is 30.

D_E (OUTPUT, OPTIONAL) real(stnd), dimension(:,2)

On exit, the first column of D_E contains the n diagonal elements of the intermediate tridiagonal form T of MAT. The n-1 first elements of the second column of D_E contains the n-1 subdiagonal elements of the intermediate tridiagonal form T of MAT. D_E(n,2) is arbitrary.

The shape of D_E must verify: size( D_E, 1 ) = n and size( D_E, 2 ) = 2

Further Details¶

This driver subroutine is adapted from the routine DSYEV in LAPACK.

`subroutine eigval_cmp3 ( mat, eigval, failure, sort, maxiter, d_e )`¶

Purpose¶

EIGVAL_CMP3 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’ * MAT * Q = T

,then the eigenvalues are computed by the implicit QR algorithm.

Arguments¶

MAT (INPUT/OUTPUT) real(stnd), dimension(:,:)

On entry, the symmetric matrix MAT. The leading n-by-n upper triangular part of MAT contains the upper triangular part of the matrix MAT. The strictly n-by-n lower triangular part of MAT is not referenced.

On exit:

If D_E is present: The leading n-by-n upper triangular part of MAT is overwritten by the matrix Q as a product of elementary reflectors.

If D_E is absent : The leading n-by-n upper triangular part of MAT is destroyed.

The shape of MAT must verify: size( MAT, 1 ) = size( MAT, 2 ) = n .

EIGVAL (OUTPUT) real(stnd), dimension(:)

On exit, the eigenvalues.

The size of EIGVAL must verify: size( EIGVAL ) = n .

FAILURE (OUTPUT) logical(lgl)

On exit:

FAILURE = false : indicates successful exit.

FAILURE = true : indicates that the algorithm did not converge and that full accuracy was not attained in the Schur decomposition of an intermediate tridiagonal form T of MAT .

SORT (INPUT, OPTIONAL) character

Sort the eigenvalues into ascending order if SORT = ‘A’ or ‘a’, or in descending order if SORT = ‘D’ or ‘d’.

MAXITER (INPUT, OPTIONAL) integer(i4b)

MAXITER controls the maximum number of QR sweeps in the Schur decomposition of an intermediate tridiagonal form T of MAT .

The algorithm fails to converge if the number of QR sweeps exceeds MAXITER * size(EIGVAL). Convergence usually occurs in about 2 * size(EIGVAL) QR sweeps.

The default is 30.

D_E (OUTPUT, OPTIONAL) real(stnd), dimension(:,2)

On exit, the first column of D_E contains the n diagonal elements of the intermediate tridiagonal form T of MAT. The n-1 first elements of the second column of D_E contains the n-1 subdiagonal elements of the intermediate tridiagonal form T of MAT. D_E(n,2) is arbitrary.

The shape of D_E must verify: size( D_E, 1 ) = n and size( D_E, 2 ) = 2

Further Details¶

This driver subroutine is adapted from the routine DSYEV in LAPACK.

`subroutine eigval_cmp3 ( matp, eigval, failure, upper, sort, maxiter, d_e )`¶

Purpose¶

EIGVAL_CMP3 computes all eigenvalues of a n-by-n real symmetric matrix mat stored in packed form in a linear array MATP.

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

Q’ * mat * Q = T

,then the eigenvalues are computed by the implicit QR algorithm.

Arguments¶

MATP (INPUT/OUTPUT) real(stnd), dimension(:)

On entry, the upper or lower triangle of the symmetric matrix mat , packed column-wise in a linear array. The j-th column of mat is stored in the array MATP as follows:

if UPPER = true, MATP(i + (j-1) * j/2) = mat(i,j) for 1<=i<=j;

if UPPER = false, MATP(i + (j-1) * (2 * n-j)/2) = mat(i,j) for j<=i<=n.

On exit:

If D_E is present : MATP is overwritten by the matrix Q as a product of elementary reflectors.

If D_E is absent : MATP is destroyed.

The size of MATP must verify: size( MATP ) = (n * (n+1)/2)

EIGVAL (OUTPUT) real(stnd), dimension(:)

On exit, the eigenvalues.

The size of EIGVAL must verify: size( EIGVAL ) = n .

FAILURE (OUTPUT) logical(lgl)

On exit:

FAILURE = false : indicates successful exit.

FAILURE = true : indicates that the algorithm did not converge and that full accuracy was not attained in the Schur decomposition of an intermediate tridiagonal form T of mat .

UPPER (INPUT) logical(lgl)

Specifies whether the upper or lower triangular part of the symmetric matrix mat is stored in the linear array MATP. If:

UPPER = true : Upper triangle of mat is stored ;

UPPER = false: Lower triangle of mat is stored .

SORT (INPUT, OPTIONAL) character

Sort the eigenvalues into ascending order if SORT = ‘A’ or ‘a’, or in descending order if SORT = ‘D’ or ‘d’.

MAXITER (INPUT, OPTIONAL) integer(i4b)

MAXITER controls the maximum number of QR sweeps in the Schur decomposition of an intermediate tridiagonal form T of mat .

The algorithm fails to converge if the number of QR sweeps exceeds MAXITER * size(EIGVAL). Convergence usually occurs in about 2 * size(EIGVAL) QR sweeps.

The default is 30.

D_E (OUTPUT, OPTIONAL) real(stnd), dimension(:,2)

On exit, the first column of D_E contains the n diagonal elements of the intermediate tridiagonal form T of mat. The n-1 first elements of the second column of D_E contains the n-1 subdiagonal elements of the intermediate tridiagonal form T of mat. D_E(n,2) is arbitrary.

The shape of D_E must verify: size( D_E, 1 ) = n and size( D_E, 2 ) = 2

Further Details¶

This driver subroutine is adapted from the routine DSYEV in LAPACK.

`subroutine eigval_cmp3 ( matp, eigval, failure, sort, maxiter, d_e )`¶

Purpose¶

EIGVAL_CMP3 computes all eigenvalues of a n-by-n real symmetric matrix mat stored in packed form in a linear array MATP.

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

Q’ * mat * Q = T

,then the eigenvalues are computed by the implicit QR algorithm.

Arguments¶

MATP (INPUT/OUTPUT) real(stnd), dimension(:)

On entry, the upper triangle of the symmetric matrix mat , packed column-wise in a linear array. The j-th column of mat is stored in the array MATP as follows:

MATP(i + (j-1) * j/2) = mat(i,j) for 1<=i<=j;

On exit:

If D_E is present : MATP is overwritten by the matrix Q as a product of elementary reflectors.

If D_E is absent : MATP is destroyed.

The size of MATP must verify: size( MATP ) = (n * (n+1)/2)

EIGVAL (OUTPUT) real(stnd), dimension(:)

On exit, the eigenvalues.

The size of EIGVAL must verify: size( EIGVAL ) = n .

FAILURE (OUTPUT) logical(lgl)

On exit:

FAILURE = false : indicates successful exit.

FAILURE = true : indicates that the algorithm did not converge and that full accuracy was not attained in the Schur decomposition of an intermediate tridiagonal form T of mat .

SORT (INPUT, OPTIONAL) character

Sort the eigenvalues into ascending order if SORT = ‘A’ or ‘a’, or in descending order if SORT = ‘D’ or ‘d’.

MAXITER (INPUT, OPTIONAL) integer(i4b)

MAXITER controls the maximum number of QR sweeps in the Schur decomposition of an intermediate tridiagonal form T of mat .

The algorithm fails to converge if the number of QR sweeps exceeds MAXITER * size(EIGVAL). Convergence usually occurs in about 2 * size(EIGVAL) QR sweeps.

The default is 30.

D_E (OUTPUT, OPTIONAL) real(stnd), dimension(:,2)

On exit, the first column of D_E contains the n diagonal elements of the intermediate tridiagonal form T of mat. The n-1 first elements of the second column of D_E contains the n-1 subdiagonal elements of the intermediate tridiagonal form T of mat. D_E(n,2) is arbitrary.

The shape of D_E must verify: size( D_E, 1 ) = n and size( D_E, 2 ) = 2

Further Details¶

This driver subroutine is adapted from the routine DSYEV in LAPACK.

`subroutine select_eigval_cmp ( mat, eigval, small, failure, upper, d_e )`¶

Purpose¶

SELECT_EIGVAL_CMP computes the m=size( EIGVAL ) largest or smallest 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’ * MAT * Q = T

,then the eigenvalues are computed by a rational QR method.

Arguments¶

MAT (INPUT/OUTPUT) real(stnd), dimension(:,:)

On entry, the symmetric matrix MAT.

If:

UPPER = true: the leading n-by-n upper triangular part of MAT contains the upper triangular part of the matrix MAT.

UPPER = false: the leading n-by-n lower triangular part of MAT contains the lower triangular part of the matrix MAT.

On exit:

If UPPER = true and D_E is present : The leading n-by-n upper triangular part of MAT is overwritten by the matrix Q as a product of elementary reflectors.

If UPPER = false and D_E is present : The leading n-by-n lower triangular part of MAT is overwritten by the matrix Q as a product of elementary reflectors.

If UPPER = true and D_E is absent : The leading n-by-n upper triangular part of MAT is destroyed.

If UPPER = false and D_E is absent : The leading n-by-n lower triangular part of MAT is destroyed.

The shape of MAT must verify: size( MAT, 1 ) = size( MAT, 2 ) = n .

EIGVAL (OUTPUT) real(stnd), dimension(:)

On exit, the m=size( EIGVAL ) largest or smallest eigenvalues of MAT in decreasing sequence.

The size of EIGVAL must verify: size( EIGVAL )<= size( MAT, 1 ) = size( MAT, 2 ) = n .

SMALL (INPUT) logical(lgl)

On entry:

SMALL = false : indicates that the m largest eigenvalues are desired.

SMALL = true : indicates that the m smallest eigenvalues are desired.

FAILURE (OUTPUT) logical(lgl)

On exit:

FAILURE = false : indicates successful exit.

FAILURE = true : indicates that the algorithm did not converge and that full accuracy was not attained in the rational QR iterations for some eigenvalues of the intermediate tridiagonal form T of MAT .

UPPER (INPUT) logical(lgl)

Specifies whether the upper or lower triangular part of the symmetric matrix MAT is stored. If:

UPPER = true : Upper triangular is stored ;

UPPER = false: Lower triangular is stored .

D_E (OUTPUT, OPTIONAL) real(stnd), dimension(:,2)

On exit, the first column of D_E contains the n diagonal elements of the intermediate tridiagonal form T of MAT. The n-1 first elements of the second column of D_E contains the n-1 subdiagonal elements of the intermediate tridiagonal form T of MAT. D_E(n,2) is arbitrary.

The shape of D_E must verify: size( D_E, 1 ) = n and size( D_E, 2 ) = 2

Further Details¶

This driver subroutine is adapted from the routine DSYEV in LAPACK.

`subroutine select_eigval_cmp ( mat, eigval, small, failure, d_e )`¶

Purpose¶

SELECT_EIGVAL_CMP computes the m=size( EIGVAL ) largest or smallest 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’ * MAT * Q = T

,then the eigenvalues are computed by a rational QR method.

Arguments¶

MAT (INPUT/OUTPUT) real(stnd), dimension(:,:)

On entry, the symmetric matrix MAT. The leading n-by-n upper triangular part of MAT contains the upper triangular part of the matrix MAT. The strictly n-by-n lower triangular part of MAT is not referenced.

On exit:

If D_E is present: The leading n-by-n upper triangular part of MAT is overwritten by the matrix Q as a product of elementary reflectors.

If D_E is absent : The leading n-by-n upper triangular part of MAT is destroyed.

The shape of MAT must verify: size( MAT, 1 ) = size( MAT, 2 ) = n .

EIGVAL (OUTPUT) real(stnd), dimension(:)

On exit, the m=size( EIGVAL ) largest or smallest eigenvalues of MAT in decreasing sequence.

The size of EIGVAL must verify: size( EIGVAL )<= size( MAT, 1 ) = size( MAT, 2 ) = n .

SMALL (INPUT) logical(lgl)

On entry:

SMALL = false : indicates that the m largest eigenvalues are desired.

SMALL = true : indicates that the m smallest eigenvalues are desired.

FAILURE (OUTPUT) logical(lgl)

On exit:

FAILURE = false : indicates successful exit.

FAILURE = true : indicates that the algorithm did not converge and that full accuracy was not attained in the rational QR iterations for some eigenvalues of the intermediate tridiagonal form T of MAT .

D_E (OUTPUT, OPTIONAL) real(stnd), dimension(:,2)

On exit, the first column of D_E contains the n diagonal elements of the intermediate tridiagonal form T of MAT. The n-1 first elements of the second column of D_E contains the n-1 subdiagonal elements of the intermediate tridiagonal form T of MAT. D_E(n,2) is arbitrary.

The shape of D_E must verify: size( D_E, 1 ) = n and size( D_E, 2 ) = 2

Further Details¶

This driver subroutine is adapted from the routine DSYEV in LAPACK.

`subroutine select_eigval_cmp ( matp, eigval, small, failure, upper, d_e )`¶

Purpose¶

SELECT_EIGVAL_CMP computes the m=size( EIGVAL ) largest or smallest eigenvalues of a n-by-n real symmetric matrix mat stored in packed form in a linear array MATP.

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

Q’ * mat * Q = T

,then the eigenvalues are computed by a rational QR method.

Arguments¶

MATP (INPUT/OUTPUT) real(stnd), dimension(:)

On entry, the upper or lower triangle of the symmetric matrix mat, packed column-wise in a linear array. The j-th column of mat is stored in the array MATP as follows:

if UPPER = true, MATP(i + (j-1) * j/2) = mat(i,j) for 1<=i<=j;

if UPPER = false, MATP(i + (j-1) * (2 * n-j)/2) = mat(i,j) for j<=i<=n.

On exit:

If D_E is present : MATP is overwritten by the matrix Q as a product of elementary reflectors.

If D_E is absent : MATP is destroyed.

The size of MATP must verify: size( MATP ) = (n * (n+1)/2).

EIGVAL (OUTPUT) real(stnd), dimension(:)

On exit, the m=size( EIGVAL ) largest or smallest eigenvalues of MAT in decreasing sequence.

The size of EIGVAL must verify: (m * (m+1))/2 <= size( MATP ) = (n * (n+1)/2) .

SMALL (INPUT) logical(lgl)

On entry:

SMALL = false : indicates that the m largest eigenvalues are desired.

SMALL = true : indicates that the m smallest eigenvalues are desired.

FAILURE (OUTPUT) logical(lgl)

On exit:

FAILURE = false : indicates successful exit.

FAILURE = true : indicates that the algorithm did not converge and that full accuracy was not attained in the rational QR iterations for some eigenvalues of the intermediate tridiagonal form T of MAT .

UPPER (INPUT) logical(lgl)

Specifies whether the upper or lower triangular part of the symmetric matrix mat is stored in the linear array MATP. If:

UPPER = true : Upper triangle of mat is stored;

UPPER = false: Lower triangle of mat is stored.

D_E (OUTPUT, OPTIONAL) real(stnd), dimension(:,2)

On exit, the first column of D_E contains the n diagonal elements of the intermediate tridiagonal form T of mat. The n-1 first elements of the second column of D_E contains the n-1 subdiagonal elements of the intermediate tridiagonal form T of mat. D_E(n,2) is arbitrary.

The shape of D_E must verify: size( D_E, 1 ) = n and size( D_E, 2 ) = 2

Further Details¶

This driver subroutine is adapted from the routine DSYEV in LAPACK.

`subroutine select_eigval_cmp ( matp, eigval, small, failure, d_e )`¶

Purpose¶

SELECT_EIGVAL_CMP computes the m=size( EIGVAL ) largest or smallest eigenvalues of a n-by-n real symmetric matrix mat stored in packed form in a linear array MATP.

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

Q’ * mat * Q = T

,then the eigenvalues are computed by a rational QR method.

Arguments¶

MATP (INPUT/OUTPUT) real(stnd), dimension(:)

On entry, the upper triangle of the symmetric matrix mat, packed column-wise in a linear array. The j-th column of mat is stored in the array MATP as follows:

MATP(i + (j-1) * j/2) = mat(i,j) for 1<=i<=j;

On exit:

If D_E is present : MATP is overwritten by the matrix Q as a product of elementary reflectors.

If D_E is absent : MATP is destroyed.

The size of MATP must verify: size( MATP ) = (n * (n+1)/2).

EIGVAL (OUTPUT) real(stnd), dimension(:)

On exit, the m=size( EIGVAL ) largest or smallest eigenvalues of MAT in decreasing sequence.

The size of EIGVAL must verify: (m * (m+1))/2 <= size( MATP ) = (n * (n+1)/2) .

SMALL (INPUT) logical(lgl)

On entry:

SMALL = false : indicates that the m largest eigenvalues are desired.

SMALL = true : indicates that the m smallest eigenvalues are desired.

FAILURE (OUTPUT) logical(lgl)

On exit:

FAILURE = false : indicates successful exit.

FAILURE = true : indicates that the algorithm did not converge and that full accuracy was not attained in the rational QR iterations for some eigenvalues of the intermediate tridiagonal form T of MAT .

D_E (OUTPUT, OPTIONAL) real(stnd), dimension(:,2)

On exit, the first column of D_E contains the n diagonal elements of the intermediate tridiagonal form T of mat. The n-1 first elements of the second column of D_E contains the n-1 subdiagonal elements of the intermediate tridiagonal form T of mat. D_E(n,2) is arbitrary.

The shape of D_E must verify: size( D_E, 1 ) = n and size( D_E, 2 ) = 2

Further Details¶

This driver subroutine is adapted from the routine DSYEV in LAPACK.

`subroutine select_eigval_cmp2 ( mat, eigval, small, val, failure, upper, d_e )`¶

Purpose¶

SELECT_EIGVAL_CMP2 computes the largest or smallest eigenvalues of a n-by-n real symmetric matrix MAT whose sum of absolute values exceeds a given (positive) threshold.

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

Q’ * MAT * Q = T

,then the eigenvalues are computed by a rational QR method.

Arguments¶

MAT (INPUT/OUTPUT) real(stnd), dimension(:,:)

On entry, the symmetric matrix MAT.

If:

UPPER = true: the leading n-by-n upper triangular part of MAT contains the upper triangular part of the matrix MAT.

UPPER = false: the leading n-by-n lower triangular part of MAT contains the lower triangular part of the matrix MAT.

On exit:

If UPPER = true and D_E is present : The leading n-by-n upper triangular part of MAT is overwritten by the matrix Q as a product of elementary reflectors.

If UPPER = false and D_E is present : The leading n-by-n lower triangular part of MAT is overwritten by the matrix Q as a product of elementary reflectors.

If UPPER = true and D_E is absent : The leading n-by-n upper triangular part of MAT is destroyed.

If UPPER = false and D_E is absent : The leading n-by-n lower triangular part of MAT is destroyed.

The shape of MAT must verify: size( MAT, 1 ) = size( MAT, 2 ) = n .

EIGVAL (OUTPUT) real(stnd), pointer, dimension(:)

On exit, the computed largest or smallest eigenvalues of MAT in decreasing sequence.

On exit, m = size( EIGVAL ) gives the number of computed eigenvalues.

SMALL (INPUT) logical(lgl)

On entry:

SMALL = false : indicates that the largest eigenvalues are desired.

SMALL = true : indicates that the smallest eigenvalues are desired.

VAL (INPUT) real(stnd)

On entry, the sum of the m eigenvalues (in absolute value) found will exceed abs(VAL) or m is equal to n.

FAILURE (OUTPUT) logical(lgl)

On exit:

FAILURE = false : indicates successful exit.

FAILURE = true : indicates that the algorithm did not converge and that full accuracy was not attained in the rational QR iterations for some eigenvalues of the intermediate tridiagonal form T of MAT .

UPPER (INPUT) logical(lgl)

Specifies whether the upper or lower triangular part of the symmetric matrix MAT is stored. If:

UPPER = true : Upper triangular is stored ;

UPPER = false: Lower triangular is stored .

D_E (OUTPUT, OPTIONAL) real(stnd), dimension(:,2)

On exit, the first column of D_E contains the n diagonal elements of the intermediate tridiagonal form T of MAT. The n-1 first elements of the second column of D_E contains the n-1 subdiagonal elements of the intermediate tridiagonal form T of MAT. D_E(n,2) is arbitrary.

The shape of D_E must verify: size( D_E, 1 ) = n and size( D_E, 2 ) = 2

Further Details¶

This driver subroutine is adapted from the routine DSYEV in LAPACK.

`subroutine select_eigval_cmp2 ( mat, eigval, small, val, failure, d_e )`¶

Purpose¶

SELECT_EIGVAL_CMP2 computes the largest or smallest eigenvalues of a n-by-n real symmetric matrix MAT whose sum of absolute values exceeds a given (positive) threshold.

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

Q’ * MAT * Q = T

,then the eigenvalues are computed by a rational QR method.

Arguments¶

MAT (INPUT/OUTPUT) real(stnd), dimension(:,:)

On entry, the symmetric matrix MAT. The leading n-by-n upper triangular part of MAT contains the upper triangular part of the matrix MAT. The strictly n-by-n lower triangular part of MAT is not referenced.

On exit:

If D_E is present: The leading n-by-n upper triangular part of MAT is overwritten by the matrix Q as a product of elementary reflectors.

If D_E is absent : The leading n-by-n upper triangular part of MAT is destroyed.

The shape of MAT must verify: size( MAT, 1 ) = size( MAT, 2 ) = n .

EIGVAL (OUTPUT) real(stnd), pointer, dimension(:)

On exit, the computed largest or smallest eigenvalues of MAT in decreasing sequence.

On exit, m = size( EIGVAL ) gives the number of computed eigenvalues.

SMALL (INPUT) logical(lgl)

On entry:

SMALL = false : indicates that the largest eigenvalues are desired.

SMALL = true : indicates that the smallest eigenvalues are desired.

VAL (INPUT) real(stnd)

On entry, the sum of the m eigenvalues (in absolute value) found will exceed abs(VAL) or m is equal to n.

FAILURE (OUTPUT) logical(lgl)

On exit:

FAILURE = false : indicates successful exit.

FAILURE = true : indicates that the algorithm did not converge and that full accuracy was not attained in the rational QR iterations for some eigenvalues of the intermediate tridiagonal form T of MAT .

D_E (OUTPUT, OPTIONAL) real(stnd), dimension(:,2)

On exit, the first column of D_E contains the n diagonal elements of the intermediate tridiagonal form T of MAT. The n-1 first elements of the second column of D_E contains the n-1 subdiagonal elements of the intermediate tridiagonal form T of MAT. D_E(n,2) is arbitrary.

The shape of D_E must verify: size( D_E, 1 ) = n and size( D_E, 2 ) = 2

Further Details¶

This driver subroutine is adapted from the routine DSYEV in LAPACK.

`subroutine select_eigval_cmp2 ( matp, eigval, small, val, failure, upper, d_e )`¶

Purpose¶

SELECT_EIGVAL_CMP2 computes the largest or smallest eigenvalues of a n-by-n real symmetric matrix mat, stored in packed form in a linear array MATP, whose sum of absolute values exceeds a given (positive) threshold.

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

Q’ * mat * Q = T

,then the eigenvalues are computed by a rational QR method.

Arguments¶

MATP (INPUT/OUTPUT) real(stnd), dimension(:)

On entry, the upper or lower triangle of the symmetric matrix mat, packed column-wise in a linear array. The j-th column of mat is stored in the array MATP as follows:

if UPPER = true, MATP(i + (j-1) * j/2) = mat(i,j) for 1<=i<=j;

if UPPER = false, MATP(i + (j-1) * (2 * n-j)/2) = mat(i,j) for j<=i<=n.

On exit:

If D_E is present : MATP is overwritten by the matrix Q as a product of elementary reflectors.

If D_E is absent : MATP is destroyed.

The size of MATP must verify: size( MATP ) = (n * (n+1)/2).

EIGVAL (OUTPUT) real(stnd), pointer, dimension(:)

On exit, the computed largest or smallest eigenvalues of MAT in decreasing sequence.

On exit, m = size( EIGVAL ) gives the number of computed eigenvalues.

SMALL (INPUT) logical(lgl)

On entry:

SMALL = false : indicates that the largest eigenvalues are desired.

SMALL = true : indicates that the smallest eigenvalues are desired.

VAL (INPUT) real(stnd)

On entry, the sum of the m eigenvalues (in absolute value) found will exceed abs(VAL) or m is equal to n.

FAILURE (OUTPUT) logical(lgl)

On exit:

FAILURE = false : indicates successful exit.

FAILURE = true : indicates that the algorithm did not converge and that full accuracy was not attained in the rational QR iterations for some eigenvalues of the intermediate tridiagonal form T of MAT .

UPPER (INPUT) logical(lgl)

Specifies whether the upper or lower triangular part of the symmetric matrix mat is stored in the linear array MATP. If:

UPPER = true : Upper triangle of mat is stored;

UPPER = false: Lower triangle of mat is stored.

D_E (OUTPUT, OPTIONAL) real(stnd), dimension(:,2)

On exit, the first column of D_E contains the n diagonal elements of the intermediate tridiagonal form T of mat. The n-1 first elements of the second column of D_E contains the n-1 subdiagonal elements of the intermediate tridiagonal form T of mat. D_E(n,2) is arbitrary.

The shape of D_E must verify: size( D_E, 1 ) = n and size( D_E, 2 ) = 2

Further Details¶

This driver subroutine is adapted from the routine DSYEV in LAPACK.

`subroutine select_eigval_cmp2 ( matp, eigval, small, val, failure, d_e )`¶

Purpose¶

SELECT_EIGVAL_CMP2 computes the largest or smallest eigenvalues of a n-by-n real symmetric matrix mat, stored in packed form in a linear array MATP, whose sum of absolute values exceeds a given (positive) threshold.

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

Q’ * mat * Q = T

,then the eigenvalues are computed by a rational QR method.

Arguments¶

MATP (INPUT/OUTPUT) real(stnd), dimension(:)

On entry, the upper triangle of the symmetric matrix mat , packed column-wise in a linear array. The j-th column of mat is stored in the array MATP as follows:

MATP(i + (j-1) * j/2) = mat(i,j) for 1<=i<=j;

On exit:

If D_E is present : MATP is overwritten by the matrix Q as a product of elementary reflectors.

If D_E is absent : MATP is destroyed.

The size of MATP must verify: size( MATP ) = (n * (n+1)/2).

EIGVAL (OUTPUT) real(stnd), pointer, dimension(:)

On exit, the computed largest or smallest eigenvalues of MAT in decreasing sequence.

On exit, m = size( EIGVAL ) gives the number of computed eigenvalues.

SMALL (INPUT) logical(lgl)

On entry:

SMALL = false : indicates that the largest eigenvalues are desired.

SMALL = true : indicates that the smallest eigenvalues are desired.

VAL (INPUT) real(stnd)

On entry, the sum of the m eigenvalues (in absolute value) found will exceed abs(VAL) or m is equal to n.

FAILURE (OUTPUT) logical(lgl)

On exit:

FAILURE = false : indicates successful exit.

FAILURE = true : indicates that the algorithm did not converge and that full accuracy was not attained in the rational QR iterations for some eigenvalues of the intermediate tridiagonal form T of MAT .

D_E (OUTPUT, OPTIONAL) real(stnd), dimension(:,2)

On exit, the first column of D_E contains the n diagonal elements of the intermediate tridiagonal form T of mat. The n-1 first elements of the second column of D_E contains the n-1 subdiagonal elements of the intermediate tridiagonal form T of mat. D_E(n,2) is arbitrary.

The shape of D_E must verify: size( D_E, 1 ) = n and size( D_E, 2 ) = 2

Further Details¶

This driver subroutine is adapted from the routine DSYEV in LAPACK.

`subroutine select_eigval_cmp3 ( mat, neig, eigval, small, failure, upper, sort, vector, scaling, init, abstol, le, theta, d_e )`¶

Purpose¶

SELECT_EIGVAL_CMP3 computes the largest or smallest 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’ * MAT * Q = T

,then the eigenvalues are computed by a bisection method.

Arguments¶

MAT (INPUT/OUTPUT) real(stnd), dimension(:,:)

On entry, the symmetric matrix MAT.

If:

UPPER = true: the leading n-by-n upper triangular part of MAT contains the upper triangular part of the matrix MAT.

UPPER = false: the leading n-by-n lower triangular part of MAT contains the lower triangular part of the matrix MAT.

On exit:

If UPPER = true and D_E is present : The leading n-by-n upper triangular part of MAT is overwritten by the matrix Q as a product of elementary reflectors.

If UPPER = false and D_E is present : The leading n-by-n lower triangular part of MAT is overwritten by the matrix Q as a product of elementary reflectors.

If UPPER = true and D_E is absent : The leading n-by-n upper triangular part of MAT is destroyed.

If UPPER = false and D_E is absent : The leading n-by-n lower triangular part of MAT is destroyed.

The shape of MAT must verify: size( MAT, 1 ) = size( MAT, 2 ) = n .

NEIG (OUTPUT) integer(i4b)

On output, NEIG specifies the number of eigenvalues which have been computed. Note that NEIG may be greater than the optional argument LE, if multiple eigenvalues at index LE make unique selection impossible.

If none of the optional arguments LE and THETA are used, NEIG is set to n and all the eigenvalues of MAT are computed.

EIGVAL (OUTPUT) real(stnd), dimension(:)

On exit, EIGVAL(1:NEIG) contains the first NEIG largest or smallest eigenvalues of MAT. The other values in EIGVAL (e.g. EIGVAL(NEIG+1:) ) are flagged by a quiet NAN.

The size of EIGVAL must verify: size( EIGVAL ) = size( MAT, 1 ) = size( MAT, 2 ) = n .

SMALL (INPUT) logical(lgl)

On entry:

SMALL = false : indicates that the largest eigenvalues are desired.

SMALL = true : indicates that the smallest eigenvalues are desired.

FAILURE (OUTPUT) logical(lgl)

On exit:

FAILURE = false : indicates successful exit and the bisection algorithm converged for all the computed eigenvalues to the desired accuracy ;

FAILURE = true : indicates that some or all of the eigenvalues failed to converge or were not computed. This is generally caused by unexpectedly inaccurate arithmetic.

UPPER (INPUT) logical(lgl)

Specifies whether the upper or lower triangular part of the symmetric matrix MAT is stored. If:

UPPER = true : Upper triangular is stored ;

UPPER = false: Lower triangular is stored .

SORT (INPUT, OPTIONAL) character

Sort the eigenvalues into ascending order if SORT = ‘A’ or ‘a’, or in descending order if SORT = ‘D’ or ‘d’. For other values of SORT nothing is done and EIGVAL(:NEIG) may not be sorted.

VECTOR (INPUT, OPTIONAL) logical(lgl)

On entry, if VECTOR is set to true, a vectorized version of the bisection algorithm is used to find the eigenvalues of the intermediate tridiagonal form T of MAT.

The default is VECTOR=false.

SCALING (INPUT, OPTIONAL) logical(lgl)

On entry, if SCALING=true the intermediate tridiagonal matrix T is scaled before computing the eigenvalues.

The default is to scale the tridiagonal matrix.

INIT (INPUT, OPTIONAL) logical(lgl)

On entry, if INIT=true the initial intervals for the bisection steps for computing the eigenvalues of the intermediate tridiagonal matrix T are estimated from the eigenvalues of the intermediate tridiagonal matrix obtained from the Pal-Walker-Kahan algorithm.

The default is not to use the Pal-Walker-Kahan algorithm.

ABSTOL (INPUT, OPTIONAL) real(stnd)

On entry, the absolute tolerance for the eigenvalues. An eigenvalue (or cluster) is considered to be located if it has been determined to lie in an interval whose width is ABSTOL or less.

If ABSTOL is less than or equal to zero, or is not specified, then ULP * | T | will be used, where | T | means the 1-norm of T (T is the intermediate tridiagonal form of MAT) and ULP is the machine precision (distance from 1 to the next larger floating point number).

Eigenvalues will be computed most accurately when ABSTOL is set to the square root of the underflow threshold, sqrt(LAMCH(‘S’)), not zero.

LE (INPUT, OPTIONAL) integer(i4b)

On entry, LE specifies the number of eigenvalues which must be computed by the subroutine. However, on output, NEIG may be different than LE if multiple eigenvalues at index LE make unique selection impossible.

If:

SMALL=false, the subroutine computes the LE largest eigenvalues of MAT,

SMALL=true, the subroutine computes the LE smallest eigenvalues of MAT.

Only one of the optional arguments LE and THETA must be specified, otherwise the subroutine will stop with an error message.

LE must be greater than 0 and less or equal to size( EIGVAL ) .

The default is LE = size( EIGVAL ), e.g. all the eigenvalues are computed.

THETA (INPUT, OPTIONAL) real(stnd)

On entry:

if SMALL=false, THETA specifies that the eigenvalues which are greater or equal to THETA must be computed. If none of the eigenvalues are greater or equal to THETA, NEIG is set to zero and EIGVAL(:) to a quiet NAN.

if SMALL=true, THETA specifies that the eigenvalues which are less or equal to THETA must be computed. If none of the eigenvalues are smaller or equal to THETA, NEIG is set to zero and EIGVAL(:) to a quiet NAN.

Only one of the optional arguments LS and THETA must be specified, otherwise the subroutine will stop with an error message.

D_E (OUTPUT, OPTIONAL) real(stnd), dimension(:,2)

On exit, the first column of D_E contains the n diagonal elements of the intermediate tridiagonal form T of MAT. The n-1 first elements of the second column of D_E contains the n-1 subdiagonal elements of the intermediate tridiagonal form T of MAT. D_E(n,2) is arbitrary.

The shape of D_E must verify: size( D_E, 1 ) = n and size( D_E, 2 ) = 2

Further Details¶

This driver subroutine is adapted from the routine DSYEV in LAPACK.

`subroutine select_eigval_cmp3 ( mat, neig, eigval, small, failure, sort, vector, scaling, init, abstol, le, theta, d_e )`¶

Purpose¶

SELECT_EIGVAL_CMP3 computes the largest or smallest 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’ * MAT * Q = T

,then the eigenvalues are computed by a bisection method.

Arguments¶

MAT (INPUT/OUTPUT) real(stnd), dimension(:,:)

On entry, the symmetric matrix MAT. The leading n-by-n upper triangular part of MAT contains the upper triangular part of the matrix MAT. The strictly n-by-n lower triangular part of MAT is not referenced.

On exit:

If D_E is present: The leading n-by-n upper triangular part of MAT is overwritten by the matrix Q as a product of elementary reflectors.

If D_E is absent : The leading n-by-n upper triangular part of MAT is destroyed.

The shape of MAT must verify: size( MAT, 1 ) = size( MAT, 2 ) = n .

NEIG (OUTPUT) integer(i4b)

On output, NEIG specifies the number of eigenvalues which have been computed. Note that NEIG may be greater than the optional argument LE, if multiple eigenvalues at index LE make unique selection impossible.

If none of the optional arguments LE and THETA are used, NEIG is set to n and all the eigenvalues of MAT are computed.

EIGVAL (OUTPUT) real(stnd), dimension(:)

On exit, EIGVAL(1:NEIG) contains the first NEIG largest or smallest eigenvalues of MAT. The other values in EIGVAL (e.g. EIGVAL(NEIG+1:) ) are flagged by a quiet NAN.

The size of EIGVAL must verify: size( EIGVAL ) = size( MAT, 1 ) = size( MAT, 2 ) = n .

SMALL (INPUT) logical(lgl)

On entry:

SMALL = false : indicates that the m largest eigenvalues are desired.

SMALL = true : indicates that the m smallest eigenvalues are desired.

FAILURE (OUTPUT) logical(lgl)

On exit:

FAILURE = false : indicates successful exit and the bisection algorithm converged for all the computed eigenvalues to the desired accuracy ;

FAILURE = true : indicates that some or all of the eigenvalues failed to converge or were not computed. This is generally caused by unexpectedly inaccurate arithmetic.

SORT (INPUT, OPTIONAL) character

Sort the eigenvalues into ascending order if SORT = ‘A’ or ‘a’, or in descending order if SORT = ‘D’ or ‘d’. For other values of SORT nothing is done and EIGVAL(:NEIG) may not be sorted.

VECTOR (INPUT, OPTIONAL) logical(lgl)

On entry, if VECTOR is set to true, a vectorized version of the bisection algorithm is used to find the eigenvalues of the intermediate tridiagonal form T of MAT.

The default is VECTOR=false.

SCALING (INPUT, OPTIONAL) logical(lgl)

On entry, if SCALING=true the intermediate tridiagonal matrix T is scaled before computing the eigenvalues.

The default is to scale the tridiagonal matrix.

INIT (INPUT, OPTIONAL) logical(lgl)

On entry, if INIT=true the initial intervals for the bisection steps for computing the eigenvalues of the intermediate tridiagonal matrix T are estimated from the eigenvalues of the intermediate tridiagonal matrix obtained from the Pal-Walker-Kahan algorithm.

The default is not to use the Pal-Walker-Kahan algorithm.

ABSTOL (INPUT, OPTIONAL) real(stnd)

On entry, the absolute tolerance for the eigenvalues. An eigenvalue (or cluster) is considered to be located if it has been determined to lie in an interval whose width is ABSTOL or less.

If ABSTOL is less than or equal to zero, or is not specified, then ULP * | T | will be used, where | T | means the 1-norm of T (T is the intermediate tridiagonal form of MAT) and ULP is the machine precision (distance from 1 to the next larger floating point number).

Eigenvalues will be computed most accurately when ABSTOL is set to the square root of the underflow threshold, sqrt(LAMCH(‘S’)), not zero.

LE (INPUT, OPTIONAL) integer(i4b)

On entry, LE specifies the number of eigenvalues which must be computed by the subroutine. However, on output, NEIG may be different than LE if multiple eigenvalues at index LE make unique selection impossible.

If:

SMALL=false, the subroutine computes the LE largest eigenvalues of MAT,

SMALL=true, the subroutine computes the LE smallest eigenvalues of MAT.

Only one of the optional arguments LE and THETA must be specified, otherwise the subroutine will stop with an error message.

LE must be greater than 0 and less or equal to size( EIGVAL ) .

The default is LE = size( EIGVAL ), e.g. all the eigenvalues are computed.

THETA (INPUT, OPTIONAL) real(stnd)

On entry,

if SMALL=false, THETA specifies that the eigenvalues which are greater or equal to THETA must be computed. If none of the eigenvalues are greater or equal to THETA, NEIG is set to zero and EIGVAL(:) to a quiet NAN.

if SMALL=true, THETA specifies that the eigenvalues which are less or equal to THETA must be computed. If none of the eigenvalues are smaller or equal to THETA, NEIG is set to zero and EIGVAL(:) to a quiet NAN.

Only one of the optional arguments LS and THETA must be specified, otherwise the subroutine will stop with an error message.

D_E (OUTPUT, OPTIONAL) real(stnd), dimension(:,2)

On exit, the first column of D_E contains the n diagonal elements of the intermediate tridiagonal form T of MAT. The n-1 first elements of the second column of D_E contains the n-1 subdiagonal elements of the intermediate tridiagonal form T of MAT. D_E(n,2) is arbitrary.

The shape of D_E must verify: size( D_E, 1 ) = n and size( D_E, 2 ) = 2

Further Details¶

This driver subroutine is adapted from the routine DSYEV in LAPACK.

`subroutine select_eigval_cmp3 ( matp, neig, eigval, small, failure, upper, sort, vector, scaling, init, abstol, le, theta, d_e )`¶

Purpose¶

SELECT_EIGVAL_CMP3 computes the largest or smallest eigenvalues of a n-by-n real symmetric matrix mat stored in packed form in a linear array MATP.

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

Q’ * mat * Q = T

,then the eigenvalues are computed by a bisection method.

Arguments¶

MATP (INPUT/OUTPUT) real(stnd), dimension(:)

On entry, the upper or lower triangle of the symmetric matrix mat, packed column-wise in a linear array. The j-th column of mat is stored in the array MATP as follows:

if UPPER = true, MATP(i + (j-1) * j/2) = mat(i,j) for 1<=i<=j;

if UPPER = false, MATP(i + (j-1) * (2 * n-j)/2) = mat(i,j) for j<=i<=n.

On exit:

If D_E is present : MATP is overwritten by the matrix Q as a product of elementary reflectors.

If D_E is absent : MATP is destroyed.

The size of MATP must verify: size( MATP ) = (n * (n+1)/2).

NEIG (OUTPUT) integer(i4b)

On output, NEIG specifies the number of eigenvalues which have been computed. Note that NEIG may be greater than the optional argument LE, if multiple eigenvalues at index LE make unique selection impossible.

If none of the optional arguments LE and THETA are used, NEIG is set to n and all the eigenvalues of MAT are computed.

EIGVAL (OUTPUT) real(stnd), dimension(:)

On exit, EIGVAL(1:NEIG) contains the first NEIG largest or smallest eigenvalues of MAT. The other values in EIGVAL (e.g. EIGVAL(NEIG+1:) ) are flagged by a quiet NAN.

The size of EIGVAL must verify: size( EIGVAL ) = n .

SMALL (INPUT) logical(lgl)

On entry:

SMALL = false : indicates that the m largest eigenvalues are desired.

SMALL = true : indicates that the m smallest eigenvalues are desired.

FAILURE (OUTPUT) logical(lgl)

On exit:

FAILURE = false : indicates successful exit and the bisection algorithm converged for all the computed eigenvalues to the desired accuracy ;

FAILURE = true : indicates that some or all of the eigenvalues failed to converge or were not computed. This is generally caused by unexpectedly inaccurate arithmetic.

UPPER (INPUT) logical(lgl)

Specifies whether the upper or lower triangular part of the symmetric matrix mat is stored in the linear array MATP. If:

UPPER = true : Upper triangle of mat is stored;

UPPER = false: Lower triangle of mat is stored.

SORT (INPUT, OPTIONAL) character

Sort the eigenvalues into ascending order if SORT = ‘A’ or ‘a’, or in descending order if SORT = ‘D’ or ‘d’. For other values of SORT nothing is done and EIGVAL(:NEIG) may not be sorted.

VECTOR (INPUT, OPTIONAL) logical(lgl)

On entry, if VECTOR is set to true, a vectorized version of the bisection algorithm is used to find the eigenvalues of the intermediate tridiagonal form T of MAT.

The default is VECTOR=false.

SCALING (INPUT, OPTIONAL) logical(lgl)

On entry, if SCALING=true the intermediate tridiagonal matrix T is scaled before computing the eigenvalues.

The default is to scale the tridiagonal matrix.

INIT (INPUT, OPTIONAL) logical(lgl)

On entry, if INIT=true the initial intervals for the bisection steps for computing the eigenvalues of the intermediate tridiagonal matrix T are estimated from the eigenvalues of the intermediate tridiagonal matrix obtained from the Pal-Walker-Kahan algorithm.

The default is not to use the Pal-Walker-Kahan algorithm.

ABSTOL (INPUT, OPTIONAL) real(stnd)

On entry, the absolute tolerance for the eigenvalues. An eigenvalue (or cluster) is considered to be located if it has been determined to lie in an interval whose width is ABSTOL or less.

If ABSTOL is less than or equal to zero, or is not specified, then ULP * | T | will be used, where | T | means the 1-norm of T (T is the intermediate tridiagonal form of mat) and ULP is the machine precision (distance from 1 to the next larger floating point number).

Eigenvalues will be computed most accurately when ABSTOL is set to the square root of the underflow threshold, sqrt(LAMCH(‘S’)), not zero.

LE (INPUT, OPTIONAL) integer(i4b)

On entry, LE specifies the number of eigenvalues which must be computed by the subroutine. However, on output, NEIG may be different than LE if multiple eigenvalues at index LE make unique selection impossible.

If:

SMALL=false, the subroutine computes the LE largest eigenvalues of MAT,

SMALL=true, the subroutine computes the LE smallest eigenvalues of MAT.

Only one of the optional arguments LE and THETA must be specified, otherwise the subroutine will stop with an error message.

LE must be greater than 0 and less or equal to size( EIGVAL ) .

The default is LE = size( EIGVAL ), e.g. all the eigenvalues are computed.

THETA (INPUT, OPTIONAL) real(stnd)

On entry:

if SMALL=false, THETA specifies that the eigenvalues which are greater or equal to THETA must be computed. If none of the eigenvalues are greater or equal to THETA, NEIG is set to zero and EIGVAL(:) to a quiet NAN.

if SMALL=true, THETA specifies that the eigenvalues which are less or equal to THETA must be computed. If none of the eigenvalues are smaller or equal to THETA, NEIG is set to zero and EIGVAL(:) to a quiet NAN.

Only one of the optional arguments LS and THETA must be specified, otherwise the subroutine will stop with an error message.

D_E (OUTPUT, OPTIONAL) real(stnd), dimension(:,2)

On exit, the first column of D_E contains the n diagonal elements of the intermediate tridiagonal form T of mat. The n-1 first elements of the second column of D_E contains the n-1 subdiagonal elements of the intermediate tridiagonal form T of mat. D_E(n,2) is arbitrary.

The shape of D_E must verify: size( D_E, 1 ) = n and size( D_E, 2 ) = 2

Further Details¶

This driver subroutine is adapted from the routine DSYEV in LAPACK.

`subroutine select_eigval_cmp3 ( matp, neig, eigval, small, failure, sort, vector, scaling, init, abstol, le, theta, d_e )`¶

Purpose¶

SELECT_EIGVAL_CMP3 computes the largest or smallest eigenvalues of a n-by-n real symmetric matrix mat stored in packed form in a linear array MATP.

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

Q’ * mat * Q = T

,then the eigenvalues are computed by a bisection method.

Arguments¶

MATP (INPUT/OUTPUT) real(stnd), dimension(:)

On entry, the upper triangle of the symmetric matrix mat, packed column-wise in a linear array. The j-th column of mat is stored in the array MATP as follows:

MATP(i + (j-1) * j/2) = mat(i,j) for 1<=i<=j;

On exit:

If D_E is present : MATP is overwritten by the matrix Q as a product of elementary reflectors.

If D_E is absent : MATP is destroyed.

The size of MATP must verify: size( MATP ) = (n * (n+1)/2).

NEIG (OUTPUT) integer(i4b)

On output, NEIG specifies the number of eigenvalues which have been computed. Note that NEIG may be greater than the optional argument LE, if multiple eigenvalues at index LE make unique selection impossible.

If none of the optional arguments LE and THETA are used, NEIG is set to n and all the eigenvalues of MAT are computed.

EIGVAL (OUTPUT) real(stnd), dimension(:)

On exit, EIGVAL(1:NEIG) contains the first NEIG largest or smallest eigenvalues of MAT. The other values in EIGVAL (e.g. EIGVAL(NEIG+1:) ) are flagged by a quiet NAN.

The size of EIGVAL must verify: size( EIGVAL ) = n .

SMALL (INPUT) logical(lgl)

On entry:

SMALL = false : indicates that the m largest eigenvalues are desired.

SMALL = true : indicates that the m smallest eigenvalues are desired.

FAILURE (OUTPUT) logical(lgl)

On exit:

FAILURE = false : indicates successful exit and the bisection algorithm converged for all the computed eigenvalues to the desired accuracy ;

FAILURE = true : indicates that some or all of the eigenvalues failed to converge or were not computed. This is generally caused by unexpectedly inaccurate arithmetic.

SORT (INPUT, OPTIONAL) character

Sort the eigenvalues into ascending order if SORT = ‘A’ or ‘a’, or in descending order if SORT = ‘D’ or ‘d’. For other values of SORT nothing is done and EIGVAL(:NEIG) may not be sorted.

VECTOR (INPUT, OPTIONAL) logical(lgl)

On entry, if VECTOR is set to true, a vectorized version of the bisection algorithm is used to find the eigenvalues of the intermediate tridiagonal form T of MAT.

The default is VECTOR=false.

SCALING (INPUT, OPTIONAL) logical(lgl)

On entry, if SCALING=true the intermediate tridiagonal matrix T is scaled before computing the eigenvalues.

The default is to scale the tridiagonal matrix.

INIT (INPUT, OPTIONAL) logical(lgl)

On entry, if INIT=true the initial intervals for the bisection steps for computing the eigenvalues of the intermediate tridiagonal matrix T are estimated from the eigenvalues of the intermediate tridiagonal matrix obtained from the Pal-Walker-Kahan algorithm.

The default is not to use the Pal-Walker-Kahan algorithm.

ABSTOL (INPUT, OPTIONAL) real(stnd)

On entry, the absolute tolerance for the eigenvalues. An eigenvalue (or cluster) is considered to be located if it has been determined to lie in an interval whose width is ABSTOL or less.

If ABSTOL is less than or equal to zero, or is not specified, then ULP * | T | will be used, where | T | means the 1-norm of T (T is the intermediate tridiagonal form of mat) and ULP is the machine precision (distance from 1 to the next larger floating point number).

Eigenvalues will be computed most accurately when ABSTOL is set to the square root of the underflow threshold, sqrt(LAMCH(‘S’)), not zero.

LE (INPUT, OPTIONAL) integer(i4b)

On entry, LE specifies the number of eigenvalues which must be computed by the subroutine. However, on output, NEIG may be different than LE if multiple eigenvalues at index LE make unique selection impossible.

If:

SMALL=false, the subroutine computes the LE largest eigenvalues of MAT,

SMALL=true, the subroutine computes the LE smallest eigenvalues of MAT.

Only one of the optional arguments LE and THETA must be specified, otherwise the subroutine will stop with an error message.

LE must be greater than 0 and less or equal to size( EIGVAL ) .

The default is LE = size( EIGVAL ), e.g. all the eigenvalues are computed.

THETA (INPUT, OPTIONAL) real(stnd)

On entry:

if SMALL=false, THETA specifies that the eigenvalues which are greater or equal to THETA must be computed. If none of the eigenvalues are greater or equal to THETA, NEIG is set to zero and EIGVAL(:) to a quiet NAN.

if SMALL=true, THETA specifies that the eigenvalues which are less or equal to THETA must be computed. If none of the eigenvalues are smaller or equal to THETA, NEIG is set to zero and EIGVAL(:) to a quiet NAN.

Only one of the optional arguments LS and THETA must be specified, otherwise the subroutine will stop with an error message.

D_E (OUTPUT, OPTIONAL) real(stnd), dimension(:,2)

On exit, the first column of D_E contains the n diagonal elements of the intermediate tridiagonal form T of mat. The n-1 first elements of the second column of D_E contains the n-1 subdiagonal elements of the intermediate tridiagonal form T of mat. D_E(n,2) is arbitrary.

The shape of D_E must verify: size( D_E, 1 ) = n and size( D_E, 2 ) = 2

Further Details¶

This driver subroutine is adapted from the routine DSYEV in LAPACK.

`subroutine reig_cmp ( mat, eigval, eigvec, failure, niter, nover, ortho, extd_samp, rng_alg, maxiter )`¶

Purpose¶

REIG_CMP computes approximations of the neig largest eigenvalues (in absolute magnitude) and associated eigenvectors of a full n-by-n real symmetric matrix MAT using randomized power, subspace or block Krylov iterations.

neig is the target rank of the partial EigenValue Decomposition (EVD), which is sought, and is equal to the size of the output real vector argument EIGVAL, i.e., neig = size( EIGVAL ).

Arguments¶

MAT (INPUT) real(stnd), dimension(:,:)

On entry, the n-by-n symmetric matrix MAT.

MAT is not modified by the routine.

EIGVAL (OUTPUT) real(stnd), dimension(:)

On exit, EIGVAL(:) contains the first top neig eigenvalues of MAT. The eigenvalues are given in decreasing order of absolute magnitude.

The size of EIGVAL must verify:

size( EIGVAL ) = neig <= size( MAT, 1 ) = size( MAT, 2 ) = n.

EIGVEC (OUTPUT) real(stnd), dimension(:,:)

On exit, the computed neig top eigenvectors. The eigenvector associated with the eigenvalue EIGVAL(j) is stored in the j-th column of EIGVEC.

The shape of EIGVEC must verify:

size( EIGVEC, 1 ) = size( MAT, 1 ) = size( MAT, 2 ) = n,

size( EIGVEC, 2 ) = size( EIGVEC ) = neig .

FAILURE (OUTPUT, OPTIONAL) logical(lgl)

On exit, if the optional logical argument FAILURE is present, a test of the accuracy of the computed partial EVD is performed and in that case:

FAILURE = false : indicates successful exit;

FAILURE = true : indicates that some of the computed eigenvalues and eigenvectors of MAT failed to converge in NITER iterations.

If FAILURE = true on exit, results are still useful, but some of the approximated eigen couplets have a poor accuracy.

NITER (INPUT, OPTIONAL) integer(i4b)

The number of randomized power, subspace or block Krylov iterations performed in the subroutine for computing the top neig eigen triplets. NITER must be positive or null.

By default, 10 randomized power, subspace or block Krylov iterations are performed.

NOVER (INPUT, OPTIONAL) integer(i4b)

The oversampling size used in the randomized power subspace or block Krylov iterations for computing the top neig eigen triplets.

NOVER must be positive or null and verifies the relationship:

NOVER + size( EIGVAL ) <= size( MAT, 1 ) = size( MAT, 2 ) = n

and is adjusted if necessary to verify this relationship in all cases.

See Further Details and the cited references for the meaning and usefulness of the oversampling size in the randomized power, subspace or block Krylov iterations.

By default, the oversampling size is set to 10.

ORTHO (INPUT, OPTIONAL) logical(lgl)

On entry, if:

ORTHO=true, orthonormalization is carried out between each step of the power or block Krylov iterations to avoid loss of accuracy due to rounding errors. This means that subspace iterations are used instead of power iterations;

ORTHO=false, orthonormalization is not performed.

The default is to use orthonormalization, e.g., ORTHO=true.

EXTD_SAMP (INPUT, OPTIONAL) logical(lgl)

The optional argument EXTD_SAMP determines if extended sampling (e.g., block Krylov iterations) is used or not for computing the top neig eigen triplets.

On entry, if:

EXTD_SAMP=true, block Krylov iterations are used;

EXTD_SAMP=false, power or subspace iterations are used.

The default is to use power or subspace iterations, e.g., EXTD_SAMP=false.

RNG_ALG (INPUT, OPTIONAL) integer(i4b)

On entry, a scalar integer to select the random (uniform) number generator used to build the random gaussian test matrix in the randomized EVD algorithm.

The possible values are:

ALG=1 : selects the Marsaglia’s KISS random number generator;

ALG=2 : selects the fast Marsaglia’s KISS random number generator;

ALG=3 : selects the L’Ecuyer’s LFSR113 random number generator;

ALG=4 : selects the Mersenne Twister random number generator;

ALG=5 : selects the maximally equidistributed Mersenne Twister random number generator;

ALG=6 : selects the extended precision of the Marsaglia’s KISS random number generator;

ALG=7 : selects the extended precision of the fast Marsaglia’s KISS random number generator;

ALG=8 : selects the extended precision of the L’Ecuyer’s LFSR113 random number generator.

ALG=9 : selects the extended precision of Mersenne Twister random number generator;

ALG=10 : selects the extended precision of maximally equidistributed Mersenne Twister random number generator;

For other values, the current random number generator and its current state are not changed. Note further, that, on exit, the current random number generator is not reset to its previous value before the call to REIG_CMP.

See the documentation of subroutine RANDOM_SEED_ in module Random for further information.

MAXITER (INPUT, OPTIONAL) integer(i4b)

MAXITER controls the maximum number of QR sweeps in the bidiagonal SVD phase of the SVD algorithm or in the QR phase of the EVD algorithm, which are used in the last phase of the randomized algorithm.

See description of suboutines SVD_CMP and EIG_CMP for further details about this optional argument.

Further Details¶

For a good introduction to randomized linear algebra, see the references (1) and (2).

The randomized subspace iteration was proposed in (3; see Algorithm 4.4) to compute an orthonormal matrix whose range approximates the range of MAT. An approximate partial EVD decomposition can then be computed using the aforementioned orthonormal matrix, see Algorithm 5.3 in (3).

The randomized block Krylov iterations for computing an approximate partial EVD was proposed in (4; see Algorithm 2). See also the reference (1).

For further details on randomized linear algebra, computing a partial EVD decomposition using randomized power, subspace or block Krylov iterations, see:

Martinsson, P.G., 2019:

Randomized methods for matrix computations. arXiv.1607.01649

Erichson, N.B., Voronin, S., Brunton, S.L., and Kutz, J.N., 2019:

Randomized matrix decompositions using R. arXiv.1608.02148

Halko, N., Martinsson, P.G., and Tropp, J.A., 2011:

Finding structure with randomness: probabilistic algorithms for constructing approximate matrix decompositions. SIAM Rev., 53, 217-288.

Musco, C., and Musco, C., 2015:

Randomized block krylov methods for stronger and faster approximate singular value decomposition. In Proceedings of the 28th International Conference on Neural Information Processing Systems, NIPS 15, pages 1396-1404, Cambridge, MA, USA, 2015. MIT Press.

Li, H.,Linderman, G.C., Szlam, A., Stanton, K.P., Kluger, Y., and Tygert, M., 2017:

Algorithm 971: An implementation of a randomized algorithm for principal component analysis. ACM Trans. Math. Softw. 43, 3, Article 28 (January 2017).

`function maxdiag_tinv_qr ( d, e, lambda )`¶

Purpose¶

This function computes the index of the element of maximum absolute value in the diagonal entries of ( T - LAMBDA * I )**(-1) where T is a symmetric tridiagonal matrix, I is the identity matrix and LAMBDA is a scalar.

Arguments¶

D (INPUT) real(stnd), dimension(:)

On entry, the diagonal elements of the tridiagonal matrix.

E (INPUT) real(stnd), dimension(:)

On entry, the n-1 subdiagonal elements of the tridiagonal matrix.

The size of E must verify: size( E ) = size( D ) - 1.

LAMBDA (INPUT) real(stnd)

On entry, the eigenvalue or shift used in the QR factorization.

Further Details¶

The diagonal entries of ( T - LAMBDA * I )**(-1) are computed by means of the QR factorization of ( T - LAMBDA * I ). For the latter computation, the semiseparable structure of ( T - LAMBDA * I )**(-1) is used, see the reference (1). Moreover, it is assumed that T is unreduced, but no check is done in the subroutine to verify this assumption.

This subroutine is adapted from the pseudo-code trace_Tinv given in the reference (1).

For further details, see:

Bini, D.A., Gemignani, L., and Tisseur, F., 2005:

The Ehrlich-Aberth method for the nonsymmetric tridiagonal eigenvalue problem. SIAM J. Matrix Anal. Appl., 27, 153-175.

Fernando, K.V., 1997:

On computing an eigenvector of a tridiagonal matrix. Part I: Basic results. Siam J. Matrix Anal. Appl., Vol. 18, 1013-1034.

Parlett, B.N., and Dhillon, I.S., 1997:

Fernando’s solution to Wilkinson’s problem: An application of double factorization. Linear Algebra and its Appl., 267, pp.247-279.

`function maxdiag_tinv_ldu ( d, e, lambda )`¶

Purpose¶

This function computes the index of the element of maximum absolute value in the diagonal entries of ( T - LAMBDA * I )**(-1) where T is a symmetric tridiagonal matrix, I is the identity matrix and LAMBDA is a scalar.

Arguments¶

D (INPUT) real(stnd), dimension(:)

On entry, the diagonal elements of the tridiagonal matrix.

E (INPUT) real(stnd), dimension(:)

On entry, the n-1 subdiagonal elements of the tridiagonal matrix.

The size of E must verify: size( E ) = size( D ) - 1.

LAMBDA (INPUT) real(stnd)

On entry, the eigenvalue or shift used.

Further Details¶

The diagonal entries of ( T - LAMBDA * I )**(-1) are computed by means of two triangular factorizations of ( T - LAMBDA * I ) of the forms L(+)D(+)U(+) and U(-)D(-)L(-) where L(+) and L(-) are unit lower bidiagonal, U(+) and U(-) are unit upper bidiagonal, and D(+) and D(-) are diagonal.

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

This subroutine is adapted from the references (1) and (2).

For further details, on Fernando’s method for computing eigenvectors of tridiagonal matrices, see:

Fernando, K.V., 1997:

On computing an eigenvector of a tridiagonal matrix. Part I: Basic results. Siam J. Matrix Anal. Appl., Vol. 18, pp. 1013-1034.

Parlett, B.N., and Dhillon, I.S., 1997:

Fernando’s solution to Wilkinson’s problem: An application of double factorization. Linear Algebra and its Appl., 267, pp.247-279.

`subroutine trid_qr_cmp ( d, e, lambda, cs, sn, diag, sup1, sup2, maxdiag_tinv )`¶

Purpose¶

TRID_QR_CMP factorizes the symmetric matrix T - LAMBDA * I, where T is an 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 assoociated with the eigenvalue LAMBDA, see the references (1), (2) and (3) for further details.

Arguments¶

D (INPUT) real(stnd), dimension(:)

On entry, the diagonal elements of the tridiagonal matrix.

E (INPUT) real(stnd), dimension(:)

On entry, the n-1 subdiagonal elements of the tridiagonal matrix.

The size of E must verify: size( E ) = size( D ) - 1.

LAMBDA (INPUT) real(stnd)

On entry, the eigenvalue or shift used in the QR factorization.

CS (OUTPUT) real(stnd), dimension(:)

On exit, the vector of the cosines coefficients of the chain of n-1 Givens rotations for the QR factorization of T - LAMBDA * I.

The size of CS must verify: size( CS ) = size( E ) = size( D ) - 1.

SN (OUTPUT) real(stnd), dimension(:)

On exit, the vector of the sines coefficients of the chain of n-1 Givens rotations for the QR factorization of T - LAMBDA * I.

The size of SN must verify: size( SN ) = size( E ) = size( D ) - 1.

DIAG (OUTPUT) real(stnd), dimension(:)

On exit, DIAG(:) contains the n diagonal elements of the upper triangular matrix R of the QR factorization of T - LAMBDA * I.

The size of DIAG must verify: size( DIAG ) = size( D ) = n .

SUP1 (OUTPUT) real(stnd), dimension(:)

On exit, SUP1(:n-1) contains the n-1 superdiagonal elements of the upper triangular matrix R of the QR factorization of T - LAMBDA * I, SUP1(n) is arbitrary.

The size of SUP1 must verify: size( SUP1 ) = size( D ) = n .

SUP2 (OUTPUT) real(stnd), dimension(:)

On exit, SUP2(:n-2) contains the n-2 second superdiagonal elements of the upper triangular matrix R of the QR factorization of T - LAMBDA * I, SUP2(n-1:n) is arbitrary.

The size of SUP2 must verify: size( SUP2 ) = size( D ) = n .

MAXDIAG_TINV (OUPTPUT) integer(i4b)

On exit, MAXDIAG_TINV is the index of the entry of maximum modulus in the main diagonal of ( T - LAMBDA * I )**(-1).

Further Details¶

The QR factorization of ( T - LAMBDA * I ) is obtained by means of n-1 unitary Givens rotations.

The diagonal entries of ( T - LAMBDA * I )**(-1) are computed by means of this QR factorization of ( T - LAMBDA * I ). For the latter computation, the semiseparable structure of ( T - LAMBDA * I )**(-1) is used, see the reference (1). Moreover, it is assumed that T is unreduced for computing the index of the entry of maximum absolute value in the diagonal of ( T - LAMBDA * I )**(-1), but no check is done in the subroutine to verify this assumption.

For further details, see:

Bini, D.A., Gemignani, L., and Tisseur, F., 2005:

The Ehrlich-Aberth method for the nonsymmetric tridiagonal eigenvalue problem. SIAM J. Matrix Anal. Appl., 27, 153-175.

Fernando, K.V., 1997:

On computing an eigenvector of a tridiagonal matrix. Part I: Basic results. Siam J. Matrix Anal. Appl., Vol. 18, pp. 1013-1034.

Parlett, B.N., and Dhillon, I.S., 1997:

Fernando’s solution to Wilkinson’s problem: An application of double factorization. Linear Algebra and its Appl., 267, pp.247-279.

`subroutine trid_qr_solve ( cs, sn, diag, sup1, sup2, y )`¶

Purpose¶

TRID_QR_SOLVE may be used to solve for x(:) the system of equations

x(:) * (T - LAMBDA * I) = scale * y(:)

, where T is an 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.

Arguments¶

CS (INPUT) real(stnd), dimension(:)

On entry, the vector of the cosines coefficients of the chain of n-1 Givens rotations for the QR factorization of T - LAMBDA * I as computed by TRID_QR_CMP or GK_QR_CMP.

The size of CS must verify: size( CS ) = size( Y ) - 1.

SN (INPUT) real(stnd), dimension(:)

On entry, the vector of the sines coefficients of the chain of n-1 Givens rotations for the QR factorization of GK - LAMBDA * I as computed by TRID_QR_CMP or GK_QR_CMP.

The size of SN must verify: size( SN ) = size( Y ) - 1.

DIAG (INPUT) real(stnd), dimension(:)

On entry, DIAG(:) contains the n diagonal elements of the upper triangular matrix R of the QR factorization of T - LAMBDA * I.

The size of DIAG must verify: size( DIAG ) = size( Y ) = n .

SUP1 (INPUT) real(stnd), dimension(:)

On entry, SUP1(:n-1) contains the n-1 superdiagonal elements of the upper triangular matrix R of the QR factorization of T - LAMBDA * I, SUP1(n) is arbitrary.

The size of SUP1 must verify: size( SUP1 ) = size( Y ) = n .

SUP2 (INPUT) real(stnd), dimension(:)

On entry, SUP2(:n-2) contains the n-2 second superdiagonal elements of the upper triangular matrix R of the QR factorization of T - LAMBDA * I, SUP2(n-1:n) is arbitrary.

The size of SUP2 must verify: size( SUP2 ) = size( Y ) = n .

Y (INPUT/OUTPUT) real(stnd), dimension(:)

On entry, the right hand side vector y. On exit, Y is overwritten the solution vector x.

The size of Y must verify: size( Y ) = n .

Further Details¶

For further details, see:

Fernando, K.V., 1997:

On computing an eigenvector of a tridiagonal matrix. Part I: Basic results. Siam J. Matrix Anal. Appl., Vol. 18, pp. 1013-1034.

Bini, D.A., Gemignani, L., and Tisseur, F., 2005:

The Ehrlich-Aberth method for the nonsymmetric tridiagonal eigenvalue problem. SIAM J. Matrix Anal. Appl., 27, 153-175.

`subroutine trid_cmp ( d, e, eigval, sub, diag, sup1, sup2, perm, tol )`¶

Purpose¶

TRID_CMP factorizes symmetric matrices of the form (T - EIGVAL(j) * I), where T is an n-by-n symmetric tridiagonal matrix and EIGVAL(j) is a scalar, as

T - EIGVAL(j) * I = P(j) * L(j) * U(j), for j=1, SIZE(EIGVAL)

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

The factorizations, for j=1, SIZE(EIGVAL), are obtained by Gaussian elimination with partial pivoting and implicit row scaling.

The parameters EIGVAL are included in the routine so that TRID_CMP may be used to obtain eigenvectors of T by inverse iteration.

Arguments¶

D (INPUT) real(stnd), dimension(:)

On entry, the diagonal elements of the symmetric tridiagonal matrix T.

E (INPUT) real(stnd), dimension(:)

On entry, the n-1 subdiagonal elements of the symmetric tridiagonal matrix T and E(n) is arbitrary .

The size of E must verify: size( E ) = size( D ) = n .

EIGVAL (INPUT) real(stnd), dimension(:)

On entry, selected eigenvalues of the symmetric tridiagonal matrix.

The size of EIGVAL must verify: size( EIGVAL ) <= size( D ) = n .

SUB (OUTPUT) real(stnd), dimension(:,:)

On exit, SUB(j,:n-1) contains the n-1 subdiagonal elements of the lower triangular matrix L(j) of the factorization of T - EIGVAL(j) * I, for j=1, SIZE(EIGVAL). SUB(:,n) is arbitrary .

The shape of SUB must verify:

size( SUB, 1 ) = size( EIGVAL ) ;

size( SUB, 2 ) = size( D ) = n .

DIAG (OUTPUT) real(stnd), dimension(:,:)

On exit, DIAG(j,:) contains the n diagonal elements of the upper triangular matrix U(j) of the factorization of T - EIGVAL(j) * I, for j=1, SIZE(EIGVAL).

The shape of DIAG must verify:

size( DIAG, 1 ) = size( EIGVAL ) ;

size( DIAG, 2 ) = size( D ) = n .

SUP1 (OUTPUT) real(stnd), dimension(:,:)

On exit, SUP1(j,:n-1) contains the n-1 superdiagonal elements of the upper triangular matrix U(j) of the factorization of T - EIGVAL(j) * I, for j=1, SIZE(EIGVAL). SUP1(:,n) is arbitrary .

The shape of SUP1 must verify:

size( SUP1, 1 ) = size( EIGVAL ) ;

size( SUP1, 2 ) = size( D ) = n .

SUP2 (OUTPUT) real(stnd), dimension(:,:)

On exit, SUP2(j,:n-2) contains the n-2 second superdiagonal elements of the upper triangular matrix U(j) of the factorization of T - EIGVAL(j) * I, for j=1, SIZE(EIGVAL). SUP2(:,n-1:n) is arbitrary .

The shape of SUP2 must verify:

size( SUP2, 1 ) = size( EIGVAL ) ;

size( SUP2, 2 ) = size( D ) = n .

PERM (OUTPUT) logical(lgl), dimension(:,:)

On exit, PERM(j,:n-1) contains details of the permutation matrix P(j). If an interchange occurred at the kth step of the elimination in the factorization of (T - EIGVAL(j) * I), then PERM(j,k) = true, otherwise PERM(j,k) = false. The element PERM(j,n) is set to true if there is an integer l such that

abs( u(j)(l,l) ).le. norm( (T - EIGVAL(j) * I)(l) ) * TOL,

where norm( A(l) ) denotes the sum of the absolute values of the lth row of the matrix A. If no such l exists then PERM(j,n) is returned as false. If PERM(j,n) is returned as true, then a diagonal element of U(j) is small, indicating that (T - EIGVAL(j) * I) is singular or nearly singular.

The shape of PERM must verify:

size( PERM, 1 ) = size( EIGVAL ) ;

size( PERM, 2 ) = size( D ) = n .

TOL (INPUT, OPTIONAL) real(stnd)

On entry, a relative tolerance used to indicate whether or not the matrices (T - EIGVAL(j) * I) are nearly singular. TOL should normally be choose as approximately the largest relative error in the elements of T. For example, if the elements of T are correct to about 4 significant figures, then TOL should be set to about 5 * 10**(-4).

If TOL is supplied as less than eps, where eps is the relative machine precision, then the value eps is used in place of TOL.

Further Details¶

This subroutine is adapted from the routine DLAGTF in LAPACK.

`subroutine trid_cmp ( d, e, eigval, sub, diag, sup1, sup2, perm, tol )`¶

Purpose¶

TRID_CMP factorizes the symmetric matrix (T - EIGVAL * I), where T is an n-by-n symmetric tridiagonal matrix and 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 implicit row scaling.

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

Arguments¶

D (INPUT) real(stnd), dimension(:)

On entry, the diagonal elements of the symmetric tridiagonal matrix T.

E (INPUT) real(stnd), dimension(:)

On entry, the n-1 subdiagonal elements of the symmetric tridiagonal matrix T and E(n) is arbitrary.

The size of E must verify: size( E ) = size( D ) = n .

EIGVAL (INPUT) real(stnd)

On entry, an eigenvalue of the symmetric tridiagonal matrix.

SUB (OUTPUT) real(stnd), dimension(:)

On exit, SUB(:n-1) contains the n-1 subdiagonal elements of the lower triangular matrix L of the factorization of T - EIGVAL * I, SUB(n) is arbitrary .

The size of SUB must verify: size( SUB ) = size( D ) = n .

DIAG (OUTPUT) real(stnd), dimension(:)

On exit, DIAG(:) contains the n diagonal elements of the upper triangular matrix U of the factorization of T - EIGVAL * I.

The size of DIAG must verify: size( DIAG ) = size( D ) = n .

SUP1 (OUTPUT) real(stnd), dimension(:)

On exit, SUP1(:n-1) contains the n-1 superdiagonal elements of the upper triangular matrix U of the factorization of T - EIGVAL * I, SUP1(n) is arbitrary.

The size of SUP1 must verify: size( SUP1 ) = size( D ) = n .

SUP2 (OUTPUT) real(stnd), dimension(:)

On exit, SUP2(:n-2) contains the n-2 second superdiagonal elements of the upper triangular matrix U of the factorization of T - EIGVAL * I, SUP2(n-1:n) is arbitrary .

The size of SUP2 must verify: size( SUP2 ) = size( D ) = n .

PERM (OUTPUT) logical(lgl), dimension(:)

On exit, PERM(:n-1) contains details of the permutation matrix P(j). If an interchange occurred at the kth step of the elimination in the factorization of (T - EIGVAL(j) * I), then PERM(k) = true, otherwise PERM(k) = false. The element PERM(n) is set to true if there is an integer l such that

abs( u(l,l) ).le. norm( (T - EIGVAL * I)(l) ) * TOL,

where norm( A(l) ) denotes the sum of the absolute values of the lth row of the matrix A. If no such l exists then PERM(n) is returned as false. If PERM(n) is returned as true, then a diagonal element of U is small, indicating that (T - EIGVAL * I) is singular or nearly singular.

The size of PERM must verify: size( PERM ) = size( D ) = n .

TOL (INPUT, OPTIONAL) real(stnd)

On entry, a relative tolerance used to indicate whether or not the matrix (T - EIGVAL * I) is nearly singular. TOL should normally be choose as approximately the largest relative error in the elements of T. For example, if the elements of T are correct to about 4 significant figures, then TOL should be set to about 5 * 10**(-4).

If TOL is supplied as less than eps, where eps is the relative machine precision, then the value eps is used in place of TOL.

Further Details¶

This subroutine is adapted from the routine DLAGTF in LAPACK.

`subroutine trid_cmp2 ( d, e, eigval, sub, diag, sup1, sup2, perm )`¶

Purpose¶

TRID_CMP2 factorizes symmetric matrices of the form (T - EIGVAL(j) * I), where T is an n-by-n symmetric tridiagonal matrix and EIGVAL(j) is a scalar, as

T - EIGVAL(j) * I = P(j) * L(j) * U(j), for j=1, SIZE(EIGVAL)

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

The factorizations, for j=1, SIZE(EIGVAL), are obtained by Gaussian elimination with partial pivoting and row interchanges.

The parameters EIGVAL are included in the routine so that TRID_CMP2 may be used to obtain eigenvectors of T by inverse iteration.

Arguments¶

D (INPUT) real(stnd), dimension(:)

On entry, the diagonal elements of the symmetric tridiagonal matrix T.

E (INPUT) real(stnd), dimension(:)

On entry, the n-1 subdiagonal elements of the symmetric tridiagonal matrix T and E(n) is arbitrary .

The size of E must verify: size( E ) = size( D ) = n .

EIGVAL (INPUT) real(stnd), dimension(:)

On entry, selected eigenvalues of the symmetric tridiagonal matrix.

The size of EIGVAL must verify: size( EIGVAL ) <= size( D ) = n .

SUB (OUTPUT) real(stnd), dimension(:,:)

On exit, SUB(j,:n-1) contains the n-1 subdiagonal elements of the lower triangular matrix L(j) of the factorization of T - EIGVAL(j) * I, for j=1, SIZE(EIGVAL). SUB(:,n) is arbitrary .

The shape of SUB must verify:

size( SUB, 1 ) = size( EIGVAL ) ;

size( SUB, 2 ) = size( D ) = n .

DIAG (OUTPUT) real(stnd), dimension(:,:)

On exit, DIAG(j,:) contains the n diagonal elements of the upper triangular matrix U(j) of the factorization of T - EIGVAL(j) * I, for j=1, SIZE(EIGVAL).

The shape of DIAG must verify:

size( DIAG, 1 ) = size( EIGVAL ) ;

size( DIAG, 2 ) = size( D ) = n .

SUP1 (OUTPUT) real(stnd), dimension(:,:)

On exit, SUP1(j,:n-1) contains the n-1 superdiagonal elements of the upper triangular matrix U(j) of the factorization of T - EIGVAL(j) * I, for j=1, SIZE(EIGVAL). SUP1(:,n) is arbitrary .

The shape of SUP1 must verify:

size( SUP1, 1 ) = size( EIGVAL ) ;

size( SUP1, 2 ) = size( D ) = n .

SUP2 (OUTPUT) real(stnd), dimension(:,:)

On exit, SUP2(j,:n-2) contains the n-2 second superdiagonal elements of the upper triangular matrix U(j) of the factorization of T - EIGVAL(j) * I, for j=1, SIZE(EIGVAL). SUP2(:,n-1:n) is arbitrary .

The shape of SUP2 must verify:

size( SUP2, 1 ) = size( EIGVAL ) ;

size( SUP2, 2 ) = size( D ) = n .

PERM (OUTPUT) logical(lgl), dimension(:,:)

On exit, PERM(j,:n-1) contains details of the permutation matrix P(j). If an interchange occurred at the kth step of the elimination in the factorization of (T - EIGVAL(j) * I), then PERM(j,k) = true, otherwise PERM(j,k) = false. PERM(:,n) is arbitrary .

The shape of PERM must verify:

size( PERM, 1 ) = size( EIGVAL ) ;

size( PERM, 2 ) = size( D ) = n .

Further Details¶

TRID_CMP2 is a simplified version of TRID_CMP. This subroutine is adapted from the routine DGTTRF in LAPACK.

`subroutine trid_cmp2 ( d, e, eigval, sub, diag, sup1, sup2, perm )`¶

Purpose¶

TRID_CMP2 factorizes the symmetric matrix (T - EIGVAL * I), where T is an n by n symmetric tridiagonal matrix and 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.

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

Arguments¶

D (INPUT) real(stnd), dimension(:)

On entry, the diagonal elements of the symmetric tridiagonal matrix T.

E (INPUT) real(stnd), dimension(:)

On entry, the n-1 subdiagonal elements of the symmetric tridiagonal matrix T and E(n) is arbitrary .

The size of E must verify: size( E ) = size( D ) = n .

EIGVAL (INPUT) real(stnd)

On entry, an eigenvalue of the symmetric tridiagonal matrix.

SUB (OUTPUT) real(stnd), dimension(:)

On exit, SUB(:n-1) contains the n-1 subdiagonal elements of the lower triangular matrix L of the factorization of T - EIGVAL * I, SUB(n) is arbitrary .

The size of SUB must verify: size( SUB ) = size( D ) = n .

DIAG (OUTPUT) real(stnd), dimension(:)

On exit, DIAG(:) contains the n diagonal elements of the upper triangular matrix U of the factorization of T - EIGVAL * I.

The size of DIAG must verify: size( DIAG ) = size( D ) = n .

SUP1 (OUTPUT) real(stnd), dimension(:)

On exit, SUP1(:n-1) contains the n-1 superdiagonal elements of the upper triangular matrix U of the factorization of T - EIGVAL * I, SUP1(n) is arbitrary .

The size of SUP1 must verify: size( SUP1 ) = size( D ) = n .

SUP2 (OUTPUT) real(stnd), dimension(:)

On exit, SUP2(:n-2) contains the n-2 second superdiagonal elements of the upper triangular matrix U of the factorization of T - EIGVAL * I, SUP2(n-1:n) is arbitrary .

The size of SUP2 must verify: size( SUP2 ) = size( D ) = n .

PERM (OUTPUT) logical(lgl), dimension(:)

On exit, PERM(:n-1) contains details of the permutation matrix P(j). If an interchange occurred at the kth step of the elimination in the factorization of (T - EIGVAL(j) * I), then PERM(k) = true, otherwise PERM(k) = false. PERM(n) is arbitrary .

The size of PERM must verify: size( PERM ) = size( D ) = n .

Further Details¶

TRID_CMP2 is a simplified version of TRID_CMP. This subroutine is adapted from the routine DGTTRF in LAPACK.

`subroutine trid_solve ( sub, diag, sup1, sup2, perm, y )`¶

Purpose¶

TRID_SOLVE may be used to solve systems of equations of the form

x(j,:) * (T - EIGVAL(j) * I) = scale(j) * y(j,:), for j=1, SIZE(EIGVAL)

, where T is an n by n symmetric tridiagonal matrix, EIGVAL(j) and scale(j) are scalars, for x(j,:) for j=1, SIZE(EIGVAL), following the factorization of (T - EIGVAL(j) * I) by TRID_CMP or TRID_CMP2 as

T - EIGVAL(j) * I = P(j) * L(j) * U(j), for j=1, SIZE(EIGVAL)

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

The matrices (T - EIGVAL(j) * I) are 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 scalars, scale(j) for j=1, SIZE(EIGVAL), are not output by this routine since this routine being intended for use in applications such as inverse iteration.

Arguments¶

SUB (INPUT) real(stnd), dimension(:,:)

On entry, SUB(j,:n-1) contains the n-1 subdiagonal elements of the lower triangular matrix L(j) of the factorization of T - EIGVAL(j) * I, for j=1, SIZE(EIGVAL). SUB(:,n) is arbitrary .

The shape of SUB must verify:

size( SUB, 1 ) = size( Y, 1 ) = size( EIGVAL ) ;

size( SUB, 2 ) = size( Y, 2 ) = n .

DIAG (INPUT) real(stnd), dimension(:,:)

On entry, DIAG(j,:) contains the n diagonal elements of the upper triangular matrix U(j) of the factorization of T - EIGVAL(j) * I, for j=1, SIZE(EIGVAL).

The shape of DIAG must verify:

size( DIAG, 1 ) = size( Y, 1 ) = size( EIGVAL ) ;

size( DIAG, 2 ) = size( Y, 2 ) = n .

SUP1 (INPUT) real(stnd), dimension(:,:)

On entry, SUP1(j,:n-1) contains the n-1 superdiagonal elements of the upper triangular matrix U(j) of the factorization of T - EIGVAL(j) * I, for j=1, SIZE(EIGVAL). SUP1(:,n) is arbitrary.

The shape of SUP1 must verify:

size( SUP1, 1 ) = size( Y, 1 ) = size( EIGVAL ) ;

size( SUP1, 2 ) = size( Y, 2 ) = n .

SUP2 (INPUT) real(stnd), dimension(:,:)

On entry, SUP2(j,:n-2) contains the n-2 second superdiagonal elements of the upper triangular matrix U(j) of the factorization of T - EIGVAL(j) * I, for j=1, SIZE(EIGVAL). SUP2(:,n-1:n) is arbitrary.

The shape of SUP2 must verify:

size( SUP2, 1 ) = size( Y, 1 ) = size( EIGVAL ) ;

size( SUP2, 2 ) = size( Y, 2 ) = n .

PERM (INPUT) logical(lgl), dimension(:,:)

On entry, PERM(j,:n-1) contains details of the permutation matrix P(j). If an interchange occurred at the kth step of the elimination in the factorization of (T - EIGVAL(j) * I), then PERM(j,k) = true, otherwise PERM(j,k) = false. PERM(:,n) is arbitrary .

The shape of PERM must verify:

size( PERM, 1 ) = size( Y, 1 ) = size( EIGVAL ) ;

size( PERM, 2 ) = size( Y, 2 ) = n .

Y (INPUT/OUTPUT) real(stnd), dimension(:,:)

On entry, the right hand side matrix y. On exit, Y is overwritten the solution matrix x.

The shape of Y must verify:

size( Y, 1 ) = size( EIGVAL ) ;

size( Y, 2 ) = n .

Further Details¶

This subroutine is adapted from the routine DLAGTS in LAPACK.

`subroutine trid_solve ( sub, diag, sup1, sup2, perm, y )`¶

Purpose¶

TRID_SOLVE may be used to solve the system of equations

x(:) * (T - EIGVAL * I) = scale * y(:)

, where T is an n by n symmetric tridiagonal matrix, EIGVAL and scale are scalars, for x(:), 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. The scalar, scale, is not output by this routine since this routine being intended for use in applications such as inverse iteration.

Arguments¶

SUB (INPUT) real(stnd), dimension(:)

On entry, SUB(:n-1) contains the n-1 subdiagonal elements of the lower triangular matrix L of the factorization of T - EIGVAL * I, SUB(n) is arbitrary .

The size of SUB must verify: size( SUB ) = size( Y ) = n .

DIAG (INPUT) real(stnd), dimension(:)

On entry, DIAG(:) contains the n diagonal elements of the upper triangular matrix U of the factorization of T - EIGVAL * I.

The shape of DIAG must verify: size( DIAG ) = size( Y ) = n .

SUP1 (INPUT) real(stnd), dimension(:)

On entry, SUP1(:n-1) contains the n-1 superdiagonal elements of the upper triangular matrix U of the factorization of T - EIGVAL * I, SUP1(n) is arbitrary.

The shape of SUP1 must verify: size( SUP1 ) = size( Y ) = n .

SUP2 (INPUT) real(stnd), dimension(:)

On entry, SUP2(:n-2) contains the n-2 second superdiagonal elements of the upper triangular matrix U of the factorization of T - EIGVAL * I, SUP2(n-1:n) is arbitrary.

The shape of SUP2 must verify: size( SUP2 ) = size( Y ) = n .

PERM (INPUT) logical(lgl), dimension(:)

On entry, PERM(:n-1) contains details of the permutation matrix P. If an interchange occurred at the kth step of the elimination in the factorization of (T - EIGVAL * I), then PERM(k) = true, otherwise PERM(k) = false. PERM(n) is arbitrary .

The shape of PERM must verify: size( PERM ) = size( Y ) = n .

Y (INPUT/OUTPUT) real(stnd), dimension(:)

On entry, the right hand side vector y. On exit, Y is overwritten the solution vector x.

The shape of Y must verify: size( Y ) = n .

Further Details¶

This subroutine is adapted from the routine DLAGTS in LAPACK.

`subroutine trid_inviter ( d, e, eigval, eigvec, failure, maxiter, scaling, initvec )`¶

Purpose¶

TRID_INVITER computes an eigenvector of a real n-by-n symmetric tridiagonal matrix T corresponding to a specified eigenvalue, by combining Fernando’s method for computing an eigenvector of a real n-by-n symmetric tridiagonal matrix and an inverse iteration algorithm.

Arguments¶

D (INPUT) real(stnd), dimension(:)

On entry, the diagonal elements of the symmetric tridiagonal matrix T.

E (INPUT) real(stnd), dimension(:)

On entry, the n-1 subdiagonal elements of the symmetric tridiagonal matrix T and E(n) is arbitrary .

The size of E must verify: size( E ) = size( D ) = n .

EIGVAL (INPUT) real(stnd)

On entry, an eigenvalue of the symmetric tridiagonal matrix.

EIGVEC (OUTPUT) real(stnd), dimension(:)

On exit, the computed eigenvector.

The shape of EIGVEC must verify: size( EIGVEC ) = size( D ) = n .

FAILURE (OUTPUT) logical(lgl)

On exit:

FAILURE = false indicates successful exit.

FAILURE = true indicates that the eigenvector failed to converge in MAXITER iterations.

MAXITER (INPUT, OPTIONAL) integer(i4b)

The number of inverse iterations performed in the subroutine.

By default, 2 inverse iterations are performed.

SCALING (INPUT, OPTIONAL) logical(lgl)

On entry, if SCALING=true the tridiagonal matrix T is scaled before computing the eigenvector.

The default is to scale the tridiagonal matrix.

INITVEC (INPUT, OPTIONAL) logical(lgl)

On entry, if INITVEC=true a Fernando vector is used to start the inverse iteration process; if INITVEC=false a random uniform starting vector is used.

For unreduced tridiagonal matrices, the default is to use a Fernando starting vector. For reduced tridiagonal matrices, the default is to use a random uniform starting vector.

Further Details¶

TRID_INVITER uses Fernando’s method for computing a first estimate of an eigenvector corresponding to an approximate eigenvalue of a real n-by-n symmetric tridiagonal matrix T (by default, only if the input tridiagonal matrix T is unreduced).

This approximate eigenvector is then refined (or computed if Fernando’s method is not used) using an inverse iteration algorithm.

For further details, on Fernando’s method for computing eigenvectors of tridiagonal matrices or inverse iteration, see:

Fernando, K.V., 1997:

On computing an eigenvector of a tridiagonal matrix. Part I: Basic results. Siam J. Matrix Anal. Appl., Vol. 18, pp. 1013-1034.

Parlett, B.N., and Dhillon, I.S., 1997:

Fernando’s solution to Wilkinson’s problem: An application of double factorization. Linear Algebra and its Appl., 267, pp.247-279.

Golub, G.H., and Van Loan, C.F., 1996:

Matrix Computations. 3rd ed. The Johns Hopkins University Press, Baltimore.

Bini, D.A., Gemignani, L., and Tisseur, F., 2005:

The Ehrlich-Aberth method for the nonsymmetric tridiagonal eigenvalue problem. SIAM J. Matrix Anal. Appl., 27, 153-175.

`subroutine trid_inviter ( d, e, eigval, eigvec, failure, maxiter, ortho, backward_sweep, scaling, initvec )`¶

Purpose¶

TRID_INVITER computes the eigenvectors of a real n-by-n symmetric tridiagonal matrix T corresponding to specified eigenvalues, by combining Fernando’s method for computing (selected) eigenvectors of a real n-by-n symmetric tridiagonal matrix and an inverse iteration algorithm.

Arguments¶

D (INPUT) real(stnd), dimension(:)

On entry, the diagonal elements of the symmetric tridiagonal matrix T.

E (INPUT) real(stnd), dimension(:)

On entry, the n-1 subdiagonal elements of the symmetric tridiagonal matrix T and E(n) is arbitrary .

The size of E must verify: size( E ) = size( D ) = n .

EIGVAL (INPUT) real(stnd), dimension(:)

On entry, selected eigenvalues of the symmetric tridiagonal matrix. The eigenvalues must be given in decreasing order.

The size of EIGVAL must verify: size( EIGVAL ) <= size( D ) = n .

EIGVEC (OUTPUT) real(stnd), dimension(:,:)

On exit, the computed eigenvectors. The eigenvector associated with the eigenvalue EIGVAL(j) is stored in the j-th column of EIGVEC.

The shape of EIGVEC must verify:

size( EIGVEC, 1 ) = size( D ) = n ,

size( EIGVEC, 2 ) = size( EIGVAL ) .

FAILURE (OUTPUT) logical(lgl)

On exit:

FAILURE = false indicates successful exit.

FAILURE = true indicates that some eigenvectors failed to converge in MAXITER iterations.

MAXITER (INPUT, OPTIONAL) integer(i4b)

The number of inverse iterations performed in the subroutine. By default, 2 inverse iterations are performed for all the eigenvectors.

ORTHO (INPUT, OPTIONAL) logical(lgl)

On entry, if:

ORTHO=true, all the eigenvectors are orthogonalized by the Modified Gram-Schmidt or QR algorithm;

ORTHO=false, the eigenvectors are not orthogonalized by the Modified Gram-Schmidt or QR algorithm.

The default is to orthogonalize the eigenvectors only for the eigenvalues, which are not well-separated.

BACKWARD_SWEEP (INPUT, OPTIONAL) logical(lgl)

On entry, if:

BACKWARD_SWEEP=true and the eigenvectors are orthogonalized by the modified Gram-Schmidt algorithm, a backward sweep of the modified Gram-Schmidt algorithm is also performed;

BACKWARD_SWEEP=false a backward sweep is not performed.

The default is not to perform a backward sweep of the modified Gram-Schmidt algorithm.

SCALING (INPUT, OPTIONAL) logical(lgl)

On entry, if:

SCALING=true, the tridiagonal matrix T is scaled before computing the eigenvectors;

SCALING=false, the tridiagonal matrix T is not scaled.

The default is to scale the tridiagonal matrix.

INITVEC (INPUT, OPTIONAL) logical(lgl)

On entry, if:

INITVEC=true, Fernando vectors are used to start the inverse iteration process;

INITVEC=false, random uniform starting vectors are used.

For unreduced tridiagonal matrices, the default is to use Fernando starting vectors if the eigenvalues are well-separated and random uniform starting vectors otherwise. For reduced tridiagonal matrices, the default is to use random uniform starting vectors.

Further Details¶

TRID_INVITER uses Fernando’s method for computing a first estimate of (selected) eigenvectors corresponding to (selected) approximate eigenvalues of a real n-by-n symmetric tridiagonal matrix T (by default, only for the eigenvalues which are well separated and if the input tridiagonal matrix T is unreduced).

These approximate eigenvectors are then refined (or computed if Fernando’s method is not used) using an inverse iteration algorithm for all the eigenvalues at one step. The eigenvectors are then orthogonalized by the Modified Gram-Schmidt or QR algorithm if the eigenvalues are not well-separated.

TRID_INVITER may fail if some the eigenvalues specified in parameter EIGVAL are nearly identical.

For further details, on Fernando’s method for computing eigenvectors of tridiagonal matrices or inverse iteration, see:

Fernando, K.V., 1997:

On computing an eigenvector of a tridiagonal matrix. Part I: Basic results. Siam J. Matrix Anal. Appl., Vol. 18, pp. 1013-1034.

Parlett, B.N., and Dhillon, I.S., 1997:

Fernando’s solution to Wilkinson’s problem: An application of double factorization. Linear Algebra and its Appl., 267, pp.247-279.

Golub, G.H., and Van Loan, C.F., 1996:

Matrix Computations. 3rd ed. The Johns Hopkins University Press, Baltimore.

Bini, D.A., Gemignani, L., and Tisseur, F., 2005:

The Ehrlich-Aberth method for the nonsymmetric tridiagonal eigenvalue problem. SIAM J. Matrix Anal. Appl., 27, 153-175.

`subroutine trid_inviter ( d, e, eigval, eigvec, failure, mat, maxiter, ortho, backward_sweep, scaling, initvec )`¶

Purpose¶

TRID_INVITER computes the eigenvectors of a full real n-by-n symmetric matrix MAT corresponding to specified eigenvalues, by combining Fernando’s method for computing (selected) eigenvectors of a real n-by-n symmetric tridiagonal matrix and inverse iteration followed by a back-transformation procedure.

It is required that the original symmetric matrix MAT has been reduced to symmetric tridiagonal form T by an orthogonal similarity transformation:

Q’ * MAT * Q = T

with a call to SYMTRID_CMP with parameter STORE_Q set to true, before calling TRID_INVITER.

Arguments¶

D (INPUT) real(stnd), dimension(:)

On entry, the diagonal elements of the symmetric tridiagonal form T of MAT.

E (INPUT) real(stnd), dimension(:)

On entry, the n-1 subdiagonal elements of the symmetric tridiagonal form T of MAT. E(n) is arbitrary .

The size of E must verify: size( E ) = size( D ) = n .

EIGVAL (INPUT) real(stnd), dimension(:)

On entry, selected eigenvalues of the symmetric matrix MAT. The eigenvalues must be given in decreasing order.

The size of EIGVAL must verify: size( EIGVAL ) <= size( D ) = n .

EIGVEC (OUTPUT) real(stnd), dimension(:,:)

On exit, the computed eigenvectors. The eigenvector associated with the eigenvalue EIGVAL(j) is stored in the j-th column of EIGVEC.

The shape of EIGVEC must verify:

size( EIGVEC, 1 ) = size( D ) = n ;

size( EIGVEC, 2 ) = size( EIGVAL ) .

FAILURE (OUTPUT) logical(lgl)

On exit:

FAILURE = false indicates successful exit.

FAILURE = true indicates that some eigenvectors failed to converge in MAXITER iterations.

MAT (INPUT) real(stnd), dimension(:,:)

On entry, the vectors and the scalars which define the elementary reflectors used to reduce the full real n-by-n symmetric matrix MAT to symmetric tridiagonal form T, as returned by SYMTRID_CMP with STORE_Q=true, in its argument MAT. MAT is not modified by the routine.

Back-transformation is used to find the selected eigenvectors of the original matrix MAT and these eigenvectors are stored in argument EIGVEC.

The shape of MAT must verify:

size( MAT, 1 ) = size( D ) = n ;

size( MAT, 2 ) = size( D ) = n .

MAXITER (INPUT, OPTIONAL) integer(i4b)

The number of inverse iterations performed in the subroutine. By default, 2 inverse iterations are performed for all the eigenvectors of the tridiagonal matrix T.

ORTHO (INPUT, OPTIONAL) logical(lgl)

On entry, if:

ORTHO=true, all the eigenvectors are orthogonalized by the Modified Gram-Schmidt or QR algorithm;

ORTHO=false, the eigenvectors are not orthogonalized by the Modified Gram-Schmidt or QR algorithm.

The default is to orthogonalize the eigenvectors only if the eigenvalues are not well-separated.

BACKWARD_SWEEP (INPUT, OPTIONAL) logical(lgl)

On entry, if:

BACKWARD_SWEEP=true and the eigenvectors are orthogonalized by the modified Gram-Schmidt algorithm, a backward sweep of the modified Gram-Schmidt algorithm is also performed;

BACKWARD_SWEEP=false a backward sweep is not performed.

The default is not to perform a backward sweep of the modified Gram-Schmidt algorithm.

SCALING (INPUT, OPTIONAL) logical(lgl)

On entry, if:

SCALING=true, the tridiagonal matrix T is scaled before computing the eigenvectors;

SCALING=false, the tridiagonal matrix T is not scaled.

The default is to scale the tridiagonal matrix.

INITVEC (INPUT, OPTIONAL) logical(lgl)

On entry, if:

INITVEC=true, Fernando vectors are used to start the inverse iteration process;

INITVEC=false, random uniform starting vectors are used.

For unreduced tridiagonal matrices, the default is to use Fernando starting vectors if the eigenvalues are well-separated and random uniform starting vectors otherwise. For reduced tridiagonal matrices, the default is to use random uniform starting vectors.

Further Details¶

TRID_INVITER uses Fernando’s method for computing a first estimate of (selected) eigenvectors corresponding to (selected) approximate eigenvalues of a real n-by-n symmetric tridiagonal matrix T (by default, only for the eigenvalues which are well separated and if the input tridiagonal matrix T is unreduced). See the references (1), (2) and (4) for details.

These approximate eigenvectors are then refined (or computed if Fernando’s method is not used) using an inverse iteration algorithm for all the eigenvalues at one step. See the reference (3) for details.

The eigenvectors are then orthogonalized by the Modified Gram-Schmidt or QR algorithm if clusters of eigenvalues are present, in a second step.

In a last step, the corresponding (selected) eigenvectors of the full real n-by-n symmetric matrix MAT are computed by a blocked back-transformation algorithm with the Householder transformations used to reduce the full real n-by-n symmetric matrix MAT to symmetric tridiagonal form T (see the references (5) and (6)).

Furthermore, the computation of the eigenvectors is parallelized if OPENMP is used.

TRID_INVITER may fail if some the eigenvalues specified in parameter EIGVAL are nearly identical.

For further details on Fernando’s method or inverse iteration for computing eigenvectors of tridiagonal matrices or the blocked back-transformation algorithm, see:

Fernando, K.V., 1997:

On computing an eigenvector of a tridiagonal matrix. Part I: Basic results. Siam J. Matrix Anal. Appl., Vol. 18, pp. 1013-1034.

Parlett, B.N., and Dhillon, I.S., 1997:

Fernando’s solution to Wilkinson’s problem: An application of double factorization. Linear Algebra and its Appl., 267, pp.247-279.

Golub, G.H., and Van Loan, C.F., 1996:

Matrix Computations. 3rd ed. The Johns Hopkins University Press, Baltimore.

Bini, D.A., Gemignani, L., and Tisseur, F., 2005:

The Ehrlich-Aberth method for the nonsymmetric tridiagonal eigenvalue problem. SIAM J. Matrix Anal. Appl., 27, 153-175.

Dongarra, J.J., Sorensen, D.C., and Hammarling, S.J., 1989:

Block reduction of matrices to condensed form for eigenvalue computations. J. of Computational and Applied Mathematics, Vol. 27, pp. 215-227.

Walker, H.F., 1988:

Implementation of the GMRES method using Householder transformations. Siam J. Sci. Stat. Comput., Vol. 9, No 1, pp. 152-163.

`subroutine trid_inviter ( d, e, eigval, eigvec, failure, matp, maxiter, ortho, backward_sweep, scaling, initvec )`¶

Purpose¶

TRID_INVITER computes the eigenvectors of a full real n-by-n symmetric matrix MAT, packed column-wise in a linear array MATP, corresponding to specified eigenvalues, using Fernando’s method for computing (selected) eigenvectors of a real n-by-n symmetric tridiagonal matrix and inverse iteration, followed by a back-transformation procedure.

It is required that the original packed symmetric matrix MAT has been reduced to symmetric tridiagonal form T by an orthogonal similarity transformation:

Q’ * MAT * Q = T

with a call to SYMTRID_CMP with parameter STORE_Q set to true, before calling TRID_INVITER.

Arguments¶

D (INPUT) real(stnd), dimension(:)

On entry, the diagonal elements of the symmetric tridiagonal form T of MAT.

E (INPUT) real(stnd), dimension(:)

On entry, the n-1 subdiagonal elements of the symmetric tridiagonal form T of MAT. E(n) is arbitrary .

The size of E must verify: size( E ) = size( D ) = n .

EIGVAL (INPUT) real(stnd), dimension(:)

On entry, selected eigenvalues of the symmetric matrix MAT. The eigenvalues must be given in decreasing order.

The size of EIGVAL must verify: size( EIGVAL ) <= size( D ) = n .

EIGVEC (OUTPUT) real(stnd), dimension(:,:)

On exit, the computed eigenvectors. The eigenvector associated with the eigenvalue EIGVAL(j) is stored in the j-th column of EIGVEC.

The shape of EIGVEC must verify:

size( EIGVEC, 1 ) = size( D ) = n .

size( EIGVEC, 2 ) = size( EIGVAL ) ,

FAILURE (OUTPUT) logical(lgl)

On exit:

FAILURE = false indicates successful exit.

FAILURE = true indicates that some eigenvectors failed to converge in MAXITER iterations.

MATP (INPUT) real(stnd), dimension(:)

On entry, the vectors and the scalars which define the elementary reflectors used to reduce the packed real n-by-n symmetric matrix MAT to symmetric tridiagonal form T, as returned by SYMTRID_CMP with STORE_Q=true, in its argument MATP. MATP is not modified by the routine.

Back-transformation is used to find the selected eigenvectors of the original matrix MAT and these eigenvectors are stored in argument EIGVEC.

The size of MATP must verify: size( MATP ) = (n * (n+1)/2)

MAXITER (INPUT, OPTIONAL) integer(i4b)

The number of inverse iterations performed in the subroutine. By default, 2 inverse iterations are performed for all the eigenvectors of the tridiagonal matrix T.

ORTHO (INPUT, OPTIONAL) logical(lgl)

On entry, if:

ORTHO=true, all the eigenvectors are orthogonalized by the Modified Gram-Schmidt or QR algorithm;

ORTHO=false, the eigenvectors are not orthogonalized by the Modified Gram-Schmidt or QR algorithm.

The default is to orthogonalize the eigenvectors only if the eigenvalues are not well-separated.

BACKWARD_SWEEP (INPUT, OPTIONAL) logical(lgl)

On entry, if:

BACKWARD_SWEEP=true and the eigenvectors are orthogonalized by the modified Gram-Schmidt algorithm, a backward sweep of the modified Gram-Schmidt algorithm is also performed;

BACKWARD_SWEEP=false a backward sweep is not performed.

The default is not to perform a backward sweep of the modified Gram-Schmidt algorithm.

SCALING (INPUT, OPTIONAL) logical(lgl)

On entry, if:

SCALING=true, the tridiagonal matrix T is scaled before computing the eigenvectors;

SCALING=false, the tridiagonal matrix T is not scaled.

The default is to scale the tridiagonal matrix.

INITVEC (INPUT, OPTIONAL) logical(lgl)

On entry, if:

INITVEC=true, Fernando vectors are used to start the inverse iteration process;

INITVEC=false, random uniform starting vectors are used.

For unreduced tridiagonal matrices, the default is to use Fernando starting vectors if the eigenvalues are well-separated and random uniform starting vectors otherwise. For reduced tridiagonal matrices, the default is to use random uniform starting vectors.

Further Details¶

TRID_INVITER uses Fernando’s method for computing a first estimate of (selected) eigenvectors corresponding to (selected) approximate eigenvalues of a real n-by-n symmetric tridiagonal matrix T (by default, only for the eigenvalues which are well separated and if the input tridiagonal matrix T is unreduced). See the references (1), (2) and (4) for details.

These approximate eigenvectors are then refined (or computed if Fernando’s method is not used) using an inverse iteration algorithm for all the eigenvalues at one step. See the reference (3) for details.

The eigenvectors are then orthogonalized by the Modified Gram-Schmidt or QR algorithm if clusters of eigenvalues are present in a second step.

In a final step, the corresponding (selected) eigenvectors of the full real n-by-n symmetric matrix MAT are computed by a blocked back-transformation algorithm with the Householder transformations used to reduce the full real n-by-n symmetric matrix MAT to symmetric tridiagonal form T (see the references (5) and (6)). These Householder transformations must be packed in the linear array MATP (as returned by SYMTRID_CMP) on entry of TRID_INVITER.

Furthermore, the computation of the eigenvectors is parallelized if OPENMP is used.

TRID_INVITER may fail if some the eigenvalues specified in parameter EIGVAL are nearly identical.

For further details on Fernando’s method or inverse iteration for computing eigenvectors of tridiagonal matrices or the blocked back-transformation algorithm, see:

Fernando, K.V., 1997:

On computing an eigenvector of a tridiagonal matrix. Part I: Basic results. Siam J. Matrix Anal. Appl., Vol. 18, pp. 1013-1034.

Parlett, B.N., and Dhillon, I.S., 1997:

Fernando’s solution to Wilkinson’s problem: An application of double factorization. Linear Algebra and its Appl., 267, pp.247-279.

Golub, G.H., and Van Loan, C.F., 1996:

Matrix Computations. 3rd ed. The Johns Hopkins University Press, Baltimore.

Bini, D.A., Gemignani, L., and Tisseur, F., 2005:

The Ehrlich-Aberth method for the nonsymmetric tridiagonal eigenvalue problem. SIAM J. Matrix Anal. Appl., 27, 153-175.

Dongarra, J.J., Sorensen, D.C., and Hammarling, S.J., 1989:

Block reduction of matrices to condensed form for eigenvalue computations. J. of Computational and Applied Mathematics, Vol. 27, pp. 215-227.

Walker, H.F., 1988:

Implementation of the GMRES method using Householder transformations. Siam J. Sci. Stat. Comput., Vol. 9, No 1, pp. 152-163.

`subroutine gen_symtrid_mat ( type, d, e, failure, known_eigval, eigval, sort, val1, val2, l0, glu0 )`¶

Purpose¶

GEN_SYMTRID_MAT generates different types of symmetric tridiagonal matrices with known eigenvalues or specific numerical properties such as clustered eigenvalues 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.

Arguments¶

TYPE (INPUT) integer(i4b)

Select the type of symmetric tridiagonal matrix TRID to be generated by the subroutine.

If TYPE is between 1 and 49, the subroutine generates a specific symmetric tridiagonal matrix as described in the comments inside the code of the subroutine. For other values of TYPE, all diagonal and off-diagonal elements of the symmetric tridiagonal matrix are generated from an uniform random numbers distribution between 0 and 1.

For TYPE between 1 and 17, the eigenvalues of the tridiagonal symetric matrix are known analytically. For other values of TYPE, the eigenvalues are estimated by a bisection algorithm with high accuracy.

In all cases, the eigenvalues may be output in the optional parameter EIGVAL.

D (OUTPUT) real(stnd), dimension(:)

On exit, D contains the diagonal elements of the tridiagonal matrix TRID.

The size of D must verify: size( D )>=2 .

E (OUTPUT) real(stnd), dimension(:)

On exit, E contains the off-diagonal elements of the tridiagonal matrix TRID. E(size(E)) is arbitrary.

The size of E must verify: size( E ) = size( D ) .

FAILURE (OUTPUT, OPTIONAL) logical(lgl)

On exit:

FAILURE = false : indicates that the eigenvalues of TRID are known analytically or have been computed with high accuracy;

FAILURE = true : indicates that the eigenvalues of TRID are not known analytically and have not been computed with maximum accuracy with the bisection algorithm.

KNOWN_EIGVAL (OUTPUT, OPTIONAL) logical(lgl)

On exit:

KNOWN_EIGVAL = true : indicates that the eigenvalues of TRID are known analytically for the selected TYPE.

KNOWN_EIGVAL = false : indicates that the eigenvalues of TRID are not known analytically for the selected TYPE.

EIGVAL (OUTPUT, OPTIONAL) real(stnd), dimension(:)

On exit, the eigenvalues of TRID computed analytically or estimated to high accuracy with a bisection algorithm.

The size of EIGVAL must verify: size( EIGVAL ) = size( D ) .

SORT (INPUT, OPTIONAL) character

Sort the eigenvalues into ascending order if SORT = ‘A’ or ‘a’, or in descending order if SORT = ‘D’ or ‘d’, if the optional argument EIGVAL is present. For other values of SORT nothing is done and EIGVAL(:) may not be sorted.

VAL1 (INPUT, OPTIONAL) real(stnd)

On entry, specifies the parameter d0 for parametrized symmetric tridiagonal matrices (e.g. TYPE= 3-4, 6-10, 12, 35-38).

If this parameter is changed for TYPE between 35 and 38, which correspond to graded (or reversely graded) matrices with an arithmetic or geometric progression, care must be taken to insure that some elements of the arithmetic or geometric progression will not underflow or overflow as no checks are done in the subroutine for such errors.

The default is 1. .

VAL2 (INPUT, OPTIONAL) real(stnd)

On entry, specifies the parameter e0 for parametrized symmetric tridiagonal matrices (e.g. TYPE= 3-4, 6-10, 12, 35-38).

If this parameter is changed for TYPE between 35 and 38, which correspond to graded (or reversely graded) matrices with an arithmetic or geometric progression, care must be taken to insure that some elements of the arithmetic or geometric progression will not underflow or overflow as no checks are done in the subroutine for such errors.

The default is 2. .

L0 (INPUT, OPTIONAL) integer(i4b)

On entry, specify the radius of the initial matrix for parametrized form of glued tridiagonal matrices (e.g. TYPE between 45 and 49).

L0 must be greater than 0 and preferably less or equal to size( D )/2 . The default is 5. .

GLU0 (INPUT, OPTIONAL) real(stnd)

On entry, specify the glue parameter for parametrized form of glued tridiagonal matrices (e.g. TYPE between 45 and 49).

The default is sqrt( epsilon(GLU0) ).

Further Details¶

This subroutine tries to take care of imprecisions in intrinsic subroutines (e.g. like the cos function in the gfortran compiler) when computing eigenvalues by analytic formulae.

For further details on the tridiagonal matrices used for testing in GEN_SYMTRID_MAT subroutine, see:

Gladwell, G.M.L., Jones, T.H., Willms N.B., 2014:

A test matrix for an inverse eigenvalue problem. Journal of Applied Mathematics, 14, 6 pages, Article ID 515082, DOI 10.1155/2014/515082.

Clement, P.A., 1959:

A class of triple-diagonal matrices for test purposes. SIAM Review, 1(1):50-52, DOI 10.1137/1001006.

Gregory, R.T., Karney, D.L., 1969:

A collection of matrices for testing computational algorithms. New York: Wiley. Reprinted with corrections by Robert E. Krieger, Huntington, New York, 1978.

Higham, N.J., 1991: Algorithm 694:

A collection of test matrices in MATLAB. ACM Transactions on Mathematical Software 17(3):289-305 DOI 10.1145/114697.116805.

Godunov, S.K., Antonov, A.G., Kirillyuk, O.P., and Kostin, V.I., 1993:

Guaranteed Accuracy in numerical linear algebra. Kluwer Academic Publishers.

Parlett, B.N., and Vomel, C., 2005:

How the MRRR algorithm can fail on tight eigenvalue clusters. Lapack Working Note 163.

Nakatsukasa, Y., Aishima, K., and Yamazaki, I., 2012:

dqds with agressive early deflation. SIAM J. Matrix Anal. Appl., 33(1): 22-51.

Fernando, K.V., and Parlett, B.N., 1994:

Accurate singular values and differenial qd algorithms. Numer. Math., 67: 191-229.

Module_Eig_Procedures¶

subroutine symtrid_cmp ( mat, d, e, store_q, upper )¶

Purpose¶

Arguments¶

Further Details¶

subroutine symtrid_cmp ( mat, d, e, store_q )¶

Purpose¶

Arguments¶

Further Details¶

subroutine symtrid_cmp ( matp, d, e, store_q, upper )¶

Purpose¶

Arguments¶

Further Details¶

subroutine symtrid_cmp ( matp, d, e, store_q )¶

Purpose¶

Arguments¶

Further Details¶

subroutine symtrid_cmp2 ( mat, d, e, store_q )¶

Purpose¶

Arguments¶

Further Details¶

subroutine ortho_gen_symtrid ( mat, upper )¶

Purpose¶

Arguments¶

Further Details¶

subroutine ortho_gen_symtrid ( mat )¶

Purpose¶

Arguments¶

Further Details¶

subroutine apply_q_symtrid ( mat, c, left, trans, upper )¶

Purpose¶

Arguments¶

Further Details¶

subroutine apply_q_symtrid ( mat, c, left, trans )¶

Purpose¶

Arguments¶

Further Details¶

subroutine apply_q_symtrid ( matp, c, left, trans, upper )¶

Purpose¶

Arguments¶

Further Details¶

subroutine apply_q_symtrid ( matp, c, left, trans )¶

Purpose¶

Arguments¶

Further Details¶

subroutine eig_sort ( sort, d, u )¶

Purpose¶

Arguments¶

Further Details¶

subroutine eig_abs_sort ( sort, d, u )¶

Purpose¶

Arguments¶

Further Details¶

subroutine eigval_sort ( sort, d )¶

Purpose¶

Arguments¶

Further Details¶

subroutine eigval_abs_sort ( sort, d )¶

Purpose¶

Arguments¶

Further Details¶

subroutine symtrid_qri ( d, e, failure, mat, init_mat, sort, maxiter )¶

Purpose¶

Arguments¶

Further Details¶

subroutine symtrid_qri ( d, e, failure, sort, maxiter )¶

Purpose¶

Arguments¶

Further Details¶

subroutine symtrid_qri2 ( d, e, failure, mat, init_mat, sort, maxiter, max_francis_steps )¶

Purpose¶

Arguments¶

Further Details¶

subroutine symtrid_qri2 ( d, e, failure, sort, maxiter )¶

Purpose¶

Arguments¶

Further Details¶

subroutine symtrid_qri3 ( d, e, failure, mat, init_mat, sort, maxiter, max_francis_steps )¶

Purpose¶

Arguments¶

`subroutine symtrid_cmp ( mat, d, e, store_q, upper )`¶

`subroutine symtrid_cmp ( mat, d, e, store_q )`¶

`subroutine symtrid_cmp ( matp, d, e, store_q, upper )`¶

`subroutine symtrid_cmp ( matp, d, e, store_q )`¶

`subroutine symtrid_cmp2 ( mat, d, e, store_q )`¶

`subroutine ortho_gen_symtrid ( mat, upper )`¶

`subroutine ortho_gen_symtrid ( mat )`¶

`subroutine apply_q_symtrid ( mat, c, left, trans, upper )`¶

`subroutine apply_q_symtrid ( mat, c, left, trans )`¶

`subroutine apply_q_symtrid ( matp, c, left, trans, upper )`¶

`subroutine apply_q_symtrid ( matp, c, left, trans )`¶

`subroutine eig_sort ( sort, d, u )`¶

`subroutine eig_abs_sort ( sort, d, u )`¶

`subroutine eigval_sort ( sort, d )`¶

`subroutine eigval_abs_sort ( sort, d )`¶

`subroutine symtrid_qri ( d, e, failure, mat, init_mat, sort, maxiter )`¶

`subroutine symtrid_qri ( d, e, failure, sort, maxiter )`¶

`subroutine symtrid_qri2 ( d, e, failure, mat, init_mat, sort, maxiter, max_francis_steps )`¶

`subroutine symtrid_qri2 ( d, e, failure, sort, maxiter )`¶

`subroutine symtrid_qri3 ( d, e, failure, mat, init_mat, sort, maxiter, max_francis_steps )`¶

`subroutine symtrid_qri3 ( d, e, failure, sort, maxiter )`¶

`subroutine lae2 ( a, b, c, rt1, rt2 )`¶

`subroutine laev2 ( a, b, c, rt1, rt2, cs1, sn1 )`¶

`subroutine symtrid_ratqri ( d, e, m, failure, small, tol )`¶

`subroutine symtrid_ratqri2 ( d, e, val, failure, m, small, tol )`¶

`subroutine symtrid_bisect ( d, e, neig, w, failure, small, sort, vector, abstol, le, theta, scaling, init )`¶

`subroutine dflgen ( d, e, lambda, cs, sn )`¶

`subroutine dflgen2 ( d, e, lambda, cs, sn, deflate )`¶

`subroutine dflapp ( d, e, cs, sn, deflate )`¶

`subroutine qrstep ( d, e, lambda, cs, sn, deflate )`¶

`subroutine prodgiv ( cs, sn, x )`¶

`function prodgiv_eigvec ( cs, sn )`¶

`subroutine symtrid_deflate ( d, e, eigval, eigvec, failure, max_qr_steps )`¶

`subroutine symtrid_deflate ( d, e, eigval, eigvec, failure, max_qr_steps )`¶

`subroutine trid_deflate ( d, e, eigval, eigvec, failure, max_qr_steps, ortho, inviter )`¶

`subroutine trid_deflate ( d, e, eigval, eigvec, failure, mat, max_qr_steps, ortho, inviter )`¶