1 SUBROUTINE DBDT01( M, N, KD, A, LDA, Q, LDQ, D, E, PT, LDPT, WORK,
2 $ RESID )
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 KD, LDA, LDPT, LDQ, M, N
10 DOUBLE PRECISION RESID
11 * ..
12 * .. Array Arguments ..
13 DOUBLE PRECISION A( LDA, * ), D( * ), E( * ), PT( LDPT, * ),
14 $ Q( LDQ, * ), WORK( * )
15 * ..
16 *
17 * Purpose
18 * =======
19 *
20 * DBDT01 reconstructs a general matrix A from its bidiagonal form
21 * A = Q * B * P'
22 * where Q (m by min(m,n)) and P' (min(m,n) by n) are orthogonal
23 * matrices and B is bidiagonal.
24 *
25 * The test ratio to test the reduction is
26 * RESID = norm( A - Q * B * PT ) / ( n * norm(A) * EPS )
27 * where PT = P' and EPS is the machine precision.
28 *
29 * Arguments
30 * =========
31 *
32 * M (input) INTEGER
33 * The number of rows of the matrices A and Q.
34 *
35 * N (input) INTEGER
36 * The number of columns of the matrices A and P'.
37 *
38 * KD (input) INTEGER
39 * If KD = 0, B is diagonal and the array E is not referenced.
40 * If KD = 1, the reduction was performed by xGEBRD; B is upper
41 * bidiagonal if M >= N, and lower bidiagonal if M < N.
42 * If KD = -1, the reduction was performed by xGBBRD; B is
43 * always upper bidiagonal.
44 *
45 * A (input) DOUBLE PRECISION array, dimension (LDA,N)
46 * The m by n matrix A.
47 *
48 * LDA (input) INTEGER
49 * The leading dimension of the array A. LDA >= max(1,M).
50 *
51 * Q (input) DOUBLE PRECISION array, dimension (LDQ,N)
52 * The m by min(m,n) orthogonal matrix Q in the reduction
53 * A = Q * B * P'.
54 *
55 * LDQ (input) INTEGER
56 * The leading dimension of the array Q. LDQ >= max(1,M).
57 *
58 * D (input) DOUBLE PRECISION array, dimension (min(M,N))
59 * The diagonal elements of the bidiagonal matrix B.
60 *
61 * E (input) DOUBLE PRECISION array, dimension (min(M,N)-1)
62 * The superdiagonal elements of the bidiagonal matrix B if
63 * m >= n, or the subdiagonal elements of B if m < n.
64 *
65 * PT (input) DOUBLE PRECISION array, dimension (LDPT,N)
66 * The min(m,n) by n orthogonal matrix P' in the reduction
67 * A = Q * B * P'.
68 *
69 * LDPT (input) INTEGER
70 * The leading dimension of the array PT.
71 * LDPT >= max(1,min(M,N)).
72 *
73 * WORK (workspace) DOUBLE PRECISION array, dimension (M+N)
74 *
75 * RESID (output) DOUBLE PRECISION
76 * The test ratio: norm(A - Q * B * P') / ( n * norm(A) * EPS )
77 *
78 * =====================================================================
79 *
80 * .. Parameters ..
81 DOUBLE PRECISION ZERO, ONE
82 PARAMETER ( ZERO = 0.0D+0, ONE = 1.0D+0 )
83 * ..
84 * .. Local Scalars ..
85 INTEGER I, J
86 DOUBLE PRECISION ANORM, EPS
87 * ..
88 * .. External Functions ..
89 DOUBLE PRECISION DASUM, DLAMCH, DLANGE
90 EXTERNAL DASUM, DLAMCH, DLANGE
91 * ..
92 * .. External Subroutines ..
93 EXTERNAL DCOPY, DGEMV
94 * ..
95 * .. Intrinsic Functions ..
96 INTRINSIC DBLE, MAX, MIN
97 * ..
98 * .. Executable Statements ..
99 *
100 * Quick return if possible
101 *
102 IF( M.LE.0 .OR. N.LE.0 ) THEN
103 RESID = ZERO
104 RETURN
105 END IF
106 *
107 * Compute A - Q * B * P' one column at a time.
108 *
109 RESID = ZERO
110 IF( KD.NE.0 ) THEN
111 *
112 * B is bidiagonal.
113 *
114 IF( KD.NE.0 .AND. M.GE.N ) THEN
115 *
116 * B is upper bidiagonal and M >= N.
117 *
118 DO 20 J = 1, N
119 CALL DCOPY( M, A( 1, J ), 1, WORK, 1 )
120 DO 10 I = 1, N - 1
121 WORK( M+I ) = D( I )*PT( I, J ) + E( I )*PT( I+1, J )
122 10 CONTINUE
123 WORK( M+N ) = D( N )*PT( N, J )
124 CALL DGEMV( 'No transpose', M, N, -ONE, Q, LDQ,
125 $ WORK( M+1 ), 1, ONE, WORK, 1 )
126 RESID = MAX( RESID, DASUM( M, WORK, 1 ) )
127 20 CONTINUE
128 ELSE IF( KD.LT.0 ) THEN
129 *
130 * B is upper bidiagonal and M < N.
131 *
132 DO 40 J = 1, N
133 CALL DCOPY( M, A( 1, J ), 1, WORK, 1 )
134 DO 30 I = 1, M - 1
135 WORK( M+I ) = D( I )*PT( I, J ) + E( I )*PT( I+1, J )
136 30 CONTINUE
137 WORK( M+M ) = D( M )*PT( M, J )
138 CALL DGEMV( 'No transpose', M, M, -ONE, Q, LDQ,
139 $ WORK( M+1 ), 1, ONE, WORK, 1 )
140 RESID = MAX( RESID, DASUM( M, WORK, 1 ) )
141 40 CONTINUE
142 ELSE
143 *
144 * B is lower bidiagonal.
145 *
146 DO 60 J = 1, N
147 CALL DCOPY( M, A( 1, J ), 1, WORK, 1 )
148 WORK( M+1 ) = D( 1 )*PT( 1, J )
149 DO 50 I = 2, M
150 WORK( M+I ) = E( I-1 )*PT( I-1, J ) +
151 $ D( I )*PT( I, J )
152 50 CONTINUE
153 CALL DGEMV( 'No transpose', M, M, -ONE, Q, LDQ,
154 $ WORK( M+1 ), 1, ONE, WORK, 1 )
155 RESID = MAX( RESID, DASUM( M, WORK, 1 ) )
156 60 CONTINUE
157 END IF
158 ELSE
159 *
160 * B is diagonal.
161 *
162 IF( M.GE.N ) THEN
163 DO 80 J = 1, N
164 CALL DCOPY( M, A( 1, J ), 1, WORK, 1 )
165 DO 70 I = 1, N
166 WORK( M+I ) = D( I )*PT( I, J )
167 70 CONTINUE
168 CALL DGEMV( 'No transpose', M, N, -ONE, Q, LDQ,
169 $ WORK( M+1 ), 1, ONE, WORK, 1 )
170 RESID = MAX( RESID, DASUM( M, WORK, 1 ) )
171 80 CONTINUE
172 ELSE
173 DO 100 J = 1, N
174 CALL DCOPY( M, A( 1, J ), 1, WORK, 1 )
175 DO 90 I = 1, M
176 WORK( M+I ) = D( I )*PT( I, J )
177 90 CONTINUE
178 CALL DGEMV( 'No transpose', M, M, -ONE, Q, LDQ,
179 $ WORK( M+1 ), 1, ONE, WORK, 1 )
180 RESID = MAX( RESID, DASUM( M, WORK, 1 ) )
181 100 CONTINUE
182 END IF
183 END IF
184 *
185 * Compute norm(A - Q * B * P') / ( n * norm(A) * EPS )
186 *
187 ANORM = DLANGE( '1', M, N, A, LDA, WORK )
188 EPS = DLAMCH( 'Precision' )
189 *
190 IF( ANORM.LE.ZERO ) THEN
191 IF( RESID.NE.ZERO )
192 $ RESID = ONE / EPS
193 ELSE
194 IF( ANORM.GE.RESID ) THEN
195 RESID = ( RESID / ANORM ) / ( DBLE( N )*EPS )
196 ELSE
197 IF( ANORM.LT.ONE ) THEN
198 RESID = ( MIN( RESID, DBLE( N )*ANORM ) / ANORM ) /
199 $ ( DBLE( N )*EPS )
200 ELSE
201 RESID = MIN( RESID / ANORM, DBLE( N ) ) /
202 $ ( DBLE( N )*EPS )
203 END IF
204 END IF
205 END IF
206 *
207 RETURN
208 *
209 * End of DBDT01
210 *
211 END
2 $ RESID )
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 KD, LDA, LDPT, LDQ, M, N
10 DOUBLE PRECISION RESID
11 * ..
12 * .. Array Arguments ..
13 DOUBLE PRECISION A( LDA, * ), D( * ), E( * ), PT( LDPT, * ),
14 $ Q( LDQ, * ), WORK( * )
15 * ..
16 *
17 * Purpose
18 * =======
19 *
20 * DBDT01 reconstructs a general matrix A from its bidiagonal form
21 * A = Q * B * P'
22 * where Q (m by min(m,n)) and P' (min(m,n) by n) are orthogonal
23 * matrices and B is bidiagonal.
24 *
25 * The test ratio to test the reduction is
26 * RESID = norm( A - Q * B * PT ) / ( n * norm(A) * EPS )
27 * where PT = P' and EPS is the machine precision.
28 *
29 * Arguments
30 * =========
31 *
32 * M (input) INTEGER
33 * The number of rows of the matrices A and Q.
34 *
35 * N (input) INTEGER
36 * The number of columns of the matrices A and P'.
37 *
38 * KD (input) INTEGER
39 * If KD = 0, B is diagonal and the array E is not referenced.
40 * If KD = 1, the reduction was performed by xGEBRD; B is upper
41 * bidiagonal if M >= N, and lower bidiagonal if M < N.
42 * If KD = -1, the reduction was performed by xGBBRD; B is
43 * always upper bidiagonal.
44 *
45 * A (input) DOUBLE PRECISION array, dimension (LDA,N)
46 * The m by n matrix A.
47 *
48 * LDA (input) INTEGER
49 * The leading dimension of the array A. LDA >= max(1,M).
50 *
51 * Q (input) DOUBLE PRECISION array, dimension (LDQ,N)
52 * The m by min(m,n) orthogonal matrix Q in the reduction
53 * A = Q * B * P'.
54 *
55 * LDQ (input) INTEGER
56 * The leading dimension of the array Q. LDQ >= max(1,M).
57 *
58 * D (input) DOUBLE PRECISION array, dimension (min(M,N))
59 * The diagonal elements of the bidiagonal matrix B.
60 *
61 * E (input) DOUBLE PRECISION array, dimension (min(M,N)-1)
62 * The superdiagonal elements of the bidiagonal matrix B if
63 * m >= n, or the subdiagonal elements of B if m < n.
64 *
65 * PT (input) DOUBLE PRECISION array, dimension (LDPT,N)
66 * The min(m,n) by n orthogonal matrix P' in the reduction
67 * A = Q * B * P'.
68 *
69 * LDPT (input) INTEGER
70 * The leading dimension of the array PT.
71 * LDPT >= max(1,min(M,N)).
72 *
73 * WORK (workspace) DOUBLE PRECISION array, dimension (M+N)
74 *
75 * RESID (output) DOUBLE PRECISION
76 * The test ratio: norm(A - Q * B * P') / ( n * norm(A) * EPS )
77 *
78 * =====================================================================
79 *
80 * .. Parameters ..
81 DOUBLE PRECISION ZERO, ONE
82 PARAMETER ( ZERO = 0.0D+0, ONE = 1.0D+0 )
83 * ..
84 * .. Local Scalars ..
85 INTEGER I, J
86 DOUBLE PRECISION ANORM, EPS
87 * ..
88 * .. External Functions ..
89 DOUBLE PRECISION DASUM, DLAMCH, DLANGE
90 EXTERNAL DASUM, DLAMCH, DLANGE
91 * ..
92 * .. External Subroutines ..
93 EXTERNAL DCOPY, DGEMV
94 * ..
95 * .. Intrinsic Functions ..
96 INTRINSIC DBLE, MAX, MIN
97 * ..
98 * .. Executable Statements ..
99 *
100 * Quick return if possible
101 *
102 IF( M.LE.0 .OR. N.LE.0 ) THEN
103 RESID = ZERO
104 RETURN
105 END IF
106 *
107 * Compute A - Q * B * P' one column at a time.
108 *
109 RESID = ZERO
110 IF( KD.NE.0 ) THEN
111 *
112 * B is bidiagonal.
113 *
114 IF( KD.NE.0 .AND. M.GE.N ) THEN
115 *
116 * B is upper bidiagonal and M >= N.
117 *
118 DO 20 J = 1, N
119 CALL DCOPY( M, A( 1, J ), 1, WORK, 1 )
120 DO 10 I = 1, N - 1
121 WORK( M+I ) = D( I )*PT( I, J ) + E( I )*PT( I+1, J )
122 10 CONTINUE
123 WORK( M+N ) = D( N )*PT( N, J )
124 CALL DGEMV( 'No transpose', M, N, -ONE, Q, LDQ,
125 $ WORK( M+1 ), 1, ONE, WORK, 1 )
126 RESID = MAX( RESID, DASUM( M, WORK, 1 ) )
127 20 CONTINUE
128 ELSE IF( KD.LT.0 ) THEN
129 *
130 * B is upper bidiagonal and M < N.
131 *
132 DO 40 J = 1, N
133 CALL DCOPY( M, A( 1, J ), 1, WORK, 1 )
134 DO 30 I = 1, M - 1
135 WORK( M+I ) = D( I )*PT( I, J ) + E( I )*PT( I+1, J )
136 30 CONTINUE
137 WORK( M+M ) = D( M )*PT( M, J )
138 CALL DGEMV( 'No transpose', M, M, -ONE, Q, LDQ,
139 $ WORK( M+1 ), 1, ONE, WORK, 1 )
140 RESID = MAX( RESID, DASUM( M, WORK, 1 ) )
141 40 CONTINUE
142 ELSE
143 *
144 * B is lower bidiagonal.
145 *
146 DO 60 J = 1, N
147 CALL DCOPY( M, A( 1, J ), 1, WORK, 1 )
148 WORK( M+1 ) = D( 1 )*PT( 1, J )
149 DO 50 I = 2, M
150 WORK( M+I ) = E( I-1 )*PT( I-1, J ) +
151 $ D( I )*PT( I, J )
152 50 CONTINUE
153 CALL DGEMV( 'No transpose', M, M, -ONE, Q, LDQ,
154 $ WORK( M+1 ), 1, ONE, WORK, 1 )
155 RESID = MAX( RESID, DASUM( M, WORK, 1 ) )
156 60 CONTINUE
157 END IF
158 ELSE
159 *
160 * B is diagonal.
161 *
162 IF( M.GE.N ) THEN
163 DO 80 J = 1, N
164 CALL DCOPY( M, A( 1, J ), 1, WORK, 1 )
165 DO 70 I = 1, N
166 WORK( M+I ) = D( I )*PT( I, J )
167 70 CONTINUE
168 CALL DGEMV( 'No transpose', M, N, -ONE, Q, LDQ,
169 $ WORK( M+1 ), 1, ONE, WORK, 1 )
170 RESID = MAX( RESID, DASUM( M, WORK, 1 ) )
171 80 CONTINUE
172 ELSE
173 DO 100 J = 1, N
174 CALL DCOPY( M, A( 1, J ), 1, WORK, 1 )
175 DO 90 I = 1, M
176 WORK( M+I ) = D( I )*PT( I, J )
177 90 CONTINUE
178 CALL DGEMV( 'No transpose', M, M, -ONE, Q, LDQ,
179 $ WORK( M+1 ), 1, ONE, WORK, 1 )
180 RESID = MAX( RESID, DASUM( M, WORK, 1 ) )
181 100 CONTINUE
182 END IF
183 END IF
184 *
185 * Compute norm(A - Q * B * P') / ( n * norm(A) * EPS )
186 *
187 ANORM = DLANGE( '1', M, N, A, LDA, WORK )
188 EPS = DLAMCH( 'Precision' )
189 *
190 IF( ANORM.LE.ZERO ) THEN
191 IF( RESID.NE.ZERO )
192 $ RESID = ONE / EPS
193 ELSE
194 IF( ANORM.GE.RESID ) THEN
195 RESID = ( RESID / ANORM ) / ( DBLE( N )*EPS )
196 ELSE
197 IF( ANORM.LT.ONE ) THEN
198 RESID = ( MIN( RESID, DBLE( N )*ANORM ) / ANORM ) /
199 $ ( DBLE( N )*EPS )
200 ELSE
201 RESID = MIN( RESID / ANORM, DBLE( N ) ) /
202 $ ( DBLE( N )*EPS )
203 END IF
204 END IF
205 END IF
206 *
207 RETURN
208 *
209 * End of DBDT01
210 *
211 END