208 SUBROUTINE cgsvts3( M, P, N, A, AF, LDA, B, BF, LDB, U, LDU, V,
209 $ ldv, q, ldq, alpha, beta, r, ldr, iwork, work,
210 $ lwork, rwork, result )
218 INTEGER LDA, LDB, LDQ, LDR, LDU, LDV, LWORK, M, N, P
222 REAL ALPHA( * ), BETA( * ), RESULT( 6 ), RWORK( * )
223 COMPLEX A( lda, * ), AF( lda, * ), B( ldb, * ),
224 $ bf( ldb, * ), q( ldq, * ), r( ldr, * ),
225 $ u( ldu, * ), v( ldv, * ), work( lwork )
232 parameter ( zero = 0.0e+0, one = 1.0e+0 )
234 parameter ( czero = ( 0.0e+0, 0.0e+0 ),
235 $ cone = ( 1.0e+0, 0.0e+0 ) )
238 INTEGER I, INFO, J, K, L
239 REAL ANORM, BNORM, RESID, TEMP, ULP, ULPINV, UNFL
242 REAL CLANGE, CLANHE, SLAMCH
243 EXTERNAL clange, clanhe, slamch
249 INTRINSIC max, min, real
253 ulp = slamch(
'Precision' )
255 unfl = slamch(
'Safe minimum' )
259 CALL clacpy(
'Full', m, n, a, lda, af, lda )
260 CALL clacpy(
'Full', p, n, b, ldb, bf, ldb )
262 anorm = max( clange(
'1', m, n, a, lda, rwork ), unfl )
263 bnorm = max( clange(
'1', p, n, b, ldb, rwork ), unfl )
267 CALL cggsvd3(
'U',
'V',
'Q', m, n, p, k, l, af, lda, bf, ldb,
268 $ alpha, beta, u, ldu, v, ldv, q, ldq, work, lwork,
269 $ rwork, iwork, info )
273 DO 20 i = 1, min( k+l, m )
275 r( i, j ) = af( i, n-k-l+j )
279 IF( m-k-l.LT.0 )
THEN
280 DO 40 i = m + 1, k + l
282 r( i, j ) = bf( i-k, n-k-l+j )
289 CALL cgemm(
'No transpose',
'No transpose', m, n, n, cone, a, lda,
290 $ q, ldq, czero, work, lda )
292 CALL cgemm(
'Conjugate transpose',
'No transpose', m, n, m, cone,
293 $ u, ldu, work, lda, czero, a, lda )
297 a( i, n-k-l+j ) = a( i, n-k-l+j ) - r( i, j )
301 DO 80 i = k + 1, min( k+l, m )
303 a( i, n-k-l+j ) = a( i, n-k-l+j ) - alpha( i )*r( i, j )
309 resid = clange(
'1', m, n, a, lda, rwork )
310 IF( anorm.GT.zero )
THEN
311 result( 1 ) = ( ( resid /
REAL( MAX( 1, M, N ) ) ) / anorm ) /
319 CALL cgemm(
'No transpose',
'No transpose', p, n, n, cone, b, ldb,
320 $ q, ldq, czero, work, ldb )
322 CALL cgemm(
'Conjugate transpose',
'No transpose', p, n, p, cone,
323 $ v, ldv, work, ldb, czero, b, ldb )
327 b( i, n-l+j ) = b( i, n-l+j ) - beta( k+i )*r( k+i, k+j )
333 resid = clange(
'1', p, n, b, ldb, rwork )
334 IF( bnorm.GT.zero )
THEN
335 result( 2 ) = ( ( resid /
REAL( MAX( 1, P, N ) ) ) / bnorm ) /
343 CALL claset(
'Full', m, m, czero, cone, work, ldq )
344 CALL cherk(
'Upper',
'Conjugate transpose', m, m, -one, u, ldu,
349 resid = clanhe(
'1',
'Upper', m, work, ldu, rwork )
350 result( 3 ) = ( resid /
REAL( MAX( 1, M ) ) ) / ulp
354 CALL claset(
'Full', p, p, czero, cone, work, ldv )
355 CALL cherk(
'Upper',
'Conjugate transpose', p, p, -one, v, ldv,
360 resid = clanhe(
'1',
'Upper', p, work, ldv, rwork )
361 result( 4 ) = ( resid /
REAL( MAX( 1, P ) ) ) / ulp
365 CALL claset(
'Full', n, n, czero, cone, work, ldq )
366 CALL cherk(
'Upper',
'Conjugate transpose', n, n, -one, q, ldq,
371 resid = clanhe(
'1',
'Upper', n, work, ldq, rwork )
372 result( 5 ) = ( resid /
REAL( MAX( 1, N ) ) ) / ulp
376 CALL scopy( n, alpha, 1, rwork, 1 )
377 DO 110 i = k + 1, min( k+l, m )
381 rwork( i ) = rwork( j )
387 DO 120 i = k + 1, min( k+l, m ) - 1
388 IF( rwork( i ).LT.rwork( i+1 ) )
389 $ result( 6 ) = ulpinv
subroutine claset(UPLO, M, N, ALPHA, BETA, A, LDA)
CLASET initializes the off-diagonal elements and the diagonal elements of a matrix to given values...
subroutine cgsvts3(M, P, N, A, AF, LDA, B, BF, LDB, U, LDU, V, LDV, Q, LDQ, ALPHA, BETA, R, LDR, IWORK, WORK, LWORK, RWORK, RESULT)
CGSVTS3
subroutine cggsvd3(JOBU, JOBV, JOBQ, M, N, P, K, L, A, LDA, B, LDB, ALPHA, BETA, U, LDU, V, LDV, Q, LDQ, WORK, LWORK, RWORK, IWORK, INFO)
CGGSVD3 computes the singular value decomposition (SVD) for OTHER matrices
subroutine cherk(UPLO, TRANS, N, K, ALPHA, A, LDA, BETA, C, LDC)
CHERK
subroutine clacpy(UPLO, M, N, A, LDA, B, LDB)
CLACPY copies all or part of one two-dimensional array to another.
subroutine cgemm(TRANSA, TRANSB, M, N, K, ALPHA, A, LDA, B, LDB, BETA, C, LDC)
CGEMM
subroutine scopy(N, SX, INCX, SY, INCY)
SCOPY