LAPACK: TESTING/EIG/dcsdts.f Source File

Go to the documentation of this file.
00001       SUBROUTINE DCSDTS( 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.1) --
00008 *     Univ. of Tennessee, Univ. of California Berkeley and NAG Ltd..
00009 *     November 2006
00010 *
00011 *     Adapted to DCSDTS
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       DOUBLE PRECISION   U1( LDU1, * ), U2( LDU2, * ), V1T( LDV1T, * ),
00021      $                   V2T( LDV2T, * ), WORK( LWORK ), X( LDX, * ),
00022      $                   XF( LDX, * )
00023 *     ..
00024 *
00025 *  Purpose
00026 *  =======
00027 *
00028 *  DCSDTS tests DORCSD, which, given an M-by-M partitioned orthogonal
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) DOUBLE PRECISION array, dimension (LDX,M)
00060 *          The M-by-M matrix X.
00061 *
00062 *  XF      (output) DOUBLE PRECISION array, dimension (LDX,M)
00063 *          Details of the CSD of X, as returned by DORCSD;
00064 *          see DORCSD 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) DOUBLE PRECISION array, dimension(LDU1,P)
00071 *          The P-by-P orthogonal matrix U1.
00072 *
00073 *  LDU1    (input) INTEGER
00074 *          The leading dimension of the array U1. LDU >= max(1,P).
00075 *
00076 *  U2      (output) DOUBLE PRECISION array, dimension(LDU2,M-P)
00077 *          The (M-P)-by-(M-P) orthogonal matrix U2.
00078 *
00079 *  LDU2    (input) INTEGER
00080 *          The leading dimension of the array U2. LDU >= max(1,M-P).
00081 *
00082 *  V1T     (output) DOUBLE PRECISION array, dimension(LDV1T,Q)
00083 *          The Q-by-Q orthogonal matrix V1T.
00084 *
00085 *  LDV1T   (input) INTEGER
00086 *          The leading dimension of the array V1T. LDV1T >=
00087 *          max(1,Q).
00088 *
00089 *  V2T     (output) DOUBLE PRECISION array, dimension(LDV2T,M-Q)
00090 *          The (M-Q)-by-(M-Q) orthogonal 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 DORCSD for
00099 *          details.
00100 *
00101 *  IWORK   (workspace) INTEGER array, dimension (M)
00102 *
00103 *  WORK    (workspace) DOUBLE PRECISION 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       DOUBLE PRECISION   ZERO, ONE
00132       PARAMETER          ( ZERO = 0.0D0, ONE = 1.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, DLANGE, DLANSY
00140       EXTERNAL           DLAMCH, DLANGE, DLANSY
00141 *     ..
00142 *     .. External Subroutines ..
00143       EXTERNAL           DGEMM, DLACPY, DLASET, DORCSD, DSYRK
00144 *     ..
00145 *     .. Intrinsic Functions ..
00146       INTRINSIC          DBLE, MAX, MIN
00147 *     ..
00148 *     .. Executable Statements ..
00149 *
00150       ULP = DLAMCH( 'Precision' )
00151       ULPINV = REALONE / ULP
00152       CALL DLASET( 'Full', M, M, ZERO, ONE, WORK, LDX )
00153       CALL DSYRK( 'Upper', 'Conjugate transpose', M, M, -ONE, X, LDX,
00154      $            ONE, WORK, LDX )
00155       EPS2 = MAX( ULP, 
00156      $            DLANGE( '1', M, M, WORK, LDX, RWORK ) / DBLE( M ) )
00157       R = MIN( P, M-P, Q, M-Q )
00158 *
00159 *     Copy the matrix X to the array XF.
00160 *
00161       CALL DLACPY( 'Full', M, M, X, LDX, XF, LDX )
00162 *
00163 *     Compute the CSD
00164 *
00165       CALL DORCSD( '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, IWORK, INFO )
00169 *
00170 *     Compute X := diag(U1,U2)'*X*diag(V1,V2) - [D11 D12; D21 D22]
00171 *
00172       CALL DGEMM( 'No transpose', 'Conjugate transpose', P, Q, Q, ONE,
00173      $            X, LDX, V1T, LDV1T, ZERO, WORK, LDX )
00174 *
00175       CALL DGEMM( '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) - COS(THETA(I))
00184       END DO
00185 *
00186       CALL DGEMM( 'No transpose', 'Conjugate transpose', P, M-Q, M-Q,
00187      $            ONE, X(1,Q+1), LDX, V2T, LDV2T, ZERO, WORK, LDX )
00188 *
00189       CALL DGEMM( 'Conjugate transpose', 'No transpose', P, M-Q, P,
00190      $            ONE, U1, LDU1, WORK, LDX, ZERO, X(1,Q+1), LDX )
00191 *
00192       DO I = 1, MIN(P,M-Q)-R
00193          X(P-I+1,M-I+1) = X(P-I+1,M-I+1) + ONE
00194       END DO
00195       DO I = 1, R
00196          X(P-(MIN(P,M-Q)-R)+1-I,M-(MIN(P,M-Q)-R)+1-I) =
00197      $      X(P-(MIN(P,M-Q)-R)+1-I,M-(MIN(P,M-Q)-R)+1-I) +
00198      $      SIN(THETA(R-I+1))
00199       END DO
00200 *
00201       CALL DGEMM( 'No transpose', 'Conjugate transpose', M-P, Q, Q, ONE,
00202      $            X(P+1,1), LDX, V1T, LDV1T, ZERO, WORK, LDX )
00203 *
00204       CALL DGEMM( 'Conjugate transpose', 'No transpose', M-P, Q, M-P,
00205      $            ONE, U2, LDU2, WORK, LDX, ZERO, X(P+1,1), LDX )
00206 *
00207       DO I = 1, MIN(M-P,Q)-R
00208          X(M-I+1,Q-I+1) = X(M-I+1,Q-I+1) - ONE
00209       END DO
00210       DO I = 1, R
00211          X(M-(MIN(M-P,Q)-R)+1-I,Q-(MIN(M-P,Q)-R)+1-I) =
00212      $             X(M-(MIN(M-P,Q)-R)+1-I,Q-(MIN(M-P,Q)-R)+1-I) -
00213      $             SIN(THETA(R-I+1))
00214       END DO
00215 *
00216       CALL DGEMM( 'No transpose', 'Conjugate transpose', M-P, M-Q, M-Q,
00217      $            ONE, X(P+1,Q+1), LDX, V2T, LDV2T, ZERO, WORK, LDX )
00218 *
00219       CALL DGEMM( 'Conjugate transpose', 'No transpose', M-P, M-Q, M-P,
00220      $            ONE, U2, LDU2, WORK, LDX, ZERO, X(P+1,Q+1), LDX )
00221 *
00222       DO I = 1, MIN(M-P,M-Q)-R
00223          X(P+I,Q+I) = X(P+I,Q+I) - ONE
00224       END DO
00225       DO I = 1, R
00226          X(P+(MIN(M-P,M-Q)-R)+I,Q+(MIN(M-P,M-Q)-R)+I) =
00227      $      X(P+(MIN(M-P,M-Q)-R)+I,Q+(MIN(M-P,M-Q)-R)+I) -
00228      $      COS(THETA(I))
00229       END DO
00230 *
00231 *     Compute norm( U1'*X11*V1 - D11 ) / ( MAX(1,P,Q)*EPS2 ) .
00232 *
00233       RESID = DLANGE( '1', P, Q, X, LDX, RWORK )
00234       RESULT( 1 ) = ( RESID / DBLE(MAX(1,P,Q)) ) / EPS2
00235 *
00236 *     Compute norm( U1'*X12*V2 - D12 ) / ( MAX(1,P,M-Q)*EPS2 ) .
00237 *
00238       RESID = DLANGE( '1', P, M-Q, X(1,Q+1), LDX, RWORK )
00239       RESULT( 2 ) = ( RESID / DBLE(MAX(1,P,M-Q)) ) / EPS2
00240 *
00241 *     Compute norm( U2'*X21*V1 - D21 ) / ( MAX(1,M-P,Q)*EPS2 ) .
00242 *
00243       RESID = DLANGE( '1', M-P, Q, X(P+1,1), LDX, RWORK )
00244       RESULT( 3 ) = ( RESID / DBLE(MAX(1,M-P,Q)) ) / EPS2
00245 *
00246 *     Compute norm( U2'*X22*V2 - D22 ) / ( MAX(1,M-P,M-Q)*EPS2 ) .
00247 *
00248       RESID = DLANGE( '1', M-P, M-Q, X(P+1,Q+1), LDX, RWORK )
00249       RESULT( 4 ) = ( RESID / DBLE(MAX(1,M-P,M-Q)) ) / EPS2
00250 *
00251 *     Compute I - U1'*U1
00252 *
00253       CALL DLASET( 'Full', P, P, ZERO, ONE, WORK, LDU1 )
00254       CALL DSYRK( 'Upper', 'Conjugate transpose', P, P, -ONE, U1, LDU1,
00255      $            ONE, WORK, LDU1 )
00256 *
00257 *     Compute norm( I - U'*U ) / ( MAX(1,P) * ULP ) .
00258 *
00259       RESID = DLANSY( '1', 'Upper', P, WORK, LDU1, RWORK )
00260       RESULT( 5 ) = ( RESID / DBLE(MAX(1,P)) ) / ULP
00261 *
00262 *     Compute I - U2'*U2
00263 *
00264       CALL DLASET( 'Full', M-P, M-P, ZERO, ONE, WORK, LDU2 )
00265       CALL DSYRK( 'Upper', 'Conjugate transpose', M-P, M-P, -ONE, U2,
00266      $            LDU2, ONE, WORK, LDU2 )
00267 *
00268 *     Compute norm( I - U2'*U2 ) / ( MAX(1,M-P) * ULP ) .
00269 *
00270       RESID = DLANSY( '1', 'Upper', M-P, WORK, LDU2, RWORK )
00271       RESULT( 6 ) = ( RESID / DBLE(MAX(1,M-P)) ) / ULP
00272 *
00273 *     Compute I - V1T*V1T'
00274 *
00275       CALL DLASET( 'Full', Q, Q, ZERO, ONE, WORK, LDV1T )
00276       CALL DSYRK( 'Upper', 'No transpose', Q, Q, -ONE, V1T, LDV1T, ONE,
00277      $            WORK, LDV1T )
00278 *
00279 *     Compute norm( I - V1T*V1T' ) / ( MAX(1,Q) * ULP ) .
00280 *
00281       RESID = DLANSY( '1', 'Upper', Q, WORK, LDV1T, RWORK )
00282       RESULT( 7 ) = ( RESID / DBLE(MAX(1,Q)) ) / ULP
00283 *
00284 *     Compute I - V2T*V2T'
00285 *
00286       CALL DLASET( 'Full', M-Q, M-Q, ZERO, ONE, WORK, LDV2T )
00287       CALL DSYRK( 'Upper', 'No transpose', M-Q, M-Q, -ONE, V2T, LDV2T,
00288      $            ONE, WORK, LDV2T )
00289 *
00290 *     Compute norm( I - V2T*V2T' ) / ( MAX(1,M-Q) * ULP ) .
00291 *
00292       RESID = DLANSY( '1', 'Upper', M-Q, WORK, LDV2T, RWORK )
00293       RESULT( 8 ) = ( RESID / DBLE(MAX(1,M-Q)) ) / ULP
00294 *
00295 *     Check sorting
00296 *
00297       RESULT(9) = REALZERO
00298       DO I = 1, R
00299          IF( THETA(I).LT.REALZERO .OR. THETA(I).GT.PIOVER2 ) THEN
00300             RESULT(9) = ULPINV
00301          END IF
00302          IF( I.GT.1) THEN
00303             IF ( THETA(I).LT.THETA(I-1) ) THEN
00304                RESULT(9) = ULPINV
00305             END IF
00306          END IF
00307       END DO
00308 *
00309       RETURN
00310 *      
00311 *     End of DCSDTS
00312 *
00313       END
00314