LAPACK  3.4.1
LAPACK: Linear Algebra PACKage
cgsvts.f
Go to the documentation of this file.
00001 *> \brief \b CGSVTS
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 CGSVTS( M, P, N, A, AF, LDA, B, BF, LDB, U, LDU, V,
00012 *                          LDV, Q, LDQ, ALPHA, BETA, R, LDR, IWORK, WORK,
00013 *                          LWORK, RWORK, RESULT )
00014 * 
00015 *       .. Scalar Arguments ..
00016 *       INTEGER            LDA, LDB, LDQ, LDR, LDU, LDV, LWORK, M, N, P
00017 *       ..
00018 *       .. Array Arguments ..
00019 *       INTEGER            IWORK( * )
00020 *       REAL               ALPHA( * ), BETA( * ), RESULT( 6 ), RWORK( * )
00021 *       COMPLEX            A( LDA, * ), AF( LDA, * ), B( LDB, * ),
00022 *      $                   BF( LDB, * ), Q( LDQ, * ), R( LDR, * ),
00023 *      $                   U( LDU, * ), V( LDV, * ), WORK( LWORK )
00024 *       ..
00025 *  
00026 *
00027 *> \par Purpose:
00028 *  =============
00029 *>
00030 *> \verbatim
00031 *>
00032 *> CGSVTS tests CGGSVD, which computes the GSVD of an M-by-N matrix A
00033 *> and a P-by-N matrix B:
00034 *>              U'*A*Q = D1*R and V'*B*Q = D2*R.
00035 *> \endverbatim
00036 *
00037 *  Arguments:
00038 *  ==========
00039 *
00040 *> \param[in] M
00041 *> \verbatim
00042 *>          M is INTEGER
00043 *>          The number of rows of the matrix A.  M >= 0.
00044 *> \endverbatim
00045 *>
00046 *> \param[in] P
00047 *> \verbatim
00048 *>          P is INTEGER
00049 *>          The number of rows of the matrix B.  P >= 0.
00050 *> \endverbatim
00051 *>
00052 *> \param[in] N
00053 *> \verbatim
00054 *>          N is INTEGER
00055 *>          The number of columns of the matrices A and B.  N >= 0.
00056 *> \endverbatim
00057 *>
00058 *> \param[in] A
00059 *> \verbatim
00060 *>          A is COMPLEX array, dimension (LDA,M)
00061 *>          The M-by-N matrix A.
00062 *> \endverbatim
00063 *>
00064 *> \param[out] AF
00065 *> \verbatim
00066 *>          AF is COMPLEX array, dimension (LDA,N)
00067 *>          Details of the GSVD of A and B, as returned by CGGSVD,
00068 *>          see CGGSVD for further details.
00069 *> \endverbatim
00070 *>
00071 *> \param[in] LDA
00072 *> \verbatim
00073 *>          LDA is INTEGER
00074 *>          The leading dimension of the arrays A and AF.
00075 *>          LDA >= max( 1,M ).
00076 *> \endverbatim
00077 *>
00078 *> \param[in] B
00079 *> \verbatim
00080 *>          B is COMPLEX array, dimension (LDB,P)
00081 *>          On entry, the P-by-N matrix B.
00082 *> \endverbatim
00083 *>
00084 *> \param[out] BF
00085 *> \verbatim
00086 *>          BF is COMPLEX array, dimension (LDB,N)
00087 *>          Details of the GSVD of A and B, as returned by CGGSVD,
00088 *>          see CGGSVD for further details.
00089 *> \endverbatim
00090 *>
00091 *> \param[in] LDB
00092 *> \verbatim
00093 *>          LDB is INTEGER
00094 *>          The leading dimension of the arrays B and BF.
00095 *>          LDB >= max(1,P).
00096 *> \endverbatim
00097 *>
00098 *> \param[out] U
00099 *> \verbatim
00100 *>          U is COMPLEX array, dimension(LDU,M)
00101 *>          The M by M unitary matrix U.
00102 *> \endverbatim
00103 *>
00104 *> \param[in] LDU
00105 *> \verbatim
00106 *>          LDU is INTEGER
00107 *>          The leading dimension of the array U. LDU >= max(1,M).
00108 *> \endverbatim
00109 *>
00110 *> \param[out] V
00111 *> \verbatim
00112 *>          V is COMPLEX array, dimension(LDV,M)
00113 *>          The P by P unitary matrix V.
00114 *> \endverbatim
00115 *>
00116 *> \param[in] LDV
00117 *> \verbatim
00118 *>          LDV is INTEGER
00119 *>          The leading dimension of the array V. LDV >= max(1,P).
00120 *> \endverbatim
00121 *>
00122 *> \param[out] Q
00123 *> \verbatim
00124 *>          Q is COMPLEX array, dimension(LDQ,N)
00125 *>          The N by N unitary matrix Q.
00126 *> \endverbatim
00127 *>
00128 *> \param[in] LDQ
00129 *> \verbatim
00130 *>          LDQ is INTEGER
00131 *>          The leading dimension of the array Q. LDQ >= max(1,N).
00132 *> \endverbatim
00133 *>
00134 *> \param[out] ALPHA
00135 *> \verbatim
00136 *>          ALPHA is REAL array, dimension (N)
00137 *> \endverbatim
00138 *>
00139 *> \param[out] BETA
00140 *> \verbatim
00141 *>          BETA is REAL array, dimension (N)
00142 *>
00143 *>          The generalized singular value pairs of A and B, the
00144 *>          ``diagonal'' matrices D1 and D2 are constructed from
00145 *>          ALPHA and BETA, see subroutine CGGSVD for details.
00146 *> \endverbatim
00147 *>
00148 *> \param[out] R
00149 *> \verbatim
00150 *>          R is COMPLEX array, dimension(LDQ,N)
00151 *>          The upper triangular matrix R.
00152 *> \endverbatim
00153 *>
00154 *> \param[in] LDR
00155 *> \verbatim
00156 *>          LDR is INTEGER
00157 *>          The leading dimension of the array R. LDR >= max(1,N).
00158 *> \endverbatim
00159 *>
00160 *> \param[out] IWORK
00161 *> \verbatim
00162 *>          IWORK is INTEGER array, dimension (N)
00163 *> \endverbatim
00164 *>
00165 *> \param[out] WORK
00166 *> \verbatim
00167 *>          WORK is COMPLEX array, dimension (LWORK)
00168 *> \endverbatim
00169 *>
00170 *> \param[in] LWORK
00171 *> \verbatim
00172 *>          LWORK is INTEGER
00173 *>          The dimension of the array WORK,
00174 *>          LWORK >= max(M,P,N)*max(M,P,N).
00175 *> \endverbatim
00176 *>
00177 *> \param[out] RWORK
00178 *> \verbatim
00179 *>          RWORK is REAL array, dimension (max(M,P,N))
00180 *> \endverbatim
00181 *>
00182 *> \param[out] RESULT
00183 *> \verbatim
00184 *>          RESULT is REAL array, dimension (5)
00185 *>          The test ratios:
00186 *>          RESULT(1) = norm( U'*A*Q - D1*R ) / ( MAX(M,N)*norm(A)*ULP)
00187 *>          RESULT(2) = norm( V'*B*Q - D2*R ) / ( MAX(P,N)*norm(B)*ULP)
00188 *>          RESULT(3) = norm( I - U'*U ) / ( M*ULP )
00189 *>          RESULT(4) = norm( I - V'*V ) / ( P*ULP )
00190 *>          RESULT(5) = norm( I - Q'*Q ) / ( N*ULP )
00191 *>          RESULT(6) = 0        if ALPHA is in decreasing order;
00192 *>                    = ULPINV   otherwise.
00193 *> \endverbatim
00194 *
00195 *  Authors:
00196 *  ========
00197 *
00198 *> \author Univ. of Tennessee 
00199 *> \author Univ. of California Berkeley 
00200 *> \author Univ. of Colorado Denver 
00201 *> \author NAG Ltd. 
00202 *
00203 *> \date November 2011
00204 *
00205 *> \ingroup complex_eig
00206 *
00207 *  =====================================================================
00208       SUBROUTINE CGSVTS( M, P, N, A, AF, LDA, B, BF, LDB, U, LDU, V,
00209      $                   LDV, Q, LDQ, ALPHA, BETA, R, LDR, IWORK, WORK,
00210      $                   LWORK, RWORK, RESULT )
00211 *
00212 *  -- LAPACK test routine (version 3.4.0) --
00213 *  -- LAPACK is a software package provided by Univ. of Tennessee,    --
00214 *  -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
00215 *     November 2011
00216 *
00217 *     .. Scalar Arguments ..
00218       INTEGER            LDA, LDB, LDQ, LDR, LDU, LDV, LWORK, M, N, P
00219 *     ..
00220 *     .. Array Arguments ..
00221       INTEGER            IWORK( * )
00222       REAL               ALPHA( * ), BETA( * ), RESULT( 6 ), RWORK( * )
00223       COMPLEX            A( LDA, * ), AF( LDA, * ), B( LDB, * ),
00224      $                   BF( LDB, * ), Q( LDQ, * ), R( LDR, * ),
00225      $                   U( LDU, * ), V( LDV, * ), WORK( LWORK )
00226 *     ..
00227 *
00228 *  =====================================================================
00229 *
00230 *     .. Parameters ..
00231       REAL               ZERO, ONE
00232       PARAMETER          ( ZERO = 0.0E+0, ONE = 1.0E+0 )
00233       COMPLEX            CZERO, CONE
00234       PARAMETER          ( CZERO = ( 0.0E+0, 0.0E+0 ),
00235      $                   CONE = ( 1.0E+0, 0.0E+0 ) )
00236 *     ..
00237 *     .. Local Scalars ..
00238       INTEGER            I, INFO, J, K, L
00239       REAL               ANORM, BNORM, RESID, TEMP, ULP, ULPINV, UNFL
00240 *     ..
00241 *     .. External Functions ..
00242       REAL               CLANGE, CLANHE, SLAMCH
00243       EXTERNAL           CLANGE, CLANHE, SLAMCH
00244 *     ..
00245 *     .. External Subroutines ..
00246       EXTERNAL           CGEMM, CGGSVD, CHERK, CLACPY, CLASET, SCOPY
00247 *     ..
00248 *     .. Intrinsic Functions ..
00249       INTRINSIC          MAX, MIN, REAL
00250 *     ..
00251 *     .. Executable Statements ..
00252 *
00253       ULP = SLAMCH( 'Precision' )
00254       ULPINV = ONE / ULP
00255       UNFL = SLAMCH( 'Safe minimum' )
00256 *
00257 *     Copy the matrix A to the array AF.
00258 *
00259       CALL CLACPY( 'Full', M, N, A, LDA, AF, LDA )
00260       CALL CLACPY( 'Full', P, N, B, LDB, BF, LDB )
00261 *
00262       ANORM = MAX( CLANGE( '1', M, N, A, LDA, RWORK ), UNFL )
00263       BNORM = MAX( CLANGE( '1', P, N, B, LDB, RWORK ), UNFL )
00264 *
00265 *     Factorize the matrices A and B in the arrays AF and BF.
00266 *
00267       CALL CGGSVD( 'U', 'V', 'Q', M, N, P, K, L, AF, LDA, BF, LDB,
00268      $             ALPHA, BETA, U, LDU, V, LDV, Q, LDQ, WORK, RWORK,
00269      $             IWORK, INFO )
00270 *
00271 *     Copy R
00272 *
00273       DO 20 I = 1, MIN( K+L, M )
00274          DO 10 J = I, K + L
00275             R( I, J ) = AF( I, N-K-L+J )
00276    10    CONTINUE
00277    20 CONTINUE
00278 *
00279       IF( M-K-L.LT.0 ) THEN
00280          DO 40 I = M + 1, K + L
00281             DO 30 J = I, K + L
00282                R( I, J ) = BF( I-K, N-K-L+J )
00283    30       CONTINUE
00284    40    CONTINUE
00285       END IF
00286 *
00287 *     Compute A:= U'*A*Q - D1*R
00288 *
00289       CALL CGEMM( 'No transpose', 'No transpose', M, N, N, CONE, A, LDA,
00290      $            Q, LDQ, CZERO, WORK, LDA )
00291 *
00292       CALL CGEMM( 'Conjugate transpose', 'No transpose', M, N, M, CONE,
00293      $            U, LDU, WORK, LDA, CZERO, A, LDA )
00294 *
00295       DO 60 I = 1, K
00296          DO 50 J = I, K + L
00297             A( I, N-K-L+J ) = A( I, N-K-L+J ) - R( I, J )
00298    50    CONTINUE
00299    60 CONTINUE
00300 *
00301       DO 80 I = K + 1, MIN( K+L, M )
00302          DO 70 J = I, K + L
00303             A( I, N-K-L+J ) = A( I, N-K-L+J ) - ALPHA( I )*R( I, J )
00304    70    CONTINUE
00305    80 CONTINUE
00306 *
00307 *     Compute norm( U'*A*Q - D1*R ) / ( MAX(1,M,N)*norm(A)*ULP ) .
00308 *
00309       RESID = CLANGE( '1', M, N, A, LDA, RWORK )
00310       IF( ANORM.GT.ZERO ) THEN
00311          RESULT( 1 ) = ( ( RESID / REAL( MAX( 1, M, N ) ) ) / ANORM ) /
00312      $                 ULP
00313       ELSE
00314          RESULT( 1 ) = ZERO
00315       END IF
00316 *
00317 *     Compute B := V'*B*Q - D2*R
00318 *
00319       CALL CGEMM( 'No transpose', 'No transpose', P, N, N, CONE, B, LDB,
00320      $            Q, LDQ, CZERO, WORK, LDB )
00321 *
00322       CALL CGEMM( 'Conjugate transpose', 'No transpose', P, N, P, CONE,
00323      $            V, LDV, WORK, LDB, CZERO, B, LDB )
00324 *
00325       DO 100 I = 1, L
00326          DO 90 J = I, L
00327             B( I, N-L+J ) = B( I, N-L+J ) - BETA( K+I )*R( K+I, K+J )
00328    90    CONTINUE
00329   100 CONTINUE
00330 *
00331 *     Compute norm( V'*B*Q - D2*R ) / ( MAX(P,N)*norm(B)*ULP ) .
00332 *
00333       RESID = CLANGE( '1', P, N, B, LDB, RWORK )
00334       IF( BNORM.GT.ZERO ) THEN
00335          RESULT( 2 ) = ( ( RESID / REAL( MAX( 1, P, N ) ) ) / BNORM ) /
00336      $                 ULP
00337       ELSE
00338          RESULT( 2 ) = ZERO
00339       END IF
00340 *
00341 *     Compute I - U'*U
00342 *
00343       CALL CLASET( 'Full', M, M, CZERO, CONE, WORK, LDQ )
00344       CALL CHERK( 'Upper', 'Conjugate transpose', M, M, -ONE, U, LDU,
00345      $            ONE, WORK, LDU )
00346 *
00347 *     Compute norm( I - U'*U ) / ( M * ULP ) .
00348 *
00349       RESID = CLANHE( '1', 'Upper', M, WORK, LDU, RWORK )
00350       RESULT( 3 ) = ( RESID / REAL( MAX( 1, M ) ) ) / ULP
00351 *
00352 *     Compute I - V'*V
00353 *
00354       CALL CLASET( 'Full', P, P, CZERO, CONE, WORK, LDV )
00355       CALL CHERK( 'Upper', 'Conjugate transpose', P, P, -ONE, V, LDV,
00356      $            ONE, WORK, LDV )
00357 *
00358 *     Compute norm( I - V'*V ) / ( P * ULP ) .
00359 *
00360       RESID = CLANHE( '1', 'Upper', P, WORK, LDV, RWORK )
00361       RESULT( 4 ) = ( RESID / REAL( MAX( 1, P ) ) ) / ULP
00362 *
00363 *     Compute I - Q'*Q
00364 *
00365       CALL CLASET( 'Full', N, N, CZERO, CONE, WORK, LDQ )
00366       CALL CHERK( 'Upper', 'Conjugate transpose', N, N, -ONE, Q, LDQ,
00367      $            ONE, WORK, LDQ )
00368 *
00369 *     Compute norm( I - Q'*Q ) / ( N * ULP ) .
00370 *
00371       RESID = CLANHE( '1', 'Upper', N, WORK, LDQ, RWORK )
00372       RESULT( 5 ) = ( RESID / REAL( MAX( 1, N ) ) ) / ULP
00373 *
00374 *     Check sorting
00375 *
00376       CALL SCOPY( N, ALPHA, 1, RWORK, 1 )
00377       DO 110 I = K + 1, MIN( K+L, M )
00378          J = IWORK( I )
00379          IF( I.NE.J ) THEN
00380             TEMP = RWORK( I )
00381             RWORK( I ) = RWORK( J )
00382             RWORK( J ) = TEMP
00383          END IF
00384   110 CONTINUE
00385 *
00386       RESULT( 6 ) = ZERO
00387       DO 120 I = K + 1, MIN( K+L, M ) - 1
00388          IF( RWORK( I ).LT.RWORK( I+1 ) )
00389      $      RESULT( 6 ) = ULPINV
00390   120 CONTINUE
00391 *
00392       RETURN
00393 *
00394 *     End of CGSVTS
00395 *
00396       END
 All Files Functions