TESTING/EIG/sget51.f

   1 *> \brief \b SGET51
   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 SGET51( ITYPE, N, A, LDA, B, LDB, U, LDU, V, LDV, WORK,
  12 *                          RESULT )
  13 *
  14 *       .. Scalar Arguments ..
  15 *       INTEGER            ITYPE, LDA, LDB, LDU, LDV, N
  16 *       REAL               RESULT
  17 *       ..
  18 *       .. Array Arguments ..
  19 *       REAL               A( LDA, * ), B( LDB, * ), U( LDU, * ),
  20 *      $                   V( LDV, * ), WORK( * )
  21 *       ..
  22 *
  23 *
  24 *> \par Purpose:
  25 *  =============
  26 *>
  27 *> \verbatim
  28 *>
  29 *>      SGET51  generally checks a decomposition of the form
  30 *>
  31 *>              A = U B V'
  32 *>
  33 *>      where ' means transpose and U and V are orthogonal.
  34 *>
  35 *>      Specifically, if ITYPE=1
  36 *>
  37 *>              RESULT = | A - U B V' | / ( |A| n ulp )
  38 *>
  39 *>      If ITYPE=2, then:
  40 *>
  41 *>              RESULT = | A - B | / ( |A| n ulp )
  42 *>
  43 *>      If ITYPE=3, then:
  44 *>
  45 *>              RESULT = | I - UU' | / ( n ulp )
  46 *> \endverbatim
  47 *
  48 *  Arguments:
  49 *  ==========
  50 *
  51 *> \param[in] ITYPE
  52 *> \verbatim
  53 *>          ITYPE is INTEGER
  54 *>          Specifies the type of tests to be performed.
  55 *>          =1: RESULT = | A - U B V' | / ( |A| n ulp )
  56 *>          =2: RESULT = | A - B | / ( |A| n ulp )
  57 *>          =3: RESULT = | I - UU' | / ( n ulp )
  58 *> \endverbatim
  59 *>
  60 *> \param[in] N
  61 *> \verbatim
  62 *>          N is INTEGER
  63 *>          The size of the matrix.  If it is zero, SGET51 does nothing.
  64 *>          It must be at least zero.
  65 *> \endverbatim
  66 *>
  67 *> \param[in] A
  68 *> \verbatim
  69 *>          A is REAL array, dimension (LDA, N)
  70 *>          The original (unfactored) matrix.
  71 *> \endverbatim
  72 *>
  73 *> \param[in] LDA
  74 *> \verbatim
  75 *>          LDA is INTEGER
  76 *>          The leading dimension of A.  It must be at least 1
  77 *>          and at least N.
  78 *> \endverbatim
  79 *>
  80 *> \param[in] B
  81 *> \verbatim
  82 *>          B is REAL array, dimension (LDB, N)
  83 *>          The factored matrix.
  84 *> \endverbatim
  85 *>
  86 *> \param[in] LDB
  87 *> \verbatim
  88 *>          LDB is INTEGER
  89 *>          The leading dimension of B.  It must be at least 1
  90 *>          and at least N.
  91 *> \endverbatim
  92 *>
  93 *> \param[in] U
  94 *> \verbatim
  95 *>          U is REAL array, dimension (LDU, N)
  96 *>          The orthogonal matrix on the left-hand side in the
  97 *>          decomposition.
  98 *>          Not referenced if ITYPE=2
  99 *> \endverbatim
 100 *>
 101 *> \param[in] LDU
 102 *> \verbatim
 103 *>          LDU is INTEGER
 104 *>          The leading dimension of U.  LDU must be at least N and
 105 *>          at least 1.
 106 *> \endverbatim
 107 *>
 108 *> \param[in] V
 109 *> \verbatim
 110 *>          V is REAL array, dimension (LDV, N)
 111 *>          The orthogonal matrix on the left-hand side in the
 112 *>          decomposition.
 113 *>          Not referenced if ITYPE=2
 114 *> \endverbatim
 115 *>
 116 *> \param[in] LDV
 117 *> \verbatim
 118 *>          LDV is INTEGER
 119 *>          The leading dimension of V.  LDV must be at least N and
 120 *>          at least 1.
 121 *> \endverbatim
 122 *>
 123 *> \param[out] WORK
 124 *> \verbatim
 125 *>          WORK is REAL array, dimension (2*N**2)
 126 *> \endverbatim
 127 *>
 128 *> \param[out] RESULT
 129 *> \verbatim
 130 *>          RESULT is REAL
 131 *>          The values computed by the test specified by ITYPE.  The
 132 *>          value is currently limited to 1/ulp, to avoid overflow.
 133 *>          Errors are flagged by RESULT=10/ulp.
 134 *> \endverbatim
 135 *
 136 *  Authors:
 137 *  ========
 138 *
 139 *> \author Univ. of Tennessee
 140 *> \author Univ. of California Berkeley
 141 *> \author Univ. of Colorado Denver
 142 *> \author NAG Ltd.
 143 *
 144 *> \date November 2011
 145 *
 146 *> \ingroup single_eig
 147 *
 148 *  =====================================================================
 149       SUBROUTINE SGET51( ITYPE, N, A, LDA, B, LDB, U, LDU, V, LDV, WORK,
 150      $                   RESULT )
 151 *
 152 *  -- LAPACK test routine (version 3.4.0) --
 153 *  -- LAPACK is a software package provided by Univ. of Tennessee,    --
 154 *  -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
 155 *     November 2011
 156 *
 157 *     .. Scalar Arguments ..
 158       INTEGER            ITYPE, LDA, LDB, LDU, LDV, N
 159       REAL               RESULT
 160 *     ..
 161 *     .. Array Arguments ..
 162       REAL               A( LDA, * ), B( LDB, * ), U( LDU, * ),
 163      $                   V( LDV, * ), WORK( * )
 164 *     ..
 165 *
 166 *  =====================================================================
 167 *
 168 *     .. Parameters ..
 169       REAL               ZERO, ONE, TEN
 170       PARAMETER          ( ZERO = 0.0, ONE = 1.0E0, TEN = 10.0E0 )
 171 *     ..
 172 *     .. Local Scalars ..
 173       INTEGER            JCOL, JDIAG, JROW
 174       REAL               ANORM, ULP, UNFL, WNORM
 175 *     ..
 176 *     .. External Functions ..
 177       REAL               SLAMCH, SLANGE
 178       EXTERNAL           SLAMCH, SLANGE
 179 *     ..
 180 *     .. External Subroutines ..
 181       EXTERNAL           SGEMM, SLACPY
 182 *     ..
 183 *     .. Intrinsic Functions ..
 184       INTRINSIC          MAX, MIN, REAL
 185 *     ..
 186 *     .. Executable Statements ..
 187 *
 188       RESULT = ZERO
 189       IF( N.LE.0 )
 190      $   RETURN
 191 *
 192 *     Constants
 193 *
 194       UNFL = SLAMCH( 'Safe minimum' )
 195       ULP = SLAMCH( 'Epsilon' )*SLAMCH( 'Base' )
 196 *
 197 *     Some Error Checks
 198 *
 199       IF( ITYPE.LT.1 .OR. ITYPE.GT.3 ) THEN
 200          RESULT = TEN / ULP
 201          RETURN
 202       END IF
 203 *
 204       IF( ITYPE.LE.2 ) THEN
 205 *
 206 *        Tests scaled by the norm(A)
 207 *
 208          ANORM = MAX( SLANGE( '1', N, N, A, LDA, WORK ), UNFL )
 209 *
 210          IF( ITYPE.EQ.1 ) THEN
 211 *
 212 *           ITYPE=1: Compute W = A - UBV'
 213 *
 214             CALL SLACPY( ' ', N, N, A, LDA, WORK, N )
 215             CALL SGEMM( 'N', 'N', N, N, N, ONE, U, LDU, B, LDB, ZERO,
 216      $                  WORK( N**2+1 ), N )
 217 *
 218             CALL SGEMM( 'N', 'C', N, N, N, -ONE, WORK( N**2+1 ), N, V,
 219      $                  LDV, ONE, WORK, N )
 220 *
 221          ELSE
 222 *
 223 *           ITYPE=2: Compute W = A - B
 224 *
 225             CALL SLACPY( ' ', N, N, B, LDB, WORK, N )
 226 *
 227             DO 20 JCOL = 1, N
 228                DO 10 JROW = 1, N
 229                   WORK( JROW+N*( JCOL-1 ) ) = WORK( JROW+N*( JCOL-1 ) )
 230      $                - A( JROW, JCOL )
 231    10          CONTINUE
 232    20       CONTINUE
 233          END IF
 234 *
 235 *        Compute norm(W)/ ( ulp*norm(A) )
 236 *
 237          WNORM = SLANGE( '1', N, N, WORK, N, WORK( N**2+1 ) )
 238 *
 239          IF( ANORM.GT.WNORM ) THEN
 240             RESULT = ( WNORM / ANORM ) / ( N*ULP )
 241          ELSE
 242             IF( ANORM.LT.ONE ) THEN
 243                RESULT = ( MIN( WNORM, N*ANORM ) / ANORM ) / ( N*ULP )
 244             ELSE
 245                RESULT = MIN( WNORM / ANORM, REAL( N ) ) / ( N*ULP )
 246             END IF
 247          END IF
 248 *
 249       ELSE
 250 *
 251 *        Tests not scaled by norm(A)
 252 *
 253 *        ITYPE=3: Compute  UU' - I
 254 *
 255          CALL SGEMM( 'N', 'C', N, N, N, ONE, U, LDU, U, LDU, ZERO, WORK,
 256      $               N )
 257 *
 258          DO 30 JDIAG = 1, N
 259             WORK( ( N+1 )*( JDIAG-1 )+1 ) = WORK( ( N+1 )*( JDIAG-1 )+
 260      $         1 ) - ONE
 261    30    CONTINUE
 262 *
 263          RESULT = MIN( SLANGE( '1', N, N, WORK, N, WORK( N**2+1 ) ),
 264      $            REAL( N ) ) / ( N*ULP )
 265       END IF
 266 *
 267       RETURN
 268 *
 269 *     End of SGET51
 270 *
 271       END