383 lines
10 KiB
Fortran
383 lines
10 KiB
Fortran
*> \brief \b SGET22
|
|
*
|
|
* =========== DOCUMENTATION ===========
|
|
*
|
|
* Online html documentation available at
|
|
* http://www.netlib.org/lapack/explore-html/
|
|
*
|
|
* Definition:
|
|
* ===========
|
|
*
|
|
* SUBROUTINE SGET22( TRANSA, TRANSE, TRANSW, N, A, LDA, E, LDE, WR,
|
|
* WI, WORK, RESULT )
|
|
*
|
|
* .. Scalar Arguments ..
|
|
* CHARACTER TRANSA, TRANSE, TRANSW
|
|
* INTEGER LDA, LDE, N
|
|
* ..
|
|
* .. Array Arguments ..
|
|
* REAL A( LDA, * ), E( LDE, * ), RESULT( 2 ), WI( * ),
|
|
* $ WORK( * ), WR( * )
|
|
* ..
|
|
*
|
|
*
|
|
*> \par Purpose:
|
|
* =============
|
|
*>
|
|
*> \verbatim
|
|
*>
|
|
*> SGET22 does an eigenvector check.
|
|
*>
|
|
*> The basic test is:
|
|
*>
|
|
*> RESULT(1) = | A E - E W | / ( |A| |E| ulp )
|
|
*>
|
|
*> using the 1-norm. It also tests the normalization of E:
|
|
*>
|
|
*> RESULT(2) = max | m-norm(E(j)) - 1 | / ( n ulp )
|
|
*> j
|
|
*>
|
|
*> where E(j) is the j-th eigenvector, and m-norm is the max-norm of a
|
|
*> vector. If an eigenvector is complex, as determined from WI(j)
|
|
*> nonzero, then the max-norm of the vector ( er + i*ei ) is the maximum
|
|
*> of
|
|
*> |er(1)| + |ei(1)|, ... , |er(n)| + |ei(n)|
|
|
*>
|
|
*> W is a block diagonal matrix, with a 1 by 1 block for each real
|
|
*> eigenvalue and a 2 by 2 block for each complex conjugate pair.
|
|
*> If eigenvalues j and j+1 are a complex conjugate pair, so that
|
|
*> WR(j) = WR(j+1) = wr and WI(j) = - WI(j+1) = wi, then the 2 by 2
|
|
*> block corresponding to the pair will be:
|
|
*>
|
|
*> ( wr wi )
|
|
*> ( -wi wr )
|
|
*>
|
|
*> Such a block multiplying an n by 2 matrix ( ur ui ) on the right
|
|
*> will be the same as multiplying ur + i*ui by wr + i*wi.
|
|
*>
|
|
*> To handle various schemes for storage of left eigenvectors, there are
|
|
*> options to use A-transpose instead of A, E-transpose instead of E,
|
|
*> and/or W-transpose instead of W.
|
|
*> \endverbatim
|
|
*
|
|
* Arguments:
|
|
* ==========
|
|
*
|
|
*> \param[in] TRANSA
|
|
*> \verbatim
|
|
*> TRANSA is CHARACTER*1
|
|
*> Specifies whether or not A is transposed.
|
|
*> = 'N': No transpose
|
|
*> = 'T': Transpose
|
|
*> = 'C': Conjugate transpose (= Transpose)
|
|
*> \endverbatim
|
|
*>
|
|
*> \param[in] TRANSE
|
|
*> \verbatim
|
|
*> TRANSE is CHARACTER*1
|
|
*> Specifies whether or not E is transposed.
|
|
*> = 'N': No transpose, eigenvectors are in columns of E
|
|
*> = 'T': Transpose, eigenvectors are in rows of E
|
|
*> = 'C': Conjugate transpose (= Transpose)
|
|
*> \endverbatim
|
|
*>
|
|
*> \param[in] TRANSW
|
|
*> \verbatim
|
|
*> TRANSW is CHARACTER*1
|
|
*> Specifies whether or not W is transposed.
|
|
*> = 'N': No transpose
|
|
*> = 'T': Transpose, use -WI(j) instead of WI(j)
|
|
*> = 'C': Conjugate transpose, use -WI(j) instead of WI(j)
|
|
*> \endverbatim
|
|
*>
|
|
*> \param[in] N
|
|
*> \verbatim
|
|
*> N is INTEGER
|
|
*> The order of the matrix A. N >= 0.
|
|
*> \endverbatim
|
|
*>
|
|
*> \param[in] A
|
|
*> \verbatim
|
|
*> A is REAL array, dimension (LDA,N)
|
|
*> The matrix whose eigenvectors are in E.
|
|
*> \endverbatim
|
|
*>
|
|
*> \param[in] LDA
|
|
*> \verbatim
|
|
*> LDA is INTEGER
|
|
*> The leading dimension of the array A. LDA >= max(1,N).
|
|
*> \endverbatim
|
|
*>
|
|
*> \param[in] E
|
|
*> \verbatim
|
|
*> E is REAL array, dimension (LDE,N)
|
|
*> The matrix of eigenvectors. If TRANSE = 'N', the eigenvectors
|
|
*> are stored in the columns of E, if TRANSE = 'T' or 'C', the
|
|
*> eigenvectors are stored in the rows of E.
|
|
*> \endverbatim
|
|
*>
|
|
*> \param[in] LDE
|
|
*> \verbatim
|
|
*> LDE is INTEGER
|
|
*> The leading dimension of the array E. LDE >= max(1,N).
|
|
*> \endverbatim
|
|
*>
|
|
*> \param[in] WR
|
|
*> \verbatim
|
|
*> WR is REAL array, dimension (N)
|
|
*> \endverbatim
|
|
*>
|
|
*> \param[in] WI
|
|
*> \verbatim
|
|
*> WI is REAL array, dimension (N)
|
|
*>
|
|
*> The real and imaginary parts of the eigenvalues of A.
|
|
*> Purely real eigenvalues are indicated by WI(j) = 0.
|
|
*> Complex conjugate pairs are indicated by WR(j)=WR(j+1) and
|
|
*> WI(j) = - WI(j+1) non-zero; the real part is assumed to be
|
|
*> stored in the j-th row/column and the imaginary part in
|
|
*> the (j+1)-th row/column.
|
|
*> \endverbatim
|
|
*>
|
|
*> \param[out] WORK
|
|
*> \verbatim
|
|
*> WORK is REAL array, dimension (N*(N+1))
|
|
*> \endverbatim
|
|
*>
|
|
*> \param[out] RESULT
|
|
*> \verbatim
|
|
*> RESULT is REAL array, dimension (2)
|
|
*> RESULT(1) = | A E - E W | / ( |A| |E| ulp )
|
|
*> RESULT(2) = max | m-norm(E(j)) - 1 | / ( n ulp )
|
|
*> j
|
|
*> \endverbatim
|
|
*
|
|
* Authors:
|
|
* ========
|
|
*
|
|
*> \author Univ. of Tennessee
|
|
*> \author Univ. of California Berkeley
|
|
*> \author Univ. of Colorado Denver
|
|
*> \author NAG Ltd.
|
|
*
|
|
*> \ingroup single_eig
|
|
*
|
|
* =====================================================================
|
|
SUBROUTINE SGET22( TRANSA, TRANSE, TRANSW, N, A, LDA, E, LDE, WR,
|
|
$ WI, WORK, RESULT )
|
|
*
|
|
* -- LAPACK test routine --
|
|
* -- LAPACK is a software package provided by Univ. of Tennessee, --
|
|
* -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
|
|
*
|
|
* .. Scalar Arguments ..
|
|
CHARACTER TRANSA, TRANSE, TRANSW
|
|
INTEGER LDA, LDE, N
|
|
* ..
|
|
* .. Array Arguments ..
|
|
REAL A( LDA, * ), E( LDE, * ), RESULT( 2 ), WI( * ),
|
|
$ WORK( * ), WR( * )
|
|
* ..
|
|
*
|
|
* =====================================================================
|
|
*
|
|
* .. Parameters ..
|
|
REAL ZERO, ONE
|
|
PARAMETER ( ZERO = 0.0, ONE = 1.0 )
|
|
* ..
|
|
* .. Local Scalars ..
|
|
CHARACTER NORMA, NORME
|
|
INTEGER IECOL, IEROW, INCE, IPAIR, ITRNSE, J, JCOL,
|
|
$ JVEC
|
|
REAL ANORM, ENORM, ENRMAX, ENRMIN, ERRNRM, TEMP1,
|
|
$ ULP, UNFL
|
|
* ..
|
|
* .. Local Arrays ..
|
|
REAL WMAT( 2, 2 )
|
|
* ..
|
|
* .. External Functions ..
|
|
LOGICAL LSAME
|
|
REAL SLAMCH, SLANGE
|
|
EXTERNAL LSAME, SLAMCH, SLANGE
|
|
* ..
|
|
* .. External Subroutines ..
|
|
EXTERNAL SAXPY, SGEMM, SLASET
|
|
* ..
|
|
* .. Intrinsic Functions ..
|
|
INTRINSIC ABS, MAX, MIN, REAL
|
|
* ..
|
|
* .. Executable Statements ..
|
|
*
|
|
* Initialize RESULT (in case N=0)
|
|
*
|
|
RESULT( 1 ) = ZERO
|
|
RESULT( 2 ) = ZERO
|
|
IF( N.LE.0 )
|
|
$ RETURN
|
|
*
|
|
UNFL = SLAMCH( 'Safe minimum' )
|
|
ULP = SLAMCH( 'Precision' )
|
|
*
|
|
ITRNSE = 0
|
|
INCE = 1
|
|
NORMA = 'O'
|
|
NORME = 'O'
|
|
*
|
|
IF( LSAME( TRANSA, 'T' ) .OR. LSAME( TRANSA, 'C' ) ) THEN
|
|
NORMA = 'I'
|
|
END IF
|
|
IF( LSAME( TRANSE, 'T' ) .OR. LSAME( TRANSE, 'C' ) ) THEN
|
|
NORME = 'I'
|
|
ITRNSE = 1
|
|
INCE = LDE
|
|
END IF
|
|
*
|
|
* Check normalization of E
|
|
*
|
|
ENRMIN = ONE / ULP
|
|
ENRMAX = ZERO
|
|
IF( ITRNSE.EQ.0 ) THEN
|
|
*
|
|
* Eigenvectors are column vectors.
|
|
*
|
|
IPAIR = 0
|
|
DO 30 JVEC = 1, N
|
|
TEMP1 = ZERO
|
|
IF( IPAIR.EQ.0 .AND. JVEC.LT.N .AND. WI( JVEC ).NE.ZERO )
|
|
$ IPAIR = 1
|
|
IF( IPAIR.EQ.1 ) THEN
|
|
*
|
|
* Complex eigenvector
|
|
*
|
|
DO 10 J = 1, N
|
|
TEMP1 = MAX( TEMP1, ABS( E( J, JVEC ) )+
|
|
$ ABS( E( J, JVEC+1 ) ) )
|
|
10 CONTINUE
|
|
ENRMIN = MIN( ENRMIN, TEMP1 )
|
|
ENRMAX = MAX( ENRMAX, TEMP1 )
|
|
IPAIR = 2
|
|
ELSE IF( IPAIR.EQ.2 ) THEN
|
|
IPAIR = 0
|
|
ELSE
|
|
*
|
|
* Real eigenvector
|
|
*
|
|
DO 20 J = 1, N
|
|
TEMP1 = MAX( TEMP1, ABS( E( J, JVEC ) ) )
|
|
20 CONTINUE
|
|
ENRMIN = MIN( ENRMIN, TEMP1 )
|
|
ENRMAX = MAX( ENRMAX, TEMP1 )
|
|
IPAIR = 0
|
|
END IF
|
|
30 CONTINUE
|
|
*
|
|
ELSE
|
|
*
|
|
* Eigenvectors are row vectors.
|
|
*
|
|
DO 40 JVEC = 1, N
|
|
WORK( JVEC ) = ZERO
|
|
40 CONTINUE
|
|
*
|
|
DO 60 J = 1, N
|
|
IPAIR = 0
|
|
DO 50 JVEC = 1, N
|
|
IF( IPAIR.EQ.0 .AND. JVEC.LT.N .AND. WI( JVEC ).NE.ZERO )
|
|
$ IPAIR = 1
|
|
IF( IPAIR.EQ.1 ) THEN
|
|
WORK( JVEC ) = MAX( WORK( JVEC ),
|
|
$ ABS( E( J, JVEC ) )+ABS( E( J,
|
|
$ JVEC+1 ) ) )
|
|
WORK( JVEC+1 ) = WORK( JVEC )
|
|
ELSE IF( IPAIR.EQ.2 ) THEN
|
|
IPAIR = 0
|
|
ELSE
|
|
WORK( JVEC ) = MAX( WORK( JVEC ),
|
|
$ ABS( E( J, JVEC ) ) )
|
|
IPAIR = 0
|
|
END IF
|
|
50 CONTINUE
|
|
60 CONTINUE
|
|
*
|
|
DO 70 JVEC = 1, N
|
|
ENRMIN = MIN( ENRMIN, WORK( JVEC ) )
|
|
ENRMAX = MAX( ENRMAX, WORK( JVEC ) )
|
|
70 CONTINUE
|
|
END IF
|
|
*
|
|
* Norm of A:
|
|
*
|
|
ANORM = MAX( SLANGE( NORMA, N, N, A, LDA, WORK ), UNFL )
|
|
*
|
|
* Norm of E:
|
|
*
|
|
ENORM = MAX( SLANGE( NORME, N, N, E, LDE, WORK ), ULP )
|
|
*
|
|
* Norm of error:
|
|
*
|
|
* Error = AE - EW
|
|
*
|
|
CALL SLASET( 'Full', N, N, ZERO, ZERO, WORK, N )
|
|
*
|
|
IPAIR = 0
|
|
IEROW = 1
|
|
IECOL = 1
|
|
*
|
|
DO 80 JCOL = 1, N
|
|
IF( ITRNSE.EQ.1 ) THEN
|
|
IEROW = JCOL
|
|
ELSE
|
|
IECOL = JCOL
|
|
END IF
|
|
*
|
|
IF( IPAIR.EQ.0 .AND. WI( JCOL ).NE.ZERO )
|
|
$ IPAIR = 1
|
|
*
|
|
IF( IPAIR.EQ.1 ) THEN
|
|
WMAT( 1, 1 ) = WR( JCOL )
|
|
WMAT( 2, 1 ) = -WI( JCOL )
|
|
WMAT( 1, 2 ) = WI( JCOL )
|
|
WMAT( 2, 2 ) = WR( JCOL )
|
|
CALL SGEMM( TRANSE, TRANSW, N, 2, 2, ONE, E( IEROW, IECOL ),
|
|
$ LDE, WMAT, 2, ZERO, WORK( N*( JCOL-1 )+1 ), N )
|
|
IPAIR = 2
|
|
ELSE IF( IPAIR.EQ.2 ) THEN
|
|
IPAIR = 0
|
|
*
|
|
ELSE
|
|
*
|
|
CALL SAXPY( N, WR( JCOL ), E( IEROW, IECOL ), INCE,
|
|
$ WORK( N*( JCOL-1 )+1 ), 1 )
|
|
IPAIR = 0
|
|
END IF
|
|
*
|
|
80 CONTINUE
|
|
*
|
|
CALL SGEMM( TRANSA, TRANSE, N, N, N, ONE, A, LDA, E, LDE, -ONE,
|
|
$ WORK, N )
|
|
*
|
|
ERRNRM = SLANGE( 'One', N, N, WORK, N, WORK( N*N+1 ) ) / ENORM
|
|
*
|
|
* Compute RESULT(1) (avoiding under/overflow)
|
|
*
|
|
IF( ANORM.GT.ERRNRM ) THEN
|
|
RESULT( 1 ) = ( ERRNRM / ANORM ) / ULP
|
|
ELSE
|
|
IF( ANORM.LT.ONE ) THEN
|
|
RESULT( 1 ) = ONE / ULP
|
|
ELSE
|
|
RESULT( 1 ) = MIN( ERRNRM / ANORM, ONE ) / ULP
|
|
END IF
|
|
END IF
|
|
*
|
|
* Compute RESULT(2) : the normalization error in E.
|
|
*
|
|
RESULT( 2 ) = MAX( ABS( ENRMAX-ONE ), ABS( ENRMIN-ONE ) ) /
|
|
$ ( REAL( N )*ULP )
|
|
*
|
|
RETURN
|
|
*
|
|
* End of SGET22
|
|
*
|
|
END
|