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
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...
real function clange(NORM, M, N, A, LDA, WORK)
CLANGE returns the value of the 1-norm, Frobenius norm, infinity-norm, or the largest absolute value ...
real function clanhe(NORM, UPLO, N, A, LDA, WORK)
CLANHE 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.
set ue cd $ADTTMP cat<< EOF > tmp f Program LinearEquations Implicit none Real b(3) integer i
subroutine scopy(N, SX, INCX, SY, INCY)
SCOPY
subroutine cgsvts(M, P, N, A, AF, LDA, B, BF, LDB, U, LDU, V, LDV, Q, LDQ, ALPHA, BETA, R, LDR, IWORK, WORK, LWORK, RWORK, RESULT)
CGSVTS
subroutine clacpy(UPLO, M, N, A, LDA, B, LDB)
CLACPY copies all or part of one two-dimensional array to another.
real function slamch(CMACH)
SLAMCH
set ue cd $ADTTMP cat<< EOF > tmp f Program LinearEquations Implicit none Real j
subroutine cggsvd(JOBU, JOBV, JOBQ, M, N, P, K, L, A, LDA, B, LDB, ALPHA, BETA, U, LDU, V, LDV, Q, LDQ, WORK, RWORK, IWORK, INFO)
CGGSVD computes the singular value decomposition (SVD) for OTHER matrices
subroutine cgemm(TRANSA, TRANSB, M, N, K, ALPHA, A, LDA, B, LDB, BETA, C, LDC)
CGEMM
subroutine cherk(UPLO, TRANS, N, K, ALPHA, A, LDA, BETA, C, LDC)
CHERK