LAPACK: TESTING/MATGEN/zlatm6.f Source File

Go to the documentation of this file.
00001       SUBROUTINE ZLATM6( TYPE, N, A, LDA, B, X, LDX, Y, LDY, ALPHA,
00002      $                   BETA, WX, WY, S, DIF )
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       INTEGER            LDA, LDX, LDY, N, TYPE
00010       COMPLEX*16         ALPHA, BETA, WX, WY
00011 *     ..
00012 *     .. Array Arguments ..
00013       DOUBLE PRECISION   DIF( * ), S( * )
00014       COMPLEX*16         A( LDA, * ), B( LDA, * ), X( LDX, * ),
00015      $                   Y( LDY, * )
00016 *     ..
00017 *
00018 *  Purpose
00019 *  =======
00020 *
00021 *  ZLATM6 generates test matrices for the generalized eigenvalue
00022 *  problem, their corresponding right and left eigenvector matrices,
00023 *  and also reciprocal condition numbers for all eigenvalues and
00024 *  the reciprocal condition numbers of eigenvectors corresponding to
00025 *  the 1th and 5th eigenvalues.
00026 *
00027 *  Test Matrices
00028 *  =============
00029 *
00030 *  Two kinds of test matrix pairs
00031 *           (A, B) = inverse(YH) * (Da, Db) * inverse(X)
00032 *  are used in the tests:
00033 *
00034 *  Type 1:
00035 *     Da = 1+a   0    0    0    0    Db = 1   0   0   0   0
00036 *           0   2+a   0    0    0         0   1   0   0   0
00037 *           0    0   3+a   0    0         0   0   1   0   0
00038 *           0    0    0   4+a   0         0   0   0   1   0
00039 *           0    0    0    0   5+a ,      0   0   0   0   1
00040 *  and Type 2:
00041 *     Da = 1+i   0    0       0       0    Db = 1   0   0   0   0
00042 *           0   1-i   0       0       0         0   1   0   0   0
00043 *           0    0    1       0       0         0   0   1   0   0
00044 *           0    0    0 (1+a)+(1+b)i  0         0   0   0   1   0
00045 *           0    0    0       0 (1+a)-(1+b)i,   0   0   0   0   1 .
00046 *
00047 *  In both cases the same inverse(YH) and inverse(X) are used to compute
00048 *  (A, B), giving the exact eigenvectors to (A,B) as (YH, X):
00049 *
00050 *  YH:  =  1    0   -y    y   -y    X =  1   0  -x  -x   x
00051 *          0    1   -y    y   -y         0   1   x  -x  -x
00052 *          0    0    1    0    0         0   0   1   0   0
00053 *          0    0    0    1    0         0   0   0   1   0
00054 *          0    0    0    0    1,        0   0   0   0   1 , where
00055 *
00056 *  a, b, x and y will have all values independently of each other.
00057 *
00058 *  Arguments
00059 *  =========
00060 *
00061 *  TYPE    (input) INTEGER
00062 *          Specifies the problem type (see futher details).
00063 *
00064 *  N       (input) INTEGER
00065 *          Size of the matrices A and B.
00066 *
00067 *  A       (output) COMPLEX*16 array, dimension (LDA, N).
00068 *          On exit A N-by-N is initialized according to TYPE.
00069 *
00070 *  LDA     (input) INTEGER
00071 *          The leading dimension of A and of B.
00072 *
00073 *  B       (output) COMPLEX*16 array, dimension (LDA, N).
00074 *          On exit B N-by-N is initialized according to TYPE.
00075 *
00076 *  X       (output) COMPLEX*16 array, dimension (LDX, N).
00077 *          On exit X is the N-by-N matrix of right eigenvectors.
00078 *
00079 *  LDX     (input) INTEGER
00080 *          The leading dimension of X.
00081 *
00082 *  Y       (output) COMPLEX*16 array, dimension (LDY, N).
00083 *          On exit Y is the N-by-N matrix of left eigenvectors.
00084 *
00085 *  LDY     (input) INTEGER
00086 *          The leading dimension of Y.
00087 *
00088 *  ALPHA   (input) COMPLEX*16
00089 *  BETA    (input) COMPLEX*16
00090 *          Weighting constants for matrix A.
00091 *
00092 *  WX      (input) COMPLEX*16
00093 *          Constant for right eigenvector matrix.
00094 *
00095 *  WY      (input) COMPLEX*16
00096 *          Constant for left eigenvector matrix.
00097 *
00098 *  S       (output) DOUBLE PRECISION array, dimension (N)
00099 *          S(i) is the reciprocal condition number for eigenvalue i.
00100 *
00101 *  DIF     (output) DOUBLE PRECISION array, dimension (N)
00102 *          DIF(i) is the reciprocal condition number for eigenvector i.
00103 *
00104 *  =====================================================================
00105 *
00106 *     .. Parameters ..
00107       DOUBLE PRECISION   RONE, TWO, THREE
00108       PARAMETER          ( RONE = 1.0D+0, TWO = 2.0D+0, THREE = 3.0D+0 )
00109       COMPLEX*16         ZERO, ONE
00110       PARAMETER          ( ZERO = ( 0.0D+0, 0.0D+0 ),
00111      $                   ONE = ( 1.0D+0, 0.0D+0 ) )
00112 *     ..
00113 *     .. Local Scalars ..
00114       INTEGER            I, INFO, J
00115 *     ..
00116 *     .. Local Arrays ..
00117       DOUBLE PRECISION   RWORK( 50 )
00118       COMPLEX*16         WORK( 26 ), Z( 8, 8 )
00119 *     ..
00120 *     .. Intrinsic Functions ..
00121       INTRINSIC          CDABS, DBLE, DCMPLX, DCONJG, SQRT
00122 *     ..
00123 *     .. External Subroutines ..
00124       EXTERNAL           ZGESVD, ZLACPY, ZLAKF2
00125 *     ..
00126 *     .. Executable Statements ..
00127 *
00128 *     Generate test problem ...
00129 *     (Da, Db) ...
00130 *
00131       DO 20 I = 1, N
00132          DO 10 J = 1, N
00133 *
00134             IF( I.EQ.J ) THEN
00135                A( I, I ) = DCMPLX( I ) + ALPHA
00136                B( I, I ) = ONE
00137             ELSE
00138                A( I, J ) = ZERO
00139                B( I, J ) = ZERO
00140             END IF
00141 *
00142    10    CONTINUE
00143    20 CONTINUE
00144       IF( TYPE.EQ.2 ) THEN
00145          A( 1, 1 ) = DCMPLX( RONE, RONE )
00146          A( 2, 2 ) = DCONJG( A( 1, 1 ) )
00147          A( 3, 3 ) = ONE
00148          A( 4, 4 ) = DCMPLX( DBLE( ONE+ALPHA ), DBLE( ONE+BETA ) )
00149          A( 5, 5 ) = DCONJG( A( 4, 4 ) )
00150       END IF
00151 *
00152 *     Form X and Y
00153 *
00154       CALL ZLACPY( 'F', N, N, B, LDA, Y, LDY )
00155       Y( 3, 1 ) = -DCONJG( WY )
00156       Y( 4, 1 ) = DCONJG( WY )
00157       Y( 5, 1 ) = -DCONJG( WY )
00158       Y( 3, 2 ) = -DCONJG( WY )
00159       Y( 4, 2 ) = DCONJG( WY )
00160       Y( 5, 2 ) = -DCONJG( WY )
00161 *
00162       CALL ZLACPY( 'F', N, N, B, LDA, X, LDX )
00163       X( 1, 3 ) = -WX
00164       X( 1, 4 ) = -WX
00165       X( 1, 5 ) = WX
00166       X( 2, 3 ) = WX
00167       X( 2, 4 ) = -WX
00168       X( 2, 5 ) = -WX
00169 *
00170 *     Form (A, B)
00171 *
00172       B( 1, 3 ) = WX + WY
00173       B( 2, 3 ) = -WX + WY
00174       B( 1, 4 ) = WX - WY
00175       B( 2, 4 ) = WX - WY
00176       B( 1, 5 ) = -WX + WY
00177       B( 2, 5 ) = WX + WY
00178       A( 1, 3 ) = WX*A( 1, 1 ) + WY*A( 3, 3 )
00179       A( 2, 3 ) = -WX*A( 2, 2 ) + WY*A( 3, 3 )
00180       A( 1, 4 ) = WX*A( 1, 1 ) - WY*A( 4, 4 )
00181       A( 2, 4 ) = WX*A( 2, 2 ) - WY*A( 4, 4 )
00182       A( 1, 5 ) = -WX*A( 1, 1 ) + WY*A( 5, 5 )
00183       A( 2, 5 ) = WX*A( 2, 2 ) + WY*A( 5, 5 )
00184 *
00185 *     Compute condition numbers
00186 *
00187       S( 1 ) = RONE / SQRT( ( RONE+THREE*CDABS( WY )*CDABS( WY ) ) /
00188      $         ( RONE+CDABS( A( 1, 1 ) )*CDABS( A( 1, 1 ) ) ) )
00189       S( 2 ) = RONE / SQRT( ( RONE+THREE*CDABS( WY )*CDABS( WY ) ) /
00190      $         ( RONE+CDABS( A( 2, 2 ) )*CDABS( A( 2, 2 ) ) ) )
00191       S( 3 ) = RONE / SQRT( ( RONE+TWO*CDABS( WX )*CDABS( WX ) ) /
00192      $         ( RONE+CDABS( A( 3, 3 ) )*CDABS( A( 3, 3 ) ) ) )
00193       S( 4 ) = RONE / SQRT( ( RONE+TWO*CDABS( WX )*CDABS( WX ) ) /
00194      $         ( RONE+CDABS( A( 4, 4 ) )*CDABS( A( 4, 4 ) ) ) )
00195       S( 5 ) = RONE / SQRT( ( RONE+TWO*CDABS( WX )*CDABS( WX ) ) /
00196      $         ( RONE+CDABS( A( 5, 5 ) )*CDABS( A( 5, 5 ) ) ) )
00197 *
00198       CALL ZLAKF2( 1, 4, A, LDA, A( 2, 2 ), B, B( 2, 2 ), Z, 8 )
00199       CALL ZGESVD( 'N', 'N', 8, 8, Z, 8, RWORK, WORK, 1, WORK( 2 ), 1,
00200      $             WORK( 3 ), 24, RWORK( 9 ), INFO )
00201       DIF( 1 ) = RWORK( 8 )
00202 *
00203       CALL ZLAKF2( 4, 1, A, LDA, A( 5, 5 ), B, B( 5, 5 ), Z, 8 )
00204       CALL ZGESVD( 'N', 'N', 8, 8, Z, 8, RWORK, WORK, 1, WORK( 2 ), 1,
00205      $             WORK( 3 ), 24, RWORK( 9 ), INFO )
00206       DIF( 5 ) = RWORK( 8 )
00207 *
00208       RETURN
00209 *
00210 *     End of ZLATM6
00211 *
00212       END