209 SUBROUTINE dgsvts( 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
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
dggsvd(
'U',
'V',
'Q', m, n, p, k, l, af, lda, bf, ldb,
267 $ alpha, beta, u, ldu, v, ldv, q, ldq, work, iwork,
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 dgemm(TRANSA, TRANSB, M, N, K, ALPHA, A, LDA, B, LDB, BETA, C, LDC)
DGEMM
subroutine dcopy(N, DX, INCX, DY, INCY)
DCOPY
set ue cd $ADTTMP cat<< EOF > tmp f Program LinearEquations Implicit none Real b(3) integer i
subroutine dggsvd(JOBU, JOBV, JOBQ, M, N, P, K, L, A, LDA, B, LDB, ALPHA, BETA, U, LDU, V, LDV, Q, LDQ, WORK, IWORK, INFO)
DGGSVD computes the singular value decomposition (SVD) for OTHER matrices
double precision function dlange(NORM, M, N, A, LDA, WORK)
DLANGE returns the value of the 1-norm, Frobenius norm, infinity-norm, or the largest absolute value ...
subroutine dsyrk(UPLO, TRANS, N, K, ALPHA, A, LDA, BETA, C, LDC)
DSYRK
subroutine dgsvts(M, P, N, A, AF, LDA, B, BF, LDB, U, LDU, V, LDV, Q, LDQ, ALPHA, BETA, R, LDR, IWORK, WORK, LWORK, RWORK, RESULT)
DGSVTS
double precision function dlamch(CMACH)
DLAMCH
set ue cd $ADTTMP cat<< EOF > tmp f Program LinearEquations Implicit none Real j
double precision function dlansy(NORM, UPLO, N, A, LDA, WORK)
DLANSY returns the value of the 1-norm, or the Frobenius norm, or the infinity norm, or the element of largest absolute value of a real symmetric matrix.
subroutine dlacpy(UPLO, M, N, A, LDA, B, LDB)
DLACPY copies all or part of one two-dimensional array to another.
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...