LAPACK  3.10.0
LAPACK: Linear Algebra PACKage
dbdt05.f
Go to the documentation of this file.
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
subroutine dgemm(TRANSA, TRANSB, M, N, K, ALPHA, A, LDA, B, LDB, BETA, C, LDC)
DGEMM
Definition: dgemm.f:187
subroutine dbdt05(M, N, A, LDA, S, NS, U, LDU, VT, LDVT, WORK, RESID)
DBDT05
Definition: dbdt05.f:127