TESTING/EIG/dsgt01.f

   1 *> \brief \b DSGT01
   2 *
   3 *  =========== DOCUMENTATION ===========
   4 *
   5 * Online html documentation available at
   6 *            http://www.netlib.org/lapack/explore-html/
   7 *
   8 *  Definition:
   9 *  ===========
  10 *
  11 *       SUBROUTINE DSGT01( ITYPE, UPLO, N, M, A, LDA, B, LDB, Z, LDZ, D,
  12 *                          WORK, RESULT )
  13 *
  14 *       .. Scalar Arguments ..
  15 *       CHARACTER          UPLO
  16 *       INTEGER            ITYPE, LDA, LDB, LDZ, M, N
  17 *       ..
  18 *       .. Array Arguments ..
  19 *       DOUBLE PRECISION   A( LDA, * ), B( LDB, * ), D( * ), RESULT( * ),
  20 *      $                   WORK( * ), Z( LDZ, * )
  21 *       ..
  22 *
  23 *
  24 *> \par Purpose:
  25 *  =============
  26 *>
  27 *> \verbatim
  28 *>
  29 *> DDGT01 checks a decomposition of the form
  30 *>
  31 *>    A Z   =  B Z D or
  32 *>    A B Z =  Z D or
  33 *>    B A Z =  Z D
  34 *>
  35 *> where A is a symmetric matrix, B is
  36 *> symmetric positive definite, Z is orthogonal, and D is diagonal.
  37 *>
  38 *> One of the following test ratios is computed:
  39 *>
  40 *> ITYPE = 1:  RESULT(1) = | A Z - B Z D | / ( |A| |Z| n ulp )
  41 *>
  42 *> ITYPE = 2:  RESULT(1) = | A B Z - Z D | / ( |A| |Z| n ulp )
  43 *>
  44 *> ITYPE = 3:  RESULT(1) = | B A Z - Z D | / ( |A| |Z| n ulp )
  45 *> \endverbatim
  46 *
  47 *  Arguments:
  48 *  ==========
  49 *
  50 *> \param[in] ITYPE
  51 *> \verbatim
  52 *>          ITYPE is INTEGER
  53 *>          The form of the symmetric generalized eigenproblem.
  54 *>          = 1:  A*z = (lambda)*B*z
  55 *>          = 2:  A*B*z = (lambda)*z
  56 *>          = 3:  B*A*z = (lambda)*z
  57 *> \endverbatim
  58 *>
  59 *> \param[in] UPLO
  60 *> \verbatim
  61 *>          UPLO is CHARACTER*1
  62 *>          Specifies whether the upper or lower triangular part of the
  63 *>          symmetric matrices A and B is stored.
  64 *>          = 'U':  Upper triangular
  65 *>          = 'L':  Lower triangular
  66 *> \endverbatim
  67 *>
  68 *> \param[in] N
  69 *> \verbatim
  70 *>          N is INTEGER
  71 *>          The order of the matrix A.  N >= 0.
  72 *> \endverbatim
  73 *>
  74 *> \param[in] M
  75 *> \verbatim
  76 *>          M is INTEGER
  77 *>          The number of eigenvalues found.  0 <= M <= N.
  78 *> \endverbatim
  79 *>
  80 *> \param[in] A
  81 *> \verbatim
  82 *>          A is DOUBLE PRECISION array, dimension (LDA, N)
  83 *>          The original symmetric matrix A.
  84 *> \endverbatim
  85 *>
  86 *> \param[in] LDA
  87 *> \verbatim
  88 *>          LDA is INTEGER
  89 *>          The leading dimension of the array A.  LDA >= max(1,N).
  90 *> \endverbatim
  91 *>
  92 *> \param[in] B
  93 *> \verbatim
  94 *>          B is DOUBLE PRECISION array, dimension (LDB, N)
  95 *>          The original symmetric positive definite matrix B.
  96 *> \endverbatim
  97 *>
  98 *> \param[in] LDB
  99 *> \verbatim
 100 *>          LDB is INTEGER
 101 *>          The leading dimension of the array B.  LDB >= max(1,N).
 102 *> \endverbatim
 103 *>
 104 *> \param[in] Z
 105 *> \verbatim
 106 *>          Z is DOUBLE PRECISION array, dimension (LDZ, M)
 107 *>          The computed eigenvectors of the generalized eigenproblem.
 108 *> \endverbatim
 109 *>
 110 *> \param[in] LDZ
 111 *> \verbatim
 112 *>          LDZ is INTEGER
 113 *>          The leading dimension of the array Z.  LDZ >= max(1,N).
 114 *> \endverbatim
 115 *>
 116 *> \param[in] D
 117 *> \verbatim
 118 *>          D is DOUBLE PRECISION array, dimension (M)
 119 *>          The computed eigenvalues of the generalized eigenproblem.
 120 *> \endverbatim
 121 *>
 122 *> \param[out] WORK
 123 *> \verbatim
 124 *>          WORK is DOUBLE PRECISION array, dimension (N*N)
 125 *> \endverbatim
 126 *>
 127 *> \param[out] RESULT
 128 *> \verbatim
 129 *>          RESULT is DOUBLE PRECISION array, dimension (1)
 130 *>          The test ratio as described above.
 131 *> \endverbatim
 132 *
 133 *  Authors:
 134 *  ========
 135 *
 136 *> \author Univ. of Tennessee
 137 *> \author Univ. of California Berkeley
 138 *> \author Univ. of Colorado Denver
 139 *> \author NAG Ltd.
 140 *
 141 *> \date November 2011
 142 *
 143 *> \ingroup double_eig
 144 *
 145 *  =====================================================================
 146       SUBROUTINE DSGT01( ITYPE, UPLO, N, M, A, LDA, B, LDB, Z, LDZ, D,
 147      $                   WORK, RESULT )
 148 *
 149 *  -- LAPACK test routine (version 3.4.0) --
 150 *  -- LAPACK is a software package provided by Univ. of Tennessee,    --
 151 *  -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
 152 *     November 2011
 153 *
 154 *     .. Scalar Arguments ..
 155       CHARACTER          UPLO
 156       INTEGER            ITYPE, LDA, LDB, LDZ, M, N
 157 *     ..
 158 *     .. Array Arguments ..
 159       DOUBLE PRECISION   A( LDA, * ), B( LDB, * ), D( * ), RESULT( * ),
 160      $                   WORK( * ), Z( LDZ, * )
 161 *     ..
 162 *
 163 *  =====================================================================
 164 *
 165 *     .. Parameters ..
 166       DOUBLE PRECISION   ZERO, ONE
 167       PARAMETER          ( ZERO = 0.0D0, ONE = 1.0D0 )
 168 *     ..
 169 *     .. Local Scalars ..
 170       INTEGER            I
 171       DOUBLE PRECISION   ANORM, ULP
 172 *     ..
 173 *     .. External Functions ..
 174       DOUBLE PRECISION   DLAMCH, DLANGE, DLANSY
 175       EXTERNAL           DLAMCH, DLANGE, DLANSY
 176 *     ..
 177 *     .. External Subroutines ..
 178       EXTERNAL           DSCAL, DSYMM
 179 *     ..
 180 *     .. Executable Statements ..
 181 *
 182       RESULT( 1 ) = ZERO
 183       IF( N.LE.0 )
 184      $   RETURN
 185 *
 186       ULP = DLAMCH( 'Epsilon' )
 187 *
 188 *     Compute product of 1-norms of A and Z.
 189 *
 190       ANORM = DLANSY( '1', UPLO, N, A, LDA, WORK )*
 191      $        DLANGE( '1', N, M, Z, LDZ, WORK )
 192       IF( ANORM.EQ.ZERO )
 193      $   ANORM = ONE
 194 *
 195       IF( ITYPE.EQ.1 ) THEN
 196 *
 197 *        Norm of AZ - BZD
 198 *
 199          CALL DSYMM( 'Left', UPLO, N, M, ONE, A, LDA, Z, LDZ, ZERO,
 200      $               WORK, N )
 201          DO 10 I = 1, M
 202             CALL DSCAL( N, D( I ), Z( 1, I ), 1 )
 203    10    CONTINUE
 204          CALL DSYMM( 'Left', UPLO, N, M, ONE, B, LDB, Z, LDZ, -ONE,
 205      $               WORK, N )
 206 *
 207          RESULT( 1 ) = ( DLANGE( '1', N, M, WORK, N, WORK ) / ANORM ) /
 208      $                 ( N*ULP )
 209 *
 210       ELSE IF( ITYPE.EQ.2 ) THEN
 211 *
 212 *        Norm of ABZ - ZD
 213 *
 214          CALL DSYMM( 'Left', UPLO, N, M, ONE, B, LDB, Z, LDZ, ZERO,
 215      $               WORK, N )
 216          DO 20 I = 1, M
 217             CALL DSCAL( N, D( I ), Z( 1, I ), 1 )
 218    20    CONTINUE
 219          CALL DSYMM( 'Left', UPLO, N, M, ONE, A, LDA, WORK, N, -ONE, Z,
 220      $               LDZ )
 221 *
 222          RESULT( 1 ) = ( DLANGE( '1', N, M, Z, LDZ, WORK ) / ANORM ) /
 223      $                 ( N*ULP )
 224 *
 225       ELSE IF( ITYPE.EQ.3 ) THEN
 226 *
 227 *        Norm of BAZ - ZD
 228 *
 229          CALL DSYMM( 'Left', UPLO, N, M, ONE, A, LDA, Z, LDZ, ZERO,
 230      $               WORK, N )
 231          DO 30 I = 1, M
 232             CALL DSCAL( N, D( I ), Z( 1, I ), 1 )
 233    30    CONTINUE
 234          CALL DSYMM( 'Left', UPLO, N, M, ONE, B, LDB, WORK, N, -ONE, Z,
 235      $               LDZ )
 236 *
 237          RESULT( 1 ) = ( DLANGE( '1', N, M, Z, LDZ, WORK ) / ANORM ) /
 238      $                 ( N*ULP )
 239       END IF
 240 *
 241       RETURN
 242 *
 243 *     End of DDGT01
 244 *
 245       END