LAPACK: TESTING/MATGEN/dlatm6.f Source File

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