208 SUBROUTINE zgsvts( 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 DOUBLE PRECISION alpha( * ), beta( * ), result( 6 ), rwork( * )
223 COMPLEX*16 a( lda, * ), af( lda, * ),
b( ldb, * ),
224 $ bf( ldb, * ), q( ldq, * ), r( ldr, * ),
225 $ u( ldu, * ), v( ldv, * ), work( lwork )
231 DOUBLE PRECISION zero, one
232 parameter( zero = 0.0d+0, one = 1.0d+0 )
233 COMPLEX*16 czero, cone
234 parameter( czero = ( 0.0d+0, 0.0d+0 ),
235 $ cone = ( 1.0d+0, 0.0d+0 ) )
238 INTEGER i, info,
j, k, l
239 DOUBLE PRECISION anorm, bnorm, resid, temp, ulp, ulpinv, unfl
249 INTRINSIC dble, max, min
253 ulp =
dlamch(
'Precision' )
255 unfl =
dlamch(
'Safe minimum' )
259 CALL
zlacpy(
'Full', m, n, a, lda, af, lda )
260 CALL
zlacpy(
'Full', p, n,
b, ldb, bf, ldb )
262 anorm = max(
zlange(
'1', m, n, a, lda, rwork ), unfl )
263 bnorm = max(
zlange(
'1', p, n,
b, ldb, rwork ), unfl )
267 CALL
zggsvd(
'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
zgemm(
'No transpose',
'No transpose', m, n, n, cone, a, lda,
290 $ q, ldq, czero, work, lda )
292 CALL
zgemm(
'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 =
zlange(
'1', m, n, a, lda, rwork )
310 IF( anorm.GT.zero )
THEN
311 result( 1 ) = ( ( resid / dble( max( 1, m, n ) ) ) / anorm ) /
319 CALL
zgemm(
'No transpose',
'No transpose', p, n, n, cone,
b, ldb,
320 $ q, ldq, czero, work, ldb )
322 CALL
zgemm(
'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 =
zlange(
'1', p, n,
b, ldb, rwork )
334 IF( bnorm.GT.zero )
THEN
335 result( 2 ) = ( ( resid / dble( max( 1, p, n ) ) ) / bnorm ) /
343 CALL
zlaset(
'Full', m, m, czero, cone, work, ldq )
344 CALL
zherk(
'Upper',
'Conjugate transpose', m, m, -one, u, ldu,
349 resid =
zlanhe(
'1',
'Upper', m, work, ldu, rwork )
350 result( 3 ) = ( resid / dble( max( 1, m ) ) ) / ulp
354 CALL
zlaset(
'Full', p, p, czero, cone, work, ldv )
355 CALL
zherk(
'Upper',
'Conjugate transpose', p, p, -one, v, ldv,
360 resid =
zlanhe(
'1',
'Upper', p, work, ldv, rwork )
361 result( 4 ) = ( resid / dble( max( 1, p ) ) ) / ulp
365 CALL
zlaset(
'Full', n, n, czero, cone, work, ldq )
366 CALL
zherk(
'Upper',
'Conjugate transpose', n, n, -one, q, ldq,
371 resid =
zlanhe(
'1',
'Upper', n, work, ldq, rwork )
372 result( 5 ) = ( resid / dble( max( 1, n ) ) ) / ulp
376 CALL
dcopy( 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 zherk(UPLO, TRANS, N, K, ALPHA, A, LDA, BETA, C, LDC)
ZHERK
subroutine zlacpy(UPLO, M, N, A, LDA, B, LDB)
ZLACPY copies all or part of one two-dimensional array to another.
subroutine dcopy(N, DX, INCX, DY, INCY)
DCOPY
double precision function zlanhe(NORM, UPLO, N, A, LDA, WORK)
ZLANHE returns the value of the 1-norm, or the Frobenius norm, or the infinity norm, or the element of largest absolute value of a complex Hermitian matrix.
subroutine zgemm(TRANSA, TRANSB, M, N, K, ALPHA, A, LDA, B, LDB, BETA, C, LDC)
ZGEMM
subroutine zggsvd(JOBU, JOBV, JOBQ, M, N, P, K, L, A, LDA, B, LDB, ALPHA, BETA, U, LDU, V, LDV, Q, LDQ, WORK, RWORK, IWORK, INFO)
ZGGSVD computes the singular value decomposition (SVD) for OTHER matrices
set ue cd $ADTTMP cat<< EOF > tmp f Program LinearEquations Implicit none Real b(3) integer i
subroutine zlaset(UPLO, M, N, ALPHA, BETA, A, LDA)
ZLASET initializes the off-diagonal elements and the diagonal elements of a matrix to given values...
subroutine zgsvts(M, P, N, A, AF, LDA, B, BF, LDB, U, LDU, V, LDV, Q, LDQ, ALPHA, BETA, R, LDR, IWORK, WORK, LWORK, RWORK, RESULT)
ZGSVTS
double precision function zlange(NORM, M, N, A, LDA, WORK)
ZLANGE returns the value of the 1-norm, Frobenius norm, infinity-norm, or the largest absolute value ...
double precision function dlamch(CMACH)
DLAMCH
set ue cd $ADTTMP cat<< EOF > tmp f Program LinearEquations Implicit none Real j