1 SUBROUTINE DSYT22( ITYPE, UPLO, N, M, KBAND, A, LDA, D, E, U, LDU,
2 $ V, LDV, TAU, 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 CHARACTER UPLO
10 INTEGER ITYPE, KBAND, LDA, LDU, LDV, M, N
11 * ..
12 * .. Array Arguments ..
13 DOUBLE PRECISION A( LDA, * ), D( * ), E( * ), RESULT( 2 ),
14 $ TAU( * ), U( LDU, * ), V( LDV, * ), WORK( * )
15 * ..
16 *
17 * Purpose
18 * =======
19 *
20 * DSYT22 generally checks a decomposition of the form
21 *
22 * A U = U S
23 *
24 * where A is symmetric, the columns of U are orthonormal, and S
25 * is diagonal (if KBAND=0) or symmetric tridiagonal (if
26 * KBAND=1). If ITYPE=1, then U is represented as a dense matrix,
27 * otherwise the U is expressed as a product of Householder
28 * transformations, whose vectors are stored in the array "V" and
29 * whose scaling constants are in "TAU"; we shall use the letter
30 * "V" to refer to the product of Householder transformations
31 * (which should be equal to U).
32 *
33 * Specifically, if ITYPE=1, then:
34 *
35 * RESULT(1) = | U' A U - S | / ( |A| m ulp ) *and*
36 * RESULT(2) = | I - U'U | / ( m ulp )
37 *
38 * Arguments
39 * =========
40 *
41 * ITYPE INTEGER
42 * Specifies the type of tests to be performed.
43 * 1: U expressed as a dense orthogonal matrix:
44 * RESULT(1) = | A - U S U' | / ( |A| n ulp ) *and*
45 * RESULT(2) = | I - UU' | / ( n ulp )
46 *
47 * UPLO CHARACTER
48 * If UPLO='U', the upper triangle of A will be used and the
49 * (strictly) lower triangle will not be referenced. If
50 * UPLO='L', the lower triangle of A will be used and the
51 * (strictly) upper triangle will not be referenced.
52 * Not modified.
53 *
54 * N INTEGER
55 * The size of the matrix. If it is zero, DSYT22 does nothing.
56 * It must be at least zero.
57 * Not modified.
58 *
59 * M INTEGER
60 * The number of columns of U. If it is zero, DSYT22 does
61 * nothing. It must be at least zero.
62 * Not modified.
63 *
64 * KBAND INTEGER
65 * The bandwidth of the matrix. It may only be zero or one.
66 * If zero, then S is diagonal, and E is not referenced. If
67 * one, then S is symmetric tri-diagonal.
68 * Not modified.
69 *
70 * A DOUBLE PRECISION array, dimension (LDA , N)
71 * The original (unfactored) matrix. It is assumed to be
72 * symmetric, and only the upper (UPLO='U') or only the lower
73 * (UPLO='L') will be referenced.
74 * Not modified.
75 *
76 * LDA INTEGER
77 * The leading dimension of A. It must be at least 1
78 * and at least N.
79 * Not modified.
80 *
81 * D DOUBLE PRECISION array, dimension (N)
82 * The diagonal of the (symmetric tri-) diagonal matrix.
83 * Not modified.
84 *
85 * E DOUBLE PRECISION array, dimension (N)
86 * The off-diagonal of the (symmetric tri-) diagonal matrix.
87 * E(1) is ignored, E(2) is the (1,2) and (2,1) element, etc.
88 * Not referenced if KBAND=0.
89 * Not modified.
90 *
91 * U DOUBLE PRECISION array, dimension (LDU, N)
92 * If ITYPE=1 or 3, this contains the orthogonal matrix in
93 * the decomposition, expressed as a dense matrix. If ITYPE=2,
94 * then it is not referenced.
95 * Not modified.
96 *
97 * LDU INTEGER
98 * The leading dimension of U. LDU must be at least N and
99 * at least 1.
100 * Not modified.
101 *
102 * V DOUBLE PRECISION array, dimension (LDV, N)
103 * If ITYPE=2 or 3, the lower triangle of this array contains
104 * the Householder vectors used to describe the orthogonal
105 * matrix in the decomposition. If ITYPE=1, then it is not
106 * referenced.
107 * Not modified.
108 *
109 * LDV INTEGER
110 * The leading dimension of V. LDV must be at least N and
111 * at least 1.
112 * Not modified.
113 *
114 * TAU DOUBLE PRECISION array, dimension (N)
115 * If ITYPE >= 2, then TAU(j) is the scalar factor of
116 * v(j) v(j)' in the Householder transformation H(j) of
117 * the product U = H(1)...H(n-2)
118 * If ITYPE < 2, then TAU is not referenced.
119 * Not modified.
120 *
121 * WORK DOUBLE PRECISION array, dimension (2*N**2)
122 * Workspace.
123 * Modified.
124 *
125 * RESULT DOUBLE PRECISION array, dimension (2)
126 * The values computed by the two tests described above. The
127 * values are currently limited to 1/ulp, to avoid overflow.
128 * RESULT(1) is always modified. RESULT(2) is modified only
129 * if LDU is at least N.
130 * Modified.
131 *
132 * =====================================================================
133 *
134 * .. Parameters ..
135 DOUBLE PRECISION ZERO, ONE
136 PARAMETER ( ZERO = 0.0D0, ONE = 1.0D0 )
137 * ..
138 * .. Local Scalars ..
139 INTEGER J, JJ, JJ1, JJ2, NN, NNP1
140 DOUBLE PRECISION ANORM, ULP, UNFL, WNORM
141 * ..
142 * .. External Functions ..
143 DOUBLE PRECISION DLAMCH, DLANSY
144 EXTERNAL DLAMCH, DLANSY
145 * ..
146 * .. External Subroutines ..
147 EXTERNAL DGEMM, DORT01, DSYMM
148 * ..
149 * .. Intrinsic Functions ..
150 INTRINSIC DBLE, MAX, MIN
151 * ..
152 * .. Executable Statements ..
153 *
154 RESULT( 1 ) = ZERO
155 RESULT( 2 ) = ZERO
156 IF( N.LE.0 .OR. M.LE.0 )
157 $ RETURN
158 *
159 UNFL = DLAMCH( 'Safe minimum' )
160 ULP = DLAMCH( 'Precision' )
161 *
162 * Do Test 1
163 *
164 * Norm of A:
165 *
166 ANORM = MAX( DLANSY( '1', UPLO, N, A, LDA, WORK ), UNFL )
167 *
168 * Compute error matrix:
169 *
170 * ITYPE=1: error = U' A U - S
171 *
172 CALL DSYMM( 'L', UPLO, N, M, ONE, A, LDA, U, LDU, ZERO, WORK, N )
173 NN = N*N
174 NNP1 = NN + 1
175 CALL DGEMM( 'T', 'N', M, M, N, ONE, U, LDU, WORK, N, ZERO,
176 $ WORK( NNP1 ), N )
177 DO 10 J = 1, M
178 JJ = NN + ( J-1 )*N + J
179 WORK( JJ ) = WORK( JJ ) - D( J )
180 10 CONTINUE
181 IF( KBAND.EQ.1 .AND. N.GT.1 ) THEN
182 DO 20 J = 2, M
183 JJ1 = NN + ( J-1 )*N + J - 1
184 JJ2 = NN + ( J-2 )*N + J
185 WORK( JJ1 ) = WORK( JJ1 ) - E( J-1 )
186 WORK( JJ2 ) = WORK( JJ2 ) - E( J-1 )
187 20 CONTINUE
188 END IF
189 WNORM = DLANSY( '1', UPLO, M, WORK( NNP1 ), N, WORK( 1 ) )
190 *
191 IF( ANORM.GT.WNORM ) THEN
192 RESULT( 1 ) = ( WNORM / ANORM ) / ( M*ULP )
193 ELSE
194 IF( ANORM.LT.ONE ) THEN
195 RESULT( 1 ) = ( MIN( WNORM, M*ANORM ) / ANORM ) / ( M*ULP )
196 ELSE
197 RESULT( 1 ) = MIN( WNORM / ANORM, DBLE( M ) ) / ( M*ULP )
198 END IF
199 END IF
200 *
201 * Do Test 2
202 *
203 * Compute U'U - I
204 *
205 IF( ITYPE.EQ.1 )
206 $ CALL DORT01( 'Columns', N, M, U, LDU, WORK, 2*N*N,
207 $ RESULT( 2 ) )
208 *
209 RETURN
210 *
211 * End of DSYT22
212 *
213 END
2 $ V, LDV, TAU, 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 CHARACTER UPLO
10 INTEGER ITYPE, KBAND, LDA, LDU, LDV, M, N
11 * ..
12 * .. Array Arguments ..
13 DOUBLE PRECISION A( LDA, * ), D( * ), E( * ), RESULT( 2 ),
14 $ TAU( * ), U( LDU, * ), V( LDV, * ), WORK( * )
15 * ..
16 *
17 * Purpose
18 * =======
19 *
20 * DSYT22 generally checks a decomposition of the form
21 *
22 * A U = U S
23 *
24 * where A is symmetric, the columns of U are orthonormal, and S
25 * is diagonal (if KBAND=0) or symmetric tridiagonal (if
26 * KBAND=1). If ITYPE=1, then U is represented as a dense matrix,
27 * otherwise the U is expressed as a product of Householder
28 * transformations, whose vectors are stored in the array "V" and
29 * whose scaling constants are in "TAU"; we shall use the letter
30 * "V" to refer to the product of Householder transformations
31 * (which should be equal to U).
32 *
33 * Specifically, if ITYPE=1, then:
34 *
35 * RESULT(1) = | U' A U - S | / ( |A| m ulp ) *and*
36 * RESULT(2) = | I - U'U | / ( m ulp )
37 *
38 * Arguments
39 * =========
40 *
41 * ITYPE INTEGER
42 * Specifies the type of tests to be performed.
43 * 1: U expressed as a dense orthogonal matrix:
44 * RESULT(1) = | A - U S U' | / ( |A| n ulp ) *and*
45 * RESULT(2) = | I - UU' | / ( n ulp )
46 *
47 * UPLO CHARACTER
48 * If UPLO='U', the upper triangle of A will be used and the
49 * (strictly) lower triangle will not be referenced. If
50 * UPLO='L', the lower triangle of A will be used and the
51 * (strictly) upper triangle will not be referenced.
52 * Not modified.
53 *
54 * N INTEGER
55 * The size of the matrix. If it is zero, DSYT22 does nothing.
56 * It must be at least zero.
57 * Not modified.
58 *
59 * M INTEGER
60 * The number of columns of U. If it is zero, DSYT22 does
61 * nothing. It must be at least zero.
62 * Not modified.
63 *
64 * KBAND INTEGER
65 * The bandwidth of the matrix. It may only be zero or one.
66 * If zero, then S is diagonal, and E is not referenced. If
67 * one, then S is symmetric tri-diagonal.
68 * Not modified.
69 *
70 * A DOUBLE PRECISION array, dimension (LDA , N)
71 * The original (unfactored) matrix. It is assumed to be
72 * symmetric, and only the upper (UPLO='U') or only the lower
73 * (UPLO='L') will be referenced.
74 * Not modified.
75 *
76 * LDA INTEGER
77 * The leading dimension of A. It must be at least 1
78 * and at least N.
79 * Not modified.
80 *
81 * D DOUBLE PRECISION array, dimension (N)
82 * The diagonal of the (symmetric tri-) diagonal matrix.
83 * Not modified.
84 *
85 * E DOUBLE PRECISION array, dimension (N)
86 * The off-diagonal of the (symmetric tri-) diagonal matrix.
87 * E(1) is ignored, E(2) is the (1,2) and (2,1) element, etc.
88 * Not referenced if KBAND=0.
89 * Not modified.
90 *
91 * U DOUBLE PRECISION array, dimension (LDU, N)
92 * If ITYPE=1 or 3, this contains the orthogonal matrix in
93 * the decomposition, expressed as a dense matrix. If ITYPE=2,
94 * then it is not referenced.
95 * Not modified.
96 *
97 * LDU INTEGER
98 * The leading dimension of U. LDU must be at least N and
99 * at least 1.
100 * Not modified.
101 *
102 * V DOUBLE PRECISION array, dimension (LDV, N)
103 * If ITYPE=2 or 3, the lower triangle of this array contains
104 * the Householder vectors used to describe the orthogonal
105 * matrix in the decomposition. If ITYPE=1, then it is not
106 * referenced.
107 * Not modified.
108 *
109 * LDV INTEGER
110 * The leading dimension of V. LDV must be at least N and
111 * at least 1.
112 * Not modified.
113 *
114 * TAU DOUBLE PRECISION array, dimension (N)
115 * If ITYPE >= 2, then TAU(j) is the scalar factor of
116 * v(j) v(j)' in the Householder transformation H(j) of
117 * the product U = H(1)...H(n-2)
118 * If ITYPE < 2, then TAU is not referenced.
119 * Not modified.
120 *
121 * WORK DOUBLE PRECISION array, dimension (2*N**2)
122 * Workspace.
123 * Modified.
124 *
125 * RESULT DOUBLE PRECISION array, dimension (2)
126 * The values computed by the two tests described above. The
127 * values are currently limited to 1/ulp, to avoid overflow.
128 * RESULT(1) is always modified. RESULT(2) is modified only
129 * if LDU is at least N.
130 * Modified.
131 *
132 * =====================================================================
133 *
134 * .. Parameters ..
135 DOUBLE PRECISION ZERO, ONE
136 PARAMETER ( ZERO = 0.0D0, ONE = 1.0D0 )
137 * ..
138 * .. Local Scalars ..
139 INTEGER J, JJ, JJ1, JJ2, NN, NNP1
140 DOUBLE PRECISION ANORM, ULP, UNFL, WNORM
141 * ..
142 * .. External Functions ..
143 DOUBLE PRECISION DLAMCH, DLANSY
144 EXTERNAL DLAMCH, DLANSY
145 * ..
146 * .. External Subroutines ..
147 EXTERNAL DGEMM, DORT01, DSYMM
148 * ..
149 * .. Intrinsic Functions ..
150 INTRINSIC DBLE, MAX, MIN
151 * ..
152 * .. Executable Statements ..
153 *
154 RESULT( 1 ) = ZERO
155 RESULT( 2 ) = ZERO
156 IF( N.LE.0 .OR. M.LE.0 )
157 $ RETURN
158 *
159 UNFL = DLAMCH( 'Safe minimum' )
160 ULP = DLAMCH( 'Precision' )
161 *
162 * Do Test 1
163 *
164 * Norm of A:
165 *
166 ANORM = MAX( DLANSY( '1', UPLO, N, A, LDA, WORK ), UNFL )
167 *
168 * Compute error matrix:
169 *
170 * ITYPE=1: error = U' A U - S
171 *
172 CALL DSYMM( 'L', UPLO, N, M, ONE, A, LDA, U, LDU, ZERO, WORK, N )
173 NN = N*N
174 NNP1 = NN + 1
175 CALL DGEMM( 'T', 'N', M, M, N, ONE, U, LDU, WORK, N, ZERO,
176 $ WORK( NNP1 ), N )
177 DO 10 J = 1, M
178 JJ = NN + ( J-1 )*N + J
179 WORK( JJ ) = WORK( JJ ) - D( J )
180 10 CONTINUE
181 IF( KBAND.EQ.1 .AND. N.GT.1 ) THEN
182 DO 20 J = 2, M
183 JJ1 = NN + ( J-1 )*N + J - 1
184 JJ2 = NN + ( J-2 )*N + J
185 WORK( JJ1 ) = WORK( JJ1 ) - E( J-1 )
186 WORK( JJ2 ) = WORK( JJ2 ) - E( J-1 )
187 20 CONTINUE
188 END IF
189 WNORM = DLANSY( '1', UPLO, M, WORK( NNP1 ), N, WORK( 1 ) )
190 *
191 IF( ANORM.GT.WNORM ) THEN
192 RESULT( 1 ) = ( WNORM / ANORM ) / ( M*ULP )
193 ELSE
194 IF( ANORM.LT.ONE ) THEN
195 RESULT( 1 ) = ( MIN( WNORM, M*ANORM ) / ANORM ) / ( M*ULP )
196 ELSE
197 RESULT( 1 ) = MIN( WNORM / ANORM, DBLE( M ) ) / ( M*ULP )
198 END IF
199 END IF
200 *
201 * Do Test 2
202 *
203 * Compute U'U - I
204 *
205 IF( ITYPE.EQ.1 )
206 $ CALL DORT01( 'Columns', N, M, U, LDU, WORK, 2*N*N,
207 $ RESULT( 2 ) )
208 *
209 RETURN
210 *
211 * End of DSYT22
212 *
213 END