Refs #324. Upgrade LAPACK to 3.5.0 version.
This commit is contained in:
Zhang Xianyi
@@ -209,14 +209,23 @@
|
||||
*
|
||||
* CSD
|
||||
*
|
||||
WRITE( IOUNIT, FMT = 9920 )1
|
||||
WRITE( IOUNIT, FMT = 9921 )2
|
||||
WRITE( IOUNIT, FMT = 9922 )3
|
||||
WRITE( IOUNIT, FMT = 9923 )4
|
||||
WRITE( IOUNIT, FMT = 9924 )5
|
||||
WRITE( IOUNIT, FMT = 9925 )6
|
||||
WRITE( IOUNIT, FMT = 9926 )7
|
||||
WRITE( IOUNIT, FMT = 9927 )8
|
||||
WRITE( IOUNIT, FMT = 9910 )
|
||||
WRITE( IOUNIT, FMT = 9911 )1
|
||||
WRITE( IOUNIT, FMT = 9912 )2
|
||||
WRITE( IOUNIT, FMT = 9913 )3
|
||||
WRITE( IOUNIT, FMT = 9914 )4
|
||||
WRITE( IOUNIT, FMT = 9915 )5
|
||||
WRITE( IOUNIT, FMT = 9916 )6
|
||||
WRITE( IOUNIT, FMT = 9917 )7
|
||||
WRITE( IOUNIT, FMT = 9918 )8
|
||||
WRITE( IOUNIT, FMT = 9919 )9
|
||||
WRITE( IOUNIT, FMT = 9920 )
|
||||
WRITE( IOUNIT, FMT = 9921 )10
|
||||
WRITE( IOUNIT, FMT = 9922 )11
|
||||
WRITE( IOUNIT, FMT = 9923 )12
|
||||
WRITE( IOUNIT, FMT = 9924 )13
|
||||
WRITE( IOUNIT, FMT = 9925 )14
|
||||
WRITE( IOUNIT, FMT = 9926 )15
|
||||
END IF
|
||||
*
|
||||
9999 FORMAT( 1X, A )
|
||||
@@ -291,18 +300,29 @@
|
||||
*
|
||||
* CSD test ratio
|
||||
*
|
||||
9920 FORMAT( 3X, I2, ': norm( U1'' * X11 * V1 - C ) / ( max( P, Q)',
|
||||
9910 FORMAT( 3X, '2-by-2 CSD' )
|
||||
9911 FORMAT( 3X, I2, ': norm( U1'' * X11 * V1 - C ) / ( max( P, Q)',
|
||||
$ ' * max(norm(I-X''*X),EPS) )' )
|
||||
9921 FORMAT( 3X, I2, ': norm( U1'' * X12 * V2-(-S)) / ( max( P,',
|
||||
9912 FORMAT( 3X, I2, ': norm( U1'' * X12 * V2-(-S)) / ( max( P,',
|
||||
$ 'M-Q) * max(norm(I-X''*X),EPS) )' )
|
||||
9922 FORMAT( 3X, I2, ': norm( U2'' * X21 * V1 - S ) / ( max(M-P,',
|
||||
9913 FORMAT( 3X, I2, ': norm( U2'' * X21 * V1 - S ) / ( max(M-P,',
|
||||
$ ' Q) * max(norm(I-X''*X),EPS) )' )
|
||||
9923 FORMAT( 3X, I2, ': norm( U2'' * X22 * V2 - C ) / ( max(M-P,',
|
||||
9914 FORMAT( 3X, I2, ': norm( U2'' * X22 * V2 - C ) / ( max(M-P,',
|
||||
$ 'M-Q) * max(norm(I-X''*X),EPS) )' )
|
||||
9924 FORMAT( 3X, I2, ': norm( I - U1''*U1 ) / ( P * EPS )' )
|
||||
9925 FORMAT( 3X, I2, ': norm( I - U2''*U2 ) / ( (M-P) * EPS )' )
|
||||
9926 FORMAT( 3X, I2, ': norm( I - V1''*V1 ) / ( Q * EPS )' )
|
||||
9927 FORMAT( 3X, I2, ': norm( I - V2''*V2 ) / ( (M-Q) * EPS )' )
|
||||
9915 FORMAT( 3X, I2, ': norm( I - U1''*U1 ) / ( P * EPS )' )
|
||||
9916 FORMAT( 3X, I2, ': norm( I - U2''*U2 ) / ( (M-P) * EPS )' )
|
||||
9917 FORMAT( 3X, I2, ': norm( I - V1''*V1 ) / ( Q * EPS )' )
|
||||
9918 FORMAT( 3X, I2, ': norm( I - V2''*V2 ) / ( (M-Q) * EPS )' )
|
||||
9919 FORMAT( 3X, I2, ': principal angle ordering ( 0 or ULP )' )
|
||||
9920 FORMAT( 3X, '2-by-1 CSD' )
|
||||
9921 FORMAT( 3X, I2, ': norm( U1'' * X11 * V1 - C ) / ( max( P, Q)',
|
||||
$ ' * max(norm(I-X''*X),EPS) )' )
|
||||
9922 FORMAT( 3X, I2, ': norm( U2'' * X21 * V1 - S ) / ( max( M-P,',
|
||||
$ 'Q) * max(norm(I-X''*X),EPS) )' )
|
||||
9923 FORMAT( 3X, I2, ': norm( I - U1''*U1 ) / ( P * EPS )' )
|
||||
9924 FORMAT( 3X, I2, ': norm( I - U2''*U2 ) / ( (M-P) * EPS )' )
|
||||
9925 FORMAT( 3X, I2, ': norm( I - V1''*V1 ) / ( Q * EPS )' )
|
||||
9926 FORMAT( 3X, I2, ': principal angle ordering ( 0 or ULP )' )
|
||||
RETURN
|
||||
*
|
||||
* End of ALAHDG
|
||||
|
||||
@@ -1026,17 +1026,17 @@
|
||||
*> \author Univ. of Colorado Denver
|
||||
*> \author NAG Ltd.
|
||||
*
|
||||
*> \date April 2012
|
||||
*> \date November 2013
|
||||
*
|
||||
*> \ingroup complex_eig
|
||||
*
|
||||
* =====================================================================
|
||||
PROGRAM CCHKEE
|
||||
*
|
||||
* -- LAPACK test routine (version 3.4.1) --
|
||||
* -- LAPACK test routine (version 3.5.0) --
|
||||
* -- LAPACK is a software package provided by Univ. of Tennessee, --
|
||||
* -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
|
||||
* April 2012
|
||||
* November 2013
|
||||
*
|
||||
* =====================================================================
|
||||
*
|
||||
@@ -1129,6 +1129,10 @@
|
||||
* ..
|
||||
* .. Executable Statements ..
|
||||
*
|
||||
A = 0.0
|
||||
B = 0.0
|
||||
C = 0.0
|
||||
DC = 0.0
|
||||
S1 = SECOND( )
|
||||
FATAL = .FALSE.
|
||||
NUNIT = NOUT
|
||||
|
||||
@@ -205,13 +205,16 @@
|
||||
*
|
||||
* .. Parameters ..
|
||||
INTEGER NTESTS
|
||||
PARAMETER ( NTESTS = 9 )
|
||||
PARAMETER ( NTESTS = 15 )
|
||||
INTEGER NTYPES
|
||||
PARAMETER ( NTYPES = 3 )
|
||||
REAL GAPDIGIT, ORTH, PIOVER2, TEN
|
||||
PARAMETER ( NTYPES = 4 )
|
||||
REAL GAPDIGIT, ORTH, PIOVER2, REALONE, REALZERO, TEN
|
||||
PARAMETER ( GAPDIGIT = 10.0E0, ORTH = 1.0E-4,
|
||||
$ PIOVER2 = 1.57079632679489662E0,
|
||||
$ TEN = 10.0D0 )
|
||||
$ REALONE = 1.0E0, REALZERO = 0.0E0,
|
||||
$ TEN = 10.0E0 )
|
||||
COMPLEX ONE, ZERO
|
||||
PARAMETER ( ONE = (1.0E0,0.0E0), ZERO = (0.0E0,0.0E0) )
|
||||
* ..
|
||||
* .. Local Scalars ..
|
||||
LOGICAL FIRSTT
|
||||
@@ -231,8 +234,8 @@
|
||||
INTRINSIC ABS, MIN
|
||||
* ..
|
||||
* .. External Functions ..
|
||||
REAL SLARND
|
||||
EXTERNAL SLARND
|
||||
REAL SLARAN, SLARND
|
||||
EXTERNAL SLARAN, SLARND
|
||||
* ..
|
||||
* .. Executable Statements ..
|
||||
*
|
||||
@@ -286,7 +289,7 @@
|
||||
$ ORTH*SLARND(2,ISEED)
|
||||
END DO
|
||||
END DO
|
||||
ELSE
|
||||
ELSE IF( IMAT.EQ.3 ) THEN
|
||||
R = MIN( P, M-P, Q, M-Q )
|
||||
DO I = 1, R+1
|
||||
THETA(I) = TEN**(-SLARND(1,ISEED)*GAPDIGIT)
|
||||
@@ -298,9 +301,18 @@
|
||||
THETA(I) = PIOVER2 * THETA(I) / THETA(R+1)
|
||||
END DO
|
||||
CALL CLACSG( M, P, Q, THETA, ISEED, X, LDX, WORK )
|
||||
ELSE
|
||||
CALL CLASET( 'F', M, M, ZERO, ONE, X, LDX )
|
||||
DO I = 1, M
|
||||
J = INT( SLARAN( ISEED ) * M ) + 1
|
||||
IF( J .NE. I ) THEN
|
||||
CALL CSROT( M, X(1+(I-1)*LDX), 1, X(1+(J-1)*LDX),
|
||||
$ 1, REALZERO, REALONE )
|
||||
END IF
|
||||
END DO
|
||||
END IF
|
||||
*
|
||||
NT = 9
|
||||
NT = 15
|
||||
*
|
||||
CALL CCSDTS( M, P, Q, X, XF, LDX, U1, LDU1, U2, LDU2, V1T,
|
||||
$ LDV1T, V2T, LDV2T, THETA, IWORK, WORK, LWORK,
|
||||
|
||||
@@ -17,7 +17,7 @@
|
||||
* ..
|
||||
* .. Array Arguments ..
|
||||
* INTEGER IWORK( * )
|
||||
* REAL RESULT( 9 ), RWORK( * ), THETA( * )
|
||||
* REAL RESULT( 15 ), RWORK( * ), THETA( * )
|
||||
* COMPLEX U1( LDU1, * ), U2( LDU2, * ), V1T( LDV1T, * ),
|
||||
* $ V2T( LDV2T, * ), WORK( LWORK ), X( LDX, * ),
|
||||
* $ XF( LDX, * )
|
||||
@@ -47,6 +47,21 @@
|
||||
*> [ 0 0 0 | I 0 0 ] [ D21 D22 ]
|
||||
*> [ 0 S 0 | 0 C 0 ]
|
||||
*> [ 0 0 I | 0 0 0 ]
|
||||
*>
|
||||
*> and also SORCSD2BY1, which, given
|
||||
*> Q
|
||||
*> [ X11 ] P ,
|
||||
*> [ X21 ] M-P
|
||||
*>
|
||||
*> computes the 2-by-1 CSD
|
||||
*>
|
||||
*> [ I 0 0 ]
|
||||
*> [ 0 C 0 ]
|
||||
*> [ 0 0 0 ]
|
||||
*> [ U1 ]**T * [ X11 ] * V1 = [----------] = [ D11 ] ,
|
||||
*> [ U2 ] [ X21 ] [ 0 0 0 ] [ D21 ]
|
||||
*> [ 0 S 0 ]
|
||||
*> [ 0 0 I ]
|
||||
*> \endverbatim
|
||||
*
|
||||
* Arguments:
|
||||
@@ -171,8 +186,9 @@
|
||||
*>
|
||||
*> \param[out] RESULT
|
||||
*> \verbatim
|
||||
*> RESULT is REAL array, dimension (9)
|
||||
*> RESULT is REAL array, dimension (15)
|
||||
*> The test ratios:
|
||||
*> First, the 2-by-2 CSD:
|
||||
*> RESULT(1) = norm( U1'*X11*V1 - D11 ) / ( MAX(1,P,Q)*EPS2 )
|
||||
*> RESULT(2) = norm( U1'*X12*V2 - D12 ) / ( MAX(1,P,M-Q)*EPS2 )
|
||||
*> RESULT(3) = norm( U2'*X21*V1 - D21 ) / ( MAX(1,M-P,Q)*EPS2 )
|
||||
@@ -184,6 +200,15 @@
|
||||
*> RESULT(9) = 0 if THETA is in increasing order and
|
||||
*> all angles are in [0,pi/2];
|
||||
*> = ULPINV otherwise.
|
||||
*> Then, the 2-by-1 CSD:
|
||||
*> RESULT(10) = norm( U1'*X11*V1 - D11 ) / ( MAX(1,P,Q)*EPS2 )
|
||||
*> RESULT(11) = norm( U2'*X21*V1 - D21 ) / ( MAX(1,M-P,Q)*EPS2 )
|
||||
*> RESULT(12) = norm( I - U1'*U1 ) / ( MAX(1,P)*ULP )
|
||||
*> RESULT(13) = norm( I - U2'*U2 ) / ( MAX(1,M-P)*ULP )
|
||||
*> RESULT(14) = norm( I - V1T'*V1T ) / ( MAX(1,Q)*ULP )
|
||||
*> RESULT(15) = 0 if THETA is in increasing order and
|
||||
*> all angles are in [0,pi/2];
|
||||
*> = ULPINV otherwise.
|
||||
*> ( EPS2 = MAX( norm( I - X'*X ) / M, ULP ). )
|
||||
*> \endverbatim
|
||||
*
|
||||
@@ -214,7 +239,7 @@
|
||||
* ..
|
||||
* .. Array Arguments ..
|
||||
INTEGER IWORK( * )
|
||||
REAL RESULT( 9 ), RWORK( * ), THETA( * )
|
||||
REAL RESULT( 15 ), RWORK( * ), THETA( * )
|
||||
COMPLEX U1( LDU1, * ), U2( LDU2, * ), V1T( LDV1T, * ),
|
||||
$ V2T( LDV2T, * ), WORK( LWORK ), X( LDX, * ),
|
||||
$ XF( LDX, * )
|
||||
@@ -238,15 +263,19 @@
|
||||
EXTERNAL SLAMCH, CLANGE, CLANHE
|
||||
* ..
|
||||
* .. External Subroutines ..
|
||||
EXTERNAL CGEMM, CLACPY, CLASET, CUNCSD, CHERK
|
||||
EXTERNAL CGEMM, CHERK, CLACPY, CLASET, CUNCSD,
|
||||
$ CUNCSD2BY1
|
||||
* ..
|
||||
* .. Intrinsic Functions ..
|
||||
INTRINSIC REAL, MAX, MIN
|
||||
INTRINSIC CMPLX, COS, MAX, MIN, REAL, SIN
|
||||
* ..
|
||||
* .. Executable Statements ..
|
||||
*
|
||||
ULP = SLAMCH( 'Precision' )
|
||||
ULPINV = REALONE / ULP
|
||||
*
|
||||
* The first half of the routine checks the 2-by-2 CSD
|
||||
*
|
||||
CALL CLASET( 'Full', M, M, ZERO, ONE, WORK, LDX )
|
||||
CALL CHERK( 'Upper', 'Conjugate transpose', M, M, -REALONE,
|
||||
$ X, LDX, REALONE, WORK, LDX )
|
||||
@@ -269,86 +298,88 @@
|
||||
$ THETA, U1, LDU1, U2, LDU2, V1T, LDV1T, V2T, LDV2T,
|
||||
$ WORK, LWORK, RWORK, 17*(R+2), IWORK, INFO )
|
||||
*
|
||||
* Compute X := diag(U1,U2)'*X*diag(V1,V2) - [D11 D12; D21 D22]
|
||||
* Compute XF := diag(U1,U2)'*X*diag(V1,V2) - [D11 D12; D21 D22]
|
||||
*
|
||||
CALL CLACPY( 'Full', M, M, X, LDX, XF, LDX )
|
||||
*
|
||||
CALL CGEMM( 'No transpose', 'Conjugate transpose', P, Q, Q, ONE,
|
||||
$ X, LDX, V1T, LDV1T, ZERO, WORK, LDX )
|
||||
$ XF, LDX, V1T, LDV1T, ZERO, WORK, LDX )
|
||||
*
|
||||
CALL CGEMM( 'Conjugate transpose', 'No transpose', P, Q, P, ONE,
|
||||
$ U1, LDU1, WORK, LDX, ZERO, X, LDX )
|
||||
$ U1, LDU1, WORK, LDX, ZERO, XF, LDX )
|
||||
*
|
||||
DO I = 1, MIN(P,Q)-R
|
||||
X(I,I) = X(I,I) - ONE
|
||||
XF(I,I) = XF(I,I) - ONE
|
||||
END DO
|
||||
DO I = 1, R
|
||||
X(MIN(P,Q)-R+I,MIN(P,Q)-R+I) =
|
||||
$ X(MIN(P,Q)-R+I,MIN(P,Q)-R+I) - CMPLX( COS(THETA(I)),
|
||||
XF(MIN(P,Q)-R+I,MIN(P,Q)-R+I) =
|
||||
$ XF(MIN(P,Q)-R+I,MIN(P,Q)-R+I) - CMPLX( COS(THETA(I)),
|
||||
$ 0.0E0 )
|
||||
END DO
|
||||
*
|
||||
CALL CGEMM( 'No transpose', 'Conjugate transpose', P, M-Q, M-Q,
|
||||
$ ONE, X(1,Q+1), LDX, V2T, LDV2T, ZERO, WORK, LDX )
|
||||
$ ONE, XF(1,Q+1), LDX, V2T, LDV2T, ZERO, WORK, LDX )
|
||||
*
|
||||
CALL CGEMM( 'Conjugate transpose', 'No transpose', P, M-Q, P,
|
||||
$ ONE, U1, LDU1, WORK, LDX, ZERO, X(1,Q+1), LDX )
|
||||
$ ONE, U1, LDU1, WORK, LDX, ZERO, XF(1,Q+1), LDX )
|
||||
*
|
||||
DO I = 1, MIN(P,M-Q)-R
|
||||
X(P-I+1,M-I+1) = X(P-I+1,M-I+1) + ONE
|
||||
XF(P-I+1,M-I+1) = XF(P-I+1,M-I+1) + ONE
|
||||
END DO
|
||||
DO I = 1, R
|
||||
X(P-(MIN(P,M-Q)-R)+1-I,M-(MIN(P,M-Q)-R)+1-I) =
|
||||
$ X(P-(MIN(P,M-Q)-R)+1-I,M-(MIN(P,M-Q)-R)+1-I) +
|
||||
XF(P-(MIN(P,M-Q)-R)+1-I,M-(MIN(P,M-Q)-R)+1-I) =
|
||||
$ XF(P-(MIN(P,M-Q)-R)+1-I,M-(MIN(P,M-Q)-R)+1-I) +
|
||||
$ CMPLX( SIN(THETA(R-I+1)), 0.0E0 )
|
||||
END DO
|
||||
*
|
||||
CALL CGEMM( 'No transpose', 'Conjugate transpose', M-P, Q, Q, ONE,
|
||||
$ X(P+1,1), LDX, V1T, LDV1T, ZERO, WORK, LDX )
|
||||
$ XF(P+1,1), LDX, V1T, LDV1T, ZERO, WORK, LDX )
|
||||
*
|
||||
CALL CGEMM( 'Conjugate transpose', 'No transpose', M-P, Q, M-P,
|
||||
$ ONE, U2, LDU2, WORK, LDX, ZERO, X(P+1,1), LDX )
|
||||
$ ONE, U2, LDU2, WORK, LDX, ZERO, XF(P+1,1), LDX )
|
||||
*
|
||||
DO I = 1, MIN(M-P,Q)-R
|
||||
X(M-I+1,Q-I+1) = X(M-I+1,Q-I+1) - ONE
|
||||
XF(M-I+1,Q-I+1) = XF(M-I+1,Q-I+1) - ONE
|
||||
END DO
|
||||
DO I = 1, R
|
||||
X(M-(MIN(M-P,Q)-R)+1-I,Q-(MIN(M-P,Q)-R)+1-I) =
|
||||
$ X(M-(MIN(M-P,Q)-R)+1-I,Q-(MIN(M-P,Q)-R)+1-I) -
|
||||
XF(M-(MIN(M-P,Q)-R)+1-I,Q-(MIN(M-P,Q)-R)+1-I) =
|
||||
$ XF(M-(MIN(M-P,Q)-R)+1-I,Q-(MIN(M-P,Q)-R)+1-I) -
|
||||
$ CMPLX( SIN(THETA(R-I+1)), 0.0E0 )
|
||||
END DO
|
||||
*
|
||||
CALL CGEMM( 'No transpose', 'Conjugate transpose', M-P, M-Q, M-Q,
|
||||
$ ONE, X(P+1,Q+1), LDX, V2T, LDV2T, ZERO, WORK, LDX )
|
||||
$ ONE, XF(P+1,Q+1), LDX, V2T, LDV2T, ZERO, WORK, LDX )
|
||||
*
|
||||
CALL CGEMM( 'Conjugate transpose', 'No transpose', M-P, M-Q, M-P,
|
||||
$ ONE, U2, LDU2, WORK, LDX, ZERO, X(P+1,Q+1), LDX )
|
||||
$ ONE, U2, LDU2, WORK, LDX, ZERO, XF(P+1,Q+1), LDX )
|
||||
*
|
||||
DO I = 1, MIN(M-P,M-Q)-R
|
||||
X(P+I,Q+I) = X(P+I,Q+I) - ONE
|
||||
XF(P+I,Q+I) = XF(P+I,Q+I) - ONE
|
||||
END DO
|
||||
DO I = 1, R
|
||||
X(P+(MIN(M-P,M-Q)-R)+I,Q+(MIN(M-P,M-Q)-R)+I) =
|
||||
$ X(P+(MIN(M-P,M-Q)-R)+I,Q+(MIN(M-P,M-Q)-R)+I) -
|
||||
XF(P+(MIN(M-P,M-Q)-R)+I,Q+(MIN(M-P,M-Q)-R)+I) =
|
||||
$ XF(P+(MIN(M-P,M-Q)-R)+I,Q+(MIN(M-P,M-Q)-R)+I) -
|
||||
$ CMPLX( COS(THETA(I)), 0.0E0 )
|
||||
END DO
|
||||
*
|
||||
* Compute norm( U1'*X11*V1 - D11 ) / ( MAX(1,P,Q)*EPS2 ) .
|
||||
*
|
||||
RESID = CLANGE( '1', P, Q, X, LDX, RWORK )
|
||||
RESID = CLANGE( '1', P, Q, XF, LDX, RWORK )
|
||||
RESULT( 1 ) = ( RESID / REAL(MAX(1,P,Q)) ) / EPS2
|
||||
*
|
||||
* Compute norm( U1'*X12*V2 - D12 ) / ( MAX(1,P,M-Q)*EPS2 ) .
|
||||
*
|
||||
RESID = CLANGE( '1', P, M-Q, X(1,Q+1), LDX, RWORK )
|
||||
RESID = CLANGE( '1', P, M-Q, XF(1,Q+1), LDX, RWORK )
|
||||
RESULT( 2 ) = ( RESID / REAL(MAX(1,P,M-Q)) ) / EPS2
|
||||
*
|
||||
* Compute norm( U2'*X21*V1 - D21 ) / ( MAX(1,M-P,Q)*EPS2 ) .
|
||||
*
|
||||
RESID = CLANGE( '1', M-P, Q, X(P+1,1), LDX, RWORK )
|
||||
RESID = CLANGE( '1', M-P, Q, XF(P+1,1), LDX, RWORK )
|
||||
RESULT( 3 ) = ( RESID / REAL(MAX(1,M-P,Q)) ) / EPS2
|
||||
*
|
||||
* Compute norm( U2'*X22*V2 - D22 ) / ( MAX(1,M-P,M-Q)*EPS2 ) .
|
||||
*
|
||||
RESID = CLANGE( '1', M-P, M-Q, X(P+1,Q+1), LDX, RWORK )
|
||||
RESID = CLANGE( '1', M-P, M-Q, XF(P+1,Q+1), LDX, RWORK )
|
||||
RESULT( 4 ) = ( RESID / REAL(MAX(1,M-P,M-Q)) ) / EPS2
|
||||
*
|
||||
* Compute I - U1'*U1
|
||||
@@ -397,14 +428,126 @@
|
||||
*
|
||||
* Check sorting
|
||||
*
|
||||
RESULT(9) = REALZERO
|
||||
RESULT( 9 ) = REALZERO
|
||||
DO I = 1, R
|
||||
IF( THETA(I).LT.REALZERO .OR. THETA(I).GT.PIOVER2 ) THEN
|
||||
RESULT(9) = ULPINV
|
||||
RESULT( 9 ) = ULPINV
|
||||
END IF
|
||||
IF( I.GT.1) THEN
|
||||
IF ( THETA(I).LT.THETA(I-1) ) THEN
|
||||
RESULT(9) = ULPINV
|
||||
RESULT( 9 ) = ULPINV
|
||||
END IF
|
||||
END IF
|
||||
END DO
|
||||
*
|
||||
* The second half of the routine checks the 2-by-1 CSD
|
||||
*
|
||||
CALL CLASET( 'Full', Q, Q, ZERO, ONE, WORK, LDX )
|
||||
CALL CHERK( 'Upper', 'Conjugate transpose', Q, M, -REALONE,
|
||||
$ X, LDX, REALONE, WORK, LDX )
|
||||
IF (M.GT.0) THEN
|
||||
EPS2 = MAX( ULP,
|
||||
$ CLANGE( '1', Q, Q, WORK, LDX, RWORK ) / REAL( M ) )
|
||||
ELSE
|
||||
EPS2 = ULP
|
||||
END IF
|
||||
R = MIN( P, M-P, Q, M-Q )
|
||||
*
|
||||
* Copy the matrix X to the array XF.
|
||||
*
|
||||
CALL CLACPY( 'Full', M, Q, X, LDX, XF, LDX )
|
||||
*
|
||||
* Compute the CSD
|
||||
*
|
||||
CALL CUNCSD2BY1( 'Y', 'Y', 'Y', M, P, Q, XF(1,1), LDX, XF(P+1,1),
|
||||
$ LDX, THETA, U1, LDU1, U2, LDU2, V1T, LDV1T, WORK,
|
||||
$ LWORK, RWORK, 17*(R+2), IWORK, INFO )
|
||||
*
|
||||
* Compute [X11;X21] := diag(U1,U2)'*[X11;X21]*V1 - [D11;D21]
|
||||
*
|
||||
CALL CGEMM( 'No transpose', 'Conjugate transpose', P, Q, Q, ONE,
|
||||
$ X, LDX, V1T, LDV1T, ZERO, WORK, LDX )
|
||||
*
|
||||
CALL CGEMM( 'Conjugate transpose', 'No transpose', P, Q, P, ONE,
|
||||
$ U1, LDU1, WORK, LDX, ZERO, X, LDX )
|
||||
*
|
||||
DO I = 1, MIN(P,Q)-R
|
||||
X(I,I) = X(I,I) - ONE
|
||||
END DO
|
||||
DO I = 1, R
|
||||
X(MIN(P,Q)-R+I,MIN(P,Q)-R+I) =
|
||||
$ X(MIN(P,Q)-R+I,MIN(P,Q)-R+I) - CMPLX( COS(THETA(I)),
|
||||
$ 0.0E0 )
|
||||
END DO
|
||||
*
|
||||
CALL CGEMM( 'No transpose', 'Conjugate transpose', M-P, Q, Q, ONE,
|
||||
$ X(P+1,1), LDX, V1T, LDV1T, ZERO, WORK, LDX )
|
||||
*
|
||||
CALL CGEMM( 'Conjugate transpose', 'No transpose', M-P, Q, M-P,
|
||||
$ ONE, U2, LDU2, WORK, LDX, ZERO, X(P+1,1), LDX )
|
||||
*
|
||||
DO I = 1, MIN(M-P,Q)-R
|
||||
X(M-I+1,Q-I+1) = X(M-I+1,Q-I+1) - ONE
|
||||
END DO
|
||||
DO I = 1, R
|
||||
X(M-(MIN(M-P,Q)-R)+1-I,Q-(MIN(M-P,Q)-R)+1-I) =
|
||||
$ X(M-(MIN(M-P,Q)-R)+1-I,Q-(MIN(M-P,Q)-R)+1-I) -
|
||||
$ CMPLX( SIN(THETA(R-I+1)), 0.0E0 )
|
||||
END DO
|
||||
*
|
||||
* Compute norm( U1'*X11*V1 - D11 ) / ( MAX(1,P,Q)*EPS2 ) .
|
||||
*
|
||||
RESID = CLANGE( '1', P, Q, X, LDX, RWORK )
|
||||
RESULT( 10 ) = ( RESID / REAL(MAX(1,P,Q)) ) / EPS2
|
||||
*
|
||||
* Compute norm( U2'*X21*V1 - D21 ) / ( MAX(1,M-P,Q)*EPS2 ) .
|
||||
*
|
||||
RESID = CLANGE( '1', M-P, Q, X(P+1,1), LDX, RWORK )
|
||||
RESULT( 11 ) = ( RESID / REAL(MAX(1,M-P,Q)) ) / EPS2
|
||||
*
|
||||
* Compute I - U1'*U1
|
||||
*
|
||||
CALL CLASET( 'Full', P, P, ZERO, ONE, WORK, LDU1 )
|
||||
CALL CHERK( 'Upper', 'Conjugate transpose', P, P, -REALONE,
|
||||
$ U1, LDU1, REALONE, WORK, LDU1 )
|
||||
*
|
||||
* Compute norm( I - U1'*U1 ) / ( MAX(1,P) * ULP ) .
|
||||
*
|
||||
RESID = CLANHE( '1', 'Upper', P, WORK, LDU1, RWORK )
|
||||
RESULT( 12 ) = ( RESID / REAL(MAX(1,P)) ) / ULP
|
||||
*
|
||||
* Compute I - U2'*U2
|
||||
*
|
||||
CALL CLASET( 'Full', M-P, M-P, ZERO, ONE, WORK, LDU2 )
|
||||
CALL CHERK( 'Upper', 'Conjugate transpose', M-P, M-P, -REALONE,
|
||||
$ U2, LDU2, REALONE, WORK, LDU2 )
|
||||
*
|
||||
* Compute norm( I - U2'*U2 ) / ( MAX(1,M-P) * ULP ) .
|
||||
*
|
||||
RESID = CLANHE( '1', 'Upper', M-P, WORK, LDU2, RWORK )
|
||||
RESULT( 13 ) = ( RESID / REAL(MAX(1,M-P)) ) / ULP
|
||||
*
|
||||
* Compute I - V1T*V1T'
|
||||
*
|
||||
CALL CLASET( 'Full', Q, Q, ZERO, ONE, WORK, LDV1T )
|
||||
CALL CHERK( 'Upper', 'No transpose', Q, Q, -REALONE,
|
||||
$ V1T, LDV1T, REALONE, WORK, LDV1T )
|
||||
*
|
||||
* Compute norm( I - V1T*V1T' ) / ( MAX(1,Q) * ULP ) .
|
||||
*
|
||||
RESID = CLANHE( '1', 'Upper', Q, WORK, LDV1T, RWORK )
|
||||
RESULT( 14 ) = ( RESID / REAL(MAX(1,Q)) ) / ULP
|
||||
*
|
||||
* Check sorting
|
||||
*
|
||||
RESULT( 15 ) = REALZERO
|
||||
DO I = 1, R
|
||||
IF( THETA(I).LT.REALZERO .OR. THETA(I).GT.PIOVER2 ) THEN
|
||||
RESULT( 15 ) = ULPINV
|
||||
END IF
|
||||
IF( I.GT.1) THEN
|
||||
IF ( THETA(I).LT.THETA(I-1) ) THEN
|
||||
RESULT( 15 ) = ULPINV
|
||||
END IF
|
||||
END IF
|
||||
END DO
|
||||
@@ -414,4 +557,3 @@
|
||||
* End of CCSDTS
|
||||
*
|
||||
END
|
||||
|
||||
|
||||
@@ -1032,17 +1032,17 @@
|
||||
*> \author Univ. of Colorado Denver
|
||||
*> \author NAG Ltd.
|
||||
*
|
||||
*> \date April 2012
|
||||
*> \date November 2013
|
||||
*
|
||||
*> \ingroup double_eig
|
||||
*
|
||||
* =====================================================================
|
||||
PROGRAM DCHKEE
|
||||
*
|
||||
* -- LAPACK test routine (version 3.4.1) --
|
||||
* -- LAPACK test routine (version 3.5.0) --
|
||||
* -- LAPACK is a software package provided by Univ. of Tennessee, --
|
||||
* -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
|
||||
* April 2012
|
||||
* November 2013
|
||||
*
|
||||
* =====================================================================
|
||||
*
|
||||
@@ -1133,6 +1133,10 @@
|
||||
* ..
|
||||
* .. Executable Statements ..
|
||||
*
|
||||
A = 0.0
|
||||
B = 0.0
|
||||
C = 0.0
|
||||
D = 0.0
|
||||
S1 = DSECND( )
|
||||
FATAL = .FALSE.
|
||||
NUNIT = NOUT
|
||||
|
||||
@@ -205,13 +205,14 @@
|
||||
*
|
||||
* .. Parameters ..
|
||||
INTEGER NTESTS
|
||||
PARAMETER ( NTESTS = 9 )
|
||||
PARAMETER ( NTESTS = 15 )
|
||||
INTEGER NTYPES
|
||||
PARAMETER ( NTYPES = 3 )
|
||||
DOUBLE PRECISION GAPDIGIT, ORTH, PIOVER2, TEN
|
||||
PARAMETER ( GAPDIGIT = 18.0D0, ORTH = 1.0D-12,
|
||||
PARAMETER ( NTYPES = 4 )
|
||||
DOUBLE PRECISION GAPDIGIT, ONE, ORTH, PIOVER2, TEN, ZERO
|
||||
PARAMETER ( GAPDIGIT = 18.0D0, ONE = 1.0D0,
|
||||
$ ORTH = 1.0D-12,
|
||||
$ PIOVER2 = 1.57079632679489662D0,
|
||||
$ TEN = 10.0D0 )
|
||||
$ TEN = 10.0D0, ZERO = 0.0D0 )
|
||||
* ..
|
||||
* .. Local Scalars ..
|
||||
LOGICAL FIRSTT
|
||||
@@ -231,8 +232,8 @@
|
||||
INTRINSIC ABS, MIN
|
||||
* ..
|
||||
* .. External Functions ..
|
||||
DOUBLE PRECISION DLARND
|
||||
EXTERNAL DLARND
|
||||
DOUBLE PRECISION DLARAN, DLARND
|
||||
EXTERNAL DLARAN, DLARND
|
||||
* ..
|
||||
* .. Executable Statements ..
|
||||
*
|
||||
@@ -286,7 +287,7 @@
|
||||
$ ORTH*DLARND(2,ISEED)
|
||||
END DO
|
||||
END DO
|
||||
ELSE
|
||||
ELSE IF( IMAT.EQ.3 ) THEN
|
||||
R = MIN( P, M-P, Q, M-Q )
|
||||
DO I = 1, R+1
|
||||
THETA(I) = TEN**(-DLARND(1,ISEED)*GAPDIGIT)
|
||||
@@ -298,9 +299,18 @@
|
||||
THETA(I) = PIOVER2 * THETA(I) / THETA(R+1)
|
||||
END DO
|
||||
CALL DLACSG( M, P, Q, THETA, ISEED, X, LDX, WORK )
|
||||
ELSE
|
||||
CALL DLASET( 'F', M, M, ZERO, ONE, X, LDX )
|
||||
DO I = 1, M
|
||||
J = INT( DLARAN( ISEED ) * M ) + 1
|
||||
IF( J .NE. I ) THEN
|
||||
CALL DROT( M, X(1+(I-1)*LDX), 1, X(1+(J-1)*LDX), 1,
|
||||
$ ZERO, ONE )
|
||||
END IF
|
||||
END DO
|
||||
END IF
|
||||
*
|
||||
NT = 9
|
||||
NT = 15
|
||||
*
|
||||
CALL DCSDTS( M, P, Q, X, XF, LDX, U1, LDU1, U2, LDU2, V1T,
|
||||
$ LDV1T, V2T, LDV2T, THETA, IWORK, WORK, LWORK,
|
||||
|
||||
@@ -17,7 +17,7 @@
|
||||
* ..
|
||||
* .. Array Arguments ..
|
||||
* INTEGER IWORK( * )
|
||||
* DOUBLE PRECISION RESULT( 9 ), RWORK( * ), THETA( * )
|
||||
* DOUBLE PRECISION RESULT( 15 ), RWORK( * ), THETA( * )
|
||||
* DOUBLE PRECISION U1( LDU1, * ), U2( LDU2, * ), V1T( LDV1T, * ),
|
||||
* $ V2T( LDV2T, * ), WORK( LWORK ), X( LDX, * ),
|
||||
* $ XF( LDX, * )
|
||||
@@ -43,10 +43,25 @@
|
||||
*> [ I 0 0 | 0 0 0 ]
|
||||
*> [ 0 C 0 | 0 -S 0 ]
|
||||
*> [ 0 0 0 | 0 0 -I ]
|
||||
*> = [---------------------] = [ D11 D12 ] .
|
||||
*> = [---------------------] = [ D11 D12 ] ,
|
||||
*> [ 0 0 0 | I 0 0 ] [ D21 D22 ]
|
||||
*> [ 0 S 0 | 0 C 0 ]
|
||||
*> [ 0 0 I | 0 0 0 ]
|
||||
*>
|
||||
*> and also DORCSD2BY1, which, given
|
||||
*> Q
|
||||
*> [ X11 ] P ,
|
||||
*> [ X21 ] M-P
|
||||
*>
|
||||
*> computes the 2-by-1 CSD
|
||||
*>
|
||||
*> [ I 0 0 ]
|
||||
*> [ 0 C 0 ]
|
||||
*> [ 0 0 0 ]
|
||||
*> [ U1 ]**T * [ X11 ] * V1 = [----------] = [ D11 ] ,
|
||||
*> [ U2 ] [ X21 ] [ 0 0 0 ] [ D21 ]
|
||||
*> [ 0 S 0 ]
|
||||
*> [ 0 0 I ]
|
||||
*> \endverbatim
|
||||
*
|
||||
* Arguments:
|
||||
@@ -171,8 +186,9 @@
|
||||
*>
|
||||
*> \param[out] RESULT
|
||||
*> \verbatim
|
||||
*> RESULT is DOUBLE PRECISION array, dimension (9)
|
||||
*> RESULT is DOUBLE PRECISION array, dimension (15)
|
||||
*> The test ratios:
|
||||
*> First, the 2-by-2 CSD:
|
||||
*> RESULT(1) = norm( U1'*X11*V1 - D11 ) / ( MAX(1,P,Q)*EPS2 )
|
||||
*> RESULT(2) = norm( U1'*X12*V2 - D12 ) / ( MAX(1,P,M-Q)*EPS2 )
|
||||
*> RESULT(3) = norm( U2'*X21*V1 - D21 ) / ( MAX(1,M-P,Q)*EPS2 )
|
||||
@@ -184,6 +200,15 @@
|
||||
*> RESULT(9) = 0 if THETA is in increasing order and
|
||||
*> all angles are in [0,pi/2];
|
||||
*> = ULPINV otherwise.
|
||||
*> Then, the 2-by-1 CSD:
|
||||
*> RESULT(10) = norm( U1'*X11*V1 - D11 ) / ( MAX(1,P,Q)*EPS2 )
|
||||
*> RESULT(11) = norm( U2'*X21*V1 - D21 ) / ( MAX(1,M-P,Q)*EPS2 )
|
||||
*> RESULT(12) = norm( I - U1'*U1 ) / ( MAX(1,P)*ULP )
|
||||
*> RESULT(13) = norm( I - U2'*U2 ) / ( MAX(1,M-P)*ULP )
|
||||
*> RESULT(14) = norm( I - V1T'*V1T ) / ( MAX(1,Q)*ULP )
|
||||
*> RESULT(15) = 0 if THETA is in increasing order and
|
||||
*> all angles are in [0,pi/2];
|
||||
*> = ULPINV otherwise.
|
||||
*> ( EPS2 = MAX( norm( I - X'*X ) / M, ULP ). )
|
||||
*> \endverbatim
|
||||
*
|
||||
@@ -214,7 +239,7 @@
|
||||
* ..
|
||||
* .. Array Arguments ..
|
||||
INTEGER IWORK( * )
|
||||
DOUBLE PRECISION RESULT( 9 ), RWORK( * ), THETA( * )
|
||||
DOUBLE PRECISION RESULT( 15 ), RWORK( * ), THETA( * )
|
||||
DOUBLE PRECISION U1( LDU1, * ), U2( LDU2, * ), V1T( LDV1T, * ),
|
||||
$ V2T( LDV2T, * ), WORK( LWORK ), X( LDX, * ),
|
||||
$ XF( LDX, * )
|
||||
@@ -238,15 +263,19 @@
|
||||
EXTERNAL DLAMCH, DLANGE, DLANSY
|
||||
* ..
|
||||
* .. External Subroutines ..
|
||||
EXTERNAL DGEMM, DLACPY, DLASET, DORCSD, DSYRK
|
||||
EXTERNAL DGEMM, DLACPY, DLASET, DORCSD, DORCSD2BY1,
|
||||
$ DSYRK
|
||||
* ..
|
||||
* .. Intrinsic Functions ..
|
||||
INTRINSIC DBLE, MAX, MIN
|
||||
INTRINSIC COS, DBLE, MAX, MIN, SIN
|
||||
* ..
|
||||
* .. Executable Statements ..
|
||||
*
|
||||
ULP = DLAMCH( 'Precision' )
|
||||
ULPINV = REALONE / ULP
|
||||
*
|
||||
* The first half of the routine checks the 2-by-2 CSD
|
||||
*
|
||||
CALL DLASET( 'Full', M, M, ZERO, ONE, WORK, LDX )
|
||||
CALL DSYRK( 'Upper', 'Conjugate transpose', M, M, -ONE, X, LDX,
|
||||
$ ONE, WORK, LDX )
|
||||
@@ -269,85 +298,87 @@
|
||||
$ THETA, U1, LDU1, U2, LDU2, V1T, LDV1T, V2T, LDV2T,
|
||||
$ WORK, LWORK, IWORK, INFO )
|
||||
*
|
||||
* Compute X := diag(U1,U2)'*X*diag(V1,V2) - [D11 D12; D21 D22]
|
||||
* Compute XF := diag(U1,U2)'*X*diag(V1,V2) - [D11 D12; D21 D22]
|
||||
*
|
||||
CALL DLACPY( 'Full', M, M, X, LDX, XF, LDX )
|
||||
*
|
||||
CALL DGEMM( 'No transpose', 'Conjugate transpose', P, Q, Q, ONE,
|
||||
$ X, LDX, V1T, LDV1T, ZERO, WORK, LDX )
|
||||
$ XF, LDX, V1T, LDV1T, ZERO, WORK, LDX )
|
||||
*
|
||||
CALL DGEMM( 'Conjugate transpose', 'No transpose', P, Q, P, ONE,
|
||||
$ U1, LDU1, WORK, LDX, ZERO, X, LDX )
|
||||
$ U1, LDU1, WORK, LDX, ZERO, XF, LDX )
|
||||
*
|
||||
DO I = 1, MIN(P,Q)-R
|
||||
X(I,I) = X(I,I) - ONE
|
||||
XF(I,I) = XF(I,I) - ONE
|
||||
END DO
|
||||
DO I = 1, R
|
||||
X(MIN(P,Q)-R+I,MIN(P,Q)-R+I) =
|
||||
$ X(MIN(P,Q)-R+I,MIN(P,Q)-R+I) - COS(THETA(I))
|
||||
XF(MIN(P,Q)-R+I,MIN(P,Q)-R+I) =
|
||||
$ XF(MIN(P,Q)-R+I,MIN(P,Q)-R+I) - COS(THETA(I))
|
||||
END DO
|
||||
*
|
||||
CALL DGEMM( 'No transpose', 'Conjugate transpose', P, M-Q, M-Q,
|
||||
$ ONE, X(1,Q+1), LDX, V2T, LDV2T, ZERO, WORK, LDX )
|
||||
$ ONE, XF(1,Q+1), LDX, V2T, LDV2T, ZERO, WORK, LDX )
|
||||
*
|
||||
CALL DGEMM( 'Conjugate transpose', 'No transpose', P, M-Q, P,
|
||||
$ ONE, U1, LDU1, WORK, LDX, ZERO, X(1,Q+1), LDX )
|
||||
$ ONE, U1, LDU1, WORK, LDX, ZERO, XF(1,Q+1), LDX )
|
||||
*
|
||||
DO I = 1, MIN(P,M-Q)-R
|
||||
X(P-I+1,M-I+1) = X(P-I+1,M-I+1) + ONE
|
||||
XF(P-I+1,M-I+1) = XF(P-I+1,M-I+1) + ONE
|
||||
END DO
|
||||
DO I = 1, R
|
||||
X(P-(MIN(P,M-Q)-R)+1-I,M-(MIN(P,M-Q)-R)+1-I) =
|
||||
$ X(P-(MIN(P,M-Q)-R)+1-I,M-(MIN(P,M-Q)-R)+1-I) +
|
||||
XF(P-(MIN(P,M-Q)-R)+1-I,M-(MIN(P,M-Q)-R)+1-I) =
|
||||
$ XF(P-(MIN(P,M-Q)-R)+1-I,M-(MIN(P,M-Q)-R)+1-I) +
|
||||
$ SIN(THETA(R-I+1))
|
||||
END DO
|
||||
*
|
||||
CALL DGEMM( 'No transpose', 'Conjugate transpose', M-P, Q, Q, ONE,
|
||||
$ X(P+1,1), LDX, V1T, LDV1T, ZERO, WORK, LDX )
|
||||
$ XF(P+1,1), LDX, V1T, LDV1T, ZERO, WORK, LDX )
|
||||
*
|
||||
CALL DGEMM( 'Conjugate transpose', 'No transpose', M-P, Q, M-P,
|
||||
$ ONE, U2, LDU2, WORK, LDX, ZERO, X(P+1,1), LDX )
|
||||
$ ONE, U2, LDU2, WORK, LDX, ZERO, XF(P+1,1), LDX )
|
||||
*
|
||||
DO I = 1, MIN(M-P,Q)-R
|
||||
X(M-I+1,Q-I+1) = X(M-I+1,Q-I+1) - ONE
|
||||
XF(M-I+1,Q-I+1) = XF(M-I+1,Q-I+1) - ONE
|
||||
END DO
|
||||
DO I = 1, R
|
||||
X(M-(MIN(M-P,Q)-R)+1-I,Q-(MIN(M-P,Q)-R)+1-I) =
|
||||
$ X(M-(MIN(M-P,Q)-R)+1-I,Q-(MIN(M-P,Q)-R)+1-I) -
|
||||
XF(M-(MIN(M-P,Q)-R)+1-I,Q-(MIN(M-P,Q)-R)+1-I) =
|
||||
$ XF(M-(MIN(M-P,Q)-R)+1-I,Q-(MIN(M-P,Q)-R)+1-I) -
|
||||
$ SIN(THETA(R-I+1))
|
||||
END DO
|
||||
*
|
||||
CALL DGEMM( 'No transpose', 'Conjugate transpose', M-P, M-Q, M-Q,
|
||||
$ ONE, X(P+1,Q+1), LDX, V2T, LDV2T, ZERO, WORK, LDX )
|
||||
$ ONE, XF(P+1,Q+1), LDX, V2T, LDV2T, ZERO, WORK, LDX )
|
||||
*
|
||||
CALL DGEMM( 'Conjugate transpose', 'No transpose', M-P, M-Q, M-P,
|
||||
$ ONE, U2, LDU2, WORK, LDX, ZERO, X(P+1,Q+1), LDX )
|
||||
$ ONE, U2, LDU2, WORK, LDX, ZERO, XF(P+1,Q+1), LDX )
|
||||
*
|
||||
DO I = 1, MIN(M-P,M-Q)-R
|
||||
X(P+I,Q+I) = X(P+I,Q+I) - ONE
|
||||
XF(P+I,Q+I) = XF(P+I,Q+I) - ONE
|
||||
END DO
|
||||
DO I = 1, R
|
||||
X(P+(MIN(M-P,M-Q)-R)+I,Q+(MIN(M-P,M-Q)-R)+I) =
|
||||
$ X(P+(MIN(M-P,M-Q)-R)+I,Q+(MIN(M-P,M-Q)-R)+I) -
|
||||
XF(P+(MIN(M-P,M-Q)-R)+I,Q+(MIN(M-P,M-Q)-R)+I) =
|
||||
$ XF(P+(MIN(M-P,M-Q)-R)+I,Q+(MIN(M-P,M-Q)-R)+I) -
|
||||
$ COS(THETA(I))
|
||||
END DO
|
||||
*
|
||||
* Compute norm( U1'*X11*V1 - D11 ) / ( MAX(1,P,Q)*EPS2 ) .
|
||||
*
|
||||
RESID = DLANGE( '1', P, Q, X, LDX, RWORK )
|
||||
RESID = DLANGE( '1', P, Q, XF, LDX, RWORK )
|
||||
RESULT( 1 ) = ( RESID / DBLE(MAX(1,P,Q)) ) / EPS2
|
||||
*
|
||||
* Compute norm( U1'*X12*V2 - D12 ) / ( MAX(1,P,M-Q)*EPS2 ) .
|
||||
*
|
||||
RESID = DLANGE( '1', P, M-Q, X(1,Q+1), LDX, RWORK )
|
||||
RESID = DLANGE( '1', P, M-Q, XF(1,Q+1), LDX, RWORK )
|
||||
RESULT( 2 ) = ( RESID / DBLE(MAX(1,P,M-Q)) ) / EPS2
|
||||
*
|
||||
* Compute norm( U2'*X21*V1 - D21 ) / ( MAX(1,M-P,Q)*EPS2 ) .
|
||||
*
|
||||
RESID = DLANGE( '1', M-P, Q, X(P+1,1), LDX, RWORK )
|
||||
RESID = DLANGE( '1', M-P, Q, XF(P+1,1), LDX, RWORK )
|
||||
RESULT( 3 ) = ( RESID / DBLE(MAX(1,M-P,Q)) ) / EPS2
|
||||
*
|
||||
* Compute norm( U2'*X22*V2 - D22 ) / ( MAX(1,M-P,M-Q)*EPS2 ) .
|
||||
*
|
||||
RESID = DLANGE( '1', M-P, M-Q, X(P+1,Q+1), LDX, RWORK )
|
||||
RESID = DLANGE( '1', M-P, M-Q, XF(P+1,Q+1), LDX, RWORK )
|
||||
RESULT( 4 ) = ( RESID / DBLE(MAX(1,M-P,M-Q)) ) / EPS2
|
||||
*
|
||||
* Compute I - U1'*U1
|
||||
@@ -396,14 +427,125 @@
|
||||
*
|
||||
* Check sorting
|
||||
*
|
||||
RESULT(9) = REALZERO
|
||||
RESULT( 9 ) = REALZERO
|
||||
DO I = 1, R
|
||||
IF( THETA(I).LT.REALZERO .OR. THETA(I).GT.PIOVER2 ) THEN
|
||||
RESULT(9) = ULPINV
|
||||
RESULT( 9 ) = ULPINV
|
||||
END IF
|
||||
IF( I.GT.1) THEN
|
||||
IF( I.GT.1 ) THEN
|
||||
IF ( THETA(I).LT.THETA(I-1) ) THEN
|
||||
RESULT(9) = ULPINV
|
||||
RESULT( 9 ) = ULPINV
|
||||
END IF
|
||||
END IF
|
||||
END DO
|
||||
*
|
||||
* The second half of the routine checks the 2-by-1 CSD
|
||||
*
|
||||
CALL DLASET( 'Full', Q, Q, ZERO, ONE, WORK, LDX )
|
||||
CALL DSYRK( 'Upper', 'Conjugate transpose', Q, M, -ONE, X, LDX,
|
||||
$ ONE, WORK, LDX )
|
||||
IF( M.GT.0 ) THEN
|
||||
EPS2 = MAX( ULP,
|
||||
$ DLANGE( '1', Q, Q, WORK, LDX, RWORK ) / DBLE( M ) )
|
||||
ELSE
|
||||
EPS2 = ULP
|
||||
END IF
|
||||
R = MIN( P, M-P, Q, M-Q )
|
||||
*
|
||||
* Copy the matrix [ X11; X21 ] to the array XF.
|
||||
*
|
||||
CALL DLACPY( 'Full', M, Q, X, LDX, XF, LDX )
|
||||
*
|
||||
* Compute the CSD
|
||||
*
|
||||
CALL DORCSD2BY1( 'Y', 'Y', 'Y', M, P, Q, XF(1,1), LDX, XF(P+1,1),
|
||||
$ LDX, THETA, U1, LDU1, U2, LDU2, V1T, LDV1T, WORK,
|
||||
$ LWORK, IWORK, INFO )
|
||||
*
|
||||
* Compute [X11;X21] := diag(U1,U2)'*[X11;X21]*V1 - [D11;D21]
|
||||
*
|
||||
CALL DGEMM( 'No transpose', 'Conjugate transpose', P, Q, Q, ONE,
|
||||
$ X, LDX, V1T, LDV1T, ZERO, WORK, LDX )
|
||||
*
|
||||
CALL DGEMM( 'Conjugate transpose', 'No transpose', P, Q, P, ONE,
|
||||
$ U1, LDU1, WORK, LDX, ZERO, X, LDX )
|
||||
*
|
||||
DO I = 1, MIN(P,Q)-R
|
||||
X(I,I) = X(I,I) - ONE
|
||||
END DO
|
||||
DO I = 1, R
|
||||
X(MIN(P,Q)-R+I,MIN(P,Q)-R+I) =
|
||||
$ X(MIN(P,Q)-R+I,MIN(P,Q)-R+I) - COS(THETA(I))
|
||||
END DO
|
||||
*
|
||||
CALL DGEMM( 'No transpose', 'Conjugate transpose', M-P, Q, Q, ONE,
|
||||
$ X(P+1,1), LDX, V1T, LDV1T, ZERO, WORK, LDX )
|
||||
*
|
||||
CALL DGEMM( 'Conjugate transpose', 'No transpose', M-P, Q, M-P,
|
||||
$ ONE, U2, LDU2, WORK, LDX, ZERO, X(P+1,1), LDX )
|
||||
*
|
||||
DO I = 1, MIN(M-P,Q)-R
|
||||
X(M-I+1,Q-I+1) = X(M-I+1,Q-I+1) - ONE
|
||||
END DO
|
||||
DO I = 1, R
|
||||
X(M-(MIN(M-P,Q)-R)+1-I,Q-(MIN(M-P,Q)-R)+1-I) =
|
||||
$ X(M-(MIN(M-P,Q)-R)+1-I,Q-(MIN(M-P,Q)-R)+1-I) -
|
||||
$ SIN(THETA(R-I+1))
|
||||
END DO
|
||||
*
|
||||
* Compute norm( U1'*X11*V1 - D11 ) / ( MAX(1,P,Q)*EPS2 ) .
|
||||
*
|
||||
RESID = DLANGE( '1', P, Q, X, LDX, RWORK )
|
||||
RESULT( 10 ) = ( RESID / DBLE(MAX(1,P,Q)) ) / EPS2
|
||||
*
|
||||
* Compute norm( U2'*X21*V1 - D21 ) / ( MAX(1,M-P,Q)*EPS2 ) .
|
||||
*
|
||||
RESID = DLANGE( '1', M-P, Q, X(P+1,1), LDX, RWORK )
|
||||
RESULT( 11 ) = ( RESID / DBLE(MAX(1,M-P,Q)) ) / EPS2
|
||||
*
|
||||
* Compute I - U1'*U1
|
||||
*
|
||||
CALL DLASET( 'Full', P, P, ZERO, ONE, WORK, LDU1 )
|
||||
CALL DSYRK( 'Upper', 'Conjugate transpose', P, P, -ONE, U1, LDU1,
|
||||
$ ONE, WORK, LDU1 )
|
||||
*
|
||||
* Compute norm( I - U1'*U1 ) / ( MAX(1,P) * ULP ) .
|
||||
*
|
||||
RESID = DLANSY( '1', 'Upper', P, WORK, LDU1, RWORK )
|
||||
RESULT( 12 ) = ( RESID / DBLE(MAX(1,P)) ) / ULP
|
||||
*
|
||||
* Compute I - U2'*U2
|
||||
*
|
||||
CALL DLASET( 'Full', M-P, M-P, ZERO, ONE, WORK, LDU2 )
|
||||
CALL DSYRK( 'Upper', 'Conjugate transpose', M-P, M-P, -ONE, U2,
|
||||
$ LDU2, ONE, WORK, LDU2 )
|
||||
*
|
||||
* Compute norm( I - U2'*U2 ) / ( MAX(1,M-P) * ULP ) .
|
||||
*
|
||||
RESID = DLANSY( '1', 'Upper', M-P, WORK, LDU2, RWORK )
|
||||
RESULT( 13 ) = ( RESID / DBLE(MAX(1,M-P)) ) / ULP
|
||||
*
|
||||
* Compute I - V1T*V1T'
|
||||
*
|
||||
CALL DLASET( 'Full', Q, Q, ZERO, ONE, WORK, LDV1T )
|
||||
CALL DSYRK( 'Upper', 'No transpose', Q, Q, -ONE, V1T, LDV1T, ONE,
|
||||
$ WORK, LDV1T )
|
||||
*
|
||||
* Compute norm( I - V1T*V1T' ) / ( MAX(1,Q) * ULP ) .
|
||||
*
|
||||
RESID = DLANSY( '1', 'Upper', Q, WORK, LDV1T, RWORK )
|
||||
RESULT( 14 ) = ( RESID / DBLE(MAX(1,Q)) ) / ULP
|
||||
*
|
||||
* Check sorting
|
||||
*
|
||||
RESULT( 15 ) = REALZERO
|
||||
DO I = 1, R
|
||||
IF( THETA(I).LT.REALZERO .OR. THETA(I).GT.PIOVER2 ) THEN
|
||||
RESULT( 15 ) = ULPINV
|
||||
END IF
|
||||
IF( I.GT.1 ) THEN
|
||||
IF ( THETA(I).LT.THETA(I-1) ) THEN
|
||||
RESULT( 15 ) = ULPINV
|
||||
END IF
|
||||
END IF
|
||||
END DO
|
||||
|
||||
@@ -58,17 +58,17 @@
|
||||
*> \author Univ. of Colorado Denver
|
||||
*> \author NAG Ltd.
|
||||
*
|
||||
*> \date November 2011
|
||||
*> \date November 2013
|
||||
*
|
||||
*> \ingroup double_eig
|
||||
*
|
||||
* =====================================================================
|
||||
SUBROUTINE DLAHD2( IOUNIT, PATH )
|
||||
*
|
||||
* -- LAPACK test routine (version 3.4.0) --
|
||||
* -- LAPACK test routine (version 3.5.0) --
|
||||
* -- LAPACK is a software package provided by Univ. of Tennessee, --
|
||||
* -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
|
||||
* November 2011
|
||||
* November 2013
|
||||
*
|
||||
* .. Scalar Arguments ..
|
||||
CHARACTER*3 PATH
|
||||
@@ -501,8 +501,8 @@
|
||||
$ / ' 2: norm( I - Q'' Q ) / ( m ulp )',
|
||||
$ / ' 3: norm( I - PT PT'' ) / ( n ulp )',
|
||||
$ / ' 4: norm( Y - Q'' C ) / ( norm(Y) max(m,nrhs) ulp )' )
|
||||
9968 FORMAT( / ' Tests performed: See sdrvst.f' )
|
||||
9967 FORMAT( / ' Tests performed: See cdrvst.f' )
|
||||
9968 FORMAT( / ' Tests performed: See ddrvst.f' )
|
||||
9967 FORMAT( / ' Tests performed: See sdrvst.f' )
|
||||
*
|
||||
* End of DLAHD2
|
||||
*
|
||||
|
||||
@@ -1032,17 +1032,17 @@
|
||||
*> \author Univ. of Colorado Denver
|
||||
*> \author NAG Ltd.
|
||||
*
|
||||
*> \date April 2012
|
||||
*> \date November 2013
|
||||
*
|
||||
*> \ingroup single_eig
|
||||
*
|
||||
* =====================================================================
|
||||
PROGRAM SCHKEE
|
||||
*
|
||||
* -- LAPACK test routine (version 3.4.1) --
|
||||
* -- LAPACK test routine (version 3.5.0) --
|
||||
* -- LAPACK is a software package provided by Univ. of Tennessee, --
|
||||
* -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
|
||||
* April 2012
|
||||
* November 2013
|
||||
*
|
||||
* =====================================================================
|
||||
*
|
||||
@@ -1133,6 +1133,10 @@
|
||||
* ..
|
||||
* .. Executable Statements ..
|
||||
*
|
||||
A = 0.0
|
||||
B = 0.0
|
||||
C = 0.0
|
||||
D = 0.0
|
||||
S1 = SECOND( )
|
||||
FATAL = .FALSE.
|
||||
NUNIT = NOUT
|
||||
|
||||
@@ -205,13 +205,14 @@
|
||||
*
|
||||
* .. Parameters ..
|
||||
INTEGER NTESTS
|
||||
PARAMETER ( NTESTS = 9 )
|
||||
PARAMETER ( NTESTS = 15 )
|
||||
INTEGER NTYPES
|
||||
PARAMETER ( NTYPES = 3 )
|
||||
REAL GAPDIGIT, ORTH, PIOVER2, TEN
|
||||
PARAMETER ( GAPDIGIT = 10.0E0, ORTH = 1.0E-4,
|
||||
PARAMETER ( NTYPES = 4 )
|
||||
REAL GAPDIGIT, ONE, ORTH, PIOVER2, TEN, ZERO
|
||||
PARAMETER ( GAPDIGIT = 10.0E0, ONE = 1.0E0,
|
||||
$ ORTH = 1.0E-4,
|
||||
$ PIOVER2 = 1.57079632679489662E0,
|
||||
$ TEN = 10.0D0 )
|
||||
$ TEN = 10.0E0, ZERO = 0.0E0 )
|
||||
* ..
|
||||
* .. Local Scalars ..
|
||||
LOGICAL FIRSTT
|
||||
@@ -231,8 +232,8 @@
|
||||
INTRINSIC ABS, MIN
|
||||
* ..
|
||||
* .. External Functions ..
|
||||
REAL SLARND
|
||||
EXTERNAL SLARND
|
||||
REAL SLARAN, SLARND
|
||||
EXTERNAL SLARAN, SLARND
|
||||
* ..
|
||||
* .. Executable Statements ..
|
||||
*
|
||||
@@ -286,7 +287,7 @@
|
||||
$ ORTH*SLARND(2,ISEED)
|
||||
END DO
|
||||
END DO
|
||||
ELSE
|
||||
ELSE IF( IMAT.EQ.3 ) THEN
|
||||
R = MIN( P, M-P, Q, M-Q )
|
||||
DO I = 1, R+1
|
||||
THETA(I) = TEN**(-SLARND(1,ISEED)*GAPDIGIT)
|
||||
@@ -298,9 +299,18 @@
|
||||
THETA(I) = PIOVER2 * THETA(I) / THETA(R+1)
|
||||
END DO
|
||||
CALL SLACSG( M, P, Q, THETA, ISEED, X, LDX, WORK )
|
||||
ELSE
|
||||
CALL SLASET( 'F', M, M, ZERO, ONE, X, LDX )
|
||||
DO I = 1, M
|
||||
J = INT( SLARAN( ISEED ) * M ) + 1
|
||||
IF( J .NE. I ) THEN
|
||||
CALL SROT( M, X(1+(I-1)*LDX), 1, X(1+(J-1)*LDX), 1,
|
||||
$ ZERO, ONE )
|
||||
END IF
|
||||
END DO
|
||||
END IF
|
||||
*
|
||||
NT = 9
|
||||
NT = 15
|
||||
*
|
||||
CALL SCSDTS( M, P, Q, X, XF, LDX, U1, LDU1, U2, LDU2, V1T,
|
||||
$ LDV1T, V2T, LDV2T, THETA, IWORK, WORK, LWORK,
|
||||
|
||||
@@ -17,7 +17,7 @@
|
||||
* ..
|
||||
* .. Array Arguments ..
|
||||
* INTEGER IWORK( * )
|
||||
* REAL RESULT( 9 ), RWORK( * ), THETA( * )
|
||||
* REAL RESULT( 15 ), RWORK( * ), THETA( * )
|
||||
* REAL U1( LDU1, * ), U2( LDU2, * ), V1T( LDV1T, * ),
|
||||
* $ V2T( LDV2T, * ), WORK( LWORK ), X( LDX, * ),
|
||||
* $ XF( LDX, * )
|
||||
@@ -47,6 +47,21 @@
|
||||
*> [ 0 0 0 | I 0 0 ] [ D21 D22 ]
|
||||
*> [ 0 S 0 | 0 C 0 ]
|
||||
*> [ 0 0 I | 0 0 0 ]
|
||||
*>
|
||||
*> and also SORCSD2BY1, which, given
|
||||
*> Q
|
||||
*> [ X11 ] P ,
|
||||
*> [ X21 ] M-P
|
||||
*>
|
||||
*> computes the 2-by-1 CSD
|
||||
*>
|
||||
*> [ I 0 0 ]
|
||||
*> [ 0 C 0 ]
|
||||
*> [ 0 0 0 ]
|
||||
*> [ U1 ]**T * [ X11 ] * V1 = [----------] = [ D11 ] ,
|
||||
*> [ U2 ] [ X21 ] [ 0 0 0 ] [ D21 ]
|
||||
*> [ 0 S 0 ]
|
||||
*> [ 0 0 I ]
|
||||
*> \endverbatim
|
||||
*
|
||||
* Arguments:
|
||||
@@ -171,8 +186,9 @@
|
||||
*>
|
||||
*> \param[out] RESULT
|
||||
*> \verbatim
|
||||
*> RESULT is REAL array, dimension (9)
|
||||
*> RESULT is REAL array, dimension (15)
|
||||
*> The test ratios:
|
||||
*> First, the 2-by-2 CSD:
|
||||
*> RESULT(1) = norm( U1'*X11*V1 - D11 ) / ( MAX(1,P,Q)*EPS2 )
|
||||
*> RESULT(2) = norm( U1'*X12*V2 - D12 ) / ( MAX(1,P,M-Q)*EPS2 )
|
||||
*> RESULT(3) = norm( U2'*X21*V1 - D21 ) / ( MAX(1,M-P,Q)*EPS2 )
|
||||
@@ -184,6 +200,15 @@
|
||||
*> RESULT(9) = 0 if THETA is in increasing order and
|
||||
*> all angles are in [0,pi/2];
|
||||
*> = ULPINV otherwise.
|
||||
*> Then, the 2-by-1 CSD:
|
||||
*> RESULT(10) = norm( U1'*X11*V1 - D11 ) / ( MAX(1,P,Q)*EPS2 )
|
||||
*> RESULT(11) = norm( U2'*X21*V1 - D21 ) / ( MAX(1,M-P,Q)*EPS2 )
|
||||
*> RESULT(12) = norm( I - U1'*U1 ) / ( MAX(1,P)*ULP )
|
||||
*> RESULT(13) = norm( I - U2'*U2 ) / ( MAX(1,M-P)*ULP )
|
||||
*> RESULT(14) = norm( I - V1T'*V1T ) / ( MAX(1,Q)*ULP )
|
||||
*> RESULT(15) = 0 if THETA is in increasing order and
|
||||
*> all angles are in [0,pi/2];
|
||||
*> = ULPINV otherwise.
|
||||
*> ( EPS2 = MAX( norm( I - X'*X ) / M, ULP ). )
|
||||
*> \endverbatim
|
||||
*
|
||||
@@ -214,7 +239,7 @@
|
||||
* ..
|
||||
* .. Array Arguments ..
|
||||
INTEGER IWORK( * )
|
||||
REAL RESULT( 9 ), RWORK( * ), THETA( * )
|
||||
REAL RESULT( 15 ), RWORK( * ), THETA( * )
|
||||
REAL U1( LDU1, * ), U2( LDU2, * ), V1T( LDV1T, * ),
|
||||
$ V2T( LDV2T, * ), WORK( LWORK ), X( LDX, * ),
|
||||
$ XF( LDX, * )
|
||||
@@ -238,15 +263,19 @@
|
||||
EXTERNAL SLAMCH, SLANGE, SLANSY
|
||||
* ..
|
||||
* .. External Subroutines ..
|
||||
EXTERNAL SGEMM, SLACPY, SLASET, SORCSD, SSYRK
|
||||
EXTERNAL SGEMM, SLACPY, SLASET, SORCSD, SORCSD2BY1,
|
||||
$ SSYRK
|
||||
* ..
|
||||
* .. Intrinsic Functions ..
|
||||
INTRINSIC REAL, MAX, MIN
|
||||
INTRINSIC COS, MAX, MIN, REAL, SIN
|
||||
* ..
|
||||
* .. Executable Statements ..
|
||||
*
|
||||
ULP = SLAMCH( 'Precision' )
|
||||
ULPINV = REALONE / ULP
|
||||
*
|
||||
* The first half of the routine checks the 2-by-2 CSD
|
||||
*
|
||||
CALL SLASET( 'Full', M, M, ZERO, ONE, WORK, LDX )
|
||||
CALL SSYRK( 'Upper', 'Conjugate transpose', M, M, -ONE, X, LDX,
|
||||
$ ONE, WORK, LDX )
|
||||
@@ -269,85 +298,87 @@
|
||||
$ THETA, U1, LDU1, U2, LDU2, V1T, LDV1T, V2T, LDV2T,
|
||||
$ WORK, LWORK, IWORK, INFO )
|
||||
*
|
||||
* Compute X := diag(U1,U2)'*X*diag(V1,V2) - [D11 D12; D21 D22]
|
||||
* Compute XF := diag(U1,U2)'*X*diag(V1,V2) - [D11 D12; D21 D22]
|
||||
*
|
||||
CALL SLACPY( 'Full', M, M, X, LDX, XF, LDX )
|
||||
*
|
||||
CALL SGEMM( 'No transpose', 'Conjugate transpose', P, Q, Q, ONE,
|
||||
$ X, LDX, V1T, LDV1T, ZERO, WORK, LDX )
|
||||
$ XF, LDX, V1T, LDV1T, ZERO, WORK, LDX )
|
||||
*
|
||||
CALL SGEMM( 'Conjugate transpose', 'No transpose', P, Q, P, ONE,
|
||||
$ U1, LDU1, WORK, LDX, ZERO, X, LDX )
|
||||
$ U1, LDU1, WORK, LDX, ZERO, XF, LDX )
|
||||
*
|
||||
DO I = 1, MIN(P,Q)-R
|
||||
X(I,I) = X(I,I) - ONE
|
||||
XF(I,I) = XF(I,I) - ONE
|
||||
END DO
|
||||
DO I = 1, R
|
||||
X(MIN(P,Q)-R+I,MIN(P,Q)-R+I) =
|
||||
$ X(MIN(P,Q)-R+I,MIN(P,Q)-R+I) - COS(THETA(I))
|
||||
XF(MIN(P,Q)-R+I,MIN(P,Q)-R+I) =
|
||||
$ XF(MIN(P,Q)-R+I,MIN(P,Q)-R+I) - COS(THETA(I))
|
||||
END DO
|
||||
*
|
||||
CALL SGEMM( 'No transpose', 'Conjugate transpose', P, M-Q, M-Q,
|
||||
$ ONE, X(1,Q+1), LDX, V2T, LDV2T, ZERO, WORK, LDX )
|
||||
$ ONE, XF(1,Q+1), LDX, V2T, LDV2T, ZERO, WORK, LDX )
|
||||
*
|
||||
CALL SGEMM( 'Conjugate transpose', 'No transpose', P, M-Q, P,
|
||||
$ ONE, U1, LDU1, WORK, LDX, ZERO, X(1,Q+1), LDX )
|
||||
$ ONE, U1, LDU1, WORK, LDX, ZERO, XF(1,Q+1), LDX )
|
||||
*
|
||||
DO I = 1, MIN(P,M-Q)-R
|
||||
X(P-I+1,M-I+1) = X(P-I+1,M-I+1) + ONE
|
||||
XF(P-I+1,M-I+1) = XF(P-I+1,M-I+1) + ONE
|
||||
END DO
|
||||
DO I = 1, R
|
||||
X(P-(MIN(P,M-Q)-R)+1-I,M-(MIN(P,M-Q)-R)+1-I) =
|
||||
$ X(P-(MIN(P,M-Q)-R)+1-I,M-(MIN(P,M-Q)-R)+1-I) +
|
||||
XF(P-(MIN(P,M-Q)-R)+1-I,M-(MIN(P,M-Q)-R)+1-I) =
|
||||
$ XF(P-(MIN(P,M-Q)-R)+1-I,M-(MIN(P,M-Q)-R)+1-I) +
|
||||
$ SIN(THETA(R-I+1))
|
||||
END DO
|
||||
*
|
||||
CALL SGEMM( 'No transpose', 'Conjugate transpose', M-P, Q, Q, ONE,
|
||||
$ X(P+1,1), LDX, V1T, LDV1T, ZERO, WORK, LDX )
|
||||
$ XF(P+1,1), LDX, V1T, LDV1T, ZERO, WORK, LDX )
|
||||
*
|
||||
CALL SGEMM( 'Conjugate transpose', 'No transpose', M-P, Q, M-P,
|
||||
$ ONE, U2, LDU2, WORK, LDX, ZERO, X(P+1,1), LDX )
|
||||
$ ONE, U2, LDU2, WORK, LDX, ZERO, XF(P+1,1), LDX )
|
||||
*
|
||||
DO I = 1, MIN(M-P,Q)-R
|
||||
X(M-I+1,Q-I+1) = X(M-I+1,Q-I+1) - ONE
|
||||
XF(M-I+1,Q-I+1) = XF(M-I+1,Q-I+1) - ONE
|
||||
END DO
|
||||
DO I = 1, R
|
||||
X(M-(MIN(M-P,Q)-R)+1-I,Q-(MIN(M-P,Q)-R)+1-I) =
|
||||
$ X(M-(MIN(M-P,Q)-R)+1-I,Q-(MIN(M-P,Q)-R)+1-I) -
|
||||
XF(M-(MIN(M-P,Q)-R)+1-I,Q-(MIN(M-P,Q)-R)+1-I) =
|
||||
$ XF(M-(MIN(M-P,Q)-R)+1-I,Q-(MIN(M-P,Q)-R)+1-I) -
|
||||
$ SIN(THETA(R-I+1))
|
||||
END DO
|
||||
*
|
||||
CALL SGEMM( 'No transpose', 'Conjugate transpose', M-P, M-Q, M-Q,
|
||||
$ ONE, X(P+1,Q+1), LDX, V2T, LDV2T, ZERO, WORK, LDX )
|
||||
$ ONE, XF(P+1,Q+1), LDX, V2T, LDV2T, ZERO, WORK, LDX )
|
||||
*
|
||||
CALL SGEMM( 'Conjugate transpose', 'No transpose', M-P, M-Q, M-P,
|
||||
$ ONE, U2, LDU2, WORK, LDX, ZERO, X(P+1,Q+1), LDX )
|
||||
$ ONE, U2, LDU2, WORK, LDX, ZERO, XF(P+1,Q+1), LDX )
|
||||
*
|
||||
DO I = 1, MIN(M-P,M-Q)-R
|
||||
X(P+I,Q+I) = X(P+I,Q+I) - ONE
|
||||
XF(P+I,Q+I) = XF(P+I,Q+I) - ONE
|
||||
END DO
|
||||
DO I = 1, R
|
||||
X(P+(MIN(M-P,M-Q)-R)+I,Q+(MIN(M-P,M-Q)-R)+I) =
|
||||
$ X(P+(MIN(M-P,M-Q)-R)+I,Q+(MIN(M-P,M-Q)-R)+I) -
|
||||
XF(P+(MIN(M-P,M-Q)-R)+I,Q+(MIN(M-P,M-Q)-R)+I) =
|
||||
$ XF(P+(MIN(M-P,M-Q)-R)+I,Q+(MIN(M-P,M-Q)-R)+I) -
|
||||
$ COS(THETA(I))
|
||||
END DO
|
||||
*
|
||||
* Compute norm( U1'*X11*V1 - D11 ) / ( MAX(1,P,Q)*EPS2 ) .
|
||||
*
|
||||
RESID = SLANGE( '1', P, Q, X, LDX, RWORK )
|
||||
RESID = SLANGE( '1', P, Q, XF, LDX, RWORK )
|
||||
RESULT( 1 ) = ( RESID / REAL(MAX(1,P,Q)) ) / EPS2
|
||||
*
|
||||
* Compute norm( U1'*X12*V2 - D12 ) / ( MAX(1,P,M-Q)*EPS2 ) .
|
||||
*
|
||||
RESID = SLANGE( '1', P, M-Q, X(1,Q+1), LDX, RWORK )
|
||||
RESID = SLANGE( '1', P, M-Q, XF(1,Q+1), LDX, RWORK )
|
||||
RESULT( 2 ) = ( RESID / REAL(MAX(1,P,M-Q)) ) / EPS2
|
||||
*
|
||||
* Compute norm( U2'*X21*V1 - D21 ) / ( MAX(1,M-P,Q)*EPS2 ) .
|
||||
*
|
||||
RESID = SLANGE( '1', M-P, Q, X(P+1,1), LDX, RWORK )
|
||||
RESID = SLANGE( '1', M-P, Q, XF(P+1,1), LDX, RWORK )
|
||||
RESULT( 3 ) = ( RESID / REAL(MAX(1,M-P,Q)) ) / EPS2
|
||||
*
|
||||
* Compute norm( U2'*X22*V2 - D22 ) / ( MAX(1,M-P,M-Q)*EPS2 ) .
|
||||
*
|
||||
RESID = SLANGE( '1', M-P, M-Q, X(P+1,Q+1), LDX, RWORK )
|
||||
RESID = SLANGE( '1', M-P, M-Q, XF(P+1,Q+1), LDX, RWORK )
|
||||
RESULT( 4 ) = ( RESID / REAL(MAX(1,M-P,M-Q)) ) / EPS2
|
||||
*
|
||||
* Compute I - U1'*U1
|
||||
@@ -396,14 +427,125 @@
|
||||
*
|
||||
* Check sorting
|
||||
*
|
||||
RESULT(9) = REALZERO
|
||||
RESULT( 9 ) = REALZERO
|
||||
DO I = 1, R
|
||||
IF( THETA(I).LT.REALZERO .OR. THETA(I).GT.PIOVER2 ) THEN
|
||||
RESULT(9) = ULPINV
|
||||
RESULT( 9 ) = ULPINV
|
||||
END IF
|
||||
IF( I.GT.1) THEN
|
||||
IF( I.GT.1 ) THEN
|
||||
IF ( THETA(I).LT.THETA(I-1) ) THEN
|
||||
RESULT(9) = ULPINV
|
||||
RESULT( 9 ) = ULPINV
|
||||
END IF
|
||||
END IF
|
||||
END DO
|
||||
*
|
||||
* The second half of the routine checks the 2-by-1 CSD
|
||||
*
|
||||
CALL SLASET( 'Full', Q, Q, ZERO, ONE, WORK, LDX )
|
||||
CALL SSYRK( 'Upper', 'Conjugate transpose', Q, M, -ONE, X, LDX,
|
||||
$ ONE, WORK, LDX )
|
||||
IF (M.GT.0) THEN
|
||||
EPS2 = MAX( ULP,
|
||||
$ SLANGE( '1', Q, Q, WORK, LDX, RWORK ) / REAL( M ) )
|
||||
ELSE
|
||||
EPS2 = ULP
|
||||
END IF
|
||||
R = MIN( P, M-P, Q, M-Q )
|
||||
*
|
||||
* Copy the matrix [X11;X21] to the array XF.
|
||||
*
|
||||
CALL SLACPY( 'Full', M, Q, X, LDX, XF, LDX )
|
||||
*
|
||||
* Compute the CSD
|
||||
*
|
||||
CALL SORCSD2BY1( 'Y', 'Y', 'Y', M, P, Q, XF(1,1), LDX, XF(P+1,1),
|
||||
$ LDX, THETA, U1, LDU1, U2, LDU2, V1T, LDV1T, WORK,
|
||||
$ LWORK, IWORK, INFO )
|
||||
*
|
||||
* Compute [X11;X21] := diag(U1,U2)'*[X11;X21]*V1 - [D11;D21]
|
||||
*
|
||||
CALL SGEMM( 'No transpose', 'Conjugate transpose', P, Q, Q, ONE,
|
||||
$ X, LDX, V1T, LDV1T, ZERO, WORK, LDX )
|
||||
*
|
||||
CALL SGEMM( 'Conjugate transpose', 'No transpose', P, Q, P, ONE,
|
||||
$ U1, LDU1, WORK, LDX, ZERO, X, LDX )
|
||||
*
|
||||
DO I = 1, MIN(P,Q)-R
|
||||
X(I,I) = X(I,I) - ONE
|
||||
END DO
|
||||
DO I = 1, R
|
||||
X(MIN(P,Q)-R+I,MIN(P,Q)-R+I) =
|
||||
$ X(MIN(P,Q)-R+I,MIN(P,Q)-R+I) - COS(THETA(I))
|
||||
END DO
|
||||
*
|
||||
CALL SGEMM( 'No transpose', 'Conjugate transpose', M-P, Q, Q, ONE,
|
||||
$ X(P+1,1), LDX, V1T, LDV1T, ZERO, WORK, LDX )
|
||||
*
|
||||
CALL SGEMM( 'Conjugate transpose', 'No transpose', M-P, Q, M-P,
|
||||
$ ONE, U2, LDU2, WORK, LDX, ZERO, X(P+1,1), LDX )
|
||||
*
|
||||
DO I = 1, MIN(M-P,Q)-R
|
||||
X(M-I+1,Q-I+1) = X(M-I+1,Q-I+1) - ONE
|
||||
END DO
|
||||
DO I = 1, R
|
||||
X(M-(MIN(M-P,Q)-R)+1-I,Q-(MIN(M-P,Q)-R)+1-I) =
|
||||
$ X(M-(MIN(M-P,Q)-R)+1-I,Q-(MIN(M-P,Q)-R)+1-I) -
|
||||
$ SIN(THETA(R-I+1))
|
||||
END DO
|
||||
*
|
||||
* Compute norm( U1'*X11*V1 - D11 ) / ( MAX(1,P,Q)*EPS2 ) .
|
||||
*
|
||||
RESID = SLANGE( '1', P, Q, X, LDX, RWORK )
|
||||
RESULT( 10 ) = ( RESID / REAL(MAX(1,P,Q)) ) / EPS2
|
||||
*
|
||||
* Compute norm( U2'*X21*V1 - D21 ) / ( MAX(1,M-P,Q)*EPS2 ) .
|
||||
*
|
||||
RESID = SLANGE( '1', M-P, Q, X(P+1,1), LDX, RWORK )
|
||||
RESULT( 11 ) = ( RESID / REAL(MAX(1,M-P,Q)) ) / EPS2
|
||||
*
|
||||
* Compute I - U1'*U1
|
||||
*
|
||||
CALL SLASET( 'Full', P, P, ZERO, ONE, WORK, LDU1 )
|
||||
CALL SSYRK( 'Upper', 'Conjugate transpose', P, P, -ONE, U1, LDU1,
|
||||
$ ONE, WORK, LDU1 )
|
||||
*
|
||||
* Compute norm( I - U1'*U1 ) / ( MAX(1,P) * ULP ) .
|
||||
*
|
||||
RESID = SLANSY( '1', 'Upper', P, WORK, LDU1, RWORK )
|
||||
RESULT( 12 ) = ( RESID / REAL(MAX(1,P)) ) / ULP
|
||||
*
|
||||
* Compute I - U2'*U2
|
||||
*
|
||||
CALL SLASET( 'Full', M-P, M-P, ZERO, ONE, WORK, LDU2 )
|
||||
CALL SSYRK( 'Upper', 'Conjugate transpose', M-P, M-P, -ONE, U2,
|
||||
$ LDU2, ONE, WORK, LDU2 )
|
||||
*
|
||||
* Compute norm( I - U2'*U2 ) / ( MAX(1,M-P) * ULP ) .
|
||||
*
|
||||
RESID = SLANSY( '1', 'Upper', M-P, WORK, LDU2, RWORK )
|
||||
RESULT( 13 ) = ( RESID / REAL(MAX(1,M-P)) ) / ULP
|
||||
*
|
||||
* Compute I - V1T*V1T'
|
||||
*
|
||||
CALL SLASET( 'Full', Q, Q, ZERO, ONE, WORK, LDV1T )
|
||||
CALL SSYRK( 'Upper', 'No transpose', Q, Q, -ONE, V1T, LDV1T, ONE,
|
||||
$ WORK, LDV1T )
|
||||
*
|
||||
* Compute norm( I - V1T*V1T' ) / ( MAX(1,Q) * ULP ) .
|
||||
*
|
||||
RESID = SLANSY( '1', 'Upper', Q, WORK, LDV1T, RWORK )
|
||||
RESULT( 14 ) = ( RESID / REAL(MAX(1,Q)) ) / ULP
|
||||
*
|
||||
* Check sorting
|
||||
*
|
||||
RESULT( 15 ) = REALZERO
|
||||
DO I = 1, R
|
||||
IF( THETA(I).LT.REALZERO .OR. THETA(I).GT.PIOVER2 ) THEN
|
||||
RESULT( 15 ) = ULPINV
|
||||
END IF
|
||||
IF( I.GT.1 ) THEN
|
||||
IF ( THETA(I).LT.THETA(I-1) ) THEN
|
||||
RESULT( 15 ) = ULPINV
|
||||
END IF
|
||||
END IF
|
||||
END DO
|
||||
|
||||
@@ -1026,17 +1026,17 @@
|
||||
*> \author Univ. of Colorado Denver
|
||||
*> \author NAG Ltd.
|
||||
*
|
||||
*> \date April 2012
|
||||
*> \date November 2013
|
||||
*
|
||||
*> \ingroup complex16_eig
|
||||
*
|
||||
* =====================================================================
|
||||
PROGRAM ZCHKEE
|
||||
*
|
||||
* -- LAPACK test routine (version 3.4.1) --
|
||||
* -- LAPACK test routine (version 3.5.0) --
|
||||
* -- LAPACK is a software package provided by Univ. of Tennessee, --
|
||||
* -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
|
||||
* April 2012
|
||||
* November 2013
|
||||
*
|
||||
* =====================================================================
|
||||
*
|
||||
@@ -1129,6 +1129,10 @@
|
||||
* ..
|
||||
* .. Executable Statements ..
|
||||
*
|
||||
A = 0.0
|
||||
B = 0.0
|
||||
C = 0.0
|
||||
DC = 0.0
|
||||
S1 = DSECND( )
|
||||
FATAL = .FALSE.
|
||||
NUNIT = NOUT
|
||||
|
||||
@@ -205,13 +205,16 @@
|
||||
*
|
||||
* .. Parameters ..
|
||||
INTEGER NTESTS
|
||||
PARAMETER ( NTESTS = 9 )
|
||||
PARAMETER ( NTESTS = 15 )
|
||||
INTEGER NTYPES
|
||||
PARAMETER ( NTYPES = 3 )
|
||||
DOUBLE PRECISION GAPDIGIT, ORTH, PIOVER2, TEN
|
||||
PARAMETER ( NTYPES = 4 )
|
||||
DOUBLE PRECISION GAPDIGIT, ORTH, PIOVER2, REALONE, REALZERO, TEN
|
||||
PARAMETER ( GAPDIGIT = 18.0D0, ORTH = 1.0D-12,
|
||||
$ PIOVER2 = 1.57079632679489662D0,
|
||||
$ REALONE = 1.0D0, REALZERO = 0.0D0,
|
||||
$ TEN = 10.0D0 )
|
||||
COMPLEX*16 ONE, ZERO
|
||||
PARAMETER ( ONE = (1.0D0,0.0D0), ZERO = (0.0D0,0.0D0) )
|
||||
* ..
|
||||
* .. Local Scalars ..
|
||||
LOGICAL FIRSTT
|
||||
@@ -231,8 +234,8 @@
|
||||
INTRINSIC ABS, MIN
|
||||
* ..
|
||||
* .. External Functions ..
|
||||
DOUBLE PRECISION DLARND
|
||||
EXTERNAL DLARND
|
||||
DOUBLE PRECISION DLARAN, DLARND
|
||||
EXTERNAL DLARAN, DLARND
|
||||
* ..
|
||||
* .. Executable Statements ..
|
||||
*
|
||||
@@ -286,7 +289,7 @@
|
||||
$ ORTH*DLARND(2,ISEED)
|
||||
END DO
|
||||
END DO
|
||||
ELSE
|
||||
ELSE IF( IMAT.EQ.3 ) THEN
|
||||
R = MIN( P, M-P, Q, M-Q )
|
||||
DO I = 1, R+1
|
||||
THETA(I) = TEN**(-DLARND(1,ISEED)*GAPDIGIT)
|
||||
@@ -298,9 +301,18 @@
|
||||
THETA(I) = PIOVER2 * THETA(I) / THETA(R+1)
|
||||
END DO
|
||||
CALL ZLACSG( M, P, Q, THETA, ISEED, X, LDX, WORK )
|
||||
ELSE
|
||||
CALL ZLASET( 'F', M, M, ZERO, ONE, X, LDX )
|
||||
DO I = 1, M
|
||||
J = INT( DLARAN( ISEED ) * M ) + 1
|
||||
IF( J .NE. I ) THEN
|
||||
CALL ZDROT( M, X(1+(I-1)*LDX), 1, X(1+(J-1)*LDX),
|
||||
$ 1, REALZERO, REALONE )
|
||||
END IF
|
||||
END DO
|
||||
END IF
|
||||
*
|
||||
NT = 9
|
||||
NT = 15
|
||||
*
|
||||
CALL ZCSDTS( M, P, Q, X, XF, LDX, U1, LDU1, U2, LDU2, V1T,
|
||||
$ LDV1T, V2T, LDV2T, THETA, IWORK, WORK, LWORK,
|
||||
|
||||
@@ -17,7 +17,7 @@
|
||||
* ..
|
||||
* .. Array Arguments ..
|
||||
* INTEGER IWORK( * )
|
||||
* DOUBLE PRECISION RESULT( 9 ), RWORK( * ), THETA( * )
|
||||
* DOUBLE PRECISION RESULT( 15 ), RWORK( * ), THETA( * )
|
||||
* COMPLEX*16 U1( LDU1, * ), U2( LDU2, * ), V1T( LDV1T, * ),
|
||||
* $ V2T( LDV2T, * ), WORK( LWORK ), X( LDX, * ),
|
||||
* $ XF( LDX, * )
|
||||
@@ -47,6 +47,21 @@
|
||||
*> [ 0 0 0 | I 0 0 ] [ D21 D22 ]
|
||||
*> [ 0 S 0 | 0 C 0 ]
|
||||
*> [ 0 0 I | 0 0 0 ]
|
||||
*>
|
||||
*> and also SORCSD2BY1, which, given
|
||||
*> Q
|
||||
*> [ X11 ] P ,
|
||||
*> [ X21 ] M-P
|
||||
*>
|
||||
*> computes the 2-by-1 CSD
|
||||
*>
|
||||
*> [ I 0 0 ]
|
||||
*> [ 0 C 0 ]
|
||||
*> [ 0 0 0 ]
|
||||
*> [ U1 ]**T * [ X11 ] * V1 = [----------] = [ D11 ] ,
|
||||
*> [ U2 ] [ X21 ] [ 0 0 0 ] [ D21 ]
|
||||
*> [ 0 S 0 ]
|
||||
*> [ 0 0 I ]
|
||||
*> \endverbatim
|
||||
*
|
||||
* Arguments:
|
||||
@@ -171,8 +186,9 @@
|
||||
*>
|
||||
*> \param[out] RESULT
|
||||
*> \verbatim
|
||||
*> RESULT is DOUBLE PRECISION array, dimension (9)
|
||||
*> RESULT is DOUBLE PRECISION array, dimension (15)
|
||||
*> The test ratios:
|
||||
*> First, the 2-by-2 CSD:
|
||||
*> RESULT(1) = norm( U1'*X11*V1 - D11 ) / ( MAX(1,P,Q)*EPS2 )
|
||||
*> RESULT(2) = norm( U1'*X12*V2 - D12 ) / ( MAX(1,P,M-Q)*EPS2 )
|
||||
*> RESULT(3) = norm( U2'*X21*V1 - D21 ) / ( MAX(1,M-P,Q)*EPS2 )
|
||||
@@ -184,6 +200,15 @@
|
||||
*> RESULT(9) = 0 if THETA is in increasing order and
|
||||
*> all angles are in [0,pi/2];
|
||||
*> = ULPINV otherwise.
|
||||
*> Then, the 2-by-1 CSD:
|
||||
*> RESULT(10) = norm( U1'*X11*V1 - D11 ) / ( MAX(1,P,Q)*EPS2 )
|
||||
*> RESULT(11) = norm( U2'*X21*V1 - D21 ) / ( MAX(1,M-P,Q)*EPS2 )
|
||||
*> RESULT(12) = norm( I - U1'*U1 ) / ( MAX(1,P)*ULP )
|
||||
*> RESULT(13) = norm( I - U2'*U2 ) / ( MAX(1,M-P)*ULP )
|
||||
*> RESULT(14) = norm( I - V1T'*V1T ) / ( MAX(1,Q)*ULP )
|
||||
*> RESULT(15) = 0 if THETA is in increasing order and
|
||||
*> all angles are in [0,pi/2];
|
||||
*> = ULPINV otherwise.
|
||||
*> ( EPS2 = MAX( norm( I - X'*X ) / M, ULP ). )
|
||||
*> \endverbatim
|
||||
*
|
||||
@@ -214,7 +239,7 @@
|
||||
* ..
|
||||
* .. Array Arguments ..
|
||||
INTEGER IWORK( * )
|
||||
DOUBLE PRECISION RESULT( 9 ), RWORK( * ), THETA( * )
|
||||
DOUBLE PRECISION RESULT( 15 ), RWORK( * ), THETA( * )
|
||||
COMPLEX*16 U1( LDU1, * ), U2( LDU2, * ), V1T( LDV1T, * ),
|
||||
$ V2T( LDV2T, * ), WORK( LWORK ), X( LDX, * ),
|
||||
$ XF( LDX, * )
|
||||
@@ -238,15 +263,19 @@
|
||||
EXTERNAL DLAMCH, ZLANGE, ZLANHE
|
||||
* ..
|
||||
* .. External Subroutines ..
|
||||
EXTERNAL ZGEMM, ZLACPY, ZLASET, ZUNCSD, ZHERK
|
||||
EXTERNAL ZGEMM, ZHERK, ZLACPY, ZLASET, ZUNCSD,
|
||||
$ ZUNCSD2BY1
|
||||
* ..
|
||||
* .. Intrinsic Functions ..
|
||||
INTRINSIC REAL, MAX, MIN
|
||||
INTRINSIC COS, DBLE, DCMPLX, MAX, MIN, REAL, SIN
|
||||
* ..
|
||||
* .. Executable Statements ..
|
||||
*
|
||||
ULP = DLAMCH( 'Precision' )
|
||||
ULPINV = REALONE / ULP
|
||||
*
|
||||
* The first half of the routine checks the 2-by-2 CSD
|
||||
*
|
||||
CALL ZLASET( 'Full', M, M, ZERO, ONE, WORK, LDX )
|
||||
CALL ZHERK( 'Upper', 'Conjugate transpose', M, M, -REALONE,
|
||||
$ X, LDX, REALONE, WORK, LDX )
|
||||
@@ -269,86 +298,88 @@
|
||||
$ THETA, U1, LDU1, U2, LDU2, V1T, LDV1T, V2T, LDV2T,
|
||||
$ WORK, LWORK, RWORK, 17*(R+2), IWORK, INFO )
|
||||
*
|
||||
* Compute X := diag(U1,U2)'*X*diag(V1,V2) - [D11 D12; D21 D22]
|
||||
* Compute XF := diag(U1,U2)'*X*diag(V1,V2) - [D11 D12; D21 D22]
|
||||
*
|
||||
CALL ZLACPY( 'Full', M, M, X, LDX, XF, LDX )
|
||||
*
|
||||
CALL ZGEMM( 'No transpose', 'Conjugate transpose', P, Q, Q, ONE,
|
||||
$ X, LDX, V1T, LDV1T, ZERO, WORK, LDX )
|
||||
$ XF, LDX, V1T, LDV1T, ZERO, WORK, LDX )
|
||||
*
|
||||
CALL ZGEMM( 'Conjugate transpose', 'No transpose', P, Q, P, ONE,
|
||||
$ U1, LDU1, WORK, LDX, ZERO, X, LDX )
|
||||
$ U1, LDU1, WORK, LDX, ZERO, XF, LDX )
|
||||
*
|
||||
DO I = 1, MIN(P,Q)-R
|
||||
X(I,I) = X(I,I) - ONE
|
||||
XF(I,I) = XF(I,I) - ONE
|
||||
END DO
|
||||
DO I = 1, R
|
||||
X(MIN(P,Q)-R+I,MIN(P,Q)-R+I) =
|
||||
$ X(MIN(P,Q)-R+I,MIN(P,Q)-R+I) - DCMPLX( COS(THETA(I)),
|
||||
XF(MIN(P,Q)-R+I,MIN(P,Q)-R+I) =
|
||||
$ XF(MIN(P,Q)-R+I,MIN(P,Q)-R+I) - DCMPLX( COS(THETA(I)),
|
||||
$ 0.0D0 )
|
||||
END DO
|
||||
*
|
||||
CALL ZGEMM( 'No transpose', 'Conjugate transpose', P, M-Q, M-Q,
|
||||
$ ONE, X(1,Q+1), LDX, V2T, LDV2T, ZERO, WORK, LDX )
|
||||
$ ONE, XF(1,Q+1), LDX, V2T, LDV2T, ZERO, WORK, LDX )
|
||||
*
|
||||
CALL ZGEMM( 'Conjugate transpose', 'No transpose', P, M-Q, P,
|
||||
$ ONE, U1, LDU1, WORK, LDX, ZERO, X(1,Q+1), LDX )
|
||||
$ ONE, U1, LDU1, WORK, LDX, ZERO, XF(1,Q+1), LDX )
|
||||
*
|
||||
DO I = 1, MIN(P,M-Q)-R
|
||||
X(P-I+1,M-I+1) = X(P-I+1,M-I+1) + ONE
|
||||
XF(P-I+1,M-I+1) = XF(P-I+1,M-I+1) + ONE
|
||||
END DO
|
||||
DO I = 1, R
|
||||
X(P-(MIN(P,M-Q)-R)+1-I,M-(MIN(P,M-Q)-R)+1-I) =
|
||||
$ X(P-(MIN(P,M-Q)-R)+1-I,M-(MIN(P,M-Q)-R)+1-I) +
|
||||
XF(P-(MIN(P,M-Q)-R)+1-I,M-(MIN(P,M-Q)-R)+1-I) =
|
||||
$ XF(P-(MIN(P,M-Q)-R)+1-I,M-(MIN(P,M-Q)-R)+1-I) +
|
||||
$ DCMPLX( SIN(THETA(R-I+1)), 0.0D0 )
|
||||
END DO
|
||||
*
|
||||
CALL ZGEMM( 'No transpose', 'Conjugate transpose', M-P, Q, Q, ONE,
|
||||
$ X(P+1,1), LDX, V1T, LDV1T, ZERO, WORK, LDX )
|
||||
$ XF(P+1,1), LDX, V1T, LDV1T, ZERO, WORK, LDX )
|
||||
*
|
||||
CALL ZGEMM( 'Conjugate transpose', 'No transpose', M-P, Q, M-P,
|
||||
$ ONE, U2, LDU2, WORK, LDX, ZERO, X(P+1,1), LDX )
|
||||
$ ONE, U2, LDU2, WORK, LDX, ZERO, XF(P+1,1), LDX )
|
||||
*
|
||||
DO I = 1, MIN(M-P,Q)-R
|
||||
X(M-I+1,Q-I+1) = X(M-I+1,Q-I+1) - ONE
|
||||
XF(M-I+1,Q-I+1) = XF(M-I+1,Q-I+1) - ONE
|
||||
END DO
|
||||
DO I = 1, R
|
||||
X(M-(MIN(M-P,Q)-R)+1-I,Q-(MIN(M-P,Q)-R)+1-I) =
|
||||
$ X(M-(MIN(M-P,Q)-R)+1-I,Q-(MIN(M-P,Q)-R)+1-I) -
|
||||
XF(M-(MIN(M-P,Q)-R)+1-I,Q-(MIN(M-P,Q)-R)+1-I) =
|
||||
$ XF(M-(MIN(M-P,Q)-R)+1-I,Q-(MIN(M-P,Q)-R)+1-I) -
|
||||
$ DCMPLX( SIN(THETA(R-I+1)), 0.0D0 )
|
||||
END DO
|
||||
*
|
||||
CALL ZGEMM( 'No transpose', 'Conjugate transpose', M-P, M-Q, M-Q,
|
||||
$ ONE, X(P+1,Q+1), LDX, V2T, LDV2T, ZERO, WORK, LDX )
|
||||
$ ONE, XF(P+1,Q+1), LDX, V2T, LDV2T, ZERO, WORK, LDX )
|
||||
*
|
||||
CALL ZGEMM( 'Conjugate transpose', 'No transpose', M-P, M-Q, M-P,
|
||||
$ ONE, U2, LDU2, WORK, LDX, ZERO, X(P+1,Q+1), LDX )
|
||||
$ ONE, U2, LDU2, WORK, LDX, ZERO, XF(P+1,Q+1), LDX )
|
||||
*
|
||||
DO I = 1, MIN(M-P,M-Q)-R
|
||||
X(P+I,Q+I) = X(P+I,Q+I) - ONE
|
||||
XF(P+I,Q+I) = XF(P+I,Q+I) - ONE
|
||||
END DO
|
||||
DO I = 1, R
|
||||
X(P+(MIN(M-P,M-Q)-R)+I,Q+(MIN(M-P,M-Q)-R)+I) =
|
||||
$ X(P+(MIN(M-P,M-Q)-R)+I,Q+(MIN(M-P,M-Q)-R)+I) -
|
||||
XF(P+(MIN(M-P,M-Q)-R)+I,Q+(MIN(M-P,M-Q)-R)+I) =
|
||||
$ XF(P+(MIN(M-P,M-Q)-R)+I,Q+(MIN(M-P,M-Q)-R)+I) -
|
||||
$ DCMPLX( COS(THETA(I)), 0.0D0 )
|
||||
END DO
|
||||
*
|
||||
* Compute norm( U1'*X11*V1 - D11 ) / ( MAX(1,P,Q)*EPS2 ) .
|
||||
*
|
||||
RESID = ZLANGE( '1', P, Q, X, LDX, RWORK )
|
||||
RESID = ZLANGE( '1', P, Q, XF, LDX, RWORK )
|
||||
RESULT( 1 ) = ( RESID / REAL(MAX(1,P,Q)) ) / EPS2
|
||||
*
|
||||
* Compute norm( U1'*X12*V2 - D12 ) / ( MAX(1,P,M-Q)*EPS2 ) .
|
||||
*
|
||||
RESID = ZLANGE( '1', P, M-Q, X(1,Q+1), LDX, RWORK )
|
||||
RESID = ZLANGE( '1', P, M-Q, XF(1,Q+1), LDX, RWORK )
|
||||
RESULT( 2 ) = ( RESID / REAL(MAX(1,P,M-Q)) ) / EPS2
|
||||
*
|
||||
* Compute norm( U2'*X21*V1 - D21 ) / ( MAX(1,M-P,Q)*EPS2 ) .
|
||||
*
|
||||
RESID = ZLANGE( '1', M-P, Q, X(P+1,1), LDX, RWORK )
|
||||
RESID = ZLANGE( '1', M-P, Q, XF(P+1,1), LDX, RWORK )
|
||||
RESULT( 3 ) = ( RESID / REAL(MAX(1,M-P,Q)) ) / EPS2
|
||||
*
|
||||
* Compute norm( U2'*X22*V2 - D22 ) / ( MAX(1,M-P,M-Q)*EPS2 ) .
|
||||
*
|
||||
RESID = ZLANGE( '1', M-P, M-Q, X(P+1,Q+1), LDX, RWORK )
|
||||
RESID = ZLANGE( '1', M-P, M-Q, XF(P+1,Q+1), LDX, RWORK )
|
||||
RESULT( 4 ) = ( RESID / REAL(MAX(1,M-P,M-Q)) ) / EPS2
|
||||
*
|
||||
* Compute I - U1'*U1
|
||||
@@ -397,14 +428,126 @@
|
||||
*
|
||||
* Check sorting
|
||||
*
|
||||
RESULT(9) = REALZERO
|
||||
RESULT( 9 ) = REALZERO
|
||||
DO I = 1, R
|
||||
IF( THETA(I).LT.REALZERO .OR. THETA(I).GT.PIOVER2 ) THEN
|
||||
RESULT(9) = ULPINV
|
||||
RESULT( 9 ) = ULPINV
|
||||
END IF
|
||||
IF( I.GT.1) THEN
|
||||
IF ( THETA(I).LT.THETA(I-1) ) THEN
|
||||
RESULT(9) = ULPINV
|
||||
RESULT( 9 ) = ULPINV
|
||||
END IF
|
||||
END IF
|
||||
END DO
|
||||
*
|
||||
* The second half of the routine checks the 2-by-1 CSD
|
||||
*
|
||||
CALL ZLASET( 'Full', Q, Q, ZERO, ONE, WORK, LDX )
|
||||
CALL ZHERK( 'Upper', 'Conjugate transpose', Q, M, -REALONE,
|
||||
$ X, LDX, REALONE, WORK, LDX )
|
||||
IF (M.GT.0) THEN
|
||||
EPS2 = MAX( ULP,
|
||||
$ ZLANGE( '1', Q, Q, WORK, LDX, RWORK ) / DBLE( M ) )
|
||||
ELSE
|
||||
EPS2 = ULP
|
||||
END IF
|
||||
R = MIN( P, M-P, Q, M-Q )
|
||||
*
|
||||
* Copy the matrix X to the array XF.
|
||||
*
|
||||
CALL ZLACPY( 'Full', M, M, X, LDX, XF, LDX )
|
||||
*
|
||||
* Compute the CSD
|
||||
*
|
||||
CALL ZUNCSD2BY1( 'Y', 'Y', 'Y', M, P, Q, XF(1,1), LDX, XF(P+1,1),
|
||||
$ LDX, THETA, U1, LDU1, U2, LDU2, V1T, LDV1T, WORK,
|
||||
$ LWORK, RWORK, 17*(R+2), IWORK, INFO )
|
||||
*
|
||||
* Compute [X11;X21] := diag(U1,U2)'*[X11;X21]*V1 - [D11;D21]
|
||||
*
|
||||
CALL ZGEMM( 'No transpose', 'Conjugate transpose', P, Q, Q, ONE,
|
||||
$ X, LDX, V1T, LDV1T, ZERO, WORK, LDX )
|
||||
*
|
||||
CALL ZGEMM( 'Conjugate transpose', 'No transpose', P, Q, P, ONE,
|
||||
$ U1, LDU1, WORK, LDX, ZERO, X, LDX )
|
||||
*
|
||||
DO I = 1, MIN(P,Q)-R
|
||||
X(I,I) = X(I,I) - ONE
|
||||
END DO
|
||||
DO I = 1, R
|
||||
X(MIN(P,Q)-R+I,MIN(P,Q)-R+I) =
|
||||
$ X(MIN(P,Q)-R+I,MIN(P,Q)-R+I) - DCMPLX( COS(THETA(I)),
|
||||
$ 0.0D0 )
|
||||
END DO
|
||||
*
|
||||
CALL ZGEMM( 'No transpose', 'Conjugate transpose', M-P, Q, Q, ONE,
|
||||
$ X(P+1,1), LDX, V1T, LDV1T, ZERO, WORK, LDX )
|
||||
*
|
||||
CALL ZGEMM( 'Conjugate transpose', 'No transpose', M-P, Q, M-P,
|
||||
$ ONE, U2, LDU2, WORK, LDX, ZERO, X(P+1,1), LDX )
|
||||
*
|
||||
DO I = 1, MIN(M-P,Q)-R
|
||||
X(M-I+1,Q-I+1) = X(M-I+1,Q-I+1) - ONE
|
||||
END DO
|
||||
DO I = 1, R
|
||||
X(M-(MIN(M-P,Q)-R)+1-I,Q-(MIN(M-P,Q)-R)+1-I) =
|
||||
$ X(M-(MIN(M-P,Q)-R)+1-I,Q-(MIN(M-P,Q)-R)+1-I) -
|
||||
$ DCMPLX( SIN(THETA(R-I+1)), 0.0D0 )
|
||||
END DO
|
||||
*
|
||||
* Compute norm( U1'*X11*V1 - D11 ) / ( MAX(1,P,Q)*EPS2 ) .
|
||||
*
|
||||
RESID = ZLANGE( '1', P, Q, X, LDX, RWORK )
|
||||
RESULT( 10 ) = ( RESID / REAL(MAX(1,P,Q)) ) / EPS2
|
||||
*
|
||||
* Compute norm( U2'*X21*V1 - D21 ) / ( MAX(1,M-P,Q)*EPS2 ) .
|
||||
*
|
||||
RESID = ZLANGE( '1', M-P, Q, X(P+1,1), LDX, RWORK )
|
||||
RESULT( 11 ) = ( RESID / REAL(MAX(1,M-P,Q)) ) / EPS2
|
||||
*
|
||||
* Compute I - U1'*U1
|
||||
*
|
||||
CALL ZLASET( 'Full', P, P, ZERO, ONE, WORK, LDU1 )
|
||||
CALL ZHERK( 'Upper', 'Conjugate transpose', P, P, -REALONE,
|
||||
$ U1, LDU1, REALONE, WORK, LDU1 )
|
||||
*
|
||||
* Compute norm( I - U'*U ) / ( MAX(1,P) * ULP ) .
|
||||
*
|
||||
RESID = ZLANHE( '1', 'Upper', P, WORK, LDU1, RWORK )
|
||||
RESULT( 12 ) = ( RESID / REAL(MAX(1,P)) ) / ULP
|
||||
*
|
||||
* Compute I - U2'*U2
|
||||
*
|
||||
CALL ZLASET( 'Full', M-P, M-P, ZERO, ONE, WORK, LDU2 )
|
||||
CALL ZHERK( 'Upper', 'Conjugate transpose', M-P, M-P, -REALONE,
|
||||
$ U2, LDU2, REALONE, WORK, LDU2 )
|
||||
*
|
||||
* Compute norm( I - U2'*U2 ) / ( MAX(1,M-P) * ULP ) .
|
||||
*
|
||||
RESID = ZLANHE( '1', 'Upper', M-P, WORK, LDU2, RWORK )
|
||||
RESULT( 13 ) = ( RESID / REAL(MAX(1,M-P)) ) / ULP
|
||||
*
|
||||
* Compute I - V1T*V1T'
|
||||
*
|
||||
CALL ZLASET( 'Full', Q, Q, ZERO, ONE, WORK, LDV1T )
|
||||
CALL ZHERK( 'Upper', 'No transpose', Q, Q, -REALONE,
|
||||
$ V1T, LDV1T, REALONE, WORK, LDV1T )
|
||||
*
|
||||
* Compute norm( I - V1T*V1T' ) / ( MAX(1,Q) * ULP ) .
|
||||
*
|
||||
RESID = ZLANHE( '1', 'Upper', Q, WORK, LDV1T, RWORK )
|
||||
RESULT( 14 ) = ( RESID / REAL(MAX(1,Q)) ) / ULP
|
||||
*
|
||||
* Check sorting
|
||||
*
|
||||
RESULT( 15 ) = REALZERO
|
||||
DO I = 1, R
|
||||
IF( THETA(I).LT.REALZERO .OR. THETA(I).GT.PIOVER2 ) THEN
|
||||
RESULT( 15 ) = ULPINV
|
||||
END IF
|
||||
IF( I.GT.1) THEN
|
||||
IF ( THETA(I).LT.THETA(I-1) ) THEN
|
||||
RESULT( 15 ) = ULPINV
|
||||
END IF
|
||||
END IF
|
||||
END DO
|
||||
|
||||
Reference in New Issue
Block a user