LAPACK  3.8.0 LAPACK: Linear Algebra PACKage
dbdt04.f
Go to the documentation of this file.
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 *> \date December 2016
127 *
128 *> \ingroup double_eig
129 *
130 * =====================================================================
131  SUBROUTINE dbdt04( UPLO, N, D, E, S, NS, U, LDU, VT, LDVT, WORK,
132  \$ RESID )
133 *
134 * -- LAPACK test routine (version 3.7.0) --
135 * -- LAPACK is a software package provided by Univ. of Tennessee, --
136 * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
137 * December 2016
138 *
139 * .. Scalar Arguments ..
140  CHARACTER UPLO
141  INTEGER LDU, LDVT, N, NS
142  DOUBLE PRECISION RESID
143 * ..
144 * .. Array Arguments ..
145  DOUBLE PRECISION D( * ), E( * ), S( * ), U( ldu, * ),
146  \$ vt( ldvt, * ), work( * )
147 * ..
148 *
149 * ======================================================================
150 *
151 * .. Parameters ..
152  DOUBLE PRECISION ZERO, ONE
153  parameter( zero = 0.0d+0, one = 1.0d+0 )
154 * ..
155 * .. Local Scalars ..
156  INTEGER I, J, K
157  DOUBLE PRECISION BNORM, EPS
158 * ..
159 * .. External Functions ..
160  LOGICAL LSAME
161  INTEGER IDAMAX
162  DOUBLE PRECISION DASUM, DLAMCH
163  EXTERNAL lsame, idamax, dasum, dlamch
164 * ..
165 * .. External Subroutines ..
166  EXTERNAL dgemm
167 * ..
168 * .. Intrinsic Functions ..
169  INTRINSIC abs, dble, max, min
170 * ..
171 * .. Executable Statements ..
172 *
173 * Quick return if possible.
174 *
175  resid = zero
176  IF( n.LE.0 .OR. ns.LE.0 )
177  \$ RETURN
178 *
179  eps = dlamch( 'Precision' )
180 *
181 * Compute S - U' * B * V.
182 *
183  bnorm = zero
184 *
185  IF( lsame( uplo, 'U' ) ) THEN
186 *
187 * B is upper bidiagonal.
188 *
189  k = 0
190  DO 20 i = 1, ns
191  DO 10 j = 1, n-1
192  k = k + 1
193  work( k ) = d( j )*vt( i, j ) + e( j )*vt( i, j+1 )
194  10 CONTINUE
195  k = k + 1
196  work( k ) = d( n )*vt( i, n )
197  20 CONTINUE
198  bnorm = abs( d( 1 ) )
199  DO 30 i = 2, n
200  bnorm = max( bnorm, abs( d( i ) )+abs( e( i-1 ) ) )
201  30 CONTINUE
202  ELSE
203 *
204 * B is lower bidiagonal.
205 *
206  k = 0
207  DO 50 i = 1, ns
208  k = k + 1
209  work( k ) = d( 1 )*vt( i, 1 )
210  DO 40 j = 1, n-1
211  k = k + 1
212  work( k ) = e( j )*vt( i, j ) + d( j+1 )*vt( i, j+1 )
213  40 CONTINUE
214  50 CONTINUE
215  bnorm = abs( d( n ) )
216  DO 60 i = 1, n-1
217  bnorm = max( bnorm, abs( d( i ) )+abs( e( i ) ) )
218  60 CONTINUE
219  END IF
220 *
221  CALL dgemm( 'T', 'N', ns, ns, n, -one, u, ldu, work( 1 ),
222  \$ n, zero, work( 1+n*ns ), ns )
223 *
224 * norm(S - U' * B * V)
225 *
226  k = n*ns
227  DO 70 i = 1, ns
228  work( k+i ) = work( k+i ) + s( i )
229  resid = max( resid, dasum( ns, work( k+1 ), 1 ) )
230  k = k + ns
231  70 CONTINUE
232 *
233  IF( bnorm.LE.zero ) THEN
234  IF( resid.NE.zero )
235  \$ resid = one / eps
236  ELSE
237  IF( bnorm.GE.resid ) THEN
238  resid = ( resid / bnorm ) / ( dble( n )*eps )
239  ELSE
240  IF( bnorm.LT.one ) THEN
241  resid = ( min( resid, dble( n )*bnorm ) / bnorm ) /
242  \$ ( dble( n )*eps )
243  ELSE
244  resid = min( resid / bnorm, dble( n ) ) /
245  \$ ( dble( n )*eps )
246  END IF
247  END IF
248  END IF
249 *
250  RETURN
251 *
252 * End of DBDT04
253 *
254  END
subroutine dgemm(TRANSA, TRANSB, M, N, K, ALPHA, A, LDA, B, LDB, BETA, C, LDC)
DGEMM
Definition: dgemm.f:189
subroutine dbdt04(UPLO, N, D, E, S, NS, U, LDU, VT, LDVT, WORK, RESID)
DBDT04
Definition: dbdt04.f:133