LAPACK  3.4.1
LAPACK: Linear Algebra PACKage
dcsdts.f
Go to the documentation of this file.
00001 *> \brief \b DCSDTS
00002 *
00003 *  =========== DOCUMENTATION ===========
00004 *
00005 * Online html documentation available at 
00006 *            http://www.netlib.org/lapack/explore-html/ 
00007 *
00008 *  Definition:
00009 *  ===========
00010 *
00011 *       SUBROUTINE DCSDTS( M, P, Q, X, XF, LDX, U1, LDU1, U2, LDU2, V1T,
00012 *                          LDV1T, V2T, LDV2T, THETA, IWORK, WORK, LWORK,
00013 *                          RWORK, RESULT )
00014 * 
00015 *       .. Scalar Arguments ..
00016 *       INTEGER            LDX, LDU1, LDU2, LDV1T, LDV2T, LWORK, M, P, Q
00017 *       ..
00018 *       .. Array Arguments ..
00019 *       INTEGER            IWORK( * )
00020 *       DOUBLE PRECISION   RESULT( 9 ), RWORK( * ), THETA( * )
00021 *       DOUBLE PRECISION   U1( LDU1, * ), U2( LDU2, * ), V1T( LDV1T, * ),
00022 *      $                   V2T( LDV2T, * ), WORK( LWORK ), X( LDX, * ),
00023 *      $                   XF( LDX, * )
00024 *       ..
00025 *  
00026 *
00027 *> \par Purpose:
00028 *  =============
00029 *>
00030 *> \verbatim
00031 *>
00032 *> DCSDTS tests DORCSD, which, given an M-by-M partitioned orthogonal
00033 *> matrix X,
00034 *>              Q  M-Q
00035 *>       X = [ X11 X12 ] P   ,
00036 *>           [ X21 X22 ] M-P
00037 *>
00038 *> computes the CSD
00039 *>
00040 *>       [ U1    ]**T * [ X11 X12 ] * [ V1    ]
00041 *>       [    U2 ]      [ X21 X22 ]   [    V2 ]
00042 *>
00043 *>                             [  I  0  0 |  0  0  0 ]
00044 *>                             [  0  C  0 |  0 -S  0 ]
00045 *>                             [  0  0  0 |  0  0 -I ]
00046 *>                           = [---------------------] = [ D11 D12 ] .
00047 *>                             [  0  0  0 |  I  0  0 ]   [ D21 D22 ]
00048 *>                             [  0  S  0 |  0  C  0 ]
00049 *>                             [  0  0  I |  0  0  0 ]
00050 *> \endverbatim
00051 *
00052 *  Arguments:
00053 *  ==========
00054 *
00055 *> \param[in] M
00056 *> \verbatim
00057 *>          M is INTEGER
00058 *>          The number of rows of the matrix X.  M >= 0.
00059 *> \endverbatim
00060 *>
00061 *> \param[in] P
00062 *> \verbatim
00063 *>          P is INTEGER
00064 *>          The number of rows of the matrix X11.  P >= 0.
00065 *> \endverbatim
00066 *>
00067 *> \param[in] Q
00068 *> \verbatim
00069 *>          Q is INTEGER
00070 *>          The number of columns of the matrix X11.  Q >= 0.
00071 *> \endverbatim
00072 *>
00073 *> \param[in] X
00074 *> \verbatim
00075 *>          X is DOUBLE PRECISION array, dimension (LDX,M)
00076 *>          The M-by-M matrix X.
00077 *> \endverbatim
00078 *>
00079 *> \param[out] XF
00080 *> \verbatim
00081 *>          XF is DOUBLE PRECISION array, dimension (LDX,M)
00082 *>          Details of the CSD of X, as returned by DORCSD;
00083 *>          see DORCSD for further details.
00084 *> \endverbatim
00085 *>
00086 *> \param[in] LDX
00087 *> \verbatim
00088 *>          LDX is INTEGER
00089 *>          The leading dimension of the arrays X and XF.
00090 *>          LDX >= max( 1,M ).
00091 *> \endverbatim
00092 *>
00093 *> \param[out] U1
00094 *> \verbatim
00095 *>          U1 is DOUBLE PRECISION array, dimension(LDU1,P)
00096 *>          The P-by-P orthogonal matrix U1.
00097 *> \endverbatim
00098 *>
00099 *> \param[in] LDU1
00100 *> \verbatim
00101 *>          LDU1 is INTEGER
00102 *>          The leading dimension of the array U1. LDU >= max(1,P).
00103 *> \endverbatim
00104 *>
00105 *> \param[out] U2
00106 *> \verbatim
00107 *>          U2 is DOUBLE PRECISION array, dimension(LDU2,M-P)
00108 *>          The (M-P)-by-(M-P) orthogonal matrix U2.
00109 *> \endverbatim
00110 *>
00111 *> \param[in] LDU2
00112 *> \verbatim
00113 *>          LDU2 is INTEGER
00114 *>          The leading dimension of the array U2. LDU >= max(1,M-P).
00115 *> \endverbatim
00116 *>
00117 *> \param[out] V1T
00118 *> \verbatim
00119 *>          V1T is DOUBLE PRECISION array, dimension(LDV1T,Q)
00120 *>          The Q-by-Q orthogonal matrix V1T.
00121 *> \endverbatim
00122 *>
00123 *> \param[in] LDV1T
00124 *> \verbatim
00125 *>          LDV1T is INTEGER
00126 *>          The leading dimension of the array V1T. LDV1T >=
00127 *>          max(1,Q).
00128 *> \endverbatim
00129 *>
00130 *> \param[out] V2T
00131 *> \verbatim
00132 *>          V2T is DOUBLE PRECISION array, dimension(LDV2T,M-Q)
00133 *>          The (M-Q)-by-(M-Q) orthogonal matrix V2T.
00134 *> \endverbatim
00135 *>
00136 *> \param[in] LDV2T
00137 *> \verbatim
00138 *>          LDV2T is INTEGER
00139 *>          The leading dimension of the array V2T. LDV2T >=
00140 *>          max(1,M-Q).
00141 *> \endverbatim
00142 *>
00143 *> \param[out] THETA
00144 *> \verbatim
00145 *>          THETA is DOUBLE PRECISION array, dimension MIN(P,M-P,Q,M-Q)
00146 *>          The CS values of X; the essentially diagonal matrices C and
00147 *>          S are constructed from THETA; see subroutine DORCSD for
00148 *>          details.
00149 *> \endverbatim
00150 *>
00151 *> \param[out] IWORK
00152 *> \verbatim
00153 *>          IWORK is INTEGER array, dimension (M)
00154 *> \endverbatim
00155 *>
00156 *> \param[out] WORK
00157 *> \verbatim
00158 *>          WORK is DOUBLE PRECISION array, dimension (LWORK)
00159 *> \endverbatim
00160 *>
00161 *> \param[in] LWORK
00162 *> \verbatim
00163 *>          LWORK is INTEGER
00164 *>          The dimension of the array WORK
00165 *> \endverbatim
00166 *>
00167 *> \param[out] RWORK
00168 *> \verbatim
00169 *>          RWORK is DOUBLE PRECISION array
00170 *> \endverbatim
00171 *>
00172 *> \param[out] RESULT
00173 *> \verbatim
00174 *>          RESULT is DOUBLE PRECISION array, dimension (9)
00175 *>          The test ratios:
00176 *>          RESULT(1) = norm( U1'*X11*V1 - D11 ) / ( MAX(1,P,Q)*EPS2 )
00177 *>          RESULT(2) = norm( U1'*X12*V2 - D12 ) / ( MAX(1,P,M-Q)*EPS2 )
00178 *>          RESULT(3) = norm( U2'*X21*V1 - D21 ) / ( MAX(1,M-P,Q)*EPS2 )
00179 *>          RESULT(4) = norm( U2'*X22*V2 - D22 ) / ( MAX(1,M-P,M-Q)*EPS2 )
00180 *>          RESULT(5) = norm( I - U1'*U1 ) / ( MAX(1,P)*ULP )
00181 *>          RESULT(6) = norm( I - U2'*U2 ) / ( MAX(1,M-P)*ULP )
00182 *>          RESULT(7) = norm( I - V1T'*V1T ) / ( MAX(1,Q)*ULP )
00183 *>          RESULT(8) = norm( I - V2T'*V2T ) / ( MAX(1,M-Q)*ULP )
00184 *>          RESULT(9) = 0        if THETA is in increasing order and
00185 *>                               all angles are in [0,pi/2];
00186 *>                    = ULPINV   otherwise.
00187 *>          ( EPS2 = MAX( norm( I - X'*X ) / M, ULP ). )
00188 *> \endverbatim
00189 *
00190 *  Authors:
00191 *  ========
00192 *
00193 *> \author Univ. of Tennessee 
00194 *> \author Univ. of California Berkeley 
00195 *> \author Univ. of Colorado Denver 
00196 *> \author NAG Ltd. 
00197 *
00198 *> \date November 2011
00199 *
00200 *> \ingroup double_eig
00201 *
00202 *  =====================================================================
00203       SUBROUTINE DCSDTS( M, P, Q, X, XF, LDX, U1, LDU1, U2, LDU2, V1T,
00204      $                   LDV1T, V2T, LDV2T, THETA, IWORK, WORK, LWORK,
00205      $                   RWORK, RESULT )
00206 *
00207 *  -- LAPACK test routine (version 3.4.0) --
00208 *  -- LAPACK is a software package provided by Univ. of Tennessee,    --
00209 *  -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
00210 *     November 2011
00211 *
00212 *     .. Scalar Arguments ..
00213       INTEGER            LDX, LDU1, LDU2, LDV1T, LDV2T, LWORK, M, P, Q
00214 *     ..
00215 *     .. Array Arguments ..
00216       INTEGER            IWORK( * )
00217       DOUBLE PRECISION   RESULT( 9 ), RWORK( * ), THETA( * )
00218       DOUBLE PRECISION   U1( LDU1, * ), U2( LDU2, * ), V1T( LDV1T, * ),
00219      $                   V2T( LDV2T, * ), WORK( LWORK ), X( LDX, * ),
00220      $                   XF( LDX, * )
00221 *     ..
00222 *
00223 *  =====================================================================
00224 *
00225 *     .. Parameters ..
00226       DOUBLE PRECISION   PIOVER2, REALONE, REALZERO
00227       PARAMETER          ( PIOVER2 = 1.57079632679489662D0,
00228      $                     REALONE = 1.0D0, REALZERO = 0.0D0 )
00229       DOUBLE PRECISION   ZERO, ONE
00230       PARAMETER          ( ZERO = 0.0D0, ONE = 1.0D0 )
00231 *     ..
00232 *     .. Local Scalars ..
00233       INTEGER            I, INFO, R
00234       DOUBLE PRECISION   EPS2, RESID, ULP, ULPINV
00235 *     ..
00236 *     .. External Functions ..
00237       DOUBLE PRECISION   DLAMCH, DLANGE, DLANSY
00238       EXTERNAL           DLAMCH, DLANGE, DLANSY
00239 *     ..
00240 *     .. External Subroutines ..
00241       EXTERNAL           DGEMM, DLACPY, DLASET, DORCSD, DSYRK
00242 *     ..
00243 *     .. Intrinsic Functions ..
00244       INTRINSIC          DBLE, MAX, MIN
00245 *     ..
00246 *     .. Executable Statements ..
00247 *
00248       ULP = DLAMCH( 'Precision' )
00249       ULPINV = REALONE / ULP
00250       CALL DLASET( 'Full', M, M, ZERO, ONE, WORK, LDX )
00251       CALL DSYRK( 'Upper', 'Conjugate transpose', M, M, -ONE, X, LDX,
00252      $            ONE, WORK, LDX )
00253       IF (M.GT.0) THEN
00254          EPS2 = MAX( ULP, 
00255      $               DLANGE( '1', M, M, WORK, LDX, RWORK ) / DBLE( M ) )
00256       ELSE
00257          EPS2 = ULP
00258       END IF
00259       R = MIN( P, M-P, Q, M-Q )
00260 *
00261 *     Copy the matrix X to the array XF.
00262 *
00263       CALL DLACPY( 'Full', M, M, X, LDX, XF, LDX )
00264 *
00265 *     Compute the CSD
00266 *
00267       CALL DORCSD( 'Y', 'Y', 'Y', 'Y', 'N', 'D', M, P, Q, XF(1,1), LDX,
00268      $             XF(1,Q+1), LDX, XF(P+1,1), LDX, XF(P+1,Q+1), LDX,
00269      $             THETA, U1, LDU1, U2, LDU2, V1T, LDV1T, V2T, LDV2T,
00270      $             WORK, LWORK, IWORK, INFO )
00271 *
00272 *     Compute X := diag(U1,U2)'*X*diag(V1,V2) - [D11 D12; D21 D22]
00273 *
00274       CALL DGEMM( 'No transpose', 'Conjugate transpose', P, Q, Q, ONE,
00275      $            X, LDX, V1T, LDV1T, ZERO, WORK, LDX )
00276 *
00277       CALL DGEMM( 'Conjugate transpose', 'No transpose', P, Q, P, ONE,
00278      $            U1, LDU1, WORK, LDX, ZERO, X, LDX )
00279 *
00280       DO I = 1, MIN(P,Q)-R
00281          X(I,I) = X(I,I) - ONE
00282       END DO
00283       DO I = 1, R
00284          X(MIN(P,Q)-R+I,MIN(P,Q)-R+I) =
00285      $           X(MIN(P,Q)-R+I,MIN(P,Q)-R+I) - COS(THETA(I))
00286       END DO
00287 *
00288       CALL DGEMM( 'No transpose', 'Conjugate transpose', P, M-Q, M-Q,
00289      $            ONE, X(1,Q+1), LDX, V2T, LDV2T, ZERO, WORK, LDX )
00290 *
00291       CALL DGEMM( 'Conjugate transpose', 'No transpose', P, M-Q, P,
00292      $            ONE, U1, LDU1, WORK, LDX, ZERO, X(1,Q+1), LDX )
00293 *
00294       DO I = 1, MIN(P,M-Q)-R
00295          X(P-I+1,M-I+1) = X(P-I+1,M-I+1) + ONE
00296       END DO
00297       DO I = 1, R
00298          X(P-(MIN(P,M-Q)-R)+1-I,M-(MIN(P,M-Q)-R)+1-I) =
00299      $      X(P-(MIN(P,M-Q)-R)+1-I,M-(MIN(P,M-Q)-R)+1-I) +
00300      $      SIN(THETA(R-I+1))
00301       END DO
00302 *
00303       CALL DGEMM( 'No transpose', 'Conjugate transpose', M-P, Q, Q, ONE,
00304      $            X(P+1,1), LDX, V1T, LDV1T, ZERO, WORK, LDX )
00305 *
00306       CALL DGEMM( 'Conjugate transpose', 'No transpose', M-P, Q, M-P,
00307      $            ONE, U2, LDU2, WORK, LDX, ZERO, X(P+1,1), LDX )
00308 *
00309       DO I = 1, MIN(M-P,Q)-R
00310          X(M-I+1,Q-I+1) = X(M-I+1,Q-I+1) - ONE
00311       END DO
00312       DO I = 1, R
00313          X(M-(MIN(M-P,Q)-R)+1-I,Q-(MIN(M-P,Q)-R)+1-I) =
00314      $             X(M-(MIN(M-P,Q)-R)+1-I,Q-(MIN(M-P,Q)-R)+1-I) -
00315      $             SIN(THETA(R-I+1))
00316       END DO
00317 *
00318       CALL DGEMM( 'No transpose', 'Conjugate transpose', M-P, M-Q, M-Q,
00319      $            ONE, X(P+1,Q+1), LDX, V2T, LDV2T, ZERO, WORK, LDX )
00320 *
00321       CALL DGEMM( 'Conjugate transpose', 'No transpose', M-P, M-Q, M-P,
00322      $            ONE, U2, LDU2, WORK, LDX, ZERO, X(P+1,Q+1), LDX )
00323 *
00324       DO I = 1, MIN(M-P,M-Q)-R
00325          X(P+I,Q+I) = X(P+I,Q+I) - ONE
00326       END DO
00327       DO I = 1, R
00328          X(P+(MIN(M-P,M-Q)-R)+I,Q+(MIN(M-P,M-Q)-R)+I) =
00329      $      X(P+(MIN(M-P,M-Q)-R)+I,Q+(MIN(M-P,M-Q)-R)+I) -
00330      $      COS(THETA(I))
00331       END DO
00332 *
00333 *     Compute norm( U1'*X11*V1 - D11 ) / ( MAX(1,P,Q)*EPS2 ) .
00334 *
00335       RESID = DLANGE( '1', P, Q, X, LDX, RWORK )
00336       RESULT( 1 ) = ( RESID / DBLE(MAX(1,P,Q)) ) / EPS2
00337 *
00338 *     Compute norm( U1'*X12*V2 - D12 ) / ( MAX(1,P,M-Q)*EPS2 ) .
00339 *
00340       RESID = DLANGE( '1', P, M-Q, X(1,Q+1), LDX, RWORK )
00341       RESULT( 2 ) = ( RESID / DBLE(MAX(1,P,M-Q)) ) / EPS2
00342 *
00343 *     Compute norm( U2'*X21*V1 - D21 ) / ( MAX(1,M-P,Q)*EPS2 ) .
00344 *
00345       RESID = DLANGE( '1', M-P, Q, X(P+1,1), LDX, RWORK )
00346       RESULT( 3 ) = ( RESID / DBLE(MAX(1,M-P,Q)) ) / EPS2
00347 *
00348 *     Compute norm( U2'*X22*V2 - D22 ) / ( MAX(1,M-P,M-Q)*EPS2 ) .
00349 *
00350       RESID = DLANGE( '1', M-P, M-Q, X(P+1,Q+1), LDX, RWORK )
00351       RESULT( 4 ) = ( RESID / DBLE(MAX(1,M-P,M-Q)) ) / EPS2
00352 *
00353 *     Compute I - U1'*U1
00354 *
00355       CALL DLASET( 'Full', P, P, ZERO, ONE, WORK, LDU1 )
00356       CALL DSYRK( 'Upper', 'Conjugate transpose', P, P, -ONE, U1, LDU1,
00357      $            ONE, WORK, LDU1 )
00358 *
00359 *     Compute norm( I - U'*U ) / ( MAX(1,P) * ULP ) .
00360 *
00361       RESID = DLANSY( '1', 'Upper', P, WORK, LDU1, RWORK )
00362       RESULT( 5 ) = ( RESID / DBLE(MAX(1,P)) ) / ULP
00363 *
00364 *     Compute I - U2'*U2
00365 *
00366       CALL DLASET( 'Full', M-P, M-P, ZERO, ONE, WORK, LDU2 )
00367       CALL DSYRK( 'Upper', 'Conjugate transpose', M-P, M-P, -ONE, U2,
00368      $            LDU2, ONE, WORK, LDU2 )
00369 *
00370 *     Compute norm( I - U2'*U2 ) / ( MAX(1,M-P) * ULP ) .
00371 *
00372       RESID = DLANSY( '1', 'Upper', M-P, WORK, LDU2, RWORK )
00373       RESULT( 6 ) = ( RESID / DBLE(MAX(1,M-P)) ) / ULP
00374 *
00375 *     Compute I - V1T*V1T'
00376 *
00377       CALL DLASET( 'Full', Q, Q, ZERO, ONE, WORK, LDV1T )
00378       CALL DSYRK( 'Upper', 'No transpose', Q, Q, -ONE, V1T, LDV1T, ONE,
00379      $            WORK, LDV1T )
00380 *
00381 *     Compute norm( I - V1T*V1T' ) / ( MAX(1,Q) * ULP ) .
00382 *
00383       RESID = DLANSY( '1', 'Upper', Q, WORK, LDV1T, RWORK )
00384       RESULT( 7 ) = ( RESID / DBLE(MAX(1,Q)) ) / ULP
00385 *
00386 *     Compute I - V2T*V2T'
00387 *
00388       CALL DLASET( 'Full', M-Q, M-Q, ZERO, ONE, WORK, LDV2T )
00389       CALL DSYRK( 'Upper', 'No transpose', M-Q, M-Q, -ONE, V2T, LDV2T,
00390      $            ONE, WORK, LDV2T )
00391 *
00392 *     Compute norm( I - V2T*V2T' ) / ( MAX(1,M-Q) * ULP ) .
00393 *
00394       RESID = DLANSY( '1', 'Upper', M-Q, WORK, LDV2T, RWORK )
00395       RESULT( 8 ) = ( RESID / DBLE(MAX(1,M-Q)) ) / ULP
00396 *
00397 *     Check sorting
00398 *
00399       RESULT(9) = REALZERO
00400       DO I = 1, R
00401          IF( THETA(I).LT.REALZERO .OR. THETA(I).GT.PIOVER2 ) THEN
00402             RESULT(9) = ULPINV
00403          END IF
00404          IF( I.GT.1) THEN
00405             IF ( THETA(I).LT.THETA(I-1) ) THEN
00406                RESULT(9) = ULPINV
00407             END IF
00408          END IF
00409       END DO
00410 *
00411       RETURN
00412 *      
00413 *     End of DCSDTS
00414 *
00415       END
00416 
 All Files Functions