LAPACK: TESTING/EIG/dgsvts.f Source File

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