TESTING/EIG/zcsdts.f

   1 *> \brief \b ZCSDTS
   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 ZCSDTS( M, P, Q, X, XF, LDX, U1, LDU1, U2, LDU2, V1T,
  12 *                          LDV1T, V2T, LDV2T, THETA, IWORK, WORK, LWORK,
  13 *                          RWORK, RESULT )
  14 *
  15 *       .. Scalar Arguments ..
  16 *       INTEGER            LDX, LDU1, LDU2, LDV1T, LDV2T, LWORK, M, P, Q
  17 *       ..
  18 *       .. Array Arguments ..
  19 *       INTEGER            IWORK( * )
  20 *       DOUBLE PRECISION   RESULT( 15 ), RWORK( * ), THETA( * )
  21 *       COMPLEX*16         U1( LDU1, * ), U2( LDU2, * ), V1T( LDV1T, * ),
  22 *      $                   V2T( LDV2T, * ), WORK( LWORK ), X( LDX, * ),
  23 *      $                   XF( LDX, * )
  24 *       ..
  25 *
  26 *
  27 *> \par Purpose:
  28 *  =============
  29 *>
  30 *> \verbatim
  31 *>
  32 *> ZCSDTS tests ZUNCSD, which, given an M-by-M partitioned unitary
  33 *> matrix X,
  34 *>              Q  M-Q
  35 *>       X = [ X11 X12 ] P   ,
  36 *>           [ X21 X22 ] M-P
  37 *>
  38 *> computes the CSD
  39 *>
  40 *>       [ U1    ]**T * [ X11 X12 ] * [ V1    ]
  41 *>       [    U2 ]      [ X21 X22 ]   [    V2 ]
  42 *>
  43 *>                             [  I  0  0 |  0  0  0 ]
  44 *>                             [  0  C  0 |  0 -S  0 ]
  45 *>                             [  0  0  0 |  0  0 -I ]
  46 *>                           = [---------------------] = [ D11 D12 ] .
  47 *>                             [  0  0  0 |  I  0  0 ]   [ D21 D22 ]
  48 *>                             [  0  S  0 |  0  C  0 ]
  49 *>                             [  0  0  I |  0  0  0 ]
  50 *>
  51 *> and also SORCSD2BY1, which, given
  52 *>          Q
  53 *>       [ X11 ] P   ,
  54 *>       [ X21 ] M-P
  55 *>
  56 *> computes the 2-by-1 CSD
  57 *>
  58 *>                                     [  I  0  0 ]
  59 *>                                     [  0  C  0 ]
  60 *>                                     [  0  0  0 ]
  61 *>       [ U1    ]**T * [ X11 ] * V1 = [----------] = [ D11 ] ,
  62 *>       [    U2 ]      [ X21 ]        [  0  0  0 ]   [ D21 ]
  63 *>                                     [  0  S  0 ]
  64 *>                                     [  0  0  I ]
  65 *> \endverbatim
  66 *
  67 *  Arguments:
  68 *  ==========
  69 *
  70 *> \param[in] M
  71 *> \verbatim
  72 *>          M is INTEGER
  73 *>          The number of rows of the matrix X.  M >= 0.
  74 *> \endverbatim
  75 *>
  76 *> \param[in] P
  77 *> \verbatim
  78 *>          P is INTEGER
  79 *>          The number of rows of the matrix X11.  P >= 0.
  80 *> \endverbatim
  81 *>
  82 *> \param[in] Q
  83 *> \verbatim
  84 *>          Q is INTEGER
  85 *>          The number of columns of the matrix X11.  Q >= 0.
  86 *> \endverbatim
  87 *>
  88 *> \param[in] X
  89 *> \verbatim
  90 *>          X is COMPLEX*16 array, dimension (LDX,M)
  91 *>          The M-by-M matrix X.
  92 *> \endverbatim
  93 *>
  94 *> \param[out] XF
  95 *> \verbatim
  96 *>          XF is COMPLEX*16 array, dimension (LDX,M)
  97 *>          Details of the CSD of X, as returned by ZUNCSD;
  98 *>          see ZUNCSD for further details.
  99 *> \endverbatim
 100 *>
 101 *> \param[in] LDX
 102 *> \verbatim
 103 *>          LDX is INTEGER
 104 *>          The leading dimension of the arrays X and XF.
 105 *>          LDX >= max( 1,M ).
 106 *> \endverbatim
 107 *>
 108 *> \param[out] U1
 109 *> \verbatim
 110 *>          U1 is COMPLEX*16 array, dimension(LDU1,P)
 111 *>          The P-by-P unitary matrix U1.
 112 *> \endverbatim
 113 *>
 114 *> \param[in] LDU1
 115 *> \verbatim
 116 *>          LDU1 is INTEGER
 117 *>          The leading dimension of the array U1. LDU >= max(1,P).
 118 *> \endverbatim
 119 *>
 120 *> \param[out] U2
 121 *> \verbatim
 122 *>          U2 is COMPLEX*16 array, dimension(LDU2,M-P)
 123 *>          The (M-P)-by-(M-P) unitary matrix U2.
 124 *> \endverbatim
 125 *>
 126 *> \param[in] LDU2
 127 *> \verbatim
 128 *>          LDU2 is INTEGER
 129 *>          The leading dimension of the array U2. LDU >= max(1,M-P).
 130 *> \endverbatim
 131 *>
 132 *> \param[out] V1T
 133 *> \verbatim
 134 *>          V1T is COMPLEX*16 array, dimension(LDV1T,Q)
 135 *>          The Q-by-Q unitary matrix V1T.
 136 *> \endverbatim
 137 *>
 138 *> \param[in] LDV1T
 139 *> \verbatim
 140 *>          LDV1T is INTEGER
 141 *>          The leading dimension of the array V1T. LDV1T >=
 142 *>          max(1,Q).
 143 *> \endverbatim
 144 *>
 145 *> \param[out] V2T
 146 *> \verbatim
 147 *>          V2T is COMPLEX*16 array, dimension(LDV2T,M-Q)
 148 *>          The (M-Q)-by-(M-Q) unitary matrix V2T.
 149 *> \endverbatim
 150 *>
 151 *> \param[in] LDV2T
 152 *> \verbatim
 153 *>          LDV2T is INTEGER
 154 *>          The leading dimension of the array V2T. LDV2T >=
 155 *>          max(1,M-Q).
 156 *> \endverbatim
 157 *>
 158 *> \param[out] THETA
 159 *> \verbatim
 160 *>          THETA is DOUBLE PRECISION array, dimension MIN(P,M-P,Q,M-Q)
 161 *>          The CS values of X; the essentially diagonal matrices C and
 162 *>          S are constructed from THETA; see subroutine ZUNCSD for
 163 *>          details.
 164 *> \endverbatim
 165 *>
 166 *> \param[out] IWORK
 167 *> \verbatim
 168 *>          IWORK is INTEGER array, dimension (M)
 169 *> \endverbatim
 170 *>
 171 *> \param[out] WORK
 172 *> \verbatim
 173 *>          WORK is COMPLEX*16 array, dimension (LWORK)
 174 *> \endverbatim
 175 *>
 176 *> \param[in] LWORK
 177 *> \verbatim
 178 *>          LWORK is INTEGER
 179 *>          The dimension of the array WORK
 180 *> \endverbatim
 181 *>
 182 *> \param[out] RWORK
 183 *> \verbatim
 184 *>          RWORK is DOUBLE PRECISION array
 185 *> \endverbatim
 186 *>
 187 *> \param[out] RESULT
 188 *> \verbatim
 189 *>          RESULT is DOUBLE PRECISION array, dimension (15)
 190 *>          The test ratios:
 191 *>          First, the 2-by-2 CSD:
 192 *>          RESULT(1) = norm( U1'*X11*V1 - D11 ) / ( MAX(1,P,Q)*EPS2 )
 193 *>          RESULT(2) = norm( U1'*X12*V2 - D12 ) / ( MAX(1,P,M-Q)*EPS2 )
 194 *>          RESULT(3) = norm( U2'*X21*V1 - D21 ) / ( MAX(1,M-P,Q)*EPS2 )
 195 *>          RESULT(4) = norm( U2'*X22*V2 - D22 ) / ( MAX(1,M-P,M-Q)*EPS2 )
 196 *>          RESULT(5) = norm( I - U1'*U1 ) / ( MAX(1,P)*ULP )
 197 *>          RESULT(6) = norm( I - U2'*U2 ) / ( MAX(1,M-P)*ULP )
 198 *>          RESULT(7) = norm( I - V1T'*V1T ) / ( MAX(1,Q)*ULP )
 199 *>          RESULT(8) = norm( I - V2T'*V2T ) / ( MAX(1,M-Q)*ULP )
 200 *>          RESULT(9) = 0        if THETA is in increasing order and
 201 *>                               all angles are in [0,pi/2];
 202 *>                    = ULPINV   otherwise.
 203 *>          Then, the 2-by-1 CSD:
 204 *>          RESULT(10) = norm( U1'*X11*V1 - D11 ) / ( MAX(1,P,Q)*EPS2 )
 205 *>          RESULT(11) = norm( U2'*X21*V1 - D21 ) / ( MAX(1,M-P,Q)*EPS2 )
 206 *>          RESULT(12) = norm( I - U1'*U1 ) / ( MAX(1,P)*ULP )
 207 *>          RESULT(13) = norm( I - U2'*U2 ) / ( MAX(1,M-P)*ULP )
 208 *>          RESULT(14) = norm( I - V1T'*V1T ) / ( MAX(1,Q)*ULP )
 209 *>          RESULT(15) = 0        if THETA is in increasing order and
 210 *>                                all angles are in [0,pi/2];
 211 *>                     = ULPINV   otherwise.
 212 *>          ( EPS2 = MAX( norm( I - X'*X ) / M, ULP ). )
 213 *> \endverbatim
 214 *
 215 *  Authors:
 216 *  ========
 217 *
 218 *> \author Univ. of Tennessee
 219 *> \author Univ. of California Berkeley
 220 *> \author Univ. of Colorado Denver
 221 *> \author NAG Ltd.
 222 *
 223 *> \date November 2015
 224 *
 225 *> \ingroup complex16_eig
 226 *
 227 *  =====================================================================
 228       SUBROUTINE ZCSDTS( M, P, Q, X, XF, LDX, U1, LDU1, U2, LDU2, V1T,
 229      $                   LDV1T, V2T, LDV2T, THETA, IWORK, WORK, LWORK,
 230      $                   RWORK, RESULT )
 231 *
 232 *  -- LAPACK test routine (version 3.6.0) --
 233 *  -- LAPACK is a software package provided by Univ. of Tennessee,    --
 234 *  -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
 235 *     November 2015
 236 *
 237 *     .. Scalar Arguments ..
 238       INTEGER            LDX, LDU1, LDU2, LDV1T, LDV2T, LWORK, M, P, Q
 239 *     ..
 240 *     .. Array Arguments ..
 241       INTEGER            IWORK( * )
 242       DOUBLE PRECISION   RESULT( 15 ), RWORK( * ), THETA( * )
 243       COMPLEX*16         U1( LDU1, * ), U2( LDU2, * ), V1T( LDV1T, * ),
 244      $                   V2T( LDV2T, * ), WORK( LWORK ), X( LDX, * ),
 245      $                   XF( LDX, * )
 246 *     ..
 247 *
 248 *  =====================================================================
 249 *
 250 *     .. Parameters ..
 251       DOUBLE PRECISION   PIOVER2, REALONE, REALZERO
 252       PARAMETER          ( PIOVER2 = 1.57079632679489662D0,
 253      $                     REALONE = 1.0D0, REALZERO = 0.0D0 )
 254       COMPLEX*16         ZERO, ONE
 255       PARAMETER          ( ZERO = (0.0D0,0.0D0), ONE = (1.0D0,0.0D0) )
 256 *     ..
 257 *     .. Local Scalars ..
 258       INTEGER            I, INFO, R
 259       DOUBLE PRECISION   EPS2, RESID, ULP, ULPINV
 260 *     ..
 261 *     .. External Functions ..
 262       DOUBLE PRECISION   DLAMCH, ZLANGE, ZLANHE
 263       EXTERNAL           DLAMCH, ZLANGE, ZLANHE
 264 *     ..
 265 *     .. External Subroutines ..
 266       EXTERNAL           ZGEMM, ZHERK, ZLACPY, ZLASET, ZUNCSD,
 267      $                   ZUNCSD2BY1
 268 *     ..
 269 *     .. Intrinsic Functions ..
 270       INTRINSIC          COS, DBLE, DCMPLX, MAX, MIN, REAL, SIN
 271 *     ..
 272 *     .. Executable Statements ..
 273 *
 274       ULP = DLAMCH( 'Precision' )
 275       ULPINV = REALONE / ULP
 276 *
 277 *     The first half of the routine checks the 2-by-2 CSD
 278 *
 279       CALL ZLASET( 'Full', M, M, ZERO, ONE, WORK, LDX )
 280       CALL ZHERK( 'Upper', 'Conjugate transpose', M, M, -REALONE,
 281      $            X, LDX, REALONE, WORK, LDX )
 282       IF (M.GT.0) THEN
 283          EPS2 = MAX( ULP,
 284      $               ZLANGE( '1', M, M, WORK, LDX, RWORK ) / DBLE( M ) )
 285       ELSE
 286          EPS2 = ULP
 287       END IF
 288       R = MIN( P, M-P, Q, M-Q )
 289 *
 290 *     Copy the matrix X to the array XF.
 291 *
 292       CALL ZLACPY( 'Full', M, M, X, LDX, XF, LDX )
 293 *
 294 *     Compute the CSD
 295 *
 296       CALL ZUNCSD( 'Y', 'Y', 'Y', 'Y', 'N', 'D', M, P, Q, XF(1,1), LDX,
 297      $             XF(1,Q+1), LDX, XF(P+1,1), LDX, XF(P+1,Q+1), LDX,
 298      $             THETA, U1, LDU1, U2, LDU2, V1T, LDV1T, V2T, LDV2T,
 299      $             WORK, LWORK, RWORK, 17*(R+2), IWORK, INFO )
 300 *
 301 *     Compute XF := diag(U1,U2)'*X*diag(V1,V2) - [D11 D12; D21 D22]
 302 *
 303       CALL ZLACPY( 'Full', M, M, X, LDX, XF, LDX )
 304 *
 305       CALL ZGEMM( 'No transpose', 'Conjugate transpose', P, Q, Q, ONE,
 306      $            XF, LDX, V1T, LDV1T, ZERO, WORK, LDX )
 307 *
 308       CALL ZGEMM( 'Conjugate transpose', 'No transpose', P, Q, P, ONE,
 309      $            U1, LDU1, WORK, LDX, ZERO, XF, LDX )
 310 *
 311       DO I = 1, MIN(P,Q)-R
 312          XF(I,I) = XF(I,I) - ONE
 313       END DO
 314       DO I = 1, R
 315          XF(MIN(P,Q)-R+I,MIN(P,Q)-R+I) =
 316      $           XF(MIN(P,Q)-R+I,MIN(P,Q)-R+I) - DCMPLX( COS(THETA(I)),
 317      $              0.0D0 )
 318       END DO
 319 *
 320       CALL ZGEMM( 'No transpose', 'Conjugate transpose', P, M-Q, M-Q,
 321      $            ONE, XF(1,Q+1), LDX, V2T, LDV2T, ZERO, WORK, LDX )
 322 *
 323       CALL ZGEMM( 'Conjugate transpose', 'No transpose', P, M-Q, P,
 324      $            ONE, U1, LDU1, WORK, LDX, ZERO, XF(1,Q+1), LDX )
 325 *
 326       DO I = 1, MIN(P,M-Q)-R
 327          XF(P-I+1,M-I+1) = XF(P-I+1,M-I+1) + ONE
 328       END DO
 329       DO I = 1, R
 330          XF(P-(MIN(P,M-Q)-R)+1-I,M-(MIN(P,M-Q)-R)+1-I) =
 331      $      XF(P-(MIN(P,M-Q)-R)+1-I,M-(MIN(P,M-Q)-R)+1-I) +
 332      $      DCMPLX( SIN(THETA(R-I+1)), 0.0D0 )
 333       END DO
 334 *
 335       CALL ZGEMM( 'No transpose', 'Conjugate transpose', M-P, Q, Q, ONE,
 336      $            XF(P+1,1), LDX, V1T, LDV1T, ZERO, WORK, LDX )
 337 *
 338       CALL ZGEMM( 'Conjugate transpose', 'No transpose', M-P, Q, M-P,
 339      $            ONE, U2, LDU2, WORK, LDX, ZERO, XF(P+1,1), LDX )
 340 *
 341       DO I = 1, MIN(M-P,Q)-R
 342          XF(M-I+1,Q-I+1) = XF(M-I+1,Q-I+1) - ONE
 343       END DO
 344       DO I = 1, R
 345          XF(M-(MIN(M-P,Q)-R)+1-I,Q-(MIN(M-P,Q)-R)+1-I) =
 346      $             XF(M-(MIN(M-P,Q)-R)+1-I,Q-(MIN(M-P,Q)-R)+1-I) -
 347      $             DCMPLX( SIN(THETA(R-I+1)), 0.0D0 )
 348       END DO
 349 *
 350       CALL ZGEMM( 'No transpose', 'Conjugate transpose', M-P, M-Q, M-Q,
 351      $            ONE, XF(P+1,Q+1), LDX, V2T, LDV2T, ZERO, WORK, LDX )
 352 *
 353       CALL ZGEMM( 'Conjugate transpose', 'No transpose', M-P, M-Q, M-P,
 354      $            ONE, U2, LDU2, WORK, LDX, ZERO, XF(P+1,Q+1), LDX )
 355 *
 356       DO I = 1, MIN(M-P,M-Q)-R
 357          XF(P+I,Q+I) = XF(P+I,Q+I) - ONE
 358       END DO
 359       DO I = 1, R
 360          XF(P+(MIN(M-P,M-Q)-R)+I,Q+(MIN(M-P,M-Q)-R)+I) =
 361      $      XF(P+(MIN(M-P,M-Q)-R)+I,Q+(MIN(M-P,M-Q)-R)+I) -
 362      $      DCMPLX( COS(THETA(I)), 0.0D0 )
 363       END DO
 364 *
 365 *     Compute norm( U1'*X11*V1 - D11 ) / ( MAX(1,P,Q)*EPS2 ) .
 366 *
 367       RESID = ZLANGE( '1', P, Q, XF, LDX, RWORK )
 368       RESULT( 1 ) = ( RESID / REAL(MAX(1,P,Q)) ) / EPS2
 369 *
 370 *     Compute norm( U1'*X12*V2 - D12 ) / ( MAX(1,P,M-Q)*EPS2 ) .
 371 *
 372       RESID = ZLANGE( '1', P, M-Q, XF(1,Q+1), LDX, RWORK )
 373       RESULT( 2 ) = ( RESID / REAL(MAX(1,P,M-Q)) ) / EPS2
 374 *
 375 *     Compute norm( U2'*X21*V1 - D21 ) / ( MAX(1,M-P,Q)*EPS2 ) .
 376 *
 377       RESID = ZLANGE( '1', M-P, Q, XF(P+1,1), LDX, RWORK )
 378       RESULT( 3 ) = ( RESID / REAL(MAX(1,M-P,Q)) ) / EPS2
 379 *
 380 *     Compute norm( U2'*X22*V2 - D22 ) / ( MAX(1,M-P,M-Q)*EPS2 ) .
 381 *
 382       RESID = ZLANGE( '1', M-P, M-Q, XF(P+1,Q+1), LDX, RWORK )
 383       RESULT( 4 ) = ( RESID / REAL(MAX(1,M-P,M-Q)) ) / EPS2
 384 *
 385 *     Compute I - U1'*U1
 386 *
 387       CALL ZLASET( 'Full', P, P, ZERO, ONE, WORK, LDU1 )
 388       CALL ZHERK( 'Upper', 'Conjugate transpose', P, P, -REALONE,
 389      $            U1, LDU1, REALONE, WORK, LDU1 )
 390 *
 391 *     Compute norm( I - U'*U ) / ( MAX(1,P) * ULP ) .
 392 *
 393       RESID = ZLANHE( '1', 'Upper', P, WORK, LDU1, RWORK )
 394       RESULT( 5 ) = ( RESID / REAL(MAX(1,P)) ) / ULP
 395 *
 396 *     Compute I - U2'*U2
 397 *
 398       CALL ZLASET( 'Full', M-P, M-P, ZERO, ONE, WORK, LDU2 )
 399       CALL ZHERK( 'Upper', 'Conjugate transpose', M-P, M-P, -REALONE,
 400      $            U2, LDU2, REALONE, WORK, LDU2 )
 401 *
 402 *     Compute norm( I - U2'*U2 ) / ( MAX(1,M-P) * ULP ) .
 403 *
 404       RESID = ZLANHE( '1', 'Upper', M-P, WORK, LDU2, RWORK )
 405       RESULT( 6 ) = ( RESID / REAL(MAX(1,M-P)) ) / ULP
 406 *
 407 *     Compute I - V1T*V1T'
 408 *
 409       CALL ZLASET( 'Full', Q, Q, ZERO, ONE, WORK, LDV1T )
 410       CALL ZHERK( 'Upper', 'No transpose', Q, Q, -REALONE,
 411      $            V1T, LDV1T, REALONE, WORK, LDV1T )
 412 *
 413 *     Compute norm( I - V1T*V1T' ) / ( MAX(1,Q) * ULP ) .
 414 *
 415       RESID = ZLANHE( '1', 'Upper', Q, WORK, LDV1T, RWORK )
 416       RESULT( 7 ) = ( RESID / REAL(MAX(1,Q)) ) / ULP
 417 *
 418 *     Compute I - V2T*V2T'
 419 *
 420       CALL ZLASET( 'Full', M-Q, M-Q, ZERO, ONE, WORK, LDV2T )
 421       CALL ZHERK( 'Upper', 'No transpose', M-Q, M-Q, -REALONE,
 422      $            V2T, LDV2T, REALONE, WORK, LDV2T )
 423 *
 424 *     Compute norm( I - V2T*V2T' ) / ( MAX(1,M-Q) * ULP ) .
 425 *
 426       RESID = ZLANHE( '1', 'Upper', M-Q, WORK, LDV2T, RWORK )
 427       RESULT( 8 ) = ( RESID / REAL(MAX(1,M-Q)) ) / ULP
 428 *
 429 *     Check sorting
 430 *
 431       RESULT( 9 ) = REALZERO
 432       DO I = 1, R
 433          IF( THETA(I).LT.REALZERO .OR. THETA(I).GT.PIOVER2 ) THEN
 434             RESULT( 9 ) = ULPINV
 435          END IF
 436          IF( I.GT.1) THEN
 437             IF ( THETA(I).LT.THETA(I-1) ) THEN
 438                RESULT( 9 ) = ULPINV
 439             END IF
 440          END IF
 441       END DO
 442 *
 443 *     The second half of the routine checks the 2-by-1 CSD
 444 *
 445       CALL ZLASET( 'Full', Q, Q, ZERO, ONE, WORK, LDX )
 446       CALL ZHERK( 'Upper', 'Conjugate transpose', Q, M, -REALONE,
 447      $            X, LDX, REALONE, WORK, LDX )
 448       IF (M.GT.0) THEN
 449          EPS2 = MAX( ULP,
 450      $               ZLANGE( '1', Q, Q, WORK, LDX, RWORK ) / DBLE( M ) )
 451       ELSE
 452          EPS2 = ULP
 453       END IF
 454       R = MIN( P, M-P, Q, M-Q )
 455 *
 456 *     Copy the matrix X to the array XF.
 457 *
 458       CALL ZLACPY( 'Full', M, M, X, LDX, XF, LDX )
 459 *
 460 *     Compute the CSD
 461 *
 462       CALL ZUNCSD2BY1( 'Y', 'Y', 'Y', M, P, Q, XF(1,1), LDX, XF(P+1,1),
 463      $             LDX, THETA, U1, LDU1, U2, LDU2, V1T, LDV1T, WORK,
 464      $             LWORK, RWORK, 17*(R+2), IWORK, INFO )
 465 *
 466 *     Compute [X11;X21] := diag(U1,U2)'*[X11;X21]*V1 - [D11;D21]
 467 *
 468       CALL ZGEMM( 'No transpose', 'Conjugate transpose', P, Q, Q, ONE,
 469      $            X, LDX, V1T, LDV1T, ZERO, WORK, LDX )
 470 *
 471       CALL ZGEMM( 'Conjugate transpose', 'No transpose', P, Q, P, ONE,
 472      $            U1, LDU1, WORK, LDX, ZERO, X, LDX )
 473 *
 474       DO I = 1, MIN(P,Q)-R
 475          X(I,I) = X(I,I) - ONE
 476       END DO
 477       DO I = 1, R
 478          X(MIN(P,Q)-R+I,MIN(P,Q)-R+I) =
 479      $           X(MIN(P,Q)-R+I,MIN(P,Q)-R+I) - DCMPLX( COS(THETA(I)),
 480      $              0.0D0 )
 481       END DO
 482 *
 483       CALL ZGEMM( 'No transpose', 'Conjugate transpose', M-P, Q, Q, ONE,
 484      $            X(P+1,1), LDX, V1T, LDV1T, ZERO, WORK, LDX )
 485 *
 486       CALL ZGEMM( 'Conjugate transpose', 'No transpose', M-P, Q, M-P,
 487      $            ONE, U2, LDU2, WORK, LDX, ZERO, X(P+1,1), LDX )
 488 *
 489       DO I = 1, MIN(M-P,Q)-R
 490          X(M-I+1,Q-I+1) = X(M-I+1,Q-I+1) - ONE
 491       END DO
 492       DO I = 1, R
 493          X(M-(MIN(M-P,Q)-R)+1-I,Q-(MIN(M-P,Q)-R)+1-I) =
 494      $             X(M-(MIN(M-P,Q)-R)+1-I,Q-(MIN(M-P,Q)-R)+1-I) -
 495      $             DCMPLX( SIN(THETA(R-I+1)), 0.0D0 )
 496       END DO
 497 *
 498 *     Compute norm( U1'*X11*V1 - D11 ) / ( MAX(1,P,Q)*EPS2 ) .
 499 *
 500       RESID = ZLANGE( '1', P, Q, X, LDX, RWORK )
 501       RESULT( 10 ) = ( RESID / REAL(MAX(1,P,Q)) ) / EPS2
 502 *
 503 *     Compute norm( U2'*X21*V1 - D21 ) / ( MAX(1,M-P,Q)*EPS2 ) .
 504 *
 505       RESID = ZLANGE( '1', M-P, Q, X(P+1,1), LDX, RWORK )
 506       RESULT( 11 ) = ( RESID / REAL(MAX(1,M-P,Q)) ) / EPS2
 507 *
 508 *     Compute I - U1'*U1
 509 *
 510       CALL ZLASET( 'Full', P, P, ZERO, ONE, WORK, LDU1 )
 511       CALL ZHERK( 'Upper', 'Conjugate transpose', P, P, -REALONE,
 512      $            U1, LDU1, REALONE, WORK, LDU1 )
 513 *
 514 *     Compute norm( I - U'*U ) / ( MAX(1,P) * ULP ) .
 515 *
 516       RESID = ZLANHE( '1', 'Upper', P, WORK, LDU1, RWORK )
 517       RESULT( 12 ) = ( RESID / REAL(MAX(1,P)) ) / ULP
 518 *
 519 *     Compute I - U2'*U2
 520 *
 521       CALL ZLASET( 'Full', M-P, M-P, ZERO, ONE, WORK, LDU2 )
 522       CALL ZHERK( 'Upper', 'Conjugate transpose', M-P, M-P, -REALONE,
 523      $            U2, LDU2, REALONE, WORK, LDU2 )
 524 *
 525 *     Compute norm( I - U2'*U2 ) / ( MAX(1,M-P) * ULP ) .
 526 *
 527       RESID = ZLANHE( '1', 'Upper', M-P, WORK, LDU2, RWORK )
 528       RESULT( 13 ) = ( RESID / REAL(MAX(1,M-P)) ) / ULP
 529 *
 530 *     Compute I - V1T*V1T'
 531 *
 532       CALL ZLASET( 'Full', Q, Q, ZERO, ONE, WORK, LDV1T )
 533       CALL ZHERK( 'Upper', 'No transpose', Q, Q, -REALONE,
 534      $            V1T, LDV1T, REALONE, WORK, LDV1T )
 535 *
 536 *     Compute norm( I - V1T*V1T' ) / ( MAX(1,Q) * ULP ) .
 537 *
 538       RESID = ZLANHE( '1', 'Upper', Q, WORK, LDV1T, RWORK )
 539       RESULT( 14 ) = ( RESID / REAL(MAX(1,Q)) ) / ULP
 540 *
 541 *     Check sorting
 542 *
 543       RESULT( 15 ) = REALZERO
 544       DO I = 1, R
 545          IF( THETA(I).LT.REALZERO .OR. THETA(I).GT.PIOVER2 ) THEN
 546             RESULT( 15 ) = ULPINV
 547          END IF
 548          IF( I.GT.1) THEN
 549             IF ( THETA(I).LT.THETA(I-1) ) THEN
 550                RESULT( 15 ) = ULPINV
 551             END IF
 552          END IF
 553       END DO
 554 *
 555       RETURN
 556 *
 557 *     End of ZCSDTS
 558 *
 559       END
 560