1*> \brief \b DBDT05
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 DBDT05( M, N, A, LDA, S, NS, U, LDU,
12*                          VT, LDVT, WORK, RESID )
13*
14*       .. Scalar Arguments ..
15*       INTEGER            LDA, 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*> DBDT05 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] M
42*> \verbatim
43*>          M is INTEGER
44*>          The number of rows of the matrices A and U.
45*> \endverbatim
46*>
47*> \param[in] N
48*> \verbatim
49*>          N is INTEGER
50*>          The number of columns of the matrices A and VT.
51*> \endverbatim
52*>
53*> \param[in] A
54*> \verbatim
55*>          A is DOUBLE PRECISION array, dimension (LDA,N)
56*>          The m by n matrix A.
57*> \endverbatim
58*>
59*> \param[in] LDA
60*> \verbatim
61*>          LDA is INTEGER
62*>          The leading dimension of the array A.  LDA >= max(1,M).
63*> \endverbatim
64*>
65*> \param[in] S
66*> \verbatim
67*>          S is DOUBLE PRECISION array, dimension (NS)
68*>          The singular values from the (partial) SVD of B, sorted in
69*>          decreasing order.
70*> \endverbatim
71*>
72*> \param[in] NS
73*> \verbatim
74*>          NS is INTEGER
75*>          The number of singular values/vectors from the (partial)
76*>          SVD of B.
77*> \endverbatim
78*>
79*> \param[in] U
80*> \verbatim
81*>          U is DOUBLE PRECISION array, dimension (LDU,NS)
82*>          The n by ns orthogonal matrix U in S = U' * B * V.
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] VT
92*> \verbatim
93*>          VT is DOUBLE PRECISION array, dimension (LDVT,N)
94*>          The n by ns orthogonal matrix V in S = U' * B * V.
95*> \endverbatim
96*>
97*> \param[in] LDVT
98*> \verbatim
99*>          LDVT is INTEGER
100*>          The leading dimension of the array VT.
101*> \endverbatim
102*>
103*> \param[out] WORK
104*> \verbatim
105*>          WORK is DOUBLE PRECISION array, dimension (M,N)
106*> \endverbatim
107*>
108*> \param[out] RESID
109*> \verbatim
110*>          RESID is DOUBLE PRECISION
111*>          The test ratio:  norm(S - U' * A * V) / ( n * norm(A) * EPS )
112*> \endverbatim
113*
114*  Authors:
115*  ========
116*
117*> \author Univ. of Tennessee
118*> \author Univ. of California Berkeley
119*> \author Univ. of Colorado Denver
120*> \author NAG Ltd.
121*
122*> \ingroup double_eig
123*
124*  =====================================================================
125      SUBROUTINE DBDT05( M, N, A, LDA, S, NS, U, LDU,
126     $                    VT, LDVT, WORK, RESID )
127*
128*  -- LAPACK test routine --
129*  -- LAPACK is a software package provided by Univ. of Tennessee,    --
130*  -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
131*
132*     .. Scalar Arguments ..
133      INTEGER            LDA, LDU, LDVT, M, N, NS
134      DOUBLE PRECISION   RESID
135*     ..
136*     .. Array Arguments ..
137      DOUBLE PRECISION   A( LDA, * ), S( * ), U( LDU, * ),
138     $                   VT( LDVT, * ), WORK( * )
139*     ..
140*
141* ======================================================================
142*
143*     .. Parameters ..
144      DOUBLE PRECISION   ZERO, ONE
145      PARAMETER          ( ZERO = 0.0D+0, ONE = 1.0D+0 )
146*     ..
147*     .. Local Scalars ..
148      INTEGER            I, J
149      DOUBLE PRECISION   ANORM, EPS
150*     ..
151*     .. External Functions ..
152      LOGICAL            LSAME
153      INTEGER            IDAMAX
154      DOUBLE PRECISION   DASUM, DLAMCH, DLANGE
155      EXTERNAL           LSAME, IDAMAX, DASUM, DLAMCH, DLANGE
156*     ..
157*     .. External Subroutines ..
158      EXTERNAL           DGEMM
159*     ..
160*     .. Intrinsic Functions ..
161      INTRINSIC          ABS, DBLE, MAX, MIN
162*     ..
163*     .. Executable Statements ..
164*
165*     Quick return if possible.
166*
167      RESID = ZERO
168      IF( MIN( M, N ).LE.0 .OR. NS.LE.0 )
169     $   RETURN
170*
171      EPS = DLAMCH( 'Precision' )
172      ANORM = DLANGE( 'M', M, N, A, LDA, WORK )
173*
174*     Compute U' * A * V.
175*
176      CALL DGEMM( 'N', 'T', M, NS, N, ONE, A, LDA, VT,
177     $            LDVT, ZERO, WORK( 1+NS*NS ), M )
178      CALL DGEMM( 'T', 'N', NS, NS, M, -ONE, U, LDU, WORK( 1+NS*NS ),
179     $            M, ZERO, WORK, NS )
180*
181*     norm(S - U' * B * V)
182*
183      J = 0
184      DO 10 I = 1, NS
185         WORK( J+I ) =  WORK( J+I ) + S( I )
186         RESID = MAX( RESID, DASUM( NS, WORK( J+1 ), 1 ) )
187         J = J + NS
188   10 CONTINUE
189*
190      IF( ANORM.LE.ZERO ) THEN
191         IF( RESID.NE.ZERO )
192     $      RESID = ONE / EPS
193      ELSE
194         IF( ANORM.GE.RESID ) THEN
195            RESID = ( RESID / ANORM ) / ( DBLE( N )*EPS )
196         ELSE
197            IF( ANORM.LT.ONE ) THEN
198               RESID = ( MIN( RESID, DBLE( N )*ANORM ) / ANORM ) /
199     $                 ( DBLE( N )*EPS )
200            ELSE
201               RESID = MIN( RESID / ANORM, DBLE( N ) ) /
202     $                 ( DBLE( N )*EPS )
203            END IF
204         END IF
205      END IF
206*
207      RETURN
208*
209*     End of DBDT05
210*
211      END
212