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