LAPACK  3.8.0
LAPACK: Linear Algebra PACKage
zbdt05.f
Go to the documentation of this file.
1 *> \brief \b ZBDT05
2 * =========== DOCUMENTATION ===========
3 *
4 * Online html documentation available at
5 * http://www.netlib.org/lapack/explore-html/
6 *
7 * Definition:
8 * ===========
9 *
10 * SUBROUTINE ZBDT05( M, N, A, LDA, S, NS, U, LDU,
11 * VT, LDVT, WORK, RESID )
12 *
13 * .. Scalar Arguments ..
14 * INTEGER LDA, LDU, LDVT, N, NS
15 * DOUBLE PRECISION RESID
16 * ..
17 * .. Array Arguments ..
18 * DOUBLE PRECISION S( * )
19 * COMPLEX*16 A( LDA, * ), U( * ), VT( LDVT, * ), WORK( * )
20 * ..
21 *
22 *> \par Purpose:
23 * =============
24 *>
25 *> \verbatim
26 *>
27 *> ZBDT05 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 COMPLEX*16 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 COMPLEX*16 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 COMPLEX*16 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 COMPLEX*16 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 December 2016
120 *
121 *> \ingroup double_eig
122 *
123 * =====================================================================
124  SUBROUTINE zbdt05( M, N, A, LDA, S, NS, U, LDU,
125  $ VT, LDVT, WORK, RESID )
126 *
127 * -- LAPACK test routine (version 3.7.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 * December 2016
131 *
132 * .. Scalar Arguments ..
133  INTEGER LDA, LDU, LDVT, M, N, NS
134  DOUBLE PRECISION RESID
135 * ..
136 * .. Array Arguments ..
137  DOUBLE PRECISION S( * )
138  COMPLEX*16 A( lda, * ), U( * ), 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  COMPLEX*16 CZERO, CONE
147  parameter( czero = ( 0.0d+0, 0.0d+0 ),
148  $ cone = ( 1.0d+0, 0.0d+0 ) )
149 * ..
150 * .. Local Scalars ..
151  INTEGER I, J
152  DOUBLE PRECISION ANORM, EPS
153 * ..
154 * .. Local Arrays ..
155  DOUBLE PRECISION DUM( 1 )
156 * ..
157 * .. External Functions ..
158  LOGICAL LSAME
159  INTEGER IDAMAX
160  DOUBLE PRECISION DASUM, DLAMCH, ZLANGE
161  EXTERNAL lsame, idamax, dasum, dlamch, zlange
162  DOUBLE PRECISION DZASUM
163 * ..
164 * .. External Subroutines ..
165  EXTERNAL zgemm
166 * ..
167 * .. Intrinsic Functions ..
168  INTRINSIC abs, dble, max, min
169 * ..
170 * .. Executable Statements ..
171 *
172 * Quick return if possible.
173 *
174  resid = zero
175  IF( min( m, n ).LE.0 .OR. ns.LE.0 )
176  $ RETURN
177 *
178  eps = dlamch( 'Precision' )
179  anorm = zlange( 'M', m, n, a, lda, dum )
180 *
181 * Compute U' * A * V.
182 *
183  CALL zgemm( 'N', 'C', m, ns, n, cone, a, lda, vt,
184  $ ldvt, czero, work( 1+ns*ns ), m )
185  CALL zgemm( 'C', 'N', ns, ns, m, -cone, u, ldu, work( 1+ns*ns ),
186  $ m, czero, work, ns )
187 *
188 * norm(S - U' * B * V)
189 *
190  j = 0
191  DO 10 i = 1, ns
192  work( j+i ) = work( j+i ) + dcmplx( s( i ), zero )
193  resid = max( resid, dzasum( ns, work( j+1 ), 1 ) )
194  j = j + ns
195  10 CONTINUE
196 *
197  IF( anorm.LE.zero ) THEN
198  IF( resid.NE.zero )
199  $ resid = one / eps
200  ELSE
201  IF( anorm.GE.resid ) THEN
202  resid = ( resid / anorm ) / ( dble( n )*eps )
203  ELSE
204  IF( anorm.LT.one ) THEN
205  resid = ( min( resid, dble( n )*anorm ) / anorm ) /
206  $ ( dble( n )*eps )
207  ELSE
208  resid = min( resid / anorm, dble( n ) ) /
209  $ ( dble( n )*eps )
210  END IF
211  END IF
212  END IF
213 *
214  RETURN
215 *
216 * End of ZBDT05
217 *
218  END
subroutine zbdt05(M, N, A, LDA, S, NS, U, LDU, VT, LDVT, WORK, RESID)
ZBDT05
Definition: zbdt05.f:126
subroutine zgemm(TRANSA, TRANSB, M, N, K, ALPHA, A, LDA, B, LDB, BETA, C, LDC)
ZGEMM
Definition: zgemm.f:189