LAPACK: TESTING/EIG/zcsdts.f Source File

Go to the documentation of this file.
00001       SUBROUTINE ZCSDTS( M, P, Q, X, XF, LDX, U1, LDU1, U2, LDU2, V1T,
00002      $                   LDV1T, V2T, LDV2T, THETA, IWORK, WORK, LWORK,
00003      $                   RWORK, RESULT )
00004       IMPLICIT NONE
00005 *
00006 *     Originally xGSVTS
00007 *  -- LAPACK test routine (version 3.3.0) --
00008 *     Univ. of Tennessee, Univ. of California Berkeley and NAG Ltd..
00009 *     November 2010
00010 *
00011 *     Adapted to ZCSDTS
00012 *     July 2010
00013 *
00014 *     .. Scalar Arguments ..
00015       INTEGER            LDX, LDU1, LDU2, LDV1T, LDV2T, LWORK, M, P, Q
00016 *     ..
00017 *     .. Array Arguments ..
00018       INTEGER            IWORK( * )
00019       DOUBLE PRECISION   RESULT( 9 ), RWORK( * ), THETA( * )
00020       COMPLEX*16         U1( LDU1, * ), U2( LDU2, * ), V1T( LDV1T, * ),
00021      $                   V2T( LDV2T, * ), WORK( LWORK ), X( LDX, * ),
00022      $                   XF( LDX, * )
00023 *     ..
00024 *
00025 *  Purpose
00026 *  =======
00027 *
00028 *  ZCSDTS tests ZUNCSD, which, given an M-by-M partitioned unitary
00029 *  matrix X,
00030 *               Q  M-Q
00031 *        X = [ X11 X12 ] P   ,
00032 *            [ X21 X22 ] M-P
00033 *
00034 *  computes the CSD
00035 *
00036 *        [ U1    ]**T * [ X11 X12 ] * [ V1    ]
00037 *        [    U2 ]      [ X21 X22 ]   [    V2 ]
00038 *
00039 *                              [  I  0  0 |  0  0  0 ]
00040 *                              [  0  C  0 |  0 -S  0 ]
00041 *                              [  0  0  0 |  0  0 -I ]
00042 *                            = [---------------------] = [ D11 D12 ] .
00043 *                              [  0  0  0 |  I  0  0 ]   [ D21 D22 ]
00044 *                              [  0  S  0 |  0  C  0 ]
00045 *                              [  0  0  I |  0  0  0 ]
00046 *
00047 *  Arguments
00048 *  =========
00049 *
00050 *  M       (input) INTEGER
00051 *          The number of rows of the matrix X.  M >= 0.
00052 *
00053 *  P       (input) INTEGER
00054 *          The number of rows of the matrix X11.  P >= 0.
00055 *
00056 *  Q       (input) INTEGER
00057 *          The number of columns of the matrix X11.  Q >= 0.
00058 *
00059 *  X       (input) COMPLEX*16 array, dimension (LDX,M)
00060 *          The M-by-M matrix X.
00061 *
00062 *  XF      (output) COMPLEX*16 array, dimension (LDX,M)
00063 *          Details of the CSD of X, as returned by ZUNCSD;
00064 *          see ZUNCSD for further details.
00065 *
00066 *  LDX     (input) INTEGER
00067 *          The leading dimension of the arrays X and XF.
00068 *          LDX >= max( 1,M ).
00069 *
00070 *  U1      (output) COMPLEX*16 array, dimension(LDU1,P)
00071 *          The P-by-P unitary matrix U1.
00072 *
00073 *  LDU1    (input) INTEGER
00074 *          The leading dimension of the array U1. LDU >= max(1,P).
00075 *
00076 *  U2      (output) COMPLEX*16 array, dimension(LDU2,M-P)
00077 *          The (M-P)-by-(M-P) unitary matrix U2.
00078 *
00079 *  LDU2    (input) INTEGER
00080 *          The leading dimension of the array U2. LDU >= max(1,M-P).
00081 *
00082 *  V1T     (output) COMPLEX*16 array, dimension(LDV1T,Q)
00083 *          The Q-by-Q unitary matrix V1T.
00084 *
00085 *  LDV1T   (input) INTEGER
00086 *          The leading dimension of the array V1T. LDV1T >=
00087 *          max(1,Q).
00088 *
00089 *  V2T     (output) COMPLEX*16 array, dimension(LDV2T,M-Q)
00090 *          The (M-Q)-by-(M-Q) unitary matrix V2T.
00091 *
00092 *  LDV2T   (input) INTEGER
00093 *          The leading dimension of the array V2T. LDV2T >=
00094 *          max(1,M-Q).
00095 *
00096 *  THETA   (output) DOUBLE PRECISION array, dimension MIN(P,M-P,Q,M-Q)
00097 *          The CS values of X; the essentially diagonal matrices C and
00098 *          S are constructed from THETA; see subroutine ZUNCSD for
00099 *          details.
00100 *
00101 *  IWORK   (workspace) INTEGER array, dimension (M)
00102 *
00103 *  WORK    (workspace) COMPLEX*16 array, dimension (LWORK)
00104 *
00105 *  LWORK   (input) INTEGER
00106 *          The dimension of the array WORK
00107 *
00108 *  RWORK   (workspace) DOUBLE PRECISION array
00109 *
00110 *  RESULT  (output) DOUBLE PRECISION array, dimension (9)
00111 *          The test ratios:
00112 *          RESULT(1) = norm( U1'*X11*V1 - D11 ) / ( MAX(1,P,Q)*EPS2 )
00113 *          RESULT(2) = norm( U1'*X12*V2 - D12 ) / ( MAX(1,P,M-Q)*EPS2 )
00114 *          RESULT(3) = norm( U2'*X21*V1 - D21 ) / ( MAX(1,M-P,Q)*EPS2 )
00115 *          RESULT(4) = norm( U2'*X22*V2 - D22 ) / ( MAX(1,M-P,M-Q)*EPS2 )
00116 *          RESULT(5) = norm( I - U1'*U1 ) / ( MAX(1,P)*ULP )
00117 *          RESULT(6) = norm( I - U2'*U2 ) / ( MAX(1,M-P)*ULP )
00118 *          RESULT(7) = norm( I - V1T'*V1T ) / ( MAX(1,Q)*ULP )
00119 *          RESULT(8) = norm( I - V2T'*V2T ) / ( MAX(1,M-Q)*ULP )
00120 *          RESULT(9) = 0        if THETA is in increasing order and
00121 *                               all angles are in [0,pi/2];
00122 *                    = ULPINV   otherwise.
00123 *          ( EPS2 = MAX( norm( I - X'*X ) / M, ULP ). )
00124 *
00125 *  =====================================================================
00126 *
00127 *     .. Parameters ..
00128       DOUBLE PRECISION   PIOVER2, REALONE, REALZERO
00129       PARAMETER          ( PIOVER2 = 1.57079632679489662D0,
00130      $                     REALONE = 1.0D0, REALZERO = 0.0D0 )
00131       COMPLEX*16         ZERO, ONE
00132       PARAMETER          ( ZERO = (0.0D0,0.0D0), ONE = (1.0D0,0.0D0) )
00133 *     ..
00134 *     .. Local Scalars ..
00135       INTEGER            I, INFO, R
00136       DOUBLE PRECISION   EPS2, RESID, ULP, ULPINV
00137 *     ..
00138 *     .. External Functions ..
00139       DOUBLE PRECISION   DLAMCH, ZLANGE, ZLANHE
00140       EXTERNAL           DLAMCH, ZLANGE, ZLANHE
00141 *     ..
00142 *     .. External Subroutines ..
00143       EXTERNAL           ZGEMM, ZLACPY, ZLASET, ZUNCSD, ZHERK
00144 *     ..
00145 *     .. Intrinsic Functions ..
00146       INTRINSIC          REAL, MAX, MIN
00147 *     ..
00148 *     .. Executable Statements ..
00149 *
00150       ULP = DLAMCH( 'Precision' )
00151       ULPINV = REALONE / ULP
00152       CALL ZLASET( 'Full', M, M, ZERO, ONE, WORK, LDX )
00153       CALL ZHERK( 'Upper', 'Conjugate transpose', M, M, -ONE, X, LDX,
00154      $            ONE, WORK, LDX )
00155       EPS2 = MAX( ULP, 
00156      $            ZLANGE( '1', M, M, WORK, LDX, RWORK ) / REAL( M ) )
00157       R = MIN( P, M-P, Q, M-Q )
00158 *
00159 *     Copy the matrix X to the array XF.
00160 *
00161       CALL ZLACPY( 'Full', M, M, X, LDX, XF, LDX )
00162 *
00163 *     Compute the CSD
00164 *
00165       CALL ZUNCSD( 'Y', 'Y', 'Y', 'Y', 'N', 'D', M, P, Q, XF(1,1), LDX,
00166      $             XF(1,Q+1), LDX, XF(P+1,1), LDX, XF(P+1,Q+1), LDX,
00167      $             THETA, U1, LDU1, U2, LDU2, V1T, LDV1T, V2T, LDV2T,
00168      $             WORK, LWORK, RWORK, 17*(R+2), IWORK, INFO )
00169 *
00170 *     Compute X := diag(U1,U2)'*X*diag(V1,V2) - [D11 D12; D21 D22]
00171 *
00172       CALL ZGEMM( 'No transpose', 'Conjugate transpose', P, Q, Q, ONE,
00173      $            X, LDX, V1T, LDV1T, ZERO, WORK, LDX )
00174 *
00175       CALL ZGEMM( 'Conjugate transpose', 'No transpose', P, Q, P, ONE,
00176      $            U1, LDU1, WORK, LDX, ZERO, X, LDX )
00177 *
00178       DO I = 1, MIN(P,Q)-R
00179          X(I,I) = X(I,I) - ONE
00180       END DO
00181       DO I = 1, R
00182          X(MIN(P,Q)-R+I,MIN(P,Q)-R+I) =
00183      $           X(MIN(P,Q)-R+I,MIN(P,Q)-R+I) - DCMPLX( COS(THETA(I)),
00184      $              0.0D0 )
00185       END DO
00186 *
00187       CALL ZGEMM( 'No transpose', 'Conjugate transpose', P, M-Q, M-Q,
00188      $            ONE, X(1,Q+1), LDX, V2T, LDV2T, ZERO, WORK, LDX )
00189 *
00190       CALL ZGEMM( 'Conjugate transpose', 'No transpose', P, M-Q, P,
00191      $            ONE, U1, LDU1, WORK, LDX, ZERO, X(1,Q+1), LDX )
00192 *
00193       DO I = 1, MIN(P,M-Q)-R
00194          X(P-I+1,M-I+1) = X(P-I+1,M-I+1) + ONE
00195       END DO
00196       DO I = 1, R
00197          X(P-(MIN(P,M-Q)-R)+1-I,M-(MIN(P,M-Q)-R)+1-I) =
00198      $      X(P-(MIN(P,M-Q)-R)+1-I,M-(MIN(P,M-Q)-R)+1-I) +
00199      $      DCMPLX( SIN(THETA(R-I+1)), 0.0D0 )
00200       END DO
00201 *
00202       CALL ZGEMM( 'No transpose', 'Conjugate transpose', M-P, Q, Q, ONE,
00203      $            X(P+1,1), LDX, V1T, LDV1T, ZERO, WORK, LDX )
00204 *
00205       CALL ZGEMM( 'Conjugate transpose', 'No transpose', M-P, Q, M-P,
00206      $            ONE, U2, LDU2, WORK, LDX, ZERO, X(P+1,1), LDX )
00207 *
00208       DO I = 1, MIN(M-P,Q)-R
00209          X(M-I+1,Q-I+1) = X(M-I+1,Q-I+1) - ONE
00210       END DO
00211       DO I = 1, R
00212          X(M-(MIN(M-P,Q)-R)+1-I,Q-(MIN(M-P,Q)-R)+1-I) =
00213      $             X(M-(MIN(M-P,Q)-R)+1-I,Q-(MIN(M-P,Q)-R)+1-I) -
00214      $             DCMPLX( SIN(THETA(R-I+1)), 0.0D0 )
00215       END DO
00216 *
00217       CALL ZGEMM( 'No transpose', 'Conjugate transpose', M-P, M-Q, M-Q,
00218      $            ONE, X(P+1,Q+1), LDX, V2T, LDV2T, ZERO, WORK, LDX )
00219 *
00220       CALL ZGEMM( 'Conjugate transpose', 'No transpose', M-P, M-Q, M-P,
00221      $            ONE, U2, LDU2, WORK, LDX, ZERO, X(P+1,Q+1), LDX )
00222 *
00223       DO I = 1, MIN(M-P,M-Q)-R
00224          X(P+I,Q+I) = X(P+I,Q+I) - ONE
00225       END DO
00226       DO I = 1, R
00227          X(P+(MIN(M-P,M-Q)-R)+I,Q+(MIN(M-P,M-Q)-R)+I) =
00228      $      X(P+(MIN(M-P,M-Q)-R)+I,Q+(MIN(M-P,M-Q)-R)+I) -
00229      $      DCMPLX( COS(THETA(I)), 0.0D0 )
00230       END DO
00231 *
00232 *     Compute norm( U1'*X11*V1 - D11 ) / ( MAX(1,P,Q)*EPS2 ) .
00233 *
00234       RESID = ZLANGE( '1', P, Q, X, LDX, RWORK )
00235       RESULT( 1 ) = ( RESID / REAL(MAX(1,P,Q)) ) / EPS2
00236 *
00237 *     Compute norm( U1'*X12*V2 - D12 ) / ( MAX(1,P,M-Q)*EPS2 ) .
00238 *
00239       RESID = ZLANGE( '1', P, M-Q, X(1,Q+1), LDX, RWORK )
00240       RESULT( 2 ) = ( RESID / REAL(MAX(1,P,M-Q)) ) / EPS2
00241 *
00242 *     Compute norm( U2'*X21*V1 - D21 ) / ( MAX(1,M-P,Q)*EPS2 ) .
00243 *
00244       RESID = ZLANGE( '1', M-P, Q, X(P+1,1), LDX, RWORK )
00245       RESULT( 3 ) = ( RESID / REAL(MAX(1,M-P,Q)) ) / EPS2
00246 *
00247 *     Compute norm( U2'*X22*V2 - D22 ) / ( MAX(1,M-P,M-Q)*EPS2 ) .
00248 *
00249       RESID = ZLANGE( '1', M-P, M-Q, X(P+1,Q+1), LDX, RWORK )
00250       RESULT( 4 ) = ( RESID / REAL(MAX(1,M-P,M-Q)) ) / EPS2
00251 *
00252 *     Compute I - U1'*U1
00253 *
00254       CALL ZLASET( 'Full', P, P, ZERO, ONE, WORK, LDU1 )
00255       CALL ZHERK( 'Upper', 'Conjugate transpose', P, P, -ONE, U1, LDU1,
00256      $            ONE, WORK, LDU1 )
00257 *
00258 *     Compute norm( I - U'*U ) / ( MAX(1,P) * ULP ) .
00259 *
00260       RESID = ZLANHE( '1', 'Upper', P, WORK, LDU1, RWORK )
00261       RESULT( 5 ) = ( RESID / REAL(MAX(1,P)) ) / ULP
00262 *
00263 *     Compute I - U2'*U2
00264 *
00265       CALL ZLASET( 'Full', M-P, M-P, ZERO, ONE, WORK, LDU2 )
00266       CALL ZHERK( 'Upper', 'Conjugate transpose', M-P, M-P, -ONE, U2,
00267      $            LDU2, ONE, WORK, LDU2 )
00268 *
00269 *     Compute norm( I - U2'*U2 ) / ( MAX(1,M-P) * ULP ) .
00270 *
00271       RESID = ZLANHE( '1', 'Upper', M-P, WORK, LDU2, RWORK )
00272       RESULT( 6 ) = ( RESID / REAL(MAX(1,M-P)) ) / ULP
00273 *
00274 *     Compute I - V1T*V1T'
00275 *
00276       CALL ZLASET( 'Full', Q, Q, ZERO, ONE, WORK, LDV1T )
00277       CALL ZHERK( 'Upper', 'No transpose', Q, Q, -ONE, V1T, LDV1T, ONE,
00278      $            WORK, LDV1T )
00279 *
00280 *     Compute norm( I - V1T*V1T' ) / ( MAX(1,Q) * ULP ) .
00281 *
00282       RESID = ZLANHE( '1', 'Upper', Q, WORK, LDV1T, RWORK )
00283       RESULT( 7 ) = ( RESID / REAL(MAX(1,Q)) ) / ULP
00284 *
00285 *     Compute I - V2T*V2T'
00286 *
00287       CALL ZLASET( 'Full', M-Q, M-Q, ZERO, ONE, WORK, LDV2T )
00288       CALL ZHERK( 'Upper', 'No transpose', M-Q, M-Q, -ONE, V2T, LDV2T,
00289      $            ONE, WORK, LDV2T )
00290 *
00291 *     Compute norm( I - V2T*V2T' ) / ( MAX(1,M-Q) * ULP ) .
00292 *
00293       RESID = ZLANHE( '1', 'Upper', M-Q, WORK, LDV2T, RWORK )
00294       RESULT( 8 ) = ( RESID / REAL(MAX(1,M-Q)) ) / ULP
00295 *
00296 *     Check sorting
00297 *
00298       RESULT(9) = REALZERO
00299       DO I = 1, R
00300          IF( THETA(I).LT.REALZERO .OR. THETA(I).GT.PIOVER2 ) THEN
00301             RESULT(9) = ULPINV
00302          END IF
00303          IF( I.GT.1) THEN
00304             IF ( THETA(I).LT.THETA(I-1) ) THEN
00305                RESULT(9) = ULPINV
00306             END IF
00307          END IF
00308       END DO
00309 *
00310       RETURN
00311 *      
00312 *     End of ZCSDTS
00313 *
00314       END
00315