1       SUBROUTINE DSGT01( ITYPE, UPLO, N, M, A, LDA, B, LDB, Z, LDZ, D,
  2      $                   WORK, RESULT )
  3 *
  4 *  -- LAPACK test routine (version 3.1) --
  5 *     Univ. of Tennessee, Univ. of California Berkeley and NAG Ltd..
  6 *     November 2006
  7 *
  8 *     modified August 1997, a new parameter M is added to the calling
  9 *     sequence.
 10 *
 11 *     .. Scalar Arguments ..
 12       CHARACTER          UPLO
 13       INTEGER            ITYPE, LDA, LDB, LDZ, M, N
 14 *     ..
 15 *     .. Array Arguments ..
 16       DOUBLE PRECISION   A( LDA, * ), B( LDB, * ), D( * ), RESULT* ),
 17      $                   WORK( * ), Z( LDZ, * )
 18 *     ..
 19 *
 20 *  Purpose
 21 *  =======
 22 *
 23 *  DDGT01 checks a decomposition of the form
 24 *
 25 *     A Z   =  B Z D or
 26 *     A B Z =  Z D or
 27 *     B A Z =  Z D
 28 *
 29 *  where A is a symmetric matrix, B is
 30 *  symmetric positive definite, Z is orthogonal, and D is diagonal.
 31 *
 32 *  One of the following test ratios is computed:
 33 *
 34 *  ITYPE = 1:  RESULT(1) = | A Z - B Z D | / ( |A| |Z| n ulp )
 35 *
 36 *  ITYPE = 2:  RESULT(1) = | A B Z - Z D | / ( |A| |Z| n ulp )
 37 *
 38 *  ITYPE = 3:  RESULT(1) = | B A Z - Z D | / ( |A| |Z| n ulp )
 39 *
 40 *  Arguments
 41 *  =========
 42 *
 43 *  ITYPE   (input) INTEGER
 44 *          The form of the symmetric generalized eigenproblem.
 45 *          = 1:  A*z = (lambda)*B*z
 46 *          = 2:  A*B*z = (lambda)*z
 47 *          = 3:  B*A*z = (lambda)*z
 48 *
 49 *  UPLO    (input) CHARACTER*1
 50 *          Specifies whether the upper or lower triangular part of the
 51 *          symmetric matrices A and B is stored.
 52 *          = 'U':  Upper triangular
 53 *          = 'L':  Lower triangular
 54 *
 55 *  N       (input) INTEGER
 56 *          The order of the matrix A.  N >= 0.
 57 *
 58 *  M       (input) INTEGER
 59 *          The number of eigenvalues found.  0 <= M <= N.
 60 *
 61 *  A       (input) DOUBLE PRECISION array, dimension (LDA, N)
 62 *          The original symmetric matrix A.
 63 *
 64 *  LDA     (input) INTEGER
 65 *          The leading dimension of the array A.  LDA >= max(1,N).
 66 *
 67 *  B       (input) DOUBLE PRECISION array, dimension (LDB, N)
 68 *          The original symmetric positive definite matrix B.
 69 *
 70 *  LDB     (input) INTEGER
 71 *          The leading dimension of the array B.  LDB >= max(1,N).
 72 *
 73 *  Z       (input) DOUBLE PRECISION array, dimension (LDZ, M)
 74 *          The computed eigenvectors of the generalized eigenproblem.
 75 *
 76 *  LDZ     (input) INTEGER
 77 *          The leading dimension of the array Z.  LDZ >= max(1,N).
 78 *
 79 *  D       (input) DOUBLE PRECISION array, dimension (M)
 80 *          The computed eigenvalues of the generalized eigenproblem.
 81 *
 82 *  WORK    (workspace) DOUBLE PRECISION array, dimension (N*N)
 83 *
 84 *  RESULT  (output) DOUBLE PRECISION array, dimension (1)
 85 *          The test ratio as described above.
 86 *
 87 *  =====================================================================
 88 *
 89 *     .. Parameters ..
 90       DOUBLE PRECISION   ZERO, ONE
 91       PARAMETER          ( ZERO = 0.0D0, ONE = 1.0D0 )
 92 *     ..
 93 *     .. Local Scalars ..
 94       INTEGER            I
 95       DOUBLE PRECISION   ANORM, ULP
 96 *     ..
 97 *     .. External Functions ..
 98       DOUBLE PRECISION   DLAMCH, DLANGE, DLANSY
 99       EXTERNAL           DLAMCH, DLANGE, DLANSY
100 *     ..
101 *     .. External Subroutines ..
102       EXTERNAL           DSCAL, DSYMM
103 *     ..
104 *     .. Executable Statements ..
105 *
106       RESULT1 ) = ZERO
107       IF( N.LE.0 )
108      $   RETURN
109 *
110       ULP = DLAMCH( 'Epsilon' )
111 *
112 *     Compute product of 1-norms of A and Z.
113 *
114       ANORM = DLANSY( '1', UPLO, N, A, LDA, WORK )*
115      $        DLANGE( '1', N, M, Z, LDZ, WORK )
116       IF( ANORM.EQ.ZERO )
117      $   ANORM = ONE
118 *
119       IF( ITYPE.EQ.1 ) THEN
120 *
121 *        Norm of AZ - BZD
122 *
123          CALL DSYMM( 'Left', UPLO, N, M, ONE, A, LDA, Z, LDZ, ZERO,
124      $               WORK, N )
125          DO 10 I = 1, M
126             CALL DSCAL( N, D( I ), Z( 1, I ), 1 )
127    10    CONTINUE
128          CALL DSYMM( 'Left', UPLO, N, M, ONE, B, LDB, Z, LDZ, -ONE,
129      $               WORK, N )
130 *
131          RESULT1 ) = ( DLANGE( '1', N, M, WORK, N, WORK ) / ANORM ) /
132      $                 ( N*ULP )
133 *
134       ELSE IF( ITYPE.EQ.2 ) THEN
135 *
136 *        Norm of ABZ - ZD
137 *
138          CALL DSYMM( 'Left', UPLO, N, M, ONE, B, LDB, Z, LDZ, ZERO,
139      $               WORK, N )
140          DO 20 I = 1, M
141             CALL DSCAL( N, D( I ), Z( 1, I ), 1 )
142    20    CONTINUE
143          CALL DSYMM( 'Left', UPLO, N, M, ONE, A, LDA, WORK, N, -ONE, Z,
144      $               LDZ )
145 *
146          RESULT1 ) = ( DLANGE( '1', N, M, Z, LDZ, WORK ) / ANORM ) /
147      $                 ( N*ULP )
148 *
149       ELSE IF( ITYPE.EQ.3 ) THEN
150 *
151 *        Norm of BAZ - ZD
152 *
153          CALL DSYMM( 'Left', UPLO, N, M, ONE, A, LDA, Z, LDZ, ZERO,
154      $               WORK, N )
155          DO 30 I = 1, M
156             CALL DSCAL( N, D( I ), Z( 1, I ), 1 )
157    30    CONTINUE
158          CALL DSYMM( 'Left', UPLO, N, M, ONE, B, LDB, WORK, N, -ONE, Z,
159      $               LDZ )
160 *
161          RESULT1 ) = ( DLANGE( '1', N, M, Z, LDZ, WORK ) / ANORM ) /
162      $                 ( N*ULP )
163       END IF
164 *
165       RETURN
166 *
167 *     End of DDGT01
168 *
169       END