LAPACK: TESTING/EIG/dsgt01.f Source File

Go to the documentation of this file.
00001       SUBROUTINE DSGT01( ITYPE, UPLO, N, M, A, LDA, B, LDB, Z, LDZ, D,
00002      $                   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 *     modified August 1997, a new parameter M is added to the calling
00009 *     sequence.
00010 *
00011 *     .. Scalar Arguments ..
00012       CHARACTER          UPLO
00013       INTEGER            ITYPE, LDA, LDB, LDZ, M, N
00014 *     ..
00015 *     .. Array Arguments ..
00016       DOUBLE PRECISION   A( LDA, * ), B( LDB, * ), D( * ), RESULT( * ),
00017      $                   WORK( * ), Z( LDZ, * )
00018 *     ..
00019 *
00020 *  Purpose
00021 *  =======
00022 *
00023 *  DDGT01 checks a decomposition of the form
00024 *
00025 *     A Z   =  B Z D or
00026 *     A B Z =  Z D or
00027 *     B A Z =  Z D
00028 *
00029 *  where A is a symmetric matrix, B is
00030 *  symmetric positive definite, Z is orthogonal, and D is diagonal.
00031 *
00032 *  One of the following test ratios is computed:
00033 *
00034 *  ITYPE = 1:  RESULT(1) = | A Z - B Z D | / ( |A| |Z| n ulp )
00035 *
00036 *  ITYPE = 2:  RESULT(1) = | A B Z - Z D | / ( |A| |Z| n ulp )
00037 *
00038 *  ITYPE = 3:  RESULT(1) = | B A Z - Z D | / ( |A| |Z| n ulp )
00039 *
00040 *  Arguments
00041 *  =========
00042 *
00043 *  ITYPE   (input) INTEGER
00044 *          The form of the symmetric generalized eigenproblem.
00045 *          = 1:  A*z = (lambda)*B*z
00046 *          = 2:  A*B*z = (lambda)*z
00047 *          = 3:  B*A*z = (lambda)*z
00048 *
00049 *  UPLO    (input) CHARACTER*1
00050 *          Specifies whether the upper or lower triangular part of the
00051 *          symmetric matrices A and B is stored.
00052 *          = 'U':  Upper triangular
00053 *          = 'L':  Lower triangular
00054 *
00055 *  N       (input) INTEGER
00056 *          The order of the matrix A.  N >= 0.
00057 *
00058 *  M       (input) INTEGER
00059 *          The number of eigenvalues found.  0 <= M <= N.
00060 *
00061 *  A       (input) DOUBLE PRECISION array, dimension (LDA, N)
00062 *          The original symmetric matrix A.
00063 *
00064 *  LDA     (input) INTEGER
00065 *          The leading dimension of the array A.  LDA >= max(1,N).
00066 *
00067 *  B       (input) DOUBLE PRECISION array, dimension (LDB, N)
00068 *          The original symmetric positive definite matrix B.
00069 *
00070 *  LDB     (input) INTEGER
00071 *          The leading dimension of the array B.  LDB >= max(1,N).
00072 *
00073 *  Z       (input) DOUBLE PRECISION array, dimension (LDZ, M)
00074 *          The computed eigenvectors of the generalized eigenproblem.
00075 *
00076 *  LDZ     (input) INTEGER
00077 *          The leading dimension of the array Z.  LDZ >= max(1,N).
00078 *
00079 *  D       (input) DOUBLE PRECISION array, dimension (M)
00080 *          The computed eigenvalues of the generalized eigenproblem.
00081 *
00082 *  WORK    (workspace) DOUBLE PRECISION array, dimension (N*N)
00083 *
00084 *  RESULT  (output) DOUBLE PRECISION array, dimension (1)
00085 *          The test ratio as described above.
00086 *
00087 *  =====================================================================
00088 *
00089 *     .. Parameters ..
00090       DOUBLE PRECISION   ZERO, ONE
00091       PARAMETER          ( ZERO = 0.0D0, ONE = 1.0D0 )
00092 *     ..
00093 *     .. Local Scalars ..
00094       INTEGER            I
00095       DOUBLE PRECISION   ANORM, ULP
00096 *     ..
00097 *     .. External Functions ..
00098       DOUBLE PRECISION   DLAMCH, DLANGE, DLANSY
00099       EXTERNAL           DLAMCH, DLANGE, DLANSY
00100 *     ..
00101 *     .. External Subroutines ..
00102       EXTERNAL           DSCAL, DSYMM
00103 *     ..
00104 *     .. Executable Statements ..
00105 *
00106       RESULT( 1 ) = ZERO
00107       IF( N.LE.0 )
00108      $   RETURN
00109 *
00110       ULP = DLAMCH( 'Epsilon' )
00111 *
00112 *     Compute product of 1-norms of A and Z.
00113 *
00114       ANORM = DLANSY( '1', UPLO, N, A, LDA, WORK )*
00115      $        DLANGE( '1', N, M, Z, LDZ, WORK )
00116       IF( ANORM.EQ.ZERO )
00117      $   ANORM = ONE
00118 *
00119       IF( ITYPE.EQ.1 ) THEN
00120 *
00121 *        Norm of AZ - BZD
00122 *
00123          CALL DSYMM( 'Left', UPLO, N, M, ONE, A, LDA, Z, LDZ, ZERO,
00124      $               WORK, N )
00125          DO 10 I = 1, M
00126             CALL DSCAL( N, D( I ), Z( 1, I ), 1 )
00127    10    CONTINUE
00128          CALL DSYMM( 'Left', UPLO, N, M, ONE, B, LDB, Z, LDZ, -ONE,
00129      $               WORK, N )
00130 *
00131          RESULT( 1 ) = ( DLANGE( '1', N, M, WORK, N, WORK ) / ANORM ) /
00132      $                 ( N*ULP )
00133 *
00134       ELSE IF( ITYPE.EQ.2 ) THEN
00135 *
00136 *        Norm of ABZ - ZD
00137 *
00138          CALL DSYMM( 'Left', UPLO, N, M, ONE, B, LDB, Z, LDZ, ZERO,
00139      $               WORK, N )
00140          DO 20 I = 1, M
00141             CALL DSCAL( N, D( I ), Z( 1, I ), 1 )
00142    20    CONTINUE
00143          CALL DSYMM( 'Left', UPLO, N, M, ONE, A, LDA, WORK, N, -ONE, Z,
00144      $               LDZ )
00145 *
00146          RESULT( 1 ) = ( DLANGE( '1', N, M, Z, LDZ, WORK ) / ANORM ) /
00147      $                 ( N*ULP )
00148 *
00149       ELSE IF( ITYPE.EQ.3 ) THEN
00150 *
00151 *        Norm of BAZ - ZD
00152 *
00153          CALL DSYMM( 'Left', UPLO, N, M, ONE, A, LDA, Z, LDZ, ZERO,
00154      $               WORK, N )
00155          DO 30 I = 1, M
00156             CALL DSCAL( N, D( I ), Z( 1, I ), 1 )
00157    30    CONTINUE
00158          CALL DSYMM( 'Left', UPLO, N, M, ONE, B, LDB, WORK, N, -ONE, Z,
00159      $               LDZ )
00160 *
00161          RESULT( 1 ) = ( DLANGE( '1', N, M, Z, LDZ, WORK ) / ANORM ) /
00162      $                 ( N*ULP )
00163       END IF
00164 *
00165       RETURN
00166 *
00167 *     End of DDGT01
00168 *
00169       END