Purpose
To compute a symplectic URV (SURV) decomposition of a real
2N-by-2N matrix H:
[ op(A) G ] T [ op(R11) R12 ] T
H = [ ] = U R V = U * [ ] * V ,
[ Q op(B) ] [ 0 op(R22) ]
where A, B, G, Q, R12 are real N-by-N matrices, op(R11) is a real
N-by-N upper triangular matrix, op(R22) is a real N-by-N lower
Hessenberg matrix and U, V are 2N-by-2N orthogonal symplectic
matrices. Unblocked version.
Specification
SUBROUTINE MB04TS( TRANA, TRANB, N, ILO, A, LDA, B, LDB, G, LDG,
$ Q, LDQ, CSL, CSR, TAUL, TAUR, DWORK, LDWORK,
$ INFO )
C .. Scalar Arguments ..
CHARACTER TRANA, TRANB
INTEGER ILO, INFO, LDA, LDB, LDG, LDQ, LDWORK, N
C .. Array Arguments ..
DOUBLE PRECISION A(LDA,*), B(LDB,*), CSL(*), CSR(*), DWORK(*),
$ G(LDG,*), Q(LDQ,*), TAUL(*), TAUR(*)
Arguments
Mode Parameters
TRANA CHARACTER*1
Specifies the form of op( A ) as follows:
= 'N': op( A ) = A;
= 'T': op( A ) = A';
= 'C': op( A ) = A'.
TRANB CHARACTER*1
Specifies the form of op( B ) as follows:
= 'N': op( B ) = B;
= 'T': op( B ) = B';
= 'C': op( B ) = B'.
Input/Output Parameters
N (input) INTEGER
The order of the matrix A. N >= 0.
ILO (input) INTEGER
It is assumed that op(A) is already upper triangular,
op(B) is lower triangular and Q is zero in rows and
columns 1:ILO-1. ILO is normally set by a previous call
to MB04DD; otherwise it should be set to 1.
1 <= ILO <= N, if N > 0; ILO=1, if N=0.
A (input/output) DOUBLE PRECISION array, dimension (LDA,N)
On entry, the leading N-by-N part of this array must
contain the matrix A.
On exit, the leading N-by-N part of this array contains
the triangular matrix R11, and in the zero part
information about the elementary reflectors used to
compute the SURV decomposition.
LDA INTEGER
The leading dimension of the array A. LDA >= MAX(1,N).
B (input/output) DOUBLE PRECISION array, dimension (LDB,N)
On entry, the leading N-by-N part of this array must
contain the matrix B.
On exit, the leading N-by-N part of this array contains
the Hessenberg matrix R22, and in the zero part
information about the elementary reflectors used to
compute the SURV decomposition.
LDB INTEGER
The leading dimension of the array B. LDB >= MAX(1,N).
G (input/output) DOUBLE PRECISION array, dimension (LDG,N)
On entry, the leading N-by-N part of this array must
contain the matrix G.
On exit, the leading N-by-N part of this array contains
the matrix R12.
LDG INTEGER
The leading dimension of the array G. LDG >= MAX(1,N).
Q (input/output) DOUBLE PRECISION array, dimension (LDQ,N)
On entry, the leading N-by-N part of this array must
contain the matrix Q.
On exit, the leading N-by-N part of this array contains
information about the elementary reflectors used to
compute the SURV decomposition.
LDQ INTEGER
The leading dimension of the array Q. LDG >= MAX(1,N).
CSL (output) DOUBLE PRECISION array, dimension (2N)
On exit, the first 2N elements of this array contain the
cosines and sines of the symplectic Givens rotations
applied from the left-hand side used to compute the SURV
decomposition.
CSR (output) DOUBLE PRECISION array, dimension (2N-2)
On exit, the first 2N-2 elements of this array contain the
cosines and sines of the symplectic Givens rotations
applied from the right-hand side used to compute the SURV
decomposition.
TAUL (output) DOUBLE PRECISION array, dimension (N)
On exit, the first N elements of this array contain the
scalar factors of some of the elementary reflectors
applied from the left-hand side.
TAUR (output) DOUBLE PRECISION array, dimension (N-1)
On exit, the first N-1 elements of this array contain the
scalar factors of some of the elementary reflectors
applied from the right-hand side.
Workspace
DWORK DOUBLE PRECISION array, dimension (LDWORK)
On exit, if INFO = 0, DWORK(1) returns the optimal
value of LDWORK.
On exit, if INFO = -16, DWORK(1) returns the minimum
value of LDWORK.
LDWORK INTEGER
The length of the array DWORK. LDWORK >= MAX(1,N).
Error Indicator
INFO INTEGER
= 0: successful exit;
< 0: if INFO = -i, the i-th argument had an illegal
value.
Method
The matrices U and V are represented as products of symplectic
reflectors and Givens rotations
U = diag( HU(1),HU(1) ) GU(1) diag( FU(1),FU(1) )
diag( HU(2),HU(2) ) GU(2) diag( FU(2),FU(2) )
....
diag( HU(n),HU(n) ) GU(n) diag( FU(n),FU(n) ),
V = diag( HV(1),HV(1) ) GV(1) diag( FV(1),FV(1) )
diag( HV(2),HV(2) ) GV(2) diag( FV(2),FV(2) )
....
diag( HV(n-1),HV(n-1) ) GV(n-1) diag( FV(n-1),FV(n-1) ).
Each HU(i) has the form
HU(i) = I - tau * v * v'
where tau is a real scalar, and v is a real vector with
v(1:i-1) = 0 and v(i) = 1; v(i+1:n) is stored on exit in
Q(i+1:n,i), and tau in Q(i,i).
Each FU(i) has the form
FU(i) = I - nu * w * w'
where nu is a real scalar, and w is a real vector with
w(1:i-1) = 0 and w(i) = 1; w(i+1:n) is stored on exit in
A(i+1:n,i), if op(A) = 'N', and in A(i,i+1:n), otherwise. The
scalar nu is stored in TAUL(i).
Each GU(i) is a Givens rotation acting on rows i and n+i,
where the cosine is stored in CSL(2*i-1) and the sine in
CSL(2*i).
Each HV(i) has the form
HV(i) = I - tau * v * v'
where tau is a real scalar, and v is a real vector with
v(1:i) = 0 and v(i+1) = 1; v(i+2:n) is stored on exit in
Q(i,i+2:n), and tau in Q(i,i+1).
Each FV(i) has the form
FV(i) = I - nu * w * w'
where nu is a real scalar, and w is a real vector with
w(1:i) = 0 and w(i+1) = 1; w(i+2:n) is stored on exit in
B(i,i+2:n), if op(B) = 'N', and in B(i+2:n,i), otherwise.
The scalar nu is stored in TAUR(i).
Each GV(i) is a Givens rotation acting on columns i+1 and n+i+1,
where the cosine is stored in CSR(2*i-1) and the sine in
CSR(2*i).
Numerical Aspects
The algorithm requires 80/3 N**3 + 20 N**2 + O(N) floating point operations and is numerically backward stable.References
[1] Benner, P., Mehrmann, V., and Xu, H.
A numerically stable, structure preserving method for
computing the eigenvalues of real Hamiltonian or symplectic
pencils. Numer. Math., Vol 78 (3), pp. 329-358, 1998.
Further Comments
NoneExample
Program Text
* MB04TS/MB04WR EXAMPLE PROGRAM TEXT
* Copyright (c) 2002-2017 NICONET e.V.
*
* .. Parameters ..
DOUBLE PRECISION ZERO, ONE
PARAMETER ( ZERO = 0.0D0, ONE = 1.0D0 )
INTEGER NIN, NOUT
PARAMETER ( NIN = 5, NOUT = 6 )
INTEGER NMAX
PARAMETER ( NMAX = 200 )
INTEGER LDA, LDB, LDG, LDQ, LDRES, LDU1, LDU2, LDV1,
$ LDV2, LDWORK
PARAMETER ( LDA = NMAX, LDB = NMAX, LDG = NMAX, LDQ = NMAX,
$ LDRES = NMAX, LDU1 = NMAX, LDU2 = NMAX,
$ LDV1 = NMAX, LDV2 = NMAX, LDWORK = NMAX )
* .. Local Scalars ..
CHARACTER*1 TRANA, TRANB, TRANV1
INTEGER I, INFO, J, N
DOUBLE PRECISION TEMP
* .. Local Arrays ..
DOUBLE PRECISION A(LDA, NMAX), B(LDB, NMAX), CSL(2*NMAX),
$ CSR(2*NMAX), DWORK(LDWORK), G(LDG,NMAX),
$ Q(LDQ,NMAX), RES(LDRES,5*NMAX), TAUL(NMAX),
$ TAUR(NMAX), U1(LDU1,NMAX), U2(LDU2, NMAX),
$ V1(LDV1, NMAX), V2(LDV2,NMAX)
* .. External Functions ..
LOGICAL LSAME
DOUBLE PRECISION DLANGE, DLAPY2, MA02JD
EXTERNAL DLANGE, DLAPY2, LSAME, MA02JD
* .. External Subroutines ..
EXTERNAL DGEMM, DLACPY, DLASET, MB04TS, MB04WR
* .. Executable Statements ..
WRITE ( NOUT, FMT = 99999 )
* Skip the heading in the data file and read the data.
READ ( NIN, FMT = '()' )
READ ( NIN, FMT = * ) N, TRANA, TRANB
IF( N.LE.0 .OR. N.GT.NMAX ) THEN
WRITE ( NOUT, FMT = 99992 ) N
ELSE
READ ( NIN, FMT = * ) ( ( A(I,J), J = 1,N ), I = 1,N )
CALL DLACPY( 'All', N, N, A, LDA, RES, LDRES )
READ ( NIN, FMT = * ) ( ( B(I,J), J = 1,N ), I = 1,N )
CALL DLACPY( 'All', N, N, B, LDB, RES(1,N+1), LDRES )
READ ( NIN, FMT = * ) ( ( G(I,J), J = 1,N ), I = 1,N )
CALL DLACPY( 'All', N, N, G, LDG, RES(1,2*N+1), LDRES )
READ ( NIN, FMT = * ) ( ( Q(I,J), J = 1,N ), I = 1,N )
CALL DLACPY( 'All', N, N, Q, LDQ, RES(1,3*N+1), LDRES )
CALL MB04TS( TRANA, TRANB, N, 1, A, LDA, B, LDB, G, LDG, Q,
$ LDQ, CSL, CSR, TAUL, TAUR, DWORK, LDWORK, INFO )
IF ( INFO.NE.0 ) THEN
WRITE ( NOUT, FMT = 99998 ) INFO
ELSE
CALL DLACPY( 'All', N, N, A, LDA, U1, LDU1 )
CALL DLACPY( 'All', N, N, Q, LDQ, U2, LDU2 )
CALL MB04WR( 'U', TRANA, N, 1, U1, LDU1, U2, LDU2, CSL,
$ TAUL, DWORK, LDWORK, INFO )
IF ( INFO.NE.0 ) THEN
WRITE ( NOUT, FMT = 99997 ) INFO
ELSE
CALL DLACPY( 'All', N, N, Q, LDQ, V2, LDV2 )
CALL DLACPY( 'All', N, N, B, LDB, V1, LDV1 )
CALL MB04WR( 'V', TRANB, N, 1, V1, LDV1, V2, LDV2,
$ CSR, TAUR, DWORK, LDWORK, INFO )
IF ( INFO.NE.0 ) THEN
WRITE ( NOUT, FMT = 99997 ) INFO
ELSE
WRITE ( NOUT, FMT = 99996 )
IF ( LSAME( TRANA, 'N' ) ) THEN
DO 10 I = 1, N
WRITE (NOUT, FMT = 99993)
$ ( U1(I,J), J = 1,N ), ( U2(I,J), J = 1,N )
10 CONTINUE
DO 20 I = 1, N
WRITE (NOUT, FMT = 99993)
$ ( -U2(I,J), J = 1,N ), ( U1(I,J), J = 1,N )
20 CONTINUE
WRITE ( NOUT, FMT = 99991 ) MA02JD( .FALSE.,
$ .FALSE., N, U1, LDU1, U2, LDU2,
$ RES(1,4*N+1), LDRES )
ELSE
DO 30 I = 1, N
WRITE (NOUT, FMT = 99993)
$ ( U1(J,I), J = 1,N ), ( U2(I,J), J = 1,N )
30 CONTINUE
DO 40 I = 1, N
WRITE (NOUT, FMT = 99993)
$ ( -U2(I,J), J = 1,N ), ( U1(J,I), J = 1,N )
40 CONTINUE
WRITE ( NOUT, FMT = 99991 ) MA02JD( .TRUE.,
$ .FALSE., N, U1, LDU1, U2, LDU2,
$ RES(1,4*N+1), LDRES )
END IF
WRITE ( NOUT, FMT = 99995 )
CALL DLASET( 'All', N, N, ZERO, ZERO, Q, LDQ )
IF ( LSAME( TRANA, 'N' ) ) THEN
CALL DLASET( 'Lower', N-1, N-1, ZERO, ZERO,
$ A(2,1), LDA )
DO 50 I = 1, N
WRITE (NOUT, FMT = 99993)
$ ( A(I,J), J = 1,N ), ( G(I,J), J = 1,N )
50 CONTINUE
ELSE
CALL DLASET( 'Upper', N-1, N-1, ZERO, ZERO,
$ A(1,2), LDA )
DO 60 I = 1, N
WRITE (NOUT, FMT = 99993)
$ ( A(J,I), J = 1,N ), ( G(I,J), J = 1,N )
60 CONTINUE
END IF
IF ( LSAME( TRANB, 'N' ) ) THEN
IF ( N.GT.1 ) THEN
CALL DLASET( 'Upper', N-2, N-2, ZERO, ZERO,
$ B(1,3), LDB )
END IF
DO 70 I = 1, N
WRITE (NOUT, FMT = 99993)
$ ( Q(I,J), J = 1,N ), ( B(I,J), J = 1,N )
70 CONTINUE
ELSE
IF ( N.GT.1 ) THEN
CALL DLASET( 'Lower', N-2, N-2, ZERO, ZERO,
$ B(3,1), LDB )
END IF
DO 80 I = 1, N
WRITE (NOUT, FMT = 99993)
$ ( Q(I,J), J = 1,N ), ( B(J,I), J = 1,N )
80 CONTINUE
END IF
C
IF ( LSAME( TRANB, 'N' ) ) THEN
TRANV1 = 'T'
ELSE
TRANV1 = 'N'
END IF
CALL DGEMM( TRANA, TRANV1, N, N, N, ONE, RES, LDRES,
$ V1, LDV1, ZERO, RES(1,4*N+1), LDRES )
CALL DGEMM( 'No Transpose', 'Transpose', N, N, N,
$ -ONE, RES(1,2*N+1), LDRES, V2, LDV2, ONE,
$ RES(1,4*N+1), LDRES )
CALL DGEMM( TRANA, TRANA, N, N, N, -ONE, U1, LDU1,
$ A, LDA, ONE, RES(1,4*N+1), LDRES )
TEMP = DLANGE( 'Frobenius', N, N, RES(1,4*N+1),
$ LDRES, DWORK )
CALL DGEMM( TRANA, 'Transpose', N, N, N, ONE, RES,
$ LDRES, V2, LDV2, ZERO, RES(1,4*N+1),
$ LDRES )
CALL DGEMM( 'No Transpose', TRANV1, N, N, N, ONE,
$ RES(1,2*N+1), LDRES, V1, LDV1, ONE,
$ RES(1,4*N+1), LDRES )
CALL DGEMM( TRANA, 'No Transpose', N, N, N, -ONE,
$ U1, LDU1, G, LDG, ONE, RES(1,4*N+1),
$ LDRES )
CALL DGEMM( 'No Transpose', TRANB, N, N, N, -ONE,
$ U2, LDU2, B, LDB, ONE, RES(1,4*N+1),
$ LDRES )
TEMP = DLAPY2( TEMP, DLANGE( 'Frobenius', N, N,
$ RES(1,4*N+1), LDRES, DWORK ) )
CALL DGEMM( 'No Transpose', TRANV1, N, N, N, ONE,
$ RES(1,3*N+1), LDRES, V1, LDV1, ZERO,
$ RES(1,4*N+1), LDRES )
CALL DGEMM( TRANB, 'Transpose', N, N, N, -ONE,
$ RES(1,N+1), LDRES, V2, LDV2, ONE,
$ RES(1,4*N+1), LDRES )
CALL DGEMM( 'No Transpose', TRANA, N, N, N, ONE,
$ U2, LDU2, A, LDA, ONE, RES(1,4*N+1),
$ LDRES )
TEMP = DLAPY2( TEMP, DLANGE( 'Frobenius', N, N,
$ RES(1,4*N+1), LDRES, DWORK ) )
CALL DGEMM( 'No Transpose', 'Transpose', N, N, N, ONE,
$ RES(1,3*N+1), LDRES, V2, LDV2, ZERO,
$ RES(1,4*N+1), LDRES )
CALL DGEMM( TRANB, TRANV1, N, N, N, ONE, RES(1,N+1),
$ LDRES, V1, LDV1, ONE, RES(1,4*N+1),
$ LDRES )
CALL DGEMM( 'No Transpose', 'No Transpose', N, N, N,
$ ONE, U2, LDU2, G, LDG, ONE, RES(1,4*N+1),
$ LDRES )
CALL DGEMM( TRANA, TRANB, N, N, N, -ONE, U1, LDU1,
$ B, LDB, ONE, RES(1,4*N+1), LDRES )
TEMP = DLAPY2( TEMP, DLANGE( 'Frobenius', N, N,
$ RES(1,4*N+1), LDRES, DWORK ) )
WRITE ( NOUT, FMT = 99990 ) TEMP
C
WRITE ( NOUT, FMT = 99994 )
IF ( LSAME( TRANB, 'N' ) ) THEN
DO 90 I = 1, N
WRITE (NOUT, FMT = 99993)
$ ( V1(J,I), J = 1,N ), ( V2(J,I), J = 1,N )
90 CONTINUE
DO 100 I = 1, N
WRITE (NOUT, FMT = 99993)
$ ( -V2(J,I), J = 1,N ), ( V1(J,I), J = 1,N )
100 CONTINUE
WRITE ( NOUT, FMT = 99989 ) MA02JD( .TRUE.,
$ .TRUE., N, V1, LDV1, V2, LDV2,
$ RES(1,4*N+1), LDRES )
ELSE
DO 110 I = 1, N
WRITE (NOUT, FMT = 99993)
$ ( V1(I,J), J = 1,N ), ( V2(J,I), J = 1,N )
110 CONTINUE
DO 120 I = 1, N
WRITE (NOUT, FMT = 99993)
$ ( -V2(J,I), J = 1,N ), ( V1(I,J), J = 1,N )
120 CONTINUE
WRITE ( NOUT, FMT = 99989 ) MA02JD( .FALSE.,
$ .TRUE., N, V1, LDV1, V2, LDV2,
$ RES(1,4*N+1), LDRES )
END IF
END IF
END IF
END IF
END IF
*
STOP
*
99999 FORMAT (' MB04TS EXAMPLE PROGRAM RESULTS',/1X)
99998 FORMAT (' INFO on exit from MB04TS = ',I2)
99997 FORMAT (' INFO on exit from MB04WR = ',I2)
99996 FORMAT (' The orthogonal symplectic factor U is ')
99995 FORMAT (/' The factor R is ')
99994 FORMAT (/' The orthogonal symplectic factor V is ')
99993 FORMAT (20(1X,F9.4))
99992 FORMAT (/' N is out of range.',/' N = ',I5)
99991 FORMAT (/' Orthogonality of U: || U^T U - I ||_F = ',G7.2)
99990 FORMAT (/' Residual: || H*V - U*R ||_F = ',G7.2)
99989 FORMAT (/' Orthogonality of V: || V^T V - I ||_F = ',G7.2)
END
Program Data
MB04TB EXAMPLE PROGRAM DATA
5 N N
0.4643 0.3655 0.6853 0.5090 0.3718
0.3688 0.6460 0.4227 0.6798 0.5135
0.7458 0.5043 0.9419 0.9717 0.9990
0.7140 0.4941 0.7802 0.5272 0.1220
0.7418 0.0339 0.7441 0.0436 0.6564
-0.4643 -0.3688 -0.7458 -0.7140 -0.7418
-0.3655 -0.6460 -0.5043 -0.4941 -0.0339
-0.6853 -0.4227 -0.9419 -0.7802 -0.7441
-0.5090 -0.6798 -0.9717 -0.5272 -0.0436
-0.3718 -0.5135 -0.9990 -0.1220 -0.6564
0.7933 1.5765 1.0711 1.0794 0.8481
1.5765 0.1167 1.5685 0.8756 0.5037
1.0711 1.5685 0.9902 0.3858 0.2109
1.0794 0.8756 0.3858 1.8834 1.4338
0.8481 0.5037 0.2109 1.4338 0.1439
1.0786 1.5264 1.1721 1.5343 0.4756
1.5264 0.8644 0.6872 1.1379 0.6499
1.1721 0.6872 1.5194 1.1197 1.0158
1.5343 1.1379 1.1197 0.6612 0.2004
0.4756 0.6499 1.0158 0.2004 1.2188
Program Results
MB04TS EXAMPLE PROGRAM RESULTS
The orthogonal symplectic factor U is
-0.1513 0.0756 -0.0027 0.1694 -0.2999 0.3515 -0.4843 0.6545 -0.1995 -0.1627
-0.1202 0.2320 0.1662 -0.2835 -0.0508 0.4975 0.3319 -0.2686 -0.4186 -0.4649
-0.2431 0.2724 0.3439 0.3954 0.0236 0.3820 -0.2863 -0.4324 0.3706 0.1984
-0.2327 -0.1509 -0.3710 -0.1240 -0.0393 0.5000 0.3659 0.1429 0.0493 0.6015
-0.2418 -0.2928 -0.0836 -0.5549 0.4824 0.1550 -0.4441 -0.0396 0.2376 -0.1702
-0.3515 0.4843 -0.6545 0.1995 0.1627 -0.1513 0.0756 -0.0027 0.1694 -0.2999
-0.4975 -0.3319 0.2686 0.4186 0.4649 -0.1202 0.2320 0.1662 -0.2835 -0.0508
-0.3820 0.2863 0.4324 -0.3706 -0.1984 -0.2431 0.2724 0.3439 0.3954 0.0236
-0.5000 -0.3659 -0.1429 -0.0493 -0.6015 -0.2327 -0.1509 -0.3710 -0.1240 -0.0393
-0.1550 0.4441 0.0396 -0.2376 0.1702 -0.2418 -0.2928 -0.0836 -0.5549 0.4824
Orthogonality of U: || U^T U - I ||_F = .24E-14
The factor R is
-3.0684 4.6724 -0.2613 -0.1996 0.0208 -0.1071 -0.1355 -0.1400 0.4652 -0.5032
0.0000 -1.8037 -0.0301 -0.1137 0.1771 0.0277 0.3929 0.5424 0.5220 -0.4843
0.0000 0.0000 -0.7617 -0.1874 0.2557 0.1244 -0.0012 0.4091 0.5123 -0.3522
0.0000 0.0000 0.0000 -0.6931 -0.4293 -0.3718 0.1542 -0.3635 0.0336 -0.9832
0.0000 0.0000 0.0000 0.0000 0.6469 0.2074 0.0266 0.2028 0.1995 0.2517
0.0000 0.0000 0.0000 0.0000 0.0000 2.6325 -4.7377 0.0000 0.0000 0.0000
0.0000 0.0000 0.0000 0.0000 0.0000 -0.2702 0.9347 -1.1210 0.0000 0.0000
0.0000 0.0000 0.0000 0.0000 0.0000 -0.3219 -0.5394 0.1748 -0.4788 0.0000
0.0000 0.0000 0.0000 0.0000 0.0000 -0.1431 -0.1021 0.4974 -0.3565 -0.6402
0.0000 0.0000 0.0000 0.0000 0.0000 -0.1622 -0.2368 0.6126 -0.7369 0.6915
Residual: || H*V - U*R ||_F = .87E-14
The orthogonal symplectic factor V is
1.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000
0.0000 -0.4740 0.6013 -0.2299 -0.4282 0.0000 0.0061 -0.1732 0.3134 0.2220
0.0000 -0.5553 -0.2623 0.6622 -0.3042 0.0000 -0.0382 0.2453 -0.1662 0.0509
0.0000 -0.5563 0.0322 -0.1431 0.4461 0.0000 -0.0665 -0.4132 -0.3100 -0.4457
0.0000 -0.3872 -0.4022 -0.4194 0.3541 0.0000 -0.0406 0.3820 0.3006 0.3861
0.0000 0.0000 0.0000 0.0000 0.0000 1.0000 0.0000 0.0000 0.0000 0.0000
0.0000 -0.0061 0.1732 -0.3134 -0.2220 0.0000 -0.4740 0.6013 -0.2299 -0.4282
0.0000 0.0382 -0.2453 0.1662 -0.0509 0.0000 -0.5553 -0.2623 0.6622 -0.3042
0.0000 0.0665 0.4132 0.3100 0.4457 0.0000 -0.5563 0.0322 -0.1431 0.4461
0.0000 0.0406 -0.3820 -0.3006 -0.3861 0.0000 -0.3872 -0.4022 -0.4194 0.3541
Orthogonality of V: || V^T V - I ||_F = .14E-14