1*> \brief \b SBDT01
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 SBDT01( M, N, KD, A, LDA, Q, LDQ, D, E, PT, LDPT, WORK,
12*                          RESID )
13*
14*       .. Scalar Arguments ..
15*       INTEGER            KD, LDA, LDPT, LDQ, M, N
16*       REAL               RESID
17*       ..
18*       .. Array Arguments ..
19*       REAL               A( LDA, * ), D( * ), E( * ), PT( LDPT, * ),
20*      $                   Q( LDQ, * ), WORK( * )
21*       ..
22*
23*
24*> \par Purpose:
25*  =============
26*>
27*> \verbatim
28*>
29*> SBDT01 reconstructs a general matrix A from its bidiagonal form
30*>    A = Q * B * P**T
31*> where Q (m by min(m,n)) and P**T (min(m,n) by n) are orthogonal
32*> matrices and B is bidiagonal.
33*>
34*> The test ratio to test the reduction is
35*>    RESID = norm(A - Q * B * P**T) / ( n * norm(A) * EPS )
36*> where EPS is the machine precision.
37*> \endverbatim
38*
39*  Arguments:
40*  ==========
41*
42*> \param[in] M
43*> \verbatim
44*>          M is INTEGER
45*>          The number of rows of the matrices A and Q.
46*> \endverbatim
47*>
48*> \param[in] N
49*> \verbatim
50*>          N is INTEGER
51*>          The number of columns of the matrices A and P**T.
52*> \endverbatim
53*>
54*> \param[in] KD
55*> \verbatim
56*>          KD is INTEGER
57*>          If KD = 0, B is diagonal and the array E is not referenced.
58*>          If KD = 1, the reduction was performed by xGEBRD; B is upper
59*>          bidiagonal if M >= N, and lower bidiagonal if M < N.
60*>          If KD = -1, the reduction was performed by xGBBRD; B is
61*>          always upper bidiagonal.
62*> \endverbatim
63*>
64*> \param[in] A
65*> \verbatim
66*>          A is REAL array, dimension (LDA,N)
67*>          The m by n matrix A.
68*> \endverbatim
69*>
70*> \param[in] LDA
71*> \verbatim
72*>          LDA is INTEGER
73*>          The leading dimension of the array A.  LDA >= max(1,M).
74*> \endverbatim
75*>
76*> \param[in] Q
77*> \verbatim
78*>          Q is REAL array, dimension (LDQ,N)
79*>          The m by min(m,n) orthogonal matrix Q in the reduction
80*>          A = Q * B * P**T.
81*> \endverbatim
82*>
83*> \param[in] LDQ
84*> \verbatim
85*>          LDQ is INTEGER
86*>          The leading dimension of the array Q.  LDQ >= max(1,M).
87*> \endverbatim
88*>
89*> \param[in] D
90*> \verbatim
91*>          D is REAL array, dimension (min(M,N))
92*>          The diagonal elements of the bidiagonal matrix B.
93*> \endverbatim
94*>
95*> \param[in] E
96*> \verbatim
97*>          E is REAL array, dimension (min(M,N)-1)
98*>          The superdiagonal elements of the bidiagonal matrix B if
99*>          m >= n, or the subdiagonal elements of B if m < n.
100*> \endverbatim
101*>
102*> \param[in] PT
103*> \verbatim
104*>          PT is REAL array, dimension (LDPT,N)
105*>          The min(m,n) by n orthogonal matrix P**T in the reduction
106*>          A = Q * B * P**T.
107*> \endverbatim
108*>
109*> \param[in] LDPT
110*> \verbatim
111*>          LDPT is INTEGER
112*>          The leading dimension of the array PT.
113*>          LDPT >= max(1,min(M,N)).
114*> \endverbatim
115*>
116*> \param[out] WORK
117*> \verbatim
118*>          WORK is REAL array, dimension (M+N)
119*> \endverbatim
120*>
121*> \param[out] RESID
122*> \verbatim
123*>          RESID is REAL
124*>          The test ratio:
125*>          norm(A - Q * B * P**T) / ( n * norm(A) * EPS )
126*> \endverbatim
127*
128*  Authors:
129*  ========
130*
131*> \author Univ. of Tennessee
132*> \author Univ. of California Berkeley
133*> \author Univ. of Colorado Denver
134*> \author NAG Ltd.
135*
136*> \ingroup single_eig
137*
138*  =====================================================================
139      SUBROUTINE SBDT01( M, N, KD, A, LDA, Q, LDQ, D, E, PT, LDPT, WORK,
140     $                   RESID )
141*
142*  -- LAPACK test routine --
143*  -- LAPACK is a software package provided by Univ. of Tennessee,    --
144*  -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
145*
146*     .. Scalar Arguments ..
147      INTEGER            KD, LDA, LDPT, LDQ, M, N
148      REAL               RESID
149*     ..
150*     .. Array Arguments ..
151      REAL               A( LDA, * ), D( * ), E( * ), PT( LDPT, * ),
152     $                   Q( LDQ, * ), WORK( * )
153*     ..
154*
155*  =====================================================================
156*
157*     .. Parameters ..
158      REAL               ZERO, ONE
159      PARAMETER          ( ZERO = 0.0E+0, ONE = 1.0E+0 )
160*     ..
161*     .. Local Scalars ..
162      INTEGER            I, J
163      REAL               ANORM, EPS
164*     ..
165*     .. External Functions ..
166      REAL               SASUM, SLAMCH, SLANGE
167      EXTERNAL           SASUM, SLAMCH, SLANGE
168*     ..
169*     .. External Subroutines ..
170      EXTERNAL           SCOPY, SGEMV
171*     ..
172*     .. Intrinsic Functions ..
173      INTRINSIC          MAX, MIN, REAL
174*     ..
175*     .. Executable Statements ..
176*
177*     Quick return if possible
178*
179      IF( M.LE.0 .OR. N.LE.0 ) THEN
180         RESID = ZERO
181         RETURN
182      END IF
183*
184*     Compute A - Q * B * P**T one column at a time.
185*
186      RESID = ZERO
187      IF( KD.NE.0 ) THEN
188*
189*        B is bidiagonal.
190*
191         IF( KD.NE.0 .AND. M.GE.N ) THEN
192*
193*           B is upper bidiagonal and M >= N.
194*
195            DO 20 J = 1, N
196               CALL SCOPY( M, A( 1, J ), 1, WORK, 1 )
197               DO 10 I = 1, N - 1
198                  WORK( M+I ) = D( I )*PT( I, J ) + E( I )*PT( I+1, J )
199   10          CONTINUE
200               WORK( M+N ) = D( N )*PT( N, J )
201               CALL SGEMV( 'No transpose', M, N, -ONE, Q, LDQ,
202     $                     WORK( M+1 ), 1, ONE, WORK, 1 )
203               RESID = MAX( RESID, SASUM( M, WORK, 1 ) )
204   20       CONTINUE
205         ELSE IF( KD.LT.0 ) THEN
206*
207*           B is upper bidiagonal and M < N.
208*
209            DO 40 J = 1, N
210               CALL SCOPY( M, A( 1, J ), 1, WORK, 1 )
211               DO 30 I = 1, M - 1
212                  WORK( M+I ) = D( I )*PT( I, J ) + E( I )*PT( I+1, J )
213   30          CONTINUE
214               WORK( M+M ) = D( M )*PT( M, J )
215               CALL SGEMV( 'No transpose', M, M, -ONE, Q, LDQ,
216     $                     WORK( M+1 ), 1, ONE, WORK, 1 )
217               RESID = MAX( RESID, SASUM( M, WORK, 1 ) )
218   40       CONTINUE
219         ELSE
220*
221*           B is lower bidiagonal.
222*
223            DO 60 J = 1, N
224               CALL SCOPY( M, A( 1, J ), 1, WORK, 1 )
225               WORK( M+1 ) = D( 1 )*PT( 1, J )
226               DO 50 I = 2, M
227                  WORK( M+I ) = E( I-1 )*PT( I-1, J ) +
228     $                          D( I )*PT( I, J )
229   50          CONTINUE
230               CALL SGEMV( 'No transpose', M, M, -ONE, Q, LDQ,
231     $                     WORK( M+1 ), 1, ONE, WORK, 1 )
232               RESID = MAX( RESID, SASUM( M, WORK, 1 ) )
233   60       CONTINUE
234         END IF
235      ELSE
236*
237*        B is diagonal.
238*
239         IF( M.GE.N ) THEN
240            DO 80 J = 1, N
241               CALL SCOPY( M, A( 1, J ), 1, WORK, 1 )
242               DO 70 I = 1, N
243                  WORK( M+I ) = D( I )*PT( I, J )
244   70          CONTINUE
245               CALL SGEMV( 'No transpose', M, N, -ONE, Q, LDQ,
246     $                     WORK( M+1 ), 1, ONE, WORK, 1 )
247               RESID = MAX( RESID, SASUM( M, WORK, 1 ) )
248   80       CONTINUE
249         ELSE
250            DO 100 J = 1, N
251               CALL SCOPY( M, A( 1, J ), 1, WORK, 1 )
252               DO 90 I = 1, M
253                  WORK( M+I ) = D( I )*PT( I, J )
254   90          CONTINUE
255               CALL SGEMV( 'No transpose', M, M, -ONE, Q, LDQ,
256     $                     WORK( M+1 ), 1, ONE, WORK, 1 )
257               RESID = MAX( RESID, SASUM( M, WORK, 1 ) )
258  100       CONTINUE
259         END IF
260      END IF
261*
262*     Compute norm(A - Q * B * P**T) / ( n * norm(A) * EPS )
263*
264      ANORM = SLANGE( '1', M, N, A, LDA, WORK )
265      EPS = SLAMCH( 'Precision' )
266*
267      IF( ANORM.LE.ZERO ) THEN
268         IF( RESID.NE.ZERO )
269     $      RESID = ONE / EPS
270      ELSE
271         IF( ANORM.GE.RESID ) THEN
272            RESID = ( RESID / ANORM ) / ( REAL( N )*EPS )
273         ELSE
274            IF( ANORM.LT.ONE ) THEN
275               RESID = ( MIN( RESID, REAL( N )*ANORM ) / ANORM ) /
276     $                 ( REAL( N )*EPS )
277            ELSE
278               RESID = MIN( RESID / ANORM, REAL( N ) ) /
279     $                 ( REAL( N )*EPS )
280            END IF
281         END IF
282      END IF
283*
284      RETURN
285*
286*     End of SBDT01
287*
288      END
289