LAPACK 3.11.0
LAPACK: Linear Algebra PACKage
Loading...
Searching...
No Matches

◆ 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 )
Author
Univ. of Tennessee
Univ. of California Berkeley
Univ. of Colorado Denver
NAG Ltd.

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: