1 SUBROUTINE ZBDT03( UPLO, N, KD, D, E, U, LDU, S, VT, LDVT, 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 CHARACTER UPLO
10 INTEGER KD, LDU, LDVT, N
11 DOUBLE PRECISION RESID
12 * ..
13 * .. Array Arguments ..
14 DOUBLE PRECISION D( * ), E( * ), S( * )
15 COMPLEX*16 U( LDU, * ), VT( LDVT, * ), WORK( * )
16 * ..
17 *
18 * Purpose
19 * =======
20 *
21 * ZBDT03 reconstructs a bidiagonal matrix B from its SVD:
22 * S = U' * B * V
23 * where U and V are orthogonal matrices and S is diagonal.
24 *
25 * The test ratio to test the singular value decomposition is
26 * RESID = norm( B - U * S * VT ) / ( n * norm(B) * EPS )
27 * where VT = V' and EPS is the machine precision.
28 *
29 * Arguments
30 * =========
31 *
32 * UPLO (input) CHARACTER*1
33 * Specifies whether the matrix B is upper or lower bidiagonal.
34 * = 'U': Upper bidiagonal
35 * = 'L': Lower bidiagonal
36 *
37 * N (input) INTEGER
38 * The order of the matrix B.
39 *
40 * KD (input) INTEGER
41 * The bandwidth of the bidiagonal matrix B. If KD = 1, the
42 * matrix B is bidiagonal, and if KD = 0, B is diagonal and E is
43 * not referenced. If KD is greater than 1, it is assumed to be
44 * 1, and if KD is less than 0, it is assumed to be 0.
45 *
46 * D (input) DOUBLE PRECISION array, dimension (N)
47 * The n diagonal elements of the bidiagonal matrix B.
48 *
49 * E (input) DOUBLE PRECISION array, dimension (N-1)
50 * The (n-1) superdiagonal elements of the bidiagonal matrix B
51 * if UPLO = 'U', or the (n-1) subdiagonal elements of B if
52 * UPLO = 'L'.
53 *
54 * U (input) COMPLEX*16 array, dimension (LDU,N)
55 * The n by n orthogonal matrix U in the reduction B = U'*A*P.
56 *
57 * LDU (input) INTEGER
58 * The leading dimension of the array U. LDU >= max(1,N)
59 *
60 * S (input) DOUBLE PRECISION array, dimension (N)
61 * The singular values from the SVD of B, sorted in decreasing
62 * order.
63 *
64 * VT (input) COMPLEX*16 array, dimension (LDVT,N)
65 * The n by n orthogonal matrix V' in the reduction
66 * B = U * S * V'.
67 *
68 * LDVT (input) INTEGER
69 * The leading dimension of the array VT.
70 *
71 * WORK (workspace) COMPLEX*16 array, dimension (2*N)
72 *
73 * RESID (output) DOUBLE PRECISION
74 * The test ratio: norm(B - U * S * V') / ( n * norm(A) * EPS )
75 *
76 * ======================================================================
77 *
78 * .. Parameters ..
79 DOUBLE PRECISION ZERO, ONE
80 PARAMETER ( ZERO = 0.0D+0, ONE = 1.0D+0 )
81 * ..
82 * .. Local Scalars ..
83 INTEGER I, J
84 DOUBLE PRECISION BNORM, EPS
85 * ..
86 * .. External Functions ..
87 LOGICAL LSAME
88 INTEGER IDAMAX
89 DOUBLE PRECISION DLAMCH, DZASUM
90 EXTERNAL LSAME, IDAMAX, DLAMCH, DZASUM
91 * ..
92 * .. External Subroutines ..
93 EXTERNAL ZGEMV
94 * ..
95 * .. Intrinsic Functions ..
96 INTRINSIC ABS, DBLE, DCMPLX, MAX, MIN
97 * ..
98 * .. Executable Statements ..
99 *
100 * Quick return if possible
101 *
102 RESID = ZERO
103 IF( N.LE.0 )
104 $ RETURN
105 *
106 * Compute B - U * S * V' one column at a time.
107 *
108 BNORM = ZERO
109 IF( KD.GE.1 ) THEN
110 *
111 * B is bidiagonal.
112 *
113 IF( LSAME( UPLO, 'U' ) ) THEN
114 *
115 * B is upper bidiagonal.
116 *
117 DO 20 J = 1, N
118 DO 10 I = 1, N
119 WORK( N+I ) = S( I )*VT( I, J )
120 10 CONTINUE
121 CALL ZGEMV( 'No transpose', N, N, -DCMPLX( ONE ), U, LDU,
122 $ WORK( N+1 ), 1, DCMPLX( ZERO ), WORK, 1 )
123 WORK( J ) = WORK( J ) + D( J )
124 IF( J.GT.1 ) THEN
125 WORK( J-1 ) = WORK( J-1 ) + E( J-1 )
126 BNORM = MAX( BNORM, ABS( D( J ) )+ABS( E( J-1 ) ) )
127 ELSE
128 BNORM = MAX( BNORM, ABS( D( J ) ) )
129 END IF
130 RESID = MAX( RESID, DZASUM( N, WORK, 1 ) )
131 20 CONTINUE
132 ELSE
133 *
134 * B is lower bidiagonal.
135 *
136 DO 40 J = 1, N
137 DO 30 I = 1, N
138 WORK( N+I ) = S( I )*VT( I, J )
139 30 CONTINUE
140 CALL ZGEMV( 'No transpose', N, N, -DCMPLX( ONE ), U, LDU,
141 $ WORK( N+1 ), 1, DCMPLX( ZERO ), WORK, 1 )
142 WORK( J ) = WORK( J ) + D( J )
143 IF( J.LT.N ) THEN
144 WORK( J+1 ) = WORK( J+1 ) + E( J )
145 BNORM = MAX( BNORM, ABS( D( J ) )+ABS( E( J ) ) )
146 ELSE
147 BNORM = MAX( BNORM, ABS( D( J ) ) )
148 END IF
149 RESID = MAX( RESID, DZASUM( N, WORK, 1 ) )
150 40 CONTINUE
151 END IF
152 ELSE
153 *
154 * B is diagonal.
155 *
156 DO 60 J = 1, N
157 DO 50 I = 1, N
158 WORK( N+I ) = S( I )*VT( I, J )
159 50 CONTINUE
160 CALL ZGEMV( 'No transpose', N, N, -DCMPLX( ONE ), U, LDU,
161 $ WORK( N+1 ), 1, DCMPLX( ZERO ), WORK, 1 )
162 WORK( J ) = WORK( J ) + D( J )
163 RESID = MAX( RESID, DZASUM( N, WORK, 1 ) )
164 60 CONTINUE
165 J = IDAMAX( N, D, 1 )
166 BNORM = ABS( D( J ) )
167 END IF
168 *
169 * Compute norm(B - U * S * V') / ( n * norm(B) * EPS )
170 *
171 EPS = DLAMCH( 'Precision' )
172 *
173 IF( BNORM.LE.ZERO ) THEN
174 IF( RESID.NE.ZERO )
175 $ RESID = ONE / EPS
176 ELSE
177 IF( BNORM.GE.RESID ) THEN
178 RESID = ( RESID / BNORM ) / ( DBLE( N )*EPS )
179 ELSE
180 IF( BNORM.LT.ONE ) THEN
181 RESID = ( MIN( RESID, DBLE( N )*BNORM ) / BNORM ) /
182 $ ( DBLE( N )*EPS )
183 ELSE
184 RESID = MIN( RESID / BNORM, DBLE( N ) ) /
185 $ ( DBLE( N )*EPS )
186 END IF
187 END IF
188 END IF
189 *
190 RETURN
191 *
192 * End of ZBDT03
193 *
194 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 CHARACTER UPLO
10 INTEGER KD, LDU, LDVT, N
11 DOUBLE PRECISION RESID
12 * ..
13 * .. Array Arguments ..
14 DOUBLE PRECISION D( * ), E( * ), S( * )
15 COMPLEX*16 U( LDU, * ), VT( LDVT, * ), WORK( * )
16 * ..
17 *
18 * Purpose
19 * =======
20 *
21 * ZBDT03 reconstructs a bidiagonal matrix B from its SVD:
22 * S = U' * B * V
23 * where U and V are orthogonal matrices and S is diagonal.
24 *
25 * The test ratio to test the singular value decomposition is
26 * RESID = norm( B - U * S * VT ) / ( n * norm(B) * EPS )
27 * where VT = V' and EPS is the machine precision.
28 *
29 * Arguments
30 * =========
31 *
32 * UPLO (input) CHARACTER*1
33 * Specifies whether the matrix B is upper or lower bidiagonal.
34 * = 'U': Upper bidiagonal
35 * = 'L': Lower bidiagonal
36 *
37 * N (input) INTEGER
38 * The order of the matrix B.
39 *
40 * KD (input) INTEGER
41 * The bandwidth of the bidiagonal matrix B. If KD = 1, the
42 * matrix B is bidiagonal, and if KD = 0, B is diagonal and E is
43 * not referenced. If KD is greater than 1, it is assumed to be
44 * 1, and if KD is less than 0, it is assumed to be 0.
45 *
46 * D (input) DOUBLE PRECISION array, dimension (N)
47 * The n diagonal elements of the bidiagonal matrix B.
48 *
49 * E (input) DOUBLE PRECISION array, dimension (N-1)
50 * The (n-1) superdiagonal elements of the bidiagonal matrix B
51 * if UPLO = 'U', or the (n-1) subdiagonal elements of B if
52 * UPLO = 'L'.
53 *
54 * U (input) COMPLEX*16 array, dimension (LDU,N)
55 * The n by n orthogonal matrix U in the reduction B = U'*A*P.
56 *
57 * LDU (input) INTEGER
58 * The leading dimension of the array U. LDU >= max(1,N)
59 *
60 * S (input) DOUBLE PRECISION array, dimension (N)
61 * The singular values from the SVD of B, sorted in decreasing
62 * order.
63 *
64 * VT (input) COMPLEX*16 array, dimension (LDVT,N)
65 * The n by n orthogonal matrix V' in the reduction
66 * B = U * S * V'.
67 *
68 * LDVT (input) INTEGER
69 * The leading dimension of the array VT.
70 *
71 * WORK (workspace) COMPLEX*16 array, dimension (2*N)
72 *
73 * RESID (output) DOUBLE PRECISION
74 * The test ratio: norm(B - U * S * V') / ( n * norm(A) * EPS )
75 *
76 * ======================================================================
77 *
78 * .. Parameters ..
79 DOUBLE PRECISION ZERO, ONE
80 PARAMETER ( ZERO = 0.0D+0, ONE = 1.0D+0 )
81 * ..
82 * .. Local Scalars ..
83 INTEGER I, J
84 DOUBLE PRECISION BNORM, EPS
85 * ..
86 * .. External Functions ..
87 LOGICAL LSAME
88 INTEGER IDAMAX
89 DOUBLE PRECISION DLAMCH, DZASUM
90 EXTERNAL LSAME, IDAMAX, DLAMCH, DZASUM
91 * ..
92 * .. External Subroutines ..
93 EXTERNAL ZGEMV
94 * ..
95 * .. Intrinsic Functions ..
96 INTRINSIC ABS, DBLE, DCMPLX, MAX, MIN
97 * ..
98 * .. Executable Statements ..
99 *
100 * Quick return if possible
101 *
102 RESID = ZERO
103 IF( N.LE.0 )
104 $ RETURN
105 *
106 * Compute B - U * S * V' one column at a time.
107 *
108 BNORM = ZERO
109 IF( KD.GE.1 ) THEN
110 *
111 * B is bidiagonal.
112 *
113 IF( LSAME( UPLO, 'U' ) ) THEN
114 *
115 * B is upper bidiagonal.
116 *
117 DO 20 J = 1, N
118 DO 10 I = 1, N
119 WORK( N+I ) = S( I )*VT( I, J )
120 10 CONTINUE
121 CALL ZGEMV( 'No transpose', N, N, -DCMPLX( ONE ), U, LDU,
122 $ WORK( N+1 ), 1, DCMPLX( ZERO ), WORK, 1 )
123 WORK( J ) = WORK( J ) + D( J )
124 IF( J.GT.1 ) THEN
125 WORK( J-1 ) = WORK( J-1 ) + E( J-1 )
126 BNORM = MAX( BNORM, ABS( D( J ) )+ABS( E( J-1 ) ) )
127 ELSE
128 BNORM = MAX( BNORM, ABS( D( J ) ) )
129 END IF
130 RESID = MAX( RESID, DZASUM( N, WORK, 1 ) )
131 20 CONTINUE
132 ELSE
133 *
134 * B is lower bidiagonal.
135 *
136 DO 40 J = 1, N
137 DO 30 I = 1, N
138 WORK( N+I ) = S( I )*VT( I, J )
139 30 CONTINUE
140 CALL ZGEMV( 'No transpose', N, N, -DCMPLX( ONE ), U, LDU,
141 $ WORK( N+1 ), 1, DCMPLX( ZERO ), WORK, 1 )
142 WORK( J ) = WORK( J ) + D( J )
143 IF( J.LT.N ) THEN
144 WORK( J+1 ) = WORK( J+1 ) + E( J )
145 BNORM = MAX( BNORM, ABS( D( J ) )+ABS( E( J ) ) )
146 ELSE
147 BNORM = MAX( BNORM, ABS( D( J ) ) )
148 END IF
149 RESID = MAX( RESID, DZASUM( N, WORK, 1 ) )
150 40 CONTINUE
151 END IF
152 ELSE
153 *
154 * B is diagonal.
155 *
156 DO 60 J = 1, N
157 DO 50 I = 1, N
158 WORK( N+I ) = S( I )*VT( I, J )
159 50 CONTINUE
160 CALL ZGEMV( 'No transpose', N, N, -DCMPLX( ONE ), U, LDU,
161 $ WORK( N+1 ), 1, DCMPLX( ZERO ), WORK, 1 )
162 WORK( J ) = WORK( J ) + D( J )
163 RESID = MAX( RESID, DZASUM( N, WORK, 1 ) )
164 60 CONTINUE
165 J = IDAMAX( N, D, 1 )
166 BNORM = ABS( D( J ) )
167 END IF
168 *
169 * Compute norm(B - U * S * V') / ( n * norm(B) * EPS )
170 *
171 EPS = DLAMCH( 'Precision' )
172 *
173 IF( BNORM.LE.ZERO ) THEN
174 IF( RESID.NE.ZERO )
175 $ RESID = ONE / EPS
176 ELSE
177 IF( BNORM.GE.RESID ) THEN
178 RESID = ( RESID / BNORM ) / ( DBLE( N )*EPS )
179 ELSE
180 IF( BNORM.LT.ONE ) THEN
181 RESID = ( MIN( RESID, DBLE( N )*BNORM ) / BNORM ) /
182 $ ( DBLE( N )*EPS )
183 ELSE
184 RESID = MIN( RESID / BNORM, DBLE( N ) ) /
185 $ ( DBLE( N )*EPS )
186 END IF
187 END IF
188 END IF
189 *
190 RETURN
191 *
192 * End of ZBDT03
193 *
194 END