1 SUBROUTINE DGET54( N, A, LDA, B, LDB, S, LDS, T, LDT, U, LDU, V,
2 $ LDV, WORK, RESULT )
3 *
4 * -- LAPACK test routine (version 3.1) --
5 * Univ. of Tennessee, Univ. of California Berkeley and NAG Ltd..
6 * November 2006
7 *
8 * .. Scalar Arguments ..
9 INTEGER LDA, LDB, LDS, LDT, LDU, LDV, N
10 DOUBLE PRECISION RESULT
11 * ..
12 * .. Array Arguments ..
13 DOUBLE PRECISION A( LDA, * ), B( LDB, * ), S( LDS, * ),
14 $ T( LDT, * ), U( LDU, * ), V( LDV, * ),
15 $ WORK( * )
16 * ..
17 *
18 * Purpose
19 * =======
20 *
21 * DGET54 checks a generalized decomposition of the form
22 *
23 * A = U*S*V' and B = U*T* V'
24 *
25 * where ' means transpose and U and V are orthogonal.
26 *
27 * Specifically,
28 *
29 * RESULT = ||( A - U*S*V', B - U*T*V' )|| / (||( A, B )||*n*ulp )
30 *
31 * Arguments
32 * =========
33 *
34 * N (input) INTEGER
35 * The size of the matrix. If it is zero, DGET54 does nothing.
36 * It must be at least zero.
37 *
38 * A (input) DOUBLE PRECISION array, dimension (LDA, N)
39 * The original (unfactored) matrix A.
40 *
41 * LDA (input) INTEGER
42 * The leading dimension of A. It must be at least 1
43 * and at least N.
44 *
45 * B (input) DOUBLE PRECISION array, dimension (LDB, N)
46 * The original (unfactored) matrix B.
47 *
48 * LDB (input) INTEGER
49 * The leading dimension of B. It must be at least 1
50 * and at least N.
51 *
52 * S (input) DOUBLE PRECISION array, dimension (LDS, N)
53 * The factored matrix S.
54 *
55 * LDS (input) INTEGER
56 * The leading dimension of S. It must be at least 1
57 * and at least N.
58 *
59 * T (input) DOUBLE PRECISION array, dimension (LDT, N)
60 * The factored matrix T.
61 *
62 * LDT (input) INTEGER
63 * The leading dimension of T. It must be at least 1
64 * and at least N.
65 *
66 * U (input) DOUBLE PRECISION array, dimension (LDU, N)
67 * The orthogonal matrix on the left-hand side in the
68 * decomposition.
69 *
70 * LDU (input) INTEGER
71 * The leading dimension of U. LDU must be at least N and
72 * at least 1.
73 *
74 * V (input) DOUBLE PRECISION array, dimension (LDV, N)
75 * The orthogonal matrix on the left-hand side in the
76 * decomposition.
77 *
78 * LDV (input) INTEGER
79 * The leading dimension of V. LDV must be at least N and
80 * at least 1.
81 *
82 * WORK (workspace) DOUBLE PRECISION array, dimension (3*N**2)
83 *
84 * RESULT (output) DOUBLE PRECISION
85 * The value RESULT, It is currently limited to 1/ulp, to
86 * avoid overflow. Errors are flagged by RESULT=10/ulp.
87 *
88 * =====================================================================
89 *
90 * .. Parameters ..
91 DOUBLE PRECISION ZERO, ONE
92 PARAMETER ( ZERO = 0.0D+0, ONE = 1.0D+0 )
93 * ..
94 * .. Local Scalars ..
95 DOUBLE PRECISION ABNORM, ULP, UNFL, WNORM
96 * ..
97 * .. Local Arrays ..
98 DOUBLE PRECISION DUM( 1 )
99 * ..
100 * .. External Functions ..
101 DOUBLE PRECISION DLAMCH, DLANGE
102 EXTERNAL DLAMCH, DLANGE
103 * ..
104 * .. External Subroutines ..
105 EXTERNAL DGEMM, DLACPY
106 * ..
107 * .. Intrinsic Functions ..
108 INTRINSIC DBLE, MAX, MIN
109 * ..
110 * .. Executable Statements ..
111 *
112 RESULT = ZERO
113 IF( N.LE.0 )
114 $ RETURN
115 *
116 * Constants
117 *
118 UNFL = DLAMCH( 'Safe minimum' )
119 ULP = DLAMCH( 'Epsilon' )*DLAMCH( 'Base' )
120 *
121 * compute the norm of (A,B)
122 *
123 CALL DLACPY( 'Full', N, N, A, LDA, WORK, N )
124 CALL DLACPY( 'Full', N, N, B, LDB, WORK( N*N+1 ), N )
125 ABNORM = MAX( DLANGE( '1', N, 2*N, WORK, N, DUM ), UNFL )
126 *
127 * Compute W1 = A - U*S*V', and put in the array WORK(1:N*N)
128 *
129 CALL DLACPY( ' ', N, N, A, LDA, WORK, N )
130 CALL DGEMM( 'N', 'N', N, N, N, ONE, U, LDU, S, LDS, ZERO,
131 $ WORK( N*N+1 ), N )
132 *
133 CALL DGEMM( 'N', 'C', N, N, N, -ONE, WORK( N*N+1 ), N, V, LDV,
134 $ ONE, WORK, N )
135 *
136 * Compute W2 = B - U*T*V', and put in the workarray W(N*N+1:2*N*N)
137 *
138 CALL DLACPY( ' ', N, N, B, LDB, WORK( N*N+1 ), N )
139 CALL DGEMM( 'N', 'N', N, N, N, ONE, U, LDU, T, LDT, ZERO,
140 $ WORK( 2*N*N+1 ), N )
141 *
142 CALL DGEMM( 'N', 'C', N, N, N, -ONE, WORK( 2*N*N+1 ), N, V, LDV,
143 $ ONE, WORK( N*N+1 ), N )
144 *
145 * Compute norm(W)/ ( ulp*norm((A,B)) )
146 *
147 WNORM = DLANGE( '1', N, 2*N, WORK, N, DUM )
148 *
149 IF( ABNORM.GT.WNORM ) THEN
150 RESULT = ( WNORM / ABNORM ) / ( 2*N*ULP )
151 ELSE
152 IF( ABNORM.LT.ONE ) THEN
153 RESULT = ( MIN( WNORM, 2*N*ABNORM ) / ABNORM ) / ( 2*N*ULP )
154 ELSE
155 RESULT = MIN( WNORM / ABNORM, DBLE( 2*N ) ) / ( 2*N*ULP )
156 END IF
157 END IF
158 *
159 RETURN
160 *
161 * End of DGET54
162 *
163 END
2 $ LDV, WORK, RESULT )
3 *
4 * -- LAPACK test routine (version 3.1) --
5 * Univ. of Tennessee, Univ. of California Berkeley and NAG Ltd..
6 * November 2006
7 *
8 * .. Scalar Arguments ..
9 INTEGER LDA, LDB, LDS, LDT, LDU, LDV, N
10 DOUBLE PRECISION RESULT
11 * ..
12 * .. Array Arguments ..
13 DOUBLE PRECISION A( LDA, * ), B( LDB, * ), S( LDS, * ),
14 $ T( LDT, * ), U( LDU, * ), V( LDV, * ),
15 $ WORK( * )
16 * ..
17 *
18 * Purpose
19 * =======
20 *
21 * DGET54 checks a generalized decomposition of the form
22 *
23 * A = U*S*V' and B = U*T* V'
24 *
25 * where ' means transpose and U and V are orthogonal.
26 *
27 * Specifically,
28 *
29 * RESULT = ||( A - U*S*V', B - U*T*V' )|| / (||( A, B )||*n*ulp )
30 *
31 * Arguments
32 * =========
33 *
34 * N (input) INTEGER
35 * The size of the matrix. If it is zero, DGET54 does nothing.
36 * It must be at least zero.
37 *
38 * A (input) DOUBLE PRECISION array, dimension (LDA, N)
39 * The original (unfactored) matrix A.
40 *
41 * LDA (input) INTEGER
42 * The leading dimension of A. It must be at least 1
43 * and at least N.
44 *
45 * B (input) DOUBLE PRECISION array, dimension (LDB, N)
46 * The original (unfactored) matrix B.
47 *
48 * LDB (input) INTEGER
49 * The leading dimension of B. It must be at least 1
50 * and at least N.
51 *
52 * S (input) DOUBLE PRECISION array, dimension (LDS, N)
53 * The factored matrix S.
54 *
55 * LDS (input) INTEGER
56 * The leading dimension of S. It must be at least 1
57 * and at least N.
58 *
59 * T (input) DOUBLE PRECISION array, dimension (LDT, N)
60 * The factored matrix T.
61 *
62 * LDT (input) INTEGER
63 * The leading dimension of T. It must be at least 1
64 * and at least N.
65 *
66 * U (input) DOUBLE PRECISION array, dimension (LDU, N)
67 * The orthogonal matrix on the left-hand side in the
68 * decomposition.
69 *
70 * LDU (input) INTEGER
71 * The leading dimension of U. LDU must be at least N and
72 * at least 1.
73 *
74 * V (input) DOUBLE PRECISION array, dimension (LDV, N)
75 * The orthogonal matrix on the left-hand side in the
76 * decomposition.
77 *
78 * LDV (input) INTEGER
79 * The leading dimension of V. LDV must be at least N and
80 * at least 1.
81 *
82 * WORK (workspace) DOUBLE PRECISION array, dimension (3*N**2)
83 *
84 * RESULT (output) DOUBLE PRECISION
85 * The value RESULT, It is currently limited to 1/ulp, to
86 * avoid overflow. Errors are flagged by RESULT=10/ulp.
87 *
88 * =====================================================================
89 *
90 * .. Parameters ..
91 DOUBLE PRECISION ZERO, ONE
92 PARAMETER ( ZERO = 0.0D+0, ONE = 1.0D+0 )
93 * ..
94 * .. Local Scalars ..
95 DOUBLE PRECISION ABNORM, ULP, UNFL, WNORM
96 * ..
97 * .. Local Arrays ..
98 DOUBLE PRECISION DUM( 1 )
99 * ..
100 * .. External Functions ..
101 DOUBLE PRECISION DLAMCH, DLANGE
102 EXTERNAL DLAMCH, DLANGE
103 * ..
104 * .. External Subroutines ..
105 EXTERNAL DGEMM, DLACPY
106 * ..
107 * .. Intrinsic Functions ..
108 INTRINSIC DBLE, MAX, MIN
109 * ..
110 * .. Executable Statements ..
111 *
112 RESULT = ZERO
113 IF( N.LE.0 )
114 $ RETURN
115 *
116 * Constants
117 *
118 UNFL = DLAMCH( 'Safe minimum' )
119 ULP = DLAMCH( 'Epsilon' )*DLAMCH( 'Base' )
120 *
121 * compute the norm of (A,B)
122 *
123 CALL DLACPY( 'Full', N, N, A, LDA, WORK, N )
124 CALL DLACPY( 'Full', N, N, B, LDB, WORK( N*N+1 ), N )
125 ABNORM = MAX( DLANGE( '1', N, 2*N, WORK, N, DUM ), UNFL )
126 *
127 * Compute W1 = A - U*S*V', and put in the array WORK(1:N*N)
128 *
129 CALL DLACPY( ' ', N, N, A, LDA, WORK, N )
130 CALL DGEMM( 'N', 'N', N, N, N, ONE, U, LDU, S, LDS, ZERO,
131 $ WORK( N*N+1 ), N )
132 *
133 CALL DGEMM( 'N', 'C', N, N, N, -ONE, WORK( N*N+1 ), N, V, LDV,
134 $ ONE, WORK, N )
135 *
136 * Compute W2 = B - U*T*V', and put in the workarray W(N*N+1:2*N*N)
137 *
138 CALL DLACPY( ' ', N, N, B, LDB, WORK( N*N+1 ), N )
139 CALL DGEMM( 'N', 'N', N, N, N, ONE, U, LDU, T, LDT, ZERO,
140 $ WORK( 2*N*N+1 ), N )
141 *
142 CALL DGEMM( 'N', 'C', N, N, N, -ONE, WORK( 2*N*N+1 ), N, V, LDV,
143 $ ONE, WORK( N*N+1 ), N )
144 *
145 * Compute norm(W)/ ( ulp*norm((A,B)) )
146 *
147 WNORM = DLANGE( '1', N, 2*N, WORK, N, DUM )
148 *
149 IF( ABNORM.GT.WNORM ) THEN
150 RESULT = ( WNORM / ABNORM ) / ( 2*N*ULP )
151 ELSE
152 IF( ABNORM.LT.ONE ) THEN
153 RESULT = ( MIN( WNORM, 2*N*ABNORM ) / ABNORM ) / ( 2*N*ULP )
154 ELSE
155 RESULT = MIN( WNORM / ABNORM, DBLE( 2*N ) ) / ( 2*N*ULP )
156 END IF
157 END IF
158 *
159 RETURN
160 *
161 * End of DGET54
162 *
163 END