LAPACK: TESTING/EIG/cgsvts.f Source File

Go to the documentation of this file.
00001       SUBROUTINE CGSVTS( M, P, N, A, AF, LDA, B, BF, LDB, U, LDU, V,
00002      $                   LDV, Q, LDQ, ALPHA, BETA, R, LDR, IWORK, WORK,
00003      $                   LWORK, RWORK, RESULT )
00004 *
00005 *  -- LAPACK test routine (version 3.1) --
00006 *     Univ. of Tennessee, Univ. of California Berkeley and NAG Ltd..
00007 *     November 2006
00008 *
00009 *     .. Scalar Arguments ..
00010       INTEGER            LDA, LDB, LDQ, LDR, LDU, LDV, LWORK, M, N, P
00011 *     ..
00012 *     .. Array Arguments ..
00013       INTEGER            IWORK( * )
00014       REAL               ALPHA( * ), BETA( * ), RESULT( 6 ), RWORK( * )
00015       COMPLEX            A( LDA, * ), AF( LDA, * ), B( LDB, * ),
00016      $                   BF( LDB, * ), Q( LDQ, * ), R( LDR, * ),
00017      $                   U( LDU, * ), V( LDV, * ), WORK( LWORK )
00018 *     ..
00019 *
00020 *  Purpose
00021 *  =======
00022 *
00023 *  CGSVTS tests CGGSVD, which computes the GSVD of an M-by-N matrix A
00024 *  and a P-by-N matrix B:
00025 *               U'*A*Q = D1*R and V'*B*Q = D2*R.
00026 *
00027 *  Arguments
00028 *  =========
00029 *
00030 *  M       (input) INTEGER
00031 *          The number of rows of the matrix A.  M >= 0.
00032 *
00033 *  P       (input) INTEGER
00034 *          The number of rows of the matrix B.  P >= 0.
00035 *
00036 *  N       (input) INTEGER
00037 *          The number of columns of the matrices A and B.  N >= 0.
00038 *
00039 *  A       (input) COMPLEX array, dimension (LDA,M)
00040 *          The M-by-N matrix A.
00041 *
00042 *  AF      (output) COMPLEX array, dimension (LDA,N)
00043 *          Details of the GSVD of A and B, as returned by CGGSVD,
00044 *          see CGGSVD for further details.
00045 *
00046 *  LDA     (input) INTEGER
00047 *          The leading dimension of the arrays A and AF.
00048 *          LDA >= max( 1,M ).
00049 *
00050 *  B       (input) COMPLEX array, dimension (LDB,P)
00051 *          On entry, the P-by-N matrix B.
00052 *
00053 *  BF      (output) COMPLEX array, dimension (LDB,N)
00054 *          Details of the GSVD of A and B, as returned by CGGSVD,
00055 *          see CGGSVD for further details.
00056 *
00057 *  LDB     (input) INTEGER
00058 *          The leading dimension of the arrays B and BF.
00059 *          LDB >= max(1,P).
00060 *
00061 *  U       (output) COMPLEX array, dimension(LDU,M)
00062 *          The M by M unitary matrix U.
00063 *
00064 *  LDU     (input) INTEGER
00065 *          The leading dimension of the array U. LDU >= max(1,M).
00066 *
00067 *  V       (output) COMPLEX array, dimension(LDV,M)
00068 *          The P by P unitary matrix V.
00069 *
00070 *  LDV     (input) INTEGER
00071 *          The leading dimension of the array V. LDV >= max(1,P).
00072 *
00073 *  Q       (output) COMPLEX array, dimension(LDQ,N)
00074 *          The N by N unitary matrix Q.
00075 *
00076 *  LDQ     (input) INTEGER
00077 *          The leading dimension of the array Q. LDQ >= max(1,N).
00078 *
00079 *  ALPHA   (output) REAL array, dimension (N)
00080 *  BETA    (output) REAL array, dimension (N)
00081 *          The generalized singular value pairs of A and B, the
00082 *          ``diagonal'' matrices D1 and D2 are constructed from
00083 *          ALPHA and BETA, see subroutine CGGSVD for details.
00084 *
00085 *  R       (output) COMPLEX array, dimension(LDQ,N)
00086 *          The upper triangular matrix R.
00087 *
00088 *  LDR     (input) INTEGER
00089 *          The leading dimension of the array R. LDR >= max(1,N).
00090 *
00091 *  IWORK   (workspace) INTEGER array, dimension (N)
00092 *
00093 *  WORK    (workspace) COMPLEX array, dimension (LWORK)
00094 *
00095 *  LWORK   (input) INTEGER
00096 *          The dimension of the array WORK,
00097 *          LWORK >= max(M,P,N)*max(M,P,N).
00098 *
00099 *  RWORK   (workspace) REAL array, dimension (max(M,P,N))
00100 *
00101 *  RESULT  (output) REAL array, dimension (5)
00102 *          The test ratios:
00103 *          RESULT(1) = norm( U'*A*Q - D1*R ) / ( MAX(M,N)*norm(A)*ULP)
00104 *          RESULT(2) = norm( V'*B*Q - D2*R ) / ( MAX(P,N)*norm(B)*ULP)
00105 *          RESULT(3) = norm( I - U'*U ) / ( M*ULP )
00106 *          RESULT(4) = norm( I - V'*V ) / ( P*ULP )
00107 *          RESULT(5) = norm( I - Q'*Q ) / ( N*ULP )
00108 *          RESULT(6) = 0        if ALPHA is in decreasing order;
00109 *                    = ULPINV   otherwise.
00110 *
00111 *  =====================================================================
00112 *
00113 *     .. Parameters ..
00114       REAL               ZERO, ONE
00115       PARAMETER          ( ZERO = 0.0E+0, ONE = 1.0E+0 )
00116       COMPLEX            CZERO, CONE
00117       PARAMETER          ( CZERO = ( 0.0E+0, 0.0E+0 ),
00118      $                   CONE = ( 1.0E+0, 0.0E+0 ) )
00119 *     ..
00120 *     .. Local Scalars ..
00121       INTEGER            I, INFO, J, K, L
00122       REAL               ANORM, BNORM, RESID, TEMP, ULP, ULPINV, UNFL
00123 *     ..
00124 *     .. External Functions ..
00125       REAL               CLANGE, CLANHE, SLAMCH
00126       EXTERNAL           CLANGE, CLANHE, SLAMCH
00127 *     ..
00128 *     .. External Subroutines ..
00129       EXTERNAL           CGEMM, CGGSVD, CHERK, CLACPY, CLASET, SCOPY
00130 *     ..
00131 *     .. Intrinsic Functions ..
00132       INTRINSIC          MAX, MIN, REAL
00133 *     ..
00134 *     .. Executable Statements ..
00135 *
00136       ULP = SLAMCH( 'Precision' )
00137       ULPINV = ONE / ULP
00138       UNFL = SLAMCH( 'Safe minimum' )
00139 *
00140 *     Copy the matrix A to the array AF.
00141 *
00142       CALL CLACPY( 'Full', M, N, A, LDA, AF, LDA )
00143       CALL CLACPY( 'Full', P, N, B, LDB, BF, LDB )
00144 *
00145       ANORM = MAX( CLANGE( '1', M, N, A, LDA, RWORK ), UNFL )
00146       BNORM = MAX( CLANGE( '1', P, N, B, LDB, RWORK ), UNFL )
00147 *
00148 *     Factorize the matrices A and B in the arrays AF and BF.
00149 *
00150       CALL CGGSVD( 'U', 'V', 'Q', M, N, P, K, L, AF, LDA, BF, LDB,
00151      $             ALPHA, BETA, U, LDU, V, LDV, Q, LDQ, WORK, RWORK,
00152      $             IWORK, INFO )
00153 *
00154 *     Copy R
00155 *
00156       DO 20 I = 1, MIN( K+L, M )
00157          DO 10 J = I, K + L
00158             R( I, J ) = AF( I, N-K-L+J )
00159    10    CONTINUE
00160    20 CONTINUE
00161 *
00162       IF( M-K-L.LT.0 ) THEN
00163          DO 40 I = M + 1, K + L
00164             DO 30 J = I, K + L
00165                R( I, J ) = BF( I-K, N-K-L+J )
00166    30       CONTINUE
00167    40    CONTINUE
00168       END IF
00169 *
00170 *     Compute A:= U'*A*Q - D1*R
00171 *
00172       CALL CGEMM( 'No transpose', 'No transpose', M, N, N, CONE, A, LDA,
00173      $            Q, LDQ, CZERO, WORK, LDA )
00174 *
00175       CALL CGEMM( 'Conjugate transpose', 'No transpose', M, N, M, CONE,
00176      $            U, LDU, WORK, LDA, CZERO, A, LDA )
00177 *
00178       DO 60 I = 1, K
00179          DO 50 J = I, K + L
00180             A( I, N-K-L+J ) = A( I, N-K-L+J ) - R( I, J )
00181    50    CONTINUE
00182    60 CONTINUE
00183 *
00184       DO 80 I = K + 1, MIN( K+L, M )
00185          DO 70 J = I, K + L
00186             A( I, N-K-L+J ) = A( I, N-K-L+J ) - ALPHA( I )*R( I, J )
00187    70    CONTINUE
00188    80 CONTINUE
00189 *
00190 *     Compute norm( U'*A*Q - D1*R ) / ( MAX(1,M,N)*norm(A)*ULP ) .
00191 *
00192       RESID = CLANGE( '1', M, N, A, LDA, RWORK )
00193       IF( ANORM.GT.ZERO ) THEN
00194          RESULT( 1 ) = ( ( RESID / REAL( MAX( 1, M, N ) ) ) / ANORM ) /
00195      $                 ULP
00196       ELSE
00197          RESULT( 1 ) = ZERO
00198       END IF
00199 *
00200 *     Compute B := V'*B*Q - D2*R
00201 *
00202       CALL CGEMM( 'No transpose', 'No transpose', P, N, N, CONE, B, LDB,
00203      $            Q, LDQ, CZERO, WORK, LDB )
00204 *
00205       CALL CGEMM( 'Conjugate transpose', 'No transpose', P, N, P, CONE,
00206      $            V, LDV, WORK, LDB, CZERO, B, LDB )
00207 *
00208       DO 100 I = 1, L
00209          DO 90 J = I, L
00210             B( I, N-L+J ) = B( I, N-L+J ) - BETA( K+I )*R( K+I, K+J )
00211    90    CONTINUE
00212   100 CONTINUE
00213 *
00214 *     Compute norm( V'*B*Q - D2*R ) / ( MAX(P,N)*norm(B)*ULP ) .
00215 *
00216       RESID = CLANGE( '1', P, N, B, LDB, RWORK )
00217       IF( BNORM.GT.ZERO ) THEN
00218          RESULT( 2 ) = ( ( RESID / REAL( MAX( 1, P, N ) ) ) / BNORM ) /
00219      $                 ULP
00220       ELSE
00221          RESULT( 2 ) = ZERO
00222       END IF
00223 *
00224 *     Compute I - U'*U
00225 *
00226       CALL CLASET( 'Full', M, M, CZERO, CONE, WORK, LDQ )
00227       CALL CHERK( 'Upper', 'Conjugate transpose', M, M, -ONE, U, LDU,
00228      $            ONE, WORK, LDU )
00229 *
00230 *     Compute norm( I - U'*U ) / ( M * ULP ) .
00231 *
00232       RESID = CLANHE( '1', 'Upper', M, WORK, LDU, RWORK )
00233       RESULT( 3 ) = ( RESID / REAL( MAX( 1, M ) ) ) / ULP
00234 *
00235 *     Compute I - V'*V
00236 *
00237       CALL CLASET( 'Full', P, P, CZERO, CONE, WORK, LDV )
00238       CALL CHERK( 'Upper', 'Conjugate transpose', P, P, -ONE, V, LDV,
00239      $            ONE, WORK, LDV )
00240 *
00241 *     Compute norm( I - V'*V ) / ( P * ULP ) .
00242 *
00243       RESID = CLANHE( '1', 'Upper', P, WORK, LDV, RWORK )
00244       RESULT( 4 ) = ( RESID / REAL( MAX( 1, P ) ) ) / ULP
00245 *
00246 *     Compute I - Q'*Q
00247 *
00248       CALL CLASET( 'Full', N, N, CZERO, CONE, WORK, LDQ )
00249       CALL CHERK( 'Upper', 'Conjugate transpose', N, N, -ONE, Q, LDQ,
00250      $            ONE, WORK, LDQ )
00251 *
00252 *     Compute norm( I - Q'*Q ) / ( N * ULP ) .
00253 *
00254       RESID = CLANHE( '1', 'Upper', N, WORK, LDQ, RWORK )
00255       RESULT( 5 ) = ( RESID / REAL( MAX( 1, N ) ) ) / ULP
00256 *
00257 *     Check sorting
00258 *
00259       CALL SCOPY( N, ALPHA, 1, RWORK, 1 )
00260       DO 110 I = K + 1, MIN( K+L, M )
00261          J = IWORK( I )
00262          IF( I.NE.J ) THEN
00263             TEMP = RWORK( I )
00264             RWORK( I ) = RWORK( J )
00265             RWORK( J ) = TEMP
00266          END IF
00267   110 CONTINUE
00268 *
00269       RESULT( 6 ) = ZERO
00270       DO 120 I = K + 1, MIN( K+L, M ) - 1
00271          IF( RWORK( I ).LT.RWORK( I+1 ) )
00272      $      RESULT( 6 ) = ULPINV
00273   120 CONTINUE
00274 *
00275       RETURN
00276 *
00277 *     End of CGSVTS
00278 *
00279       END