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