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*> \date December 2016
131*
132*> \ingroup double_eig
133*
134*  =====================================================================
135      SUBROUTINE DBDT03( UPLO, N, KD, D, E, U, LDU, S, VT, LDVT, WORK,
136     $                   RESID )
137*
138*  -- LAPACK test routine (version 3.7.0) --
139*  -- LAPACK is a software package provided by Univ. of Tennessee,    --
140*  -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
141*     December 2016
142*
143*     .. Scalar Arguments ..
144      CHARACTER          UPLO
145      INTEGER            KD, LDU, LDVT, N
146      DOUBLE PRECISION   RESID
147*     ..
148*     .. Array Arguments ..
149      DOUBLE PRECISION   D( * ), E( * ), S( * ), U( LDU, * ),
150     $                   VT( LDVT, * ), WORK( * )
151*     ..
152*
153* ======================================================================
154*
155*     .. Parameters ..
156      DOUBLE PRECISION   ZERO, ONE
157      PARAMETER          ( ZERO = 0.0D+0, ONE = 1.0D+0 )
158*     ..
159*     .. Local Scalars ..
160      INTEGER            I, J
161      DOUBLE PRECISION   BNORM, EPS
162*     ..
163*     .. External Functions ..
164      LOGICAL            LSAME
165      INTEGER            IDAMAX
166      DOUBLE PRECISION   DASUM, DLAMCH
167      EXTERNAL           LSAME, IDAMAX, DASUM, DLAMCH
168*     ..
169*     .. External Subroutines ..
170      EXTERNAL           DGEMV
171*     ..
172*     .. Intrinsic Functions ..
173      INTRINSIC          ABS, DBLE, MAX, MIN
174*     ..
175*     .. Executable Statements ..
176*
177*     Quick return if possible
178*
179      RESID = ZERO
180      IF( N.LE.0 )
181     $   RETURN
182*
183*     Compute B - U * S * V' one column at a time.
184*
185      BNORM = ZERO
186      IF( KD.GE.1 ) THEN
187*
188*        B is bidiagonal.
189*
190         IF( LSAME( UPLO, 'U' ) ) THEN
191*
192*           B is upper bidiagonal.
193*
194            DO 20 J = 1, N
195               DO 10 I = 1, N
196                  WORK( N+I ) = S( I )*VT( I, J )
197   10          CONTINUE
198               CALL DGEMV( 'No transpose', N, N, -ONE, U, LDU,
199     $                     WORK( N+1 ), 1, ZERO, WORK, 1 )
200               WORK( J ) = WORK( J ) + D( J )
201               IF( J.GT.1 ) THEN
202                  WORK( J-1 ) = WORK( J-1 ) + E( J-1 )
203                  BNORM = MAX( BNORM, ABS( D( J ) )+ABS( E( J-1 ) ) )
204               ELSE
205                  BNORM = MAX( BNORM, ABS( D( J ) ) )
206               END IF
207               RESID = MAX( RESID, DASUM( N, WORK, 1 ) )
208   20       CONTINUE
209         ELSE
210*
211*           B is lower bidiagonal.
212*
213            DO 40 J = 1, N
214               DO 30 I = 1, N
215                  WORK( N+I ) = S( I )*VT( I, J )
216   30          CONTINUE
217               CALL DGEMV( 'No transpose', N, N, -ONE, U, LDU,
218     $                     WORK( N+1 ), 1, ZERO, WORK, 1 )
219               WORK( J ) = WORK( J ) + D( J )
220               IF( J.LT.N ) THEN
221                  WORK( J+1 ) = WORK( J+1 ) + E( J )
222                  BNORM = MAX( BNORM, ABS( D( J ) )+ABS( E( J ) ) )
223               ELSE
224                  BNORM = MAX( BNORM, ABS( D( J ) ) )
225               END IF
226               RESID = MAX( RESID, DASUM( N, WORK, 1 ) )
227   40       CONTINUE
228         END IF
229      ELSE
230*
231*        B is diagonal.
232*
233         DO 60 J = 1, N
234            DO 50 I = 1, N
235               WORK( N+I ) = S( I )*VT( I, J )
236   50       CONTINUE
237            CALL DGEMV( 'No transpose', N, N, -ONE, U, LDU, WORK( N+1 ),
238     $                  1, ZERO, WORK, 1 )
239            WORK( J ) = WORK( J ) + D( J )
240            RESID = MAX( RESID, DASUM( N, WORK, 1 ) )
241   60    CONTINUE
242         J = IDAMAX( N, D, 1 )
243         BNORM = ABS( D( J ) )
244      END IF
245*
246*     Compute norm(B - U * S * V') / ( n * norm(B) * EPS )
247*
248      EPS = DLAMCH( 'Precision' )
249*
250      IF( BNORM.LE.ZERO ) THEN
251         IF( RESID.NE.ZERO )
252     $      RESID = ONE / EPS
253      ELSE
254         IF( BNORM.GE.RESID ) THEN
255            RESID = ( RESID / BNORM ) / ( DBLE( N )*EPS )
256         ELSE
257            IF( BNORM.LT.ONE ) THEN
258               RESID = ( MIN( RESID, DBLE( N )*BNORM ) / BNORM ) /
259     $                 ( DBLE( N )*EPS )
260            ELSE
261               RESID = MIN( RESID / BNORM, DBLE( N ) ) /
262     $                 ( DBLE( N )*EPS )
263            END IF
264         END IF
265      END IF
266*
267      RETURN
268*
269*     End of DBDT03
270*
271      END
272