209 SUBROUTINE dgsvts3( 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 DOUBLE PRECISION A( lda, * ), AF( lda, * ), ALPHA( * ),
224 $ b( ldb, * ), beta( * ), bf( ldb, * ),
225 $ q( ldq, * ), r( ldr, * ), result( 6 ),
226 $ rwork( * ), u( ldu, * ), v( ldv, * ),
233 DOUBLE PRECISION ZERO, ONE
234 parameter ( zero = 0.0d+0, one = 1.0d+0 )
237 INTEGER I, INFO, J, K, L
238 DOUBLE PRECISION ANORM, BNORM, RESID, TEMP, ULP, ULPINV, UNFL
241 DOUBLE PRECISION DLAMCH, DLANGE, DLANSY
242 EXTERNAL dlamch, dlange, dlansy
248 INTRINSIC dble, max, min
252 ulp = dlamch(
'Precision' )
254 unfl = dlamch(
'Safe minimum' )
258 CALL dlacpy(
'Full', m, n, a, lda, af, lda )
259 CALL dlacpy(
'Full', p, n, b, ldb, bf, ldb )
261 anorm = max( dlange(
'1', m, n, a, lda, rwork ), unfl )
262 bnorm = max( dlange(
'1', p, n, b, ldb, rwork ), unfl )
266 CALL dggsvd3(
'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 dgemm(
'No transpose',
'No transpose', m, n, n, one, a, lda,
289 $ q, ldq, zero, work, lda )
291 CALL dgemm(
'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 = dlange(
'1', m, n, a, lda, rwork )
310 IF( anorm.GT.zero )
THEN
311 result( 1 ) = ( ( resid / dble( max( 1, m, n ) ) ) / anorm ) /
319 CALL dgemm(
'No transpose',
'No transpose', p, n, n, one, b, ldb,
320 $ q, ldq, zero, work, ldb )
322 CALL dgemm(
'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 = dlange(
'1', p, n, b, ldb, rwork )
334 IF( bnorm.GT.zero )
THEN
335 result( 2 ) = ( ( resid / dble( max( 1, p, n ) ) ) / bnorm ) /
343 CALL dlaset(
'Full', m, m, zero, one, work, ldq )
344 CALL dsyrk(
'Upper',
'Transpose', m, m, -one, u, ldu, one, work,
349 resid = dlansy(
'1',
'Upper', m, work, ldu, rwork )
350 result( 3 ) = ( resid / dble( max( 1, m ) ) ) / ulp
354 CALL dlaset(
'Full', p, p, zero, one, work, ldv )
355 CALL dsyrk(
'Upper',
'Transpose', p, p, -one, v, ldv, one, work,
360 resid = dlansy(
'1',
'Upper', p, work, ldv, rwork )
361 result( 4 ) = ( resid / dble( max( 1, p ) ) ) / ulp
365 CALL dlaset(
'Full', n, n, zero, one, work, ldq )
366 CALL dsyrk(
'Upper',
'Transpose', n, n, -one, q, ldq, one, work,
371 resid = dlansy(
'1',
'Upper', n, work, ldq, rwork )
372 result( 5 ) = ( resid / dble( max( 1, n ) ) ) / ulp
376 CALL dcopy( 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 dlacpy(UPLO, M, N, A, LDA, B, LDB)
DLACPY copies all or part of one two-dimensional array to another.
subroutine dcopy(N, DX, INCX, DY, INCY)
DCOPY
subroutine dgsvts3(M, P, N, A, AF, LDA, B, BF, LDB, U, LDU, V, LDV, Q, LDQ, ALPHA, BETA, R, LDR, IWORK, WORK, LWORK, RWORK, RESULT)
DGSVTS3
subroutine dlaset(UPLO, M, N, ALPHA, BETA, A, LDA)
DLASET initializes the off-diagonal elements and the diagonal elements of a matrix to given values...
subroutine dgemm(TRANSA, TRANSB, M, N, K, ALPHA, A, LDA, B, LDB, BETA, C, LDC)
DGEMM
subroutine dsyrk(UPLO, TRANS, N, K, ALPHA, A, LDA, BETA, C, LDC)
DSYRK
subroutine dggsvd3(JOBU, JOBV, JOBQ, M, N, P, K, L, A, LDA, B, LDB, ALPHA, BETA, U, LDU, V, LDV, Q, LDQ, WORK, LWORK, IWORK, INFO)
DGGSVD3 computes the singular value decomposition (SVD) for OTHER matrices