208 SUBROUTINE cgsvts( 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
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
cggsvd(
'U',
'V',
'Q', m, n, p, k, l, af, lda, bf, ldb,
268 $ alpha, beta, u, ldu, v, ldv, q, ldq, work, rwork,
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