LAPACK: TESTING/EIG/sget22.f Source File

Go to the documentation of this file.
00001       SUBROUTINE SGET22( TRANSA, TRANSE, TRANSW, N, A, LDA, E, LDE, WR,
00002      $                   WI, WORK, RESULT )
00003 *
00004 *  -- LAPACK test routine (version 3.1) --
00005 *     Univ. of Tennessee, Univ. of California Berkeley and NAG Ltd..
00006 *     November 2006
00007 *
00008 *     .. Scalar Arguments ..
00009       CHARACTER          TRANSA, TRANSE, TRANSW
00010       INTEGER            LDA, LDE, N
00011 *     ..
00012 *     .. Array Arguments ..
00013       REAL               A( LDA, * ), E( LDE, * ), RESULT( 2 ), WI( * ),
00014      $                   WORK( * ), WR( * )
00015 *     ..
00016 *
00017 *  Purpose
00018 *  =======
00019 *
00020 *  SGET22 does an eigenvector check.
00021 *
00022 *  The basic test is:
00023 *
00024 *     RESULT(1) = | A E  -  E W | / ( |A| |E| ulp )
00025 *
00026 *  using the 1-norm.  It also tests the normalization of E:
00027 *
00028 *     RESULT(2) = max | m-norm(E(j)) - 1 | / ( n ulp )
00029 *                  j
00030 *
00031 *  where E(j) is the j-th eigenvector, and m-norm is the max-norm of a
00032 *  vector.  If an eigenvector is complex, as determined from WI(j)
00033 *  nonzero, then the max-norm of the vector ( er + i*ei ) is the maximum
00034 *  of
00035 *     |er(1)| + |ei(1)|, ... , |er(n)| + |ei(n)|
00036 *
00037 *  W is a block diagonal matrix, with a 1 by 1 block for each real
00038 *  eigenvalue and a 2 by 2 block for each complex conjugate pair.
00039 *  If eigenvalues j and j+1 are a complex conjugate pair, so that
00040 *  WR(j) = WR(j+1) = wr and WI(j) = - WI(j+1) = wi, then the 2 by 2
00041 *  block corresponding to the pair will be:
00042 *
00043 *     (  wr  wi  )
00044 *     ( -wi  wr  )
00045 *
00046 *  Such a block multiplying an n by 2 matrix ( ur ui ) on the right
00047 *  will be the same as multiplying  ur + i*ui  by  wr + i*wi.
00048 *
00049 *  To handle various schemes for storage of left eigenvectors, there are
00050 *  options to use A-transpose instead of A, E-transpose instead of E,
00051 *  and/or W-transpose instead of W.
00052 *
00053 *  Arguments
00054 *  ==========
00055 *
00056 *  TRANSA  (input) CHARACTER*1
00057 *          Specifies whether or not A is transposed.
00058 *          = 'N':  No transpose
00059 *          = 'T':  Transpose
00060 *          = 'C':  Conjugate transpose (= Transpose)
00061 *
00062 *  TRANSE  (input) CHARACTER*1
00063 *          Specifies whether or not E is transposed.
00064 *          = 'N':  No transpose, eigenvectors are in columns of E
00065 *          = 'T':  Transpose, eigenvectors are in rows of E
00066 *          = 'C':  Conjugate transpose (= Transpose)
00067 *
00068 *  TRANSW  (input) CHARACTER*1
00069 *          Specifies whether or not W is transposed.
00070 *          = 'N':  No transpose
00071 *          = 'T':  Transpose, use -WI(j) instead of WI(j)
00072 *          = 'C':  Conjugate transpose, use -WI(j) instead of WI(j)
00073 *
00074 *  N       (input) INTEGER
00075 *          The order of the matrix A.  N >= 0.
00076 *
00077 *  A       (input) REAL array, dimension (LDA,N)
00078 *          The matrix whose eigenvectors are in E.
00079 *
00080 *  LDA     (input) INTEGER
00081 *          The leading dimension of the array A.  LDA >= max(1,N).
00082 *
00083 *  E       (input) REAL array, dimension (LDE,N)
00084 *          The matrix of eigenvectors. If TRANSE = 'N', the eigenvectors
00085 *          are stored in the columns of E, if TRANSE = 'T' or 'C', the
00086 *          eigenvectors are stored in the rows of E.
00087 *
00088 *  LDE     (input) INTEGER
00089 *          The leading dimension of the array E.  LDE >= max(1,N).
00090 *
00091 *  WR      (input) REAL array, dimension (N)
00092 *  WI      (input) REAL array, dimension (N)
00093 *          The real and imaginary parts of the eigenvalues of A.
00094 *          Purely real eigenvalues are indicated by WI(j) = 0.
00095 *          Complex conjugate pairs are indicated by WR(j)=WR(j+1) and
00096 *          WI(j) = - WI(j+1) non-zero; the real part is assumed to be
00097 *          stored in the j-th row/column and the imaginary part in
00098 *          the (j+1)-th row/column.
00099 *
00100 *  WORK    (workspace) REAL array, dimension (N*(N+1))
00101 *
00102 *  RESULT  (output) REAL array, dimension (2)
00103 *          RESULT(1) = | A E  -  E W | / ( |A| |E| ulp )
00104 *          RESULT(2) = max | m-norm(E(j)) - 1 | / ( n ulp )
00105 *
00106 *  =====================================================================
00107 *
00108 *     .. Parameters ..
00109       REAL               ZERO, ONE
00110       PARAMETER          ( ZERO = 0.0, ONE = 1.0 )
00111 *     ..
00112 *     .. Local Scalars ..
00113       CHARACTER          NORMA, NORME
00114       INTEGER            IECOL, IEROW, INCE, IPAIR, ITRNSE, J, JCOL,
00115      $                   JVEC
00116       REAL               ANORM, ENORM, ENRMAX, ENRMIN, ERRNRM, TEMP1,
00117      $                   ULP, UNFL
00118 *     ..
00119 *     .. Local Arrays ..
00120       REAL               WMAT( 2, 2 )
00121 *     ..
00122 *     .. External Functions ..
00123       LOGICAL            LSAME
00124       REAL               SLAMCH, SLANGE
00125       EXTERNAL           LSAME, SLAMCH, SLANGE
00126 *     ..
00127 *     .. External Subroutines ..
00128       EXTERNAL           SAXPY, SGEMM, SLASET
00129 *     ..
00130 *     .. Intrinsic Functions ..
00131       INTRINSIC          ABS, MAX, MIN, REAL
00132 *     ..
00133 *     .. Executable Statements ..
00134 *
00135 *     Initialize RESULT (in case N=0)
00136 *
00137       RESULT( 1 ) = ZERO
00138       RESULT( 2 ) = ZERO
00139       IF( N.LE.0 )
00140      $   RETURN
00141 *
00142       UNFL = SLAMCH( 'Safe minimum' )
00143       ULP = SLAMCH( 'Precision' )
00144 *
00145       ITRNSE = 0
00146       INCE = 1
00147       NORMA = 'O'
00148       NORME = 'O'
00149 *
00150       IF( LSAME( TRANSA, 'T' ) .OR. LSAME( TRANSA, 'C' ) ) THEN
00151          NORMA = 'I'
00152       END IF
00153       IF( LSAME( TRANSE, 'T' ) .OR. LSAME( TRANSE, 'C' ) ) THEN
00154          NORME = 'I'
00155          ITRNSE = 1
00156          INCE = LDE
00157       END IF
00158 *
00159 *     Check normalization of E
00160 *
00161       ENRMIN = ONE / ULP
00162       ENRMAX = ZERO
00163       IF( ITRNSE.EQ.0 ) THEN
00164 *
00165 *        Eigenvectors are column vectors.
00166 *
00167          IPAIR = 0
00168          DO 30 JVEC = 1, N
00169             TEMP1 = ZERO
00170             IF( IPAIR.EQ.0 .AND. JVEC.LT.N .AND. WI( JVEC ).NE.ZERO )
00171      $         IPAIR = 1
00172             IF( IPAIR.EQ.1 ) THEN
00173 *
00174 *              Complex eigenvector
00175 *
00176                DO 10 J = 1, N
00177                   TEMP1 = MAX( TEMP1, ABS( E( J, JVEC ) )+
00178      $                    ABS( E( J, JVEC+1 ) ) )
00179    10          CONTINUE
00180                ENRMIN = MIN( ENRMIN, TEMP1 )
00181                ENRMAX = MAX( ENRMAX, TEMP1 )
00182                IPAIR = 2
00183             ELSE IF( IPAIR.EQ.2 ) THEN
00184                IPAIR = 0
00185             ELSE
00186 *
00187 *              Real eigenvector
00188 *
00189                DO 20 J = 1, N
00190                   TEMP1 = MAX( TEMP1, ABS( E( J, JVEC ) ) )
00191    20          CONTINUE
00192                ENRMIN = MIN( ENRMIN, TEMP1 )
00193                ENRMAX = MAX( ENRMAX, TEMP1 )
00194                IPAIR = 0
00195             END IF
00196    30    CONTINUE
00197 *
00198       ELSE
00199 *
00200 *        Eigenvectors are row vectors.
00201 *
00202          DO 40 JVEC = 1, N
00203             WORK( JVEC ) = ZERO
00204    40    CONTINUE
00205 *
00206          DO 60 J = 1, N
00207             IPAIR = 0
00208             DO 50 JVEC = 1, N
00209                IF( IPAIR.EQ.0 .AND. JVEC.LT.N .AND. WI( JVEC ).NE.ZERO )
00210      $            IPAIR = 1
00211                IF( IPAIR.EQ.1 ) THEN
00212                   WORK( JVEC ) = MAX( WORK( JVEC ),
00213      $                           ABS( E( J, JVEC ) )+ABS( E( J,
00214      $                           JVEC+1 ) ) )
00215                   WORK( JVEC+1 ) = WORK( JVEC )
00216                ELSE IF( IPAIR.EQ.2 ) THEN
00217                   IPAIR = 0
00218                ELSE
00219                   WORK( JVEC ) = MAX( WORK( JVEC ),
00220      $                           ABS( E( J, JVEC ) ) )
00221                   IPAIR = 0
00222                END IF
00223    50       CONTINUE
00224    60    CONTINUE
00225 *
00226          DO 70 JVEC = 1, N
00227             ENRMIN = MIN( ENRMIN, WORK( JVEC ) )
00228             ENRMAX = MAX( ENRMAX, WORK( JVEC ) )
00229    70    CONTINUE
00230       END IF
00231 *
00232 *     Norm of A:
00233 *
00234       ANORM = MAX( SLANGE( NORMA, N, N, A, LDA, WORK ), UNFL )
00235 *
00236 *     Norm of E:
00237 *
00238       ENORM = MAX( SLANGE( NORME, N, N, E, LDE, WORK ), ULP )
00239 *
00240 *     Norm of error:
00241 *
00242 *     Error =  AE - EW
00243 *
00244       CALL SLASET( 'Full', N, N, ZERO, ZERO, WORK, N )
00245 *
00246       IPAIR = 0
00247       IEROW = 1
00248       IECOL = 1
00249 *
00250       DO 80 JCOL = 1, N
00251          IF( ITRNSE.EQ.1 ) THEN
00252             IEROW = JCOL
00253          ELSE
00254             IECOL = JCOL
00255          END IF
00256 *
00257          IF( IPAIR.EQ.0 .AND. WI( JCOL ).NE.ZERO )
00258      $      IPAIR = 1
00259 *
00260          IF( IPAIR.EQ.1 ) THEN
00261             WMAT( 1, 1 ) = WR( JCOL )
00262             WMAT( 2, 1 ) = -WI( JCOL )
00263             WMAT( 1, 2 ) = WI( JCOL )
00264             WMAT( 2, 2 ) = WR( JCOL )
00265             CALL SGEMM( TRANSE, TRANSW, N, 2, 2, ONE, E( IEROW, IECOL ),
00266      $                  LDE, WMAT, 2, ZERO, WORK( N*( JCOL-1 )+1 ), N )
00267             IPAIR = 2
00268          ELSE IF( IPAIR.EQ.2 ) THEN
00269             IPAIR = 0
00270 *
00271          ELSE
00272 *
00273             CALL SAXPY( N, WR( JCOL ), E( IEROW, IECOL ), INCE,
00274      $                  WORK( N*( JCOL-1 )+1 ), 1 )
00275             IPAIR = 0
00276          END IF
00277 *
00278    80 CONTINUE
00279 *
00280       CALL SGEMM( TRANSA, TRANSE, N, N, N, ONE, A, LDA, E, LDE, -ONE,
00281      $            WORK, N )
00282 *
00283       ERRNRM = SLANGE( 'One', N, N, WORK, N, WORK( N*N+1 ) ) / ENORM
00284 *
00285 *     Compute RESULT(1) (avoiding under/overflow)
00286 *
00287       IF( ANORM.GT.ERRNRM ) THEN
00288          RESULT( 1 ) = ( ERRNRM / ANORM ) / ULP
00289       ELSE
00290          IF( ANORM.LT.ONE ) THEN
00291             RESULT( 1 ) = ( MIN( ERRNRM, ANORM ) / ANORM ) / ULP
00292          ELSE
00293             RESULT( 1 ) = MIN( ERRNRM / ANORM, ONE ) / ULP
00294          END IF
00295       END IF
00296 *
00297 *     Compute RESULT(2) : the normalization error in E.
00298 *
00299       RESULT( 2 ) = MAX( ABS( ENRMAX-ONE ), ABS( ENRMIN-ONE ) ) /
00300      $              ( REAL( N )*ULP )
00301 *
00302       RETURN
00303 *
00304 *     End of SGET22
00305 *
00306       END