LAPACK  3.6.1
LAPACK: Linear Algebra PACKage
dbdt05.f
Go to the documentation of this file.
1 * =========== DOCUMENTATION ===========
2 *
3 * Online html documentation available at
4 * http://www.netlib.org/lapack/explore-html/
5 *
6 * Definition:
7 * ===========
8 *
9 * SUBROUTINE DBDT05( M, N, A, LDA, S, NS, U, LDU,
10 * VT, LDVT, WORK, RESID )
11 *
12 * .. Scalar Arguments ..
13 * INTEGER LDA, LDU, LDVT, N, NS
14 * DOUBLE PRECISION RESID
15 * ..
16 * .. Array Arguments ..
17 * DOUBLE PRECISION D( * ), E( * ), S( * ), U( LDU, * ),
18 * $ VT( LDVT, * ), WORK( * )
19 * ..
20 *
21 *
22 *> \par Purpose:
23 * =============
24 *>
25 *> \verbatim
26 *>
27 *> DBDT05 reconstructs a bidiagonal matrix B from its (partial) SVD:
28 *> S = U' * B * V
29 *> where U and V are orthogonal matrices and S is diagonal.
30 *>
31 *> The test ratio to test the singular value decomposition is
32 *> RESID = norm( S - U' * B * V ) / ( n * norm(B) * EPS )
33 *> where VT = V' and EPS is the machine precision.
34 *> \endverbatim
35 *
36 * Arguments:
37 * ==========
38 *
39 *> \param[in] M
40 *> \verbatim
41 *> M is INTEGER
42 *> The number of rows of the matrices A and U.
43 *> \endverbatim
44 *>
45 *> \param[in] N
46 *> \verbatim
47 *> N is INTEGER
48 *> The number of columns of the matrices A and VT.
49 *> \endverbatim
50 *>
51 *> \param[in] A
52 *> \verbatim
53 *> A is DOUBLE PRECISION array, dimension (LDA,N)
54 *> The m by n matrix A.
55 *>
56 *> \param[in] LDA
57 *> \verbatim
58 *> LDA is INTEGER
59 *> The leading dimension of the array A. LDA >= max(1,M).
60 *> \endverbatim
61 *>
62 *> \param[in] S
63 *> \verbatim
64 *> S is DOUBLE PRECISION array, dimension (NS)
65 *> The singular values from the (partial) SVD of B, sorted in
66 *> decreasing order.
67 *> \endverbatim
68 *>
69 *> \param[in] NS
70 *> \verbatim
71 *> NS is INTEGER
72 *> The number of singular values/vectors from the (partial)
73 *> SVD of B.
74 *> \endverbatim
75 *>
76 *> \param[in] U
77 *> \verbatim
78 *> U is DOUBLE PRECISION array, dimension (LDU,NS)
79 *> The n by ns orthogonal matrix U in S = U' * B * V.
80 *> \endverbatim
81 *>
82 *> \param[in] LDU
83 *> \verbatim
84 *> LDU is INTEGER
85 *> The leading dimension of the array U. LDU >= max(1,N)
86 *> \endverbatim
87 *>
88 *> \param[in] VT
89 *> \verbatim
90 *> VT is DOUBLE PRECISION array, dimension (LDVT,N)
91 *> The n by ns orthogonal matrix V in S = U' * B * V.
92 *> \endverbatim
93 *>
94 *> \param[in] LDVT
95 *> \verbatim
96 *> LDVT is INTEGER
97 *> The leading dimension of the array VT.
98 *> \endverbatim
99 *>
100 *> \param[out] WORK
101 *> \verbatim
102 *> WORK is DOUBLE PRECISION array, dimension (M,N)
103 *> \endverbatim
104 *>
105 *> \param[out] RESID
106 *> \verbatim
107 *> RESID is DOUBLE PRECISION
108 *> The test ratio: norm(S - U' * A * V) / ( n * norm(A) * EPS )
109 *> \endverbatim
110 *
111 * Authors:
112 * ========
113 *
114 *> \author Univ. of Tennessee
115 *> \author Univ. of California Berkeley
116 *> \author Univ. of Colorado Denver
117 *> \author NAG Ltd.
118 *
119 *> \date November 2011
120 *
121 *> \ingroup double_eig
122 *
123 * =====================================================================
124  SUBROUTINE dbdt05( M, N, A, LDA, S, NS, U, LDU,
125  $ vt, ldvt, work, resid )
126 *
127 * -- LAPACK test routine (version 3.4.0) --
128 * -- LAPACK is a software package provided by Univ. of Tennessee, --
129 * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
130 * November 2011
131 *
132 * .. Scalar Arguments ..
133  CHARACTER UPLO
134  INTEGER LDA, LDU, LDVT, M, N, NS
135  DOUBLE PRECISION RESID
136 * ..
137 * .. Array Arguments ..
138  DOUBLE PRECISION A( lda, * ), S( * ), U( ldu, * ),
139  $ vt( ldvt, * ), work( * )
140 * ..
141 *
142 * ======================================================================
143 *
144 * .. Parameters ..
145  DOUBLE PRECISION ZERO, ONE
146  parameter ( zero = 0.0d+0, one = 1.0d+0 )
147 * ..
148 * .. Local Scalars ..
149  INTEGER I, J
150  DOUBLE PRECISION ANORM, EPS
151 * ..
152 * .. External Functions ..
153  LOGICAL LSAME
154  INTEGER IDAMAX
155  DOUBLE PRECISION DASUM, DLAMCH, DLANGE
156  EXTERNAL lsame, idamax, dasum, dlamch, dlange
157 * ..
158 * .. External Subroutines ..
159  EXTERNAL dgemm
160 * ..
161 * .. Intrinsic Functions ..
162  INTRINSIC abs, dble, max, min
163 * ..
164 * .. Executable Statements ..
165 *
166 * Quick return if possible.
167 *
168  resid = zero
169  IF( min( m, n ).LE.0 .OR. ns.LE.0 )
170  $ RETURN
171 *
172  eps = dlamch( 'Precision' )
173  anorm = dlange( 'M', m, n, a, lda, work )
174 *
175 * Compute U' * A * V.
176 *
177  CALL dgemm( 'N', 'T', m, ns, n, one, a, lda, vt,
178  $ ldvt, zero, work( 1+ns*ns ), m )
179  CALL dgemm( 'T', 'N', ns, ns, m, -one, u, ldu, work( 1+ns*ns ),
180  $ m, zero, work, ns )
181 *
182 * norm(S - U' * B * V)
183 *
184  j = 0
185  DO 10 i = 1, ns
186  work( j+i ) = work( j+i ) + s( i )
187  resid = max( resid, dasum( ns, work( j+1 ), 1 ) )
188  j = j + ns
189  10 CONTINUE
190 *
191  IF( anorm.LE.zero ) THEN
192  IF( resid.NE.zero )
193  $ resid = one / eps
194  ELSE
195  IF( anorm.GE.resid ) THEN
196  resid = ( resid / anorm ) / ( dble( n )*eps )
197  ELSE
198  IF( anorm.LT.one ) THEN
199  resid = ( min( resid, dble( n )*anorm ) / anorm ) /
200  $ ( dble( n )*eps )
201  ELSE
202  resid = min( resid / anorm, dble( n ) ) /
203  $ ( dble( n )*eps )
204  END IF
205  END IF
206  END IF
207 *
208  RETURN
209 *
210 * End of DBDT05
211 *
212  END
subroutine dgemm(TRANSA, TRANSB, M, N, K, ALPHA, A, LDA, B, LDB, BETA, C, LDC)
DGEMM
Definition: dgemm.f:189
subroutine dbdt05(M, N, A, LDA, S, NS, U, LDU, VT, LDVT, WORK, RESID)
Definition: dbdt05.f:126