LAPACK: TESTING/EIG/zget52.f Source File

Go to the documentation of this file.
00001       SUBROUTINE ZGET52( LEFT, N, A, LDA, B, LDB, E, LDE, ALPHA, BETA,
00002      $                   WORK, RWORK, 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       LOGICAL            LEFT
00010       INTEGER            LDA, LDB, LDE, N
00011 *     ..
00012 *     .. Array Arguments ..
00013       DOUBLE PRECISION   RESULT( 2 ), RWORK( * )
00014       COMPLEX*16         A( LDA, * ), ALPHA( * ), B( LDB, * ),
00015      $                   BETA( * ), E( LDE, * ), WORK( * )
00016 *     ..
00017 *
00018 *  Purpose
00019 *  =======
00020 *
00021 *  ZGET52  does an eigenvector check for the generalized eigenvalue
00022 *  problem.
00023 *
00024 *  The basic test for right eigenvectors is:
00025 *
00026 *                            | b(i) A E(i) -  a(i) B E(i) |
00027 *          RESULT(1) = max   -------------------------------
00028 *                       i    n ulp max( |b(i) A|, |a(i) B| )
00029 *
00030 *  using the 1-norm.  Here, a(i)/b(i) = w is the i-th generalized
00031 *  eigenvalue of A - w B, or, equivalently, b(i)/a(i) = m is the i-th
00032 *  generalized eigenvalue of m A - B.
00033 *
00034 *                          H   H  _      _
00035 *  For left eigenvectors, A , B , a, and b  are used.
00036 *
00037 *  ZGET52 also tests the normalization of E.  Each eigenvector is
00038 *  supposed to be normalized so that the maximum "absolute value"
00039 *  of its elements is 1, where in this case, "absolute value"
00040 *  of a complex value x is  |Re(x)| + |Im(x)| ; let us call this
00041 *  maximum "absolute value" norm of a vector v  M(v).  
00042 *  If a(i)=b(i)=0, then the eigenvector is set to be the jth coordinate
00043 *  vector. The normalization test is:
00044 *
00045 *          RESULT(2) =      max       | M(v(i)) - 1 | / ( n ulp )
00046 *                     eigenvectors v(i)
00047 *
00048 *
00049 *  Arguments
00050 *  =========
00051 *
00052 *  LEFT    (input) LOGICAL
00053 *          =.TRUE.:  The eigenvectors in the columns of E are assumed
00054 *                    to be *left* eigenvectors.
00055 *          =.FALSE.: The eigenvectors in the columns of E are assumed
00056 *                    to be *right* eigenvectors.
00057 *
00058 *  N       (input) INTEGER
00059 *          The size of the matrices.  If it is zero, ZGET52 does
00060 *          nothing.  It must be at least zero.
00061 *
00062 *  A       (input) COMPLEX*16 array, dimension (LDA, N)
00063 *          The matrix A.
00064 *
00065 *  LDA     (input) INTEGER
00066 *          The leading dimension of A.  It must be at least 1
00067 *          and at least N.
00068 *
00069 *  B       (input) COMPLEX*16 array, dimension (LDB, N)
00070 *          The matrix B.
00071 *
00072 *  LDB     (input) INTEGER
00073 *          The leading dimension of B.  It must be at least 1
00074 *          and at least N.
00075 *
00076 *  E       (input) COMPLEX*16 array, dimension (LDE, N)
00077 *          The matrix of eigenvectors.  It must be O( 1 ).
00078 *
00079 *  LDE     (input) INTEGER
00080 *          The leading dimension of E.  It must be at least 1 and at
00081 *          least N.
00082 *
00083 *  ALPHA   (input) COMPLEX*16 array, dimension (N)
00084 *          The values a(i) as described above, which, along with b(i),
00085 *          define the generalized eigenvalues.
00086 *
00087 *  BETA    (input) COMPLEX*16 array, dimension (N)
00088 *          The values b(i) as described above, which, along with a(i),
00089 *          define the generalized eigenvalues.
00090 *
00091 *  WORK    (workspace) COMPLEX*16 array, dimension (N**2)
00092 *
00093 *  RWORK   (workspace) DOUBLE PRECISION array, dimension (N)
00094 *
00095 *  RESULT  (output) DOUBLE PRECISION array, dimension (2)
00096 *          The values computed by the test described above.  If A E or
00097 *          B E is likely to overflow, then RESULT(1:2) is set to
00098 *          10 / ulp.
00099 *
00100 *  =====================================================================
00101 *
00102 *     .. Parameters ..
00103       DOUBLE PRECISION   ZERO, ONE
00104       PARAMETER          ( ZERO = 0.0D+0, ONE = 1.0D+0 )
00105       COMPLEX*16         CZERO, CONE
00106       PARAMETER          ( CZERO = ( 0.0D+0, 0.0D+0 ),
00107      $                   CONE = ( 1.0D+0, 0.0D+0 ) )
00108 *     ..
00109 *     .. Local Scalars ..
00110       CHARACTER          NORMAB, TRANS
00111       INTEGER            J, JVEC
00112       DOUBLE PRECISION   ABMAX, ALFMAX, ANORM, BETMAX, BNORM, ENORM,
00113      $                   ENRMER, ERRNRM, SAFMAX, SAFMIN, SCALE, TEMP1,
00114      $                   ULP
00115       COMPLEX*16         ACOEFF, ALPHAI, BCOEFF, BETAI, X
00116 *     ..
00117 *     .. External Functions ..
00118       DOUBLE PRECISION   DLAMCH, ZLANGE
00119       EXTERNAL           DLAMCH, ZLANGE
00120 *     ..
00121 *     .. External Subroutines ..
00122       EXTERNAL           ZGEMV
00123 *     ..
00124 *     .. Intrinsic Functions ..
00125       INTRINSIC          ABS, DBLE, DCONJG, DIMAG, MAX
00126 *     ..
00127 *     .. Statement Functions ..
00128       DOUBLE PRECISION   ABS1
00129 *     ..
00130 *     .. Statement Function definitions ..
00131       ABS1( X ) = ABS( DBLE( X ) ) + ABS( DIMAG( X ) )
00132 *     ..
00133 *     .. Executable Statements ..
00134 *
00135       RESULT( 1 ) = ZERO
00136       RESULT( 2 ) = ZERO
00137       IF( N.LE.0 )
00138      $   RETURN
00139 *
00140       SAFMIN = DLAMCH( 'Safe minimum' )
00141       SAFMAX = ONE / SAFMIN
00142       ULP = DLAMCH( 'Epsilon' )*DLAMCH( 'Base' )
00143 *
00144       IF( LEFT ) THEN
00145          TRANS = 'C'
00146          NORMAB = 'I'
00147       ELSE
00148          TRANS = 'N'
00149          NORMAB = 'O'
00150       END IF
00151 *
00152 *     Norm of A, B, and E:
00153 *
00154       ANORM = MAX( ZLANGE( NORMAB, N, N, A, LDA, RWORK ), SAFMIN )
00155       BNORM = MAX( ZLANGE( NORMAB, N, N, B, LDB, RWORK ), SAFMIN )
00156       ENORM = MAX( ZLANGE( 'O', N, N, E, LDE, RWORK ), ULP )
00157       ALFMAX = SAFMAX / MAX( ONE, BNORM )
00158       BETMAX = SAFMAX / MAX( ONE, ANORM )
00159 *
00160 *     Compute error matrix.
00161 *     Column i = ( b(i) A - a(i) B ) E(i) / max( |a(i) B| |b(i) A| )
00162 *
00163       DO 10 JVEC = 1, N
00164          ALPHAI = ALPHA( JVEC )
00165          BETAI = BETA( JVEC )
00166          ABMAX = MAX( ABS1( ALPHAI ), ABS1( BETAI ) )
00167          IF( ABS1( ALPHAI ).GT.ALFMAX .OR. ABS1( BETAI ).GT.BETMAX .OR.
00168      $       ABMAX.LT.ONE ) THEN
00169             SCALE = ONE / MAX( ABMAX, SAFMIN )
00170             ALPHAI = SCALE*ALPHAI
00171             BETAI = SCALE*BETAI
00172          END IF
00173          SCALE = ONE / MAX( ABS1( ALPHAI )*BNORM, ABS1( BETAI )*ANORM,
00174      $           SAFMIN )
00175          ACOEFF = SCALE*BETAI
00176          BCOEFF = SCALE*ALPHAI
00177          IF( LEFT ) THEN
00178             ACOEFF = DCONJG( ACOEFF )
00179             BCOEFF = DCONJG( BCOEFF )
00180          END IF
00181          CALL ZGEMV( TRANS, N, N, ACOEFF, A, LDA, E( 1, JVEC ), 1,
00182      $               CZERO, WORK( N*( JVEC-1 )+1 ), 1 )
00183          CALL ZGEMV( TRANS, N, N, -BCOEFF, B, LDA, E( 1, JVEC ), 1,
00184      $               CONE, WORK( N*( JVEC-1 )+1 ), 1 )
00185    10 CONTINUE
00186 *
00187       ERRNRM = ZLANGE( 'One', N, N, WORK, N, RWORK ) / ENORM
00188 *
00189 *     Compute RESULT(1)
00190 *
00191       RESULT( 1 ) = ERRNRM / ULP
00192 *
00193 *     Normalization of E:
00194 *
00195       ENRMER = ZERO
00196       DO 30 JVEC = 1, N
00197          TEMP1 = ZERO
00198          DO 20 J = 1, N
00199             TEMP1 = MAX( TEMP1, ABS1( E( J, JVEC ) ) )
00200    20    CONTINUE
00201          ENRMER = MAX( ENRMER, TEMP1-ONE )
00202    30 CONTINUE
00203 *
00204 *     Compute RESULT(2) : the normalization error in E.
00205 *
00206       RESULT( 2 ) = ENRMER / ( DBLE( N )*ULP )
00207 *
00208       RETURN
00209 *
00210 *     End of ZGET52
00211 *
00212       END