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