LAPACK  3.4.1
LAPACK: Linear Algebra PACKage
zget22.f
Go to the documentation of this file.
00001 *> \brief \b ZGET22
00002 *
00003 *  =========== DOCUMENTATION ===========
00004 *
00005 * Online html documentation available at 
00006 *            http://www.netlib.org/lapack/explore-html/ 
00007 *
00008 *  Definition:
00009 *  ===========
00010 *
00011 *       SUBROUTINE ZGET22( TRANSA, TRANSE, TRANSW, N, A, LDA, E, LDE, W,
00012 *                          WORK, RWORK, RESULT )
00013 * 
00014 *       .. Scalar Arguments ..
00015 *       CHARACTER          TRANSA, TRANSE, TRANSW
00016 *       INTEGER            LDA, LDE, N
00017 *       ..
00018 *       .. Array Arguments ..
00019 *       DOUBLE PRECISION   RESULT( 2 ), RWORK( * )
00020 *       COMPLEX*16         A( LDA, * ), E( LDE, * ), W( * ), WORK( * )
00021 *       ..
00022 *  
00023 *
00024 *> \par Purpose:
00025 *  =============
00026 *>
00027 *> \verbatim
00028 *>
00029 *> ZGET22 does an eigenvector check.
00030 *>
00031 *> The basic test is:
00032 *>
00033 *>    RESULT(1) = | A E  -  E W | / ( |A| |E| ulp )
00034 *>
00035 *> using the 1-norm.  It also tests the normalization of E:
00036 *>
00037 *>    RESULT(2) = max | m-norm(E(j)) - 1 | / ( n ulp )
00038 *>                 j
00039 *>
00040 *> where E(j) is the j-th eigenvector, and m-norm is the max-norm of a
00041 *> vector.  The max-norm of a complex n-vector x in this case is the
00042 *> maximum of |re(x(i)| + |im(x(i)| over i = 1, ..., n.
00043 *> \endverbatim
00044 *
00045 *  Arguments:
00046 *  ==========
00047 *
00048 *> \param[in] TRANSA
00049 *> \verbatim
00050 *>          TRANSA is CHARACTER*1
00051 *>          Specifies whether or not A is transposed.
00052 *>          = 'N':  No transpose
00053 *>          = 'T':  Transpose
00054 *>          = 'C':  Conjugate transpose
00055 *> \endverbatim
00056 *>
00057 *> \param[in] TRANSE
00058 *> \verbatim
00059 *>          TRANSE is CHARACTER*1
00060 *>          Specifies whether or not E is transposed.
00061 *>          = 'N':  No transpose, eigenvectors are in columns of E
00062 *>          = 'T':  Transpose, eigenvectors are in rows of E
00063 *>          = 'C':  Conjugate transpose, eigenvectors are in rows of E
00064 *> \endverbatim
00065 *>
00066 *> \param[in] TRANSW
00067 *> \verbatim
00068 *>          TRANSW is CHARACTER*1
00069 *>          Specifies whether or not W is transposed.
00070 *>          = 'N':  No transpose
00071 *>          = 'T':  Transpose, same as TRANSW = 'N'
00072 *>          = 'C':  Conjugate transpose, use -WI(j) instead of WI(j)
00073 *> \endverbatim
00074 *>
00075 *> \param[in] N
00076 *> \verbatim
00077 *>          N is INTEGER
00078 *>          The order of the matrix A.  N >= 0.
00079 *> \endverbatim
00080 *>
00081 *> \param[in] A
00082 *> \verbatim
00083 *>          A is COMPLEX*16 array, dimension (LDA,N)
00084 *>          The matrix whose eigenvectors are in E.
00085 *> \endverbatim
00086 *>
00087 *> \param[in] LDA
00088 *> \verbatim
00089 *>          LDA is INTEGER
00090 *>          The leading dimension of the array A.  LDA >= max(1,N).
00091 *> \endverbatim
00092 *>
00093 *> \param[in] E
00094 *> \verbatim
00095 *>          E is COMPLEX*16 array, dimension (LDE,N)
00096 *>          The matrix of eigenvectors. If TRANSE = 'N', the eigenvectors
00097 *>          are stored in the columns of E, if TRANSE = 'T' or 'C', the
00098 *>          eigenvectors are stored in the rows of E.
00099 *> \endverbatim
00100 *>
00101 *> \param[in] LDE
00102 *> \verbatim
00103 *>          LDE is INTEGER
00104 *>          The leading dimension of the array E.  LDE >= max(1,N).
00105 *> \endverbatim
00106 *>
00107 *> \param[in] W
00108 *> \verbatim
00109 *>          W is COMPLEX*16 array, dimension (N)
00110 *>          The eigenvalues of A.
00111 *> \endverbatim
00112 *>
00113 *> \param[out] WORK
00114 *> \verbatim
00115 *>          WORK is COMPLEX*16 array, dimension (N*N)
00116 *> \endverbatim
00117 *>
00118 *> \param[out] RWORK
00119 *> \verbatim
00120 *>          RWORK is DOUBLE PRECISION array, dimension (N)
00121 *> \endverbatim
00122 *>
00123 *> \param[out] RESULT
00124 *> \verbatim
00125 *>          RESULT is DOUBLE PRECISION array, dimension (2)
00126 *>          RESULT(1) = | A E  -  E W | / ( |A| |E| ulp )
00127 *>          RESULT(2) = max | m-norm(E(j)) - 1 | / ( n ulp )
00128 *> \endverbatim
00129 *
00130 *  Authors:
00131 *  ========
00132 *
00133 *> \author Univ. of Tennessee 
00134 *> \author Univ. of California Berkeley 
00135 *> \author Univ. of Colorado Denver 
00136 *> \author NAG Ltd. 
00137 *
00138 *> \date November 2011
00139 *
00140 *> \ingroup complex16_eig
00141 *
00142 *  =====================================================================
00143       SUBROUTINE ZGET22( TRANSA, TRANSE, TRANSW, N, A, LDA, E, LDE, W,
00144      $                   WORK, RWORK, RESULT )
00145 *
00146 *  -- LAPACK test routine (version 3.4.0) --
00147 *  -- LAPACK is a software package provided by Univ. of Tennessee,    --
00148 *  -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
00149 *     November 2011
00150 *
00151 *     .. Scalar Arguments ..
00152       CHARACTER          TRANSA, TRANSE, TRANSW
00153       INTEGER            LDA, LDE, N
00154 *     ..
00155 *     .. Array Arguments ..
00156       DOUBLE PRECISION   RESULT( 2 ), RWORK( * )
00157       COMPLEX*16         A( LDA, * ), E( LDE, * ), W( * ), WORK( * )
00158 *     ..
00159 *
00160 *  =====================================================================
00161 *
00162 *     .. Parameters ..
00163       DOUBLE PRECISION   ZERO, ONE
00164       PARAMETER          ( ZERO = 0.0D+0, ONE = 1.0D+0 )
00165       COMPLEX*16         CZERO, CONE
00166       PARAMETER          ( CZERO = ( 0.0D+0, 0.0D+0 ),
00167      $                   CONE = ( 1.0D+0, 0.0D+0 ) )
00168 *     ..
00169 *     .. Local Scalars ..
00170       CHARACTER          NORMA, NORME
00171       INTEGER            ITRNSE, ITRNSW, J, JCOL, JOFF, JROW, JVEC
00172       DOUBLE PRECISION   ANORM, ENORM, ENRMAX, ENRMIN, ERRNRM, TEMP1,
00173      $                   ULP, UNFL
00174       COMPLEX*16         WTEMP
00175 *     ..
00176 *     .. External Functions ..
00177       LOGICAL            LSAME
00178       DOUBLE PRECISION   DLAMCH, ZLANGE
00179       EXTERNAL           LSAME, DLAMCH, ZLANGE
00180 *     ..
00181 *     .. External Subroutines ..
00182       EXTERNAL           ZGEMM, ZLASET
00183 *     ..
00184 *     .. Intrinsic Functions ..
00185       INTRINSIC          ABS, DBLE, DCONJG, DIMAG, MAX, MIN
00186 *     ..
00187 *     .. Executable Statements ..
00188 *
00189 *     Initialize RESULT (in case N=0)
00190 *
00191       RESULT( 1 ) = ZERO
00192       RESULT( 2 ) = ZERO
00193       IF( N.LE.0 )
00194      $   RETURN
00195 *
00196       UNFL = DLAMCH( 'Safe minimum' )
00197       ULP = DLAMCH( 'Precision' )
00198 *
00199       ITRNSE = 0
00200       ITRNSW = 0
00201       NORMA = 'O'
00202       NORME = 'O'
00203 *
00204       IF( LSAME( TRANSA, 'T' ) .OR. LSAME( TRANSA, 'C' ) ) THEN
00205          NORMA = 'I'
00206       END IF
00207 *
00208       IF( LSAME( TRANSE, 'T' ) ) THEN
00209          ITRNSE = 1
00210          NORME = 'I'
00211       ELSE IF( LSAME( TRANSE, 'C' ) ) THEN
00212          ITRNSE = 2
00213          NORME = 'I'
00214       END IF
00215 *
00216       IF( LSAME( TRANSW, 'C' ) ) THEN
00217          ITRNSW = 1
00218       END IF
00219 *
00220 *     Normalization of E:
00221 *
00222       ENRMIN = ONE / ULP
00223       ENRMAX = ZERO
00224       IF( ITRNSE.EQ.0 ) THEN
00225          DO 20 JVEC = 1, N
00226             TEMP1 = ZERO
00227             DO 10 J = 1, N
00228                TEMP1 = MAX( TEMP1, ABS( DBLE( E( J, JVEC ) ) )+
00229      $                 ABS( DIMAG( E( J, JVEC ) ) ) )
00230    10       CONTINUE
00231             ENRMIN = MIN( ENRMIN, TEMP1 )
00232             ENRMAX = MAX( ENRMAX, TEMP1 )
00233    20    CONTINUE
00234       ELSE
00235          DO 30 JVEC = 1, N
00236             RWORK( JVEC ) = ZERO
00237    30    CONTINUE
00238 *
00239          DO 50 J = 1, N
00240             DO 40 JVEC = 1, N
00241                RWORK( JVEC ) = MAX( RWORK( JVEC ),
00242      $                         ABS( DBLE( E( JVEC, J ) ) )+
00243      $                         ABS( DIMAG( E( JVEC, J ) ) ) )
00244    40       CONTINUE
00245    50    CONTINUE
00246 *
00247          DO 60 JVEC = 1, N
00248             ENRMIN = MIN( ENRMIN, RWORK( JVEC ) )
00249             ENRMAX = MAX( ENRMAX, RWORK( JVEC ) )
00250    60    CONTINUE
00251       END IF
00252 *
00253 *     Norm of A:
00254 *
00255       ANORM = MAX( ZLANGE( NORMA, N, N, A, LDA, RWORK ), UNFL )
00256 *
00257 *     Norm of E:
00258 *
00259       ENORM = MAX( ZLANGE( NORME, N, N, E, LDE, RWORK ), ULP )
00260 *
00261 *     Norm of error:
00262 *
00263 *     Error =  AE - EW
00264 *
00265       CALL ZLASET( 'Full', N, N, CZERO, CZERO, WORK, N )
00266 *
00267       JOFF = 0
00268       DO 100 JCOL = 1, N
00269          IF( ITRNSW.EQ.0 ) THEN
00270             WTEMP = W( JCOL )
00271          ELSE
00272             WTEMP = DCONJG( W( JCOL ) )
00273          END IF
00274 *
00275          IF( ITRNSE.EQ.0 ) THEN
00276             DO 70 JROW = 1, N
00277                WORK( JOFF+JROW ) = E( JROW, JCOL )*WTEMP
00278    70       CONTINUE
00279          ELSE IF( ITRNSE.EQ.1 ) THEN
00280             DO 80 JROW = 1, N
00281                WORK( JOFF+JROW ) = E( JCOL, JROW )*WTEMP
00282    80       CONTINUE
00283          ELSE
00284             DO 90 JROW = 1, N
00285                WORK( JOFF+JROW ) = DCONJG( E( JCOL, JROW ) )*WTEMP
00286    90       CONTINUE
00287          END IF
00288          JOFF = JOFF + N
00289   100 CONTINUE
00290 *
00291       CALL ZGEMM( TRANSA, TRANSE, N, N, N, CONE, A, LDA, E, LDE, -CONE,
00292      $            WORK, N )
00293 *
00294       ERRNRM = ZLANGE( 'One', N, N, WORK, N, RWORK ) / ENORM
00295 *
00296 *     Compute RESULT(1) (avoiding under/overflow)
00297 *
00298       IF( ANORM.GT.ERRNRM ) THEN
00299          RESULT( 1 ) = ( ERRNRM / ANORM ) / ULP
00300       ELSE
00301          IF( ANORM.LT.ONE ) THEN
00302             RESULT( 1 ) = ( MIN( ERRNRM, ANORM ) / ANORM ) / ULP
00303          ELSE
00304             RESULT( 1 ) = MIN( ERRNRM / ANORM, ONE ) / ULP
00305          END IF
00306       END IF
00307 *
00308 *     Compute RESULT(2) : the normalization error in E.
00309 *
00310       RESULT( 2 ) = MAX( ABS( ENRMAX-ONE ), ABS( ENRMIN-ONE ) ) /
00311      $              ( DBLE( N )*ULP )
00312 *
00313       RETURN
00314 *
00315 *     End of ZGET22
00316 *
00317       END
 All Files Functions