1 SUBROUTINE DGET51( ITYPE, N, A, LDA, B, LDB, U, LDU, V, LDV, WORK,
2 $ 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 ITYPE, LDA, LDB, LDU, LDV, N
10 DOUBLE PRECISION RESULT
11 * ..
12 * .. Array Arguments ..
13 DOUBLE PRECISION A( LDA, * ), B( LDB, * ), U( LDU, * ),
14 $ V( LDV, * ), WORK( * )
15 * ..
16 *
17 * Purpose
18 * =======
19 *
20 * DGET51 generally checks a decomposition of the form
21 *
22 * A = U B V'
23 *
24 * where ' means transpose and U and V are orthogonal.
25 *
26 * Specifically, if ITYPE=1
27 *
28 * RESULT = | A - U B V' | / ( |A| n ulp )
29 *
30 * If ITYPE=2, then:
31 *
32 * RESULT = | A - B | / ( |A| n ulp )
33 *
34 * If ITYPE=3, then:
35 *
36 * RESULT = | I - UU' | / ( n ulp )
37 *
38 * Arguments
39 * =========
40 *
41 * ITYPE (input) INTEGER
42 * Specifies the type of tests to be performed.
43 * =1: RESULT = | A - U B V' | / ( |A| n ulp )
44 * =2: RESULT = | A - B | / ( |A| n ulp )
45 * =3: RESULT = | I - UU' | / ( n ulp )
46 *
47 * N (input) INTEGER
48 * The size of the matrix. If it is zero, DGET51 does nothing.
49 * It must be at least zero.
50 *
51 * A (input) DOUBLE PRECISION array, dimension (LDA, N)
52 * The original (unfactored) matrix.
53 *
54 * LDA (input) INTEGER
55 * The leading dimension of A. It must be at least 1
56 * and at least N.
57 *
58 * B (input) DOUBLE PRECISION array, dimension (LDB, N)
59 * The factored matrix.
60 *
61 * LDB (input) INTEGER
62 * The leading dimension of B. It must be at least 1
63 * and at least N.
64 *
65 * U (input) DOUBLE PRECISION array, dimension (LDU, N)
66 * The orthogonal matrix on the left-hand side in the
67 * decomposition.
68 * Not referenced if ITYPE=2
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 * Not referenced if ITYPE=2
78 *
79 * LDV (input) INTEGER
80 * The leading dimension of V. LDV must be at least N and
81 * at least 1.
82 *
83 * WORK (workspace) DOUBLE PRECISION array, dimension (2*N**2)
84 *
85 * RESULT (output) DOUBLE PRECISION
86 * The values computed by the test specified by ITYPE. The
87 * value is currently limited to 1/ulp, to avoid overflow.
88 * Errors are flagged by RESULT=10/ulp.
89 *
90 * =====================================================================
91 *
92 * .. Parameters ..
93 DOUBLE PRECISION ZERO, ONE, TEN
94 PARAMETER ( ZERO = 0.0D0, ONE = 1.0D0, TEN = 10.0D0 )
95 * ..
96 * .. Local Scalars ..
97 INTEGER JCOL, JDIAG, JROW
98 DOUBLE PRECISION ANORM, ULP, UNFL, WNORM
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 * Some Error Checks
122 *
123 IF( ITYPE.LT.1 .OR. ITYPE.GT.3 ) THEN
124 RESULT = TEN / ULP
125 RETURN
126 END IF
127 *
128 IF( ITYPE.LE.2 ) THEN
129 *
130 * Tests scaled by the norm(A)
131 *
132 ANORM = MAX( DLANGE( '1', N, N, A, LDA, WORK ), UNFL )
133 *
134 IF( ITYPE.EQ.1 ) THEN
135 *
136 * ITYPE=1: Compute W = A - UBV'
137 *
138 CALL DLACPY( ' ', N, N, A, LDA, WORK, N )
139 CALL DGEMM( 'N', 'N', N, N, N, ONE, U, LDU, B, LDB, ZERO,
140 $ WORK( N**2+1 ), N )
141 *
142 CALL DGEMM( 'N', 'C', N, N, N, -ONE, WORK( N**2+1 ), N, V,
143 $ LDV, ONE, WORK, N )
144 *
145 ELSE
146 *
147 * ITYPE=2: Compute W = A - B
148 *
149 CALL DLACPY( ' ', N, N, B, LDB, WORK, N )
150 *
151 DO 20 JCOL = 1, N
152 DO 10 JROW = 1, N
153 WORK( JROW+N*( JCOL-1 ) ) = WORK( JROW+N*( JCOL-1 ) )
154 $ - A( JROW, JCOL )
155 10 CONTINUE
156 20 CONTINUE
157 END IF
158 *
159 * Compute norm(W)/ ( ulp*norm(A) )
160 *
161 WNORM = DLANGE( '1', N, N, WORK, N, WORK( N**2+1 ) )
162 *
163 IF( ANORM.GT.WNORM ) THEN
164 RESULT = ( WNORM / ANORM ) / ( N*ULP )
165 ELSE
166 IF( ANORM.LT.ONE ) THEN
167 RESULT = ( MIN( WNORM, N*ANORM ) / ANORM ) / ( N*ULP )
168 ELSE
169 RESULT = MIN( WNORM / ANORM, DBLE( N ) ) / ( N*ULP )
170 END IF
171 END IF
172 *
173 ELSE
174 *
175 * Tests not scaled by norm(A)
176 *
177 * ITYPE=3: Compute UU' - I
178 *
179 CALL DGEMM( 'N', 'C', N, N, N, ONE, U, LDU, U, LDU, ZERO, WORK,
180 $ N )
181 *
182 DO 30 JDIAG = 1, N
183 WORK( ( N+1 )*( JDIAG-1 )+1 ) = WORK( ( N+1 )*( JDIAG-1 )+
184 $ 1 ) - ONE
185 30 CONTINUE
186 *
187 RESULT = MIN( DLANGE( '1', N, N, WORK, N, WORK( N**2+1 ) ),
188 $ DBLE( N ) ) / ( N*ULP )
189 END IF
190 *
191 RETURN
192 *
193 * End of DGET51
194 *
195 END
2 $ 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 ITYPE, LDA, LDB, LDU, LDV, N
10 DOUBLE PRECISION RESULT
11 * ..
12 * .. Array Arguments ..
13 DOUBLE PRECISION A( LDA, * ), B( LDB, * ), U( LDU, * ),
14 $ V( LDV, * ), WORK( * )
15 * ..
16 *
17 * Purpose
18 * =======
19 *
20 * DGET51 generally checks a decomposition of the form
21 *
22 * A = U B V'
23 *
24 * where ' means transpose and U and V are orthogonal.
25 *
26 * Specifically, if ITYPE=1
27 *
28 * RESULT = | A - U B V' | / ( |A| n ulp )
29 *
30 * If ITYPE=2, then:
31 *
32 * RESULT = | A - B | / ( |A| n ulp )
33 *
34 * If ITYPE=3, then:
35 *
36 * RESULT = | I - UU' | / ( n ulp )
37 *
38 * Arguments
39 * =========
40 *
41 * ITYPE (input) INTEGER
42 * Specifies the type of tests to be performed.
43 * =1: RESULT = | A - U B V' | / ( |A| n ulp )
44 * =2: RESULT = | A - B | / ( |A| n ulp )
45 * =3: RESULT = | I - UU' | / ( n ulp )
46 *
47 * N (input) INTEGER
48 * The size of the matrix. If it is zero, DGET51 does nothing.
49 * It must be at least zero.
50 *
51 * A (input) DOUBLE PRECISION array, dimension (LDA, N)
52 * The original (unfactored) matrix.
53 *
54 * LDA (input) INTEGER
55 * The leading dimension of A. It must be at least 1
56 * and at least N.
57 *
58 * B (input) DOUBLE PRECISION array, dimension (LDB, N)
59 * The factored matrix.
60 *
61 * LDB (input) INTEGER
62 * The leading dimension of B. It must be at least 1
63 * and at least N.
64 *
65 * U (input) DOUBLE PRECISION array, dimension (LDU, N)
66 * The orthogonal matrix on the left-hand side in the
67 * decomposition.
68 * Not referenced if ITYPE=2
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 * Not referenced if ITYPE=2
78 *
79 * LDV (input) INTEGER
80 * The leading dimension of V. LDV must be at least N and
81 * at least 1.
82 *
83 * WORK (workspace) DOUBLE PRECISION array, dimension (2*N**2)
84 *
85 * RESULT (output) DOUBLE PRECISION
86 * The values computed by the test specified by ITYPE. The
87 * value is currently limited to 1/ulp, to avoid overflow.
88 * Errors are flagged by RESULT=10/ulp.
89 *
90 * =====================================================================
91 *
92 * .. Parameters ..
93 DOUBLE PRECISION ZERO, ONE, TEN
94 PARAMETER ( ZERO = 0.0D0, ONE = 1.0D0, TEN = 10.0D0 )
95 * ..
96 * .. Local Scalars ..
97 INTEGER JCOL, JDIAG, JROW
98 DOUBLE PRECISION ANORM, ULP, UNFL, WNORM
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 * Some Error Checks
122 *
123 IF( ITYPE.LT.1 .OR. ITYPE.GT.3 ) THEN
124 RESULT = TEN / ULP
125 RETURN
126 END IF
127 *
128 IF( ITYPE.LE.2 ) THEN
129 *
130 * Tests scaled by the norm(A)
131 *
132 ANORM = MAX( DLANGE( '1', N, N, A, LDA, WORK ), UNFL )
133 *
134 IF( ITYPE.EQ.1 ) THEN
135 *
136 * ITYPE=1: Compute W = A - UBV'
137 *
138 CALL DLACPY( ' ', N, N, A, LDA, WORK, N )
139 CALL DGEMM( 'N', 'N', N, N, N, ONE, U, LDU, B, LDB, ZERO,
140 $ WORK( N**2+1 ), N )
141 *
142 CALL DGEMM( 'N', 'C', N, N, N, -ONE, WORK( N**2+1 ), N, V,
143 $ LDV, ONE, WORK, N )
144 *
145 ELSE
146 *
147 * ITYPE=2: Compute W = A - B
148 *
149 CALL DLACPY( ' ', N, N, B, LDB, WORK, N )
150 *
151 DO 20 JCOL = 1, N
152 DO 10 JROW = 1, N
153 WORK( JROW+N*( JCOL-1 ) ) = WORK( JROW+N*( JCOL-1 ) )
154 $ - A( JROW, JCOL )
155 10 CONTINUE
156 20 CONTINUE
157 END IF
158 *
159 * Compute norm(W)/ ( ulp*norm(A) )
160 *
161 WNORM = DLANGE( '1', N, N, WORK, N, WORK( N**2+1 ) )
162 *
163 IF( ANORM.GT.WNORM ) THEN
164 RESULT = ( WNORM / ANORM ) / ( N*ULP )
165 ELSE
166 IF( ANORM.LT.ONE ) THEN
167 RESULT = ( MIN( WNORM, N*ANORM ) / ANORM ) / ( N*ULP )
168 ELSE
169 RESULT = MIN( WNORM / ANORM, DBLE( N ) ) / ( N*ULP )
170 END IF
171 END IF
172 *
173 ELSE
174 *
175 * Tests not scaled by norm(A)
176 *
177 * ITYPE=3: Compute UU' - I
178 *
179 CALL DGEMM( 'N', 'C', N, N, N, ONE, U, LDU, U, LDU, ZERO, WORK,
180 $ N )
181 *
182 DO 30 JDIAG = 1, N
183 WORK( ( N+1 )*( JDIAG-1 )+1 ) = WORK( ( N+1 )*( JDIAG-1 )+
184 $ 1 ) - ONE
185 30 CONTINUE
186 *
187 RESULT = MIN( DLANGE( '1', N, N, WORK, N, WORK( N**2+1 ) ),
188 $ DBLE( N ) ) / ( N*ULP )
189 END IF
190 *
191 RETURN
192 *
193 * End of DGET51
194 *
195 END