 LAPACK  3.10.0 LAPACK: Linear Algebra PACKage

◆ sbdt04()

 subroutine sbdt04 ( character UPLO, integer N, real, dimension( * ) D, real, dimension( * ) E, real, dimension( * ) S, integer NS, real, dimension( ldu, * ) U, integer LDU, real, dimension( ldvt, * ) VT, integer LDVT, real, dimension( * ) WORK, real RESID )

SBDT04

Purpose:
SBDT04 reconstructs a bidiagonal matrix B from its (partial) SVD:
S = U' * B * V
where U and V are orthogonal matrices and S is diagonal.

The test ratio to test the singular value decomposition is
RESID = norm( S - U' * B * V ) / ( n * norm(B) * EPS )
where VT = V' and EPS is the machine precision.
Parameters
 [in] UPLO UPLO is CHARACTER*1 Specifies whether the matrix B is upper or lower bidiagonal. = 'U': Upper bidiagonal = 'L': Lower bidiagonal [in] N N is INTEGER The order of the matrix B. [in] D D is REAL array, dimension (N) The n diagonal elements of the bidiagonal matrix B. [in] E E is REAL array, dimension (N-1) The (n-1) superdiagonal elements of the bidiagonal matrix B if UPLO = 'U', or the (n-1) subdiagonal elements of B if UPLO = 'L'. [in] S S is REAL array, dimension (NS) The singular values from the (partial) SVD of B, sorted in decreasing order. [in] NS NS is INTEGER The number of singular values/vectors from the (partial) SVD of B. [in] U U is REAL array, dimension (LDU,NS) The n by ns orthogonal matrix U in S = U' * B * V. [in] LDU LDU is INTEGER The leading dimension of the array U. LDU >= max(1,N) [in] VT VT is REAL array, dimension (LDVT,N) The n by ns orthogonal matrix V in S = U' * B * V. [in] LDVT LDVT is INTEGER The leading dimension of the array VT. [out] WORK WORK is REAL array, dimension (2*N) [out] RESID RESID is REAL The test ratio: norm(S - U' * B * V) / ( n * norm(B) * EPS )

Definition at line 129 of file sbdt04.f.

131 *
132 * -- LAPACK test routine --
133 * -- LAPACK is a software package provided by Univ. of Tennessee, --
134 * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
135 *
136 * .. Scalar Arguments ..
137  CHARACTER UPLO
138  INTEGER LDU, LDVT, N, NS
139  REAL RESID
140 * ..
141 * .. Array Arguments ..
142  REAL D( * ), E( * ), S( * ), U( LDU, * ),
143  \$ VT( LDVT, * ), WORK( * )
144 * ..
145 *
146 * ======================================================================
147 *
148 * .. Parameters ..
149  REAL ZERO, ONE
150  parameter( zero = 0.0e+0, one = 1.0e+0 )
151 * ..
152 * .. Local Scalars ..
153  INTEGER I, J, K
154  REAL BNORM, EPS
155 * ..
156 * .. External Functions ..
157  LOGICAL LSAME
158  INTEGER ISAMAX
159  REAL SASUM, SLAMCH
160  EXTERNAL lsame, isamax, sasum, slamch
161 * ..
162 * .. External Subroutines ..
163  EXTERNAL sgemm
164 * ..
165 * .. Intrinsic Functions ..
166  INTRINSIC abs, real, max, min
167 * ..
168 * .. Executable Statements ..
169 *
170 * Quick return if possible.
171 *
172  resid = zero
173  IF( n.LE.0 .OR. ns.LE.0 )
174  \$ RETURN
175 *
176  eps = slamch( 'Precision' )
177 *
178 * Compute S - U' * B * V.
179 *
180  bnorm = zero
181 *
182  IF( lsame( uplo, 'U' ) ) THEN
183 *
184 * B is upper bidiagonal.
185 *
186  k = 0
187  DO 20 i = 1, ns
188  DO 10 j = 1, n-1
189  k = k + 1
190  work( k ) = d( j )*vt( i, j ) + e( j )*vt( i, j+1 )
191  10 CONTINUE
192  k = k + 1
193  work( k ) = d( n )*vt( i, n )
194  20 CONTINUE
195  bnorm = abs( d( 1 ) )
196  DO 30 i = 2, n
197  bnorm = max( bnorm, abs( d( i ) )+abs( e( i-1 ) ) )
198  30 CONTINUE
199  ELSE
200 *
201 * B is lower bidiagonal.
202 *
203  k = 0
204  DO 50 i = 1, ns
205  k = k + 1
206  work( k ) = d( 1 )*vt( i, 1 )
207  DO 40 j = 1, n-1
208  k = k + 1
209  work( k ) = e( j )*vt( i, j ) + d( j+1 )*vt( i, j+1 )
210  40 CONTINUE
211  50 CONTINUE
212  bnorm = abs( d( n ) )
213  DO 60 i = 1, n-1
214  bnorm = max( bnorm, abs( d( i ) )+abs( e( i ) ) )
215  60 CONTINUE
216  END IF
217 *
218  CALL sgemm( 'T', 'N', ns, ns, n, -one, u, ldu, work( 1 ),
219  \$ n, zero, work( 1+n*ns ), ns )
220 *
221 * norm(S - U' * B * V)
222 *
223  k = n*ns
224  DO 70 i = 1, ns
225  work( k+i ) = work( k+i ) + s( i )
226  resid = max( resid, sasum( ns, work( k+1 ), 1 ) )
227  k = k + ns
228  70 CONTINUE
229 *
230  IF( bnorm.LE.zero ) THEN
231  IF( resid.NE.zero )
232  \$ resid = one / eps
233  ELSE
234  IF( bnorm.GE.resid ) THEN
235  resid = ( resid / bnorm ) / ( real( n )*eps )
236  ELSE
237  IF( bnorm.LT.one ) THEN
238  resid = ( min( resid, real( n )*bnorm ) / bnorm ) /
239  \$ ( real( n )*eps )
240  ELSE
241  resid = min( resid / bnorm, real( n ) ) /
242  \$ ( real( n )*eps )
243  END IF
244  END IF
245  END IF
246 *
247  RETURN
248 *
249 * End of SBDT04
250 *
integer function isamax(N, SX, INCX)
ISAMAX
Definition: isamax.f:71
logical function lsame(CA, CB)
LSAME
Definition: lsame.f:53
real function sasum(N, SX, INCX)
SASUM
Definition: sasum.f:72
subroutine sgemm(TRANSA, TRANSB, M, N, K, ALPHA, A, LDA, B, LDB, BETA, C, LDC)
SGEMM
Definition: sgemm.f:187
real function slamch(CMACH)
SLAMCH
Definition: slamch.f:68
Here is the call graph for this function:
Here is the caller graph for this function: