1*> \brief \b DBDT03
2*
3*  =========== DOCUMENTATION ===========
4*
5* Online html documentation available at
6*            http://www.netlib.org/lapack/explore-html/
7*
8*  Definition:
9*  ===========
10*
11*       SUBROUTINE DBDT03( UPLO, N, KD, D, E, U, LDU, S, VT, LDVT, WORK,
12*                          RESID )
13*
14*       .. Scalar Arguments ..
15*       CHARACTER          UPLO
16*       INTEGER            KD, LDU, LDVT, N
17*       DOUBLE PRECISION   RESID
18*       ..
19*       .. Array Arguments ..
20*       DOUBLE PRECISION   D( * ), E( * ), S( * ), U( LDU, * ),
21*      $                   VT( LDVT, * ), WORK( * )
22*       ..
23*
24*
25*> \par Purpose:
26*  =============
27*>
28*> \verbatim
29*>
30*> DBDT03 reconstructs a bidiagonal matrix B from its SVD:
31*>    S = U' * B * V
32*> where U and V are orthogonal matrices and S is diagonal.
33*>
34*> The test ratio to test the singular value decomposition is
35*>    RESID = norm( B - U * S * VT ) / ( n * norm(B) * EPS )
36*> where VT = V' and EPS is the machine precision.
37*> \endverbatim
38*
39*  Arguments:
40*  ==========
41*
42*> \param[in] UPLO
43*> \verbatim
44*>          UPLO is CHARACTER*1
45*>          Specifies whether the matrix B is upper or lower bidiagonal.
46*>          = 'U':  Upper bidiagonal
47*>          = 'L':  Lower bidiagonal
48*> \endverbatim
49*>
50*> \param[in] N
51*> \verbatim
52*>          N is INTEGER
53*>          The order of the matrix B.
54*> \endverbatim
55*>
56*> \param[in] KD
57*> \verbatim
58*>          KD is INTEGER
59*>          The bandwidth of the bidiagonal matrix B.  If KD = 1, the
60*>          matrix B is bidiagonal, and if KD = 0, B is diagonal and E is
61*>          not referenced.  If KD is greater than 1, it is assumed to be
62*>          1, and if KD is less than 0, it is assumed to be 0.
63*> \endverbatim
64*>
65*> \param[in] D
66*> \verbatim
67*>          D is DOUBLE PRECISION array, dimension (N)
68*>          The n diagonal elements of the bidiagonal matrix B.
69*> \endverbatim
70*>
71*> \param[in] E
72*> \verbatim
73*>          E is DOUBLE PRECISION array, dimension (N-1)
74*>          The (n-1) superdiagonal elements of the bidiagonal matrix B
75*>          if UPLO = 'U', or the (n-1) subdiagonal elements of B if
76*>          UPLO = 'L'.
77*> \endverbatim
78*>
79*> \param[in] U
80*> \verbatim
81*>          U is DOUBLE PRECISION array, dimension (LDU,N)
82*>          The n by n orthogonal matrix U in the reduction B = U'*A*P.
83*> \endverbatim
84*>
85*> \param[in] LDU
86*> \verbatim
87*>          LDU is INTEGER
88*>          The leading dimension of the array U.  LDU >= max(1,N)
89*> \endverbatim
90*>
91*> \param[in] S
92*> \verbatim
93*>          S is DOUBLE PRECISION array, dimension (N)
94*>          The singular values from the SVD of B, sorted in decreasing
95*>          order.
96*> \endverbatim
97*>
98*> \param[in] VT
99*> \verbatim
100*>          VT is DOUBLE PRECISION array, dimension (LDVT,N)
101*>          The n by n orthogonal matrix V' in the reduction
102*>          B = U * S * V'.
103*> \endverbatim
104*>
105*> \param[in] LDVT
106*> \verbatim
107*>          LDVT is INTEGER
108*>          The leading dimension of the array VT.
109*> \endverbatim
110*>
111*> \param[out] WORK
112*> \verbatim
113*>          WORK is DOUBLE PRECISION array, dimension (2*N)
114*> \endverbatim
115*>
116*> \param[out] RESID
117*> \verbatim
118*>          RESID is DOUBLE PRECISION
119*>          The test ratio:  norm(B - U * S * V') / ( n * norm(A) * EPS )
120*> \endverbatim
121*
122*  Authors:
123*  ========
124*
125*> \author Univ. of Tennessee
126*> \author Univ. of California Berkeley
127*> \author Univ. of Colorado Denver
128*> \author NAG Ltd.
129*
130*> \ingroup double_eig
131*
132*  =====================================================================
133      SUBROUTINE DBDT03( UPLO, N, KD, D, E, U, LDU, S, VT, LDVT, WORK,
134     $                   RESID )
135*
136*  -- LAPACK test routine --
137*  -- LAPACK is a software package provided by Univ. of Tennessee,    --
138*  -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
139*
140*     .. Scalar Arguments ..
141      CHARACTER          UPLO
142      INTEGER            KD, LDU, LDVT, N
143      DOUBLE PRECISION   RESID
144*     ..
145*     .. Array Arguments ..
146      DOUBLE PRECISION   D( * ), E( * ), S( * ), U( LDU, * ),
147     $                   VT( LDVT, * ), WORK( * )
148*     ..
149*
150* ======================================================================
151*
152*     .. Parameters ..
153      DOUBLE PRECISION   ZERO, ONE
154      PARAMETER          ( ZERO = 0.0D+0, ONE = 1.0D+0 )
155*     ..
156*     .. Local Scalars ..
157      INTEGER            I, J
158      DOUBLE PRECISION   BNORM, EPS
159*     ..
160*     .. External Functions ..
161      LOGICAL            LSAME
162      INTEGER            IDAMAX
163      DOUBLE PRECISION   DASUM, DLAMCH
164      EXTERNAL           LSAME, IDAMAX, DASUM, DLAMCH
165*     ..
166*     .. External Subroutines ..
167      EXTERNAL           DGEMV
168*     ..
169*     .. Intrinsic Functions ..
170      INTRINSIC          ABS, DBLE, MAX, MIN
171*     ..
172*     .. Executable Statements ..
173*
174*     Quick return if possible
175*
176      RESID = ZERO
177      IF( N.LE.0 )
178     $   RETURN
179*
180*     Compute B - U * S * V' one column at a time.
181*
182      BNORM = ZERO
183      IF( KD.GE.1 ) THEN
184*
185*        B is bidiagonal.
186*
187         IF( LSAME( UPLO, 'U' ) ) THEN
188*
189*           B is upper bidiagonal.
190*
191            DO 20 J = 1, N
192               DO 10 I = 1, N
193                  WORK( N+I ) = S( I )*VT( I, J )
194   10          CONTINUE
195               CALL DGEMV( 'No transpose', N, N, -ONE, U, LDU,
196     $                     WORK( N+1 ), 1, ZERO, WORK, 1 )
197               WORK( J ) = WORK( J ) + D( J )
198               IF( J.GT.1 ) THEN
199                  WORK( J-1 ) = WORK( J-1 ) + E( J-1 )
200                  BNORM = MAX( BNORM, ABS( D( J ) )+ABS( E( J-1 ) ) )
201               ELSE
202                  BNORM = MAX( BNORM, ABS( D( J ) ) )
203               END IF
204               RESID = MAX( RESID, DASUM( N, WORK, 1 ) )
205   20       CONTINUE
206         ELSE
207*
208*           B is lower bidiagonal.
209*
210            DO 40 J = 1, N
211               DO 30 I = 1, N
212                  WORK( N+I ) = S( I )*VT( I, J )
213   30          CONTINUE
214               CALL DGEMV( 'No transpose', N, N, -ONE, U, LDU,
215     $                     WORK( N+1 ), 1, ZERO, WORK, 1 )
216               WORK( J ) = WORK( J ) + D( J )
217               IF( J.LT.N ) THEN
218                  WORK( J+1 ) = WORK( J+1 ) + E( J )
219                  BNORM = MAX( BNORM, ABS( D( J ) )+ABS( E( J ) ) )
220               ELSE
221                  BNORM = MAX( BNORM, ABS( D( J ) ) )
222               END IF
223               RESID = MAX( RESID, DASUM( N, WORK, 1 ) )
224   40       CONTINUE
225         END IF
226      ELSE
227*
228*        B is diagonal.
229*
230         DO 60 J = 1, N
231            DO 50 I = 1, N
232               WORK( N+I ) = S( I )*VT( I, J )
233   50       CONTINUE
234            CALL DGEMV( 'No transpose', N, N, -ONE, U, LDU, WORK( N+1 ),
235     $                  1, ZERO, WORK, 1 )
236            WORK( J ) = WORK( J ) + D( J )
237            RESID = MAX( RESID, DASUM( N, WORK, 1 ) )
238   60    CONTINUE
239         J = IDAMAX( N, D, 1 )
240         BNORM = ABS( D( J ) )
241      END IF
242*
243*     Compute norm(B - U * S * V') / ( n * norm(B) * EPS )
244*
245      EPS = DLAMCH( 'Precision' )
246*
247      IF( BNORM.LE.ZERO ) THEN
248         IF( RESID.NE.ZERO )
249     $      RESID = ONE / EPS
250      ELSE
251         IF( BNORM.GE.RESID ) THEN
252            RESID = ( RESID / BNORM ) / ( DBLE( N )*EPS )
253         ELSE
254            IF( BNORM.LT.ONE ) THEN
255               RESID = ( MIN( RESID, DBLE( N )*BNORM ) / BNORM ) /
256     $                 ( DBLE( N )*EPS )
257            ELSE
258               RESID = MIN( RESID / BNORM, DBLE( N ) ) /
259     $                 ( DBLE( N )*EPS )
260            END IF
261         END IF
262      END IF
263*
264      RETURN
265*
266*     End of DBDT03
267*
268      END
269