209 SUBROUTINE sgsvts( 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 REAL a( lda, * ), af( lda, * ), alpha( * ),
224 $
b( ldb, * ), beta( * ), bf( ldb, * ),
225 $ q( ldq, * ), r( ldr, * ), result( 6 ),
226 $ rwork( * ), u( ldu, * ), v( ldv, * ),
234 parameter( zero = 0.0e+0, one = 1.0e+0 )
237 INTEGER i, info,
j, k, l
238 REAL anorm, bnorm, resid, temp, ulp, ulpinv, unfl
248 INTRINSIC max, min, real
252 ulp =
slamch(
'Precision' )
254 unfl =
slamch(
'Safe minimum' )
258 CALL
slacpy(
'Full', m, n, a, lda, af, lda )
259 CALL
slacpy(
'Full', p, n,
b, ldb, bf, ldb )
261 anorm = max(
slange(
'1', m, n, a, lda, rwork ), unfl )
262 bnorm = max(
slange(
'1', p, n,
b, ldb, rwork ), unfl )
266 CALL
sggsvd(
'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
sgemm(
'No transpose',
'No transpose', m, n, n, one, a, lda,
289 $ q, ldq, zero, work, lda )
291 CALL
sgemm(
'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 =
slange(
'1', m, n, a, lda, rwork )
310 IF( anorm.GT.zero )
THEN
311 result( 1 ) = ( ( resid /
REAL( MAX( 1, M, N ) ) ) / anorm ) /
319 CALL
sgemm(
'No transpose',
'No transpose', p, n, n, one,
b, ldb,
320 $ q, ldq, zero, work, ldb )
322 CALL
sgemm(
'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 =
slange(
'1', p, n,
b, ldb, rwork )
334 IF( bnorm.GT.zero )
THEN
335 result( 2 ) = ( ( resid /
REAL( MAX( 1, P, N ) ) ) / bnorm ) /
343 CALL
slaset(
'Full', m, m, zero, one, work, ldq )
344 CALL
ssyrk(
'Upper',
'Transpose', m, m, -one, u, ldu, one, work,
349 resid =
slansy(
'1',
'Upper', m, work, ldu, rwork )
350 result( 3 ) = ( resid /
REAL( MAX( 1, M ) ) ) / ulp
354 CALL
slaset(
'Full', p, p, zero, one, work, ldv )
355 CALL
ssyrk(
'Upper',
'Transpose', p, p, -one, v, ldv, one, work,
360 resid =
slansy(
'1',
'Upper', p, work, ldv, rwork )
361 result( 4 ) = ( resid /
REAL( MAX( 1, P ) ) ) / ulp
365 CALL
slaset(
'Full', n, n, zero, one, work, ldq )
366 CALL
ssyrk(
'Upper',
'Transpose', n, n, -one, q, ldq, one, work,
371 resid =
slansy(
'1',
'Upper', n, work, ldq, rwork )
372 result( 5 ) = ( resid /
REAL( MAX( 1, N ) ) ) / ulp
376 CALL
scopy( 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 slaset(UPLO, M, N, ALPHA, BETA, A, LDA)
SLASET initializes the off-diagonal elements and the diagonal elements of a matrix to given values...
REAL function slansy(NORM, UPLO, N, A, LDA, WORK)
SLANSY 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.
REAL function slamch(CMACH)
SLAMCH
subroutine ssyrk(UPLO, TRANS, N, K, ALPHA, A, LDA, BETA, C, LDC)
SSYRK
subroutine sgsvts(M, P, N, A, AF, LDA, B, BF, LDB, U, LDU, V, LDV, Q, LDQ, ALPHA, BETA, R, LDR, IWORK, WORK, LWORK, RWORK, RESULT)
SGSVTS
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 slacpy(UPLO, M, N, A, LDA, B, LDB)
SLACPY copies all or part of one two-dimensional array to another.
REAL function slange(NORM, M, N, A, LDA, WORK)
SLANGE returns the value of the 1-norm, Frobenius norm, infinity-norm, or the largest absolute value ...
subroutine sgemm(TRANSA, TRANSB, M, N, K, ALPHA, A, LDA, B, LDB, BETA, C, LDC)
SGEMM
set ue cd $ADTTMP cat<< EOF > tmp f Program LinearEquations Implicit none Real j
subroutine sggsvd(JOBU, JOBV, JOBQ, M, N, P, K, L, A, LDA, B, LDB, ALPHA, BETA, U, LDU, V, LDV, Q, LDQ, WORK, IWORK, INFO)
SGGSVD computes the singular value decomposition (SVD) for OTHER matrices