209 SUBROUTINE sgsvts( M, P, N, A, AF, LDA, B, BF, LDB, U, LDU, V,
210 $ ldv, q, ldq, alpha, beta, r, ldr, iwork, work,
211 $ lwork, rwork, result )
219 INTEGER LDA, LDB, LDQ, LDR, LDU, LDV, LWORK, M, N, P
223 REAL A( lda, * ), AF( lda, * ), ALPHA( * ),
224 $ b( ldb, * ), beta( * ), bf( ldb, * ),
225 $ q( ldq, * ), r( ldr, * ), result( 6 ),
226 $ rwork( * ), u( ldu, * ), v( ldv, * ),
234 parameter( zero = 0.0e+0, one = 1.0e+0 )
237 INTEGER I, INFO, J, K, L
238 REAL ANORM, BNORM, RESID, TEMP, ULP, ULPINV, UNFL
241 REAL SLAMCH, SLANGE, SLANSY
242 EXTERNAL slamch, slange, slansy
248 INTRINSIC max, min, real
252 ulp = slamch(
'Precision' )
254 unfl = slamch(
'Safe minimum' )
258 CALL
slacpy(
'Full', m, n, a, lda, af, lda )
259 CALL
slacpy(
'Full', p, n, b, ldb, bf, ldb )
261 anorm = max( slange(
'1', m, n, a, lda, rwork ), unfl )
262 bnorm = max( slange(
'1', p, n, b, ldb, rwork ), unfl )
266 CALL
sggsvd(
'U',
'V',
'Q', m, n, p, k, l, af, lda, bf, ldb,
267 $ alpha, beta, u, ldu, v, ldv, q, ldq, work, iwork,
272 DO 20 i = 1, min( k+l, m )
274 r( i, j ) = af( i, n-k-l+j )
278 IF( m-k-l.LT.0 )
THEN
279 DO 40 i = m + 1, k + l
281 r( i, j ) = bf( i-k, n-k-l+j )
288 CALL
sgemm(
'No transpose',
'No transpose', m, n, n, one, a, lda,
289 $ q, ldq, zero, work, lda )
291 CALL
sgemm(
'Transpose',
'No transpose', m, n, m, one, u, ldu,
292 $ work, lda, zero, a, lda )
296 a( i, n-k-l+j ) = a( i, n-k-l+j ) - r( i, j )
300 DO 80 i = k + 1, min( k+l, m )
302 a( i, n-k-l+j ) = a( i, n-k-l+j ) - alpha( i )*r( i, j )
308 resid = slange(
'1', m, n, a, lda, rwork )
310 IF( anorm.GT.zero )
THEN
311 result( 1 ) = ( ( resid /
REAL( MAX( 1, M, N ) ) ) / anorm ) /
319 CALL
sgemm(
'No transpose',
'No transpose', p, n, n, one, b, ldb,
320 $ q, ldq, zero, work, ldb )
322 CALL
sgemm(
'Transpose',
'No transpose', p, n, p, one, v, ldv,
323 $ work, ldb, zero, b, ldb )
327 b( i, n-l+j ) = b( i, n-l+j ) - beta( k+i )*r( k+i, k+j )
333 resid = slange(
'1', p, n, b, ldb, rwork )
334 IF( bnorm.GT.zero )
THEN
335 result( 2 ) = ( ( resid /
REAL( MAX( 1, P, N ) ) ) / bnorm ) /
343 CALL
slaset(
'Full', m, m, zero, one, work, ldq )
344 CALL
ssyrk(
'Upper',
'Transpose', m, m, -one, u, ldu, one, work,
349 resid = slansy(
'1',
'Upper', m, work, ldu, rwork )
350 result( 3 ) = ( resid /
REAL( MAX( 1, M ) ) ) / ulp
354 CALL
slaset(
'Full', p, p, zero, one, work, ldv )
355 CALL
ssyrk(
'Upper',
'Transpose', p, p, -one, v, ldv, one, work,
360 resid = slansy(
'1',
'Upper', p, work, ldv, rwork )
361 result( 4 ) = ( resid /
REAL( MAX( 1, P ) ) ) / ulp
365 CALL
slaset(
'Full', n, n, zero, one, work, ldq )
366 CALL
ssyrk(
'Upper',
'Transpose', n, n, -one, q, ldq, one, work,
371 resid = slansy(
'1',
'Upper', n, work, ldq, rwork )
372 result( 5 ) = ( resid /
REAL( MAX( 1, N ) ) ) / ulp
376 CALL
scopy( n, alpha, 1, work, 1 )
377 DO 110 i = k + 1, min( k+l, m )
381 work( i ) = work( j )
387 DO 120 i = k + 1, min( k+l, m ) - 1
388 IF( work( i ).LT.work( i+1 ) )
389 $ result( 6 ) = ulpinv
subroutine slaset(UPLO, M, N, ALPHA, BETA, A, LDA)
SLASET initializes the off-diagonal elements and the diagonal elements of a matrix to given values...
subroutine ssyrk(UPLO, TRANS, N, K, ALPHA, A, LDA, BETA, C, LDC)
SSYRK
subroutine sgsvts(M, P, N, A, AF, LDA, B, BF, LDB, U, LDU, V, LDV, Q, LDQ, ALPHA, BETA, R, LDR, IWORK, WORK, LWORK, RWORK, RESULT)
SGSVTS
subroutine scopy(N, SX, INCX, SY, INCY)
SCOPY
subroutine slacpy(UPLO, M, N, A, LDA, B, LDB)
SLACPY copies all or part of one two-dimensional array to another.
subroutine sgemm(TRANSA, TRANSB, M, N, K, ALPHA, A, LDA, B, LDB, BETA, C, LDC)
SGEMM
subroutine sggsvd(JOBU, JOBV, JOBQ, M, N, P, K, L, A, LDA, B, LDB, ALPHA, BETA, U, LDU, V, LDV, Q, LDQ, WORK, IWORK, INFO)
SGGSVD computes the singular value decomposition (SVD) for OTHER matrices