209 SUBROUTINE sgsvts3( 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 sggsvd3(
'U',
'V',
'Q', m, n, p, k, l, af, lda, bf, ldb,
267 $ alpha, beta, u, ldu, v, ldv, q, ldq, work, lwork,
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 sggsvd3(JOBU, JOBV, JOBQ, M, N, P, K, L, A, LDA, B, LDB, ALPHA, BETA, U, LDU, V, LDV, Q, LDQ, WORK, LWORK, IWORK, INFO)
SGGSVD3 computes the singular value decomposition (SVD) for OTHER matrices
subroutine ssyrk(UPLO, TRANS, N, K, ALPHA, A, LDA, BETA, C, LDC)
SSYRK
subroutine sgemm(TRANSA, TRANSB, M, N, K, ALPHA, A, LDA, B, LDB, BETA, C, LDC)
SGEMM
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 sgsvts3(M, P, N, A, AF, LDA, B, BF, LDB, U, LDU, V, LDV, Q, LDQ, ALPHA, BETA, R, LDR, IWORK, WORK, LWORK, RWORK, RESULT)
SGSVTS3
subroutine slacpy(UPLO, M, N, A, LDA, B, LDB)
SLACPY copies all or part of one two-dimensional array to another.
subroutine scopy(N, SX, INCX, SY, INCY)
SCOPY