LAPACK  3.10.1
LAPACK: Linear Algebra PACKage

◆ dgsvts3()

subroutine dgsvts3 ( integer  M,
integer  P,
integer  N,
double precision, dimension( lda, * )  A,
double precision, dimension( lda, * )  AF,
integer  LDA,
double precision, dimension( ldb, * )  B,
double precision, dimension( ldb, * )  BF,
integer  LDB,
double precision, dimension( ldu, * )  U,
integer  LDU,
double precision, dimension( ldv, * )  V,
integer  LDV,
double precision, dimension( ldq, * )  Q,
integer  LDQ,
double precision, dimension( * )  ALPHA,
double precision, dimension( * )  BETA,
double precision, dimension( ldr, * )  R,
integer  LDR,
integer, dimension( * )  IWORK,
double precision, dimension( lwork )  WORK,
integer  LWORK,
double precision, dimension( * )  RWORK,
double precision, dimension( 6 )  RESULT 
)

DGSVTS3

Purpose:
 DGSVTS3 tests DGGSVD3, which computes the GSVD of an M-by-N matrix A
 and a P-by-N matrix B:
              U'*A*Q = D1*R and V'*B*Q = D2*R.
Parameters
[in]M
          M is INTEGER
          The number of rows of the matrix A.  M >= 0.
[in]P
          P is INTEGER
          The number of rows of the matrix B.  P >= 0.
[in]N
          N is INTEGER
          The number of columns of the matrices A and B.  N >= 0.
[in]A
          A is DOUBLE PRECISION array, dimension (LDA,M)
          The M-by-N matrix A.
[out]AF
          AF is DOUBLE PRECISION array, dimension (LDA,N)
          Details of the GSVD of A and B, as returned by DGGSVD3,
          see DGGSVD3 for further details.
[in]LDA
          LDA is INTEGER
          The leading dimension of the arrays A and AF.
          LDA >= max( 1,M ).
[in]B
          B is DOUBLE PRECISION array, dimension (LDB,P)
          On entry, the P-by-N matrix B.
[out]BF
          BF is DOUBLE PRECISION array, dimension (LDB,N)
          Details of the GSVD of A and B, as returned by DGGSVD3,
          see DGGSVD3 for further details.
[in]LDB
          LDB is INTEGER
          The leading dimension of the arrays B and BF.
          LDB >= max(1,P).
[out]U
          U is DOUBLE PRECISION array, dimension(LDU,M)
          The M by M orthogonal matrix U.
[in]LDU
          LDU is INTEGER
          The leading dimension of the array U. LDU >= max(1,M).
[out]V
          V is DOUBLE PRECISION array, dimension(LDV,M)
          The P by P orthogonal matrix V.
[in]LDV
          LDV is INTEGER
          The leading dimension of the array V. LDV >= max(1,P).
[out]Q
          Q is DOUBLE PRECISION array, dimension(LDQ,N)
          The N by N orthogonal matrix Q.
[in]LDQ
          LDQ is INTEGER
          The leading dimension of the array Q. LDQ >= max(1,N).
[out]ALPHA
          ALPHA is DOUBLE PRECISION array, dimension (N)
[out]BETA
          BETA is DOUBLE PRECISION array, dimension (N)

          The generalized singular value pairs of A and B, the
          ``diagonal'' matrices D1 and D2 are constructed from
          ALPHA and BETA, see subroutine DGGSVD3 for details.
[out]R
          R is DOUBLE PRECISION array, dimension(LDQ,N)
          The upper triangular matrix R.
[in]LDR
          LDR is INTEGER
          The leading dimension of the array R. LDR >= max(1,N).
[out]IWORK
          IWORK is INTEGER array, dimension (N)
[out]WORK
          WORK is DOUBLE PRECISION array, dimension (LWORK)
[in]LWORK
          LWORK is INTEGER
          The dimension of the array WORK,
          LWORK >= max(M,P,N)*max(M,P,N).
[out]RWORK
          RWORK is DOUBLE PRECISION array, dimension (max(M,P,N))
[out]RESULT
          RESULT is DOUBLE PRECISION array, dimension (6)
          The test ratios:
          RESULT(1) = norm( U'*A*Q - D1*R ) / ( MAX(M,N)*norm(A)*ULP)
          RESULT(2) = norm( V'*B*Q - D2*R ) / ( MAX(P,N)*norm(B)*ULP)
          RESULT(3) = norm( I - U'*U ) / ( M*ULP )
          RESULT(4) = norm( I - V'*V ) / ( P*ULP )
          RESULT(5) = norm( I - Q'*Q ) / ( N*ULP )
          RESULT(6) = 0        if ALPHA is in decreasing order;
                    = ULPINV   otherwise.
Author
Univ. of Tennessee
Univ. of California Berkeley
Univ. of Colorado Denver
NAG Ltd.

Definition at line 207 of file dgsvts3.f.

210 *
211 * -- LAPACK test routine --
212 * -- LAPACK is a software package provided by Univ. of Tennessee, --
213 * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
214 *
215 * .. Scalar Arguments ..
216  INTEGER LDA, LDB, LDQ, LDR, LDU, LDV, LWORK, M, N, P
217 * ..
218 * .. Array Arguments ..
219  INTEGER IWORK( * )
220  DOUBLE PRECISION A( LDA, * ), AF( LDA, * ), ALPHA( * ),
221  $ B( LDB, * ), BETA( * ), BF( LDB, * ),
222  $ Q( LDQ, * ), R( LDR, * ), RESULT( 6 ),
223  $ RWORK( * ), U( LDU, * ), V( LDV, * ),
224  $ WORK( LWORK )
225 * ..
226 *
227 * =====================================================================
228 *
229 * .. Parameters ..
230  DOUBLE PRECISION ZERO, ONE
231  parameter( zero = 0.0d+0, one = 1.0d+0 )
232 * ..
233 * .. Local Scalars ..
234  INTEGER I, INFO, J, K, L
235  DOUBLE PRECISION ANORM, BNORM, RESID, TEMP, ULP, ULPINV, UNFL
236 * ..
237 * .. External Functions ..
238  DOUBLE PRECISION DLAMCH, DLANGE, DLANSY
239  EXTERNAL dlamch, dlange, dlansy
240 * ..
241 * .. External Subroutines ..
242  EXTERNAL dcopy, dgemm, dggsvd3, dlacpy, dlaset, dsyrk
243 * ..
244 * .. Intrinsic Functions ..
245  INTRINSIC dble, max, min
246 * ..
247 * .. Executable Statements ..
248 *
249  ulp = dlamch( 'Precision' )
250  ulpinv = one / ulp
251  unfl = dlamch( 'Safe minimum' )
252 *
253 * Copy the matrix A to the array AF.
254 *
255  CALL dlacpy( 'Full', m, n, a, lda, af, lda )
256  CALL dlacpy( 'Full', p, n, b, ldb, bf, ldb )
257 *
258  anorm = max( dlange( '1', m, n, a, lda, rwork ), unfl )
259  bnorm = max( dlange( '1', p, n, b, ldb, rwork ), unfl )
260 *
261 * Factorize the matrices A and B in the arrays AF and BF.
262 *
263  CALL dggsvd3( 'U', 'V', 'Q', m, n, p, k, l, af, lda, bf, ldb,
264  $ alpha, beta, u, ldu, v, ldv, q, ldq, work, lwork,
265  $ iwork, info )
266 *
267 * Copy R
268 *
269  DO 20 i = 1, min( k+l, m )
270  DO 10 j = i, k + l
271  r( i, j ) = af( i, n-k-l+j )
272  10 CONTINUE
273  20 CONTINUE
274 *
275  IF( m-k-l.LT.0 ) THEN
276  DO 40 i = m + 1, k + l
277  DO 30 j = i, k + l
278  r( i, j ) = bf( i-k, n-k-l+j )
279  30 CONTINUE
280  40 CONTINUE
281  END IF
282 *
283 * Compute A:= U'*A*Q - D1*R
284 *
285  CALL dgemm( 'No transpose', 'No transpose', m, n, n, one, a, lda,
286  $ q, ldq, zero, work, lda )
287 *
288  CALL dgemm( 'Transpose', 'No transpose', m, n, m, one, u, ldu,
289  $ work, lda, zero, a, lda )
290 *
291  DO 60 i = 1, k
292  DO 50 j = i, k + l
293  a( i, n-k-l+j ) = a( i, n-k-l+j ) - r( i, j )
294  50 CONTINUE
295  60 CONTINUE
296 *
297  DO 80 i = k + 1, min( k+l, m )
298  DO 70 j = i, k + l
299  a( i, n-k-l+j ) = a( i, n-k-l+j ) - alpha( i )*r( i, j )
300  70 CONTINUE
301  80 CONTINUE
302 *
303 * Compute norm( U'*A*Q - D1*R ) / ( MAX(1,M,N)*norm(A)*ULP ) .
304 *
305  resid = dlange( '1', m, n, a, lda, rwork )
306 *
307  IF( anorm.GT.zero ) THEN
308  result( 1 ) = ( ( resid / dble( max( 1, m, n ) ) ) / anorm ) /
309  $ ulp
310  ELSE
311  result( 1 ) = zero
312  END IF
313 *
314 * Compute B := V'*B*Q - D2*R
315 *
316  CALL dgemm( 'No transpose', 'No transpose', p, n, n, one, b, ldb,
317  $ q, ldq, zero, work, ldb )
318 *
319  CALL dgemm( 'Transpose', 'No transpose', p, n, p, one, v, ldv,
320  $ work, ldb, zero, b, ldb )
321 *
322  DO 100 i = 1, l
323  DO 90 j = i, l
324  b( i, n-l+j ) = b( i, n-l+j ) - beta( k+i )*r( k+i, k+j )
325  90 CONTINUE
326  100 CONTINUE
327 *
328 * Compute norm( V'*B*Q - D2*R ) / ( MAX(P,N)*norm(B)*ULP ) .
329 *
330  resid = dlange( '1', p, n, b, ldb, rwork )
331  IF( bnorm.GT.zero ) THEN
332  result( 2 ) = ( ( resid / dble( max( 1, p, n ) ) ) / bnorm ) /
333  $ ulp
334  ELSE
335  result( 2 ) = zero
336  END IF
337 *
338 * Compute I - U'*U
339 *
340  CALL dlaset( 'Full', m, m, zero, one, work, ldq )
341  CALL dsyrk( 'Upper', 'Transpose', m, m, -one, u, ldu, one, work,
342  $ ldu )
343 *
344 * Compute norm( I - U'*U ) / ( M * ULP ) .
345 *
346  resid = dlansy( '1', 'Upper', m, work, ldu, rwork )
347  result( 3 ) = ( resid / dble( max( 1, m ) ) ) / ulp
348 *
349 * Compute I - V'*V
350 *
351  CALL dlaset( 'Full', p, p, zero, one, work, ldv )
352  CALL dsyrk( 'Upper', 'Transpose', p, p, -one, v, ldv, one, work,
353  $ ldv )
354 *
355 * Compute norm( I - V'*V ) / ( P * ULP ) .
356 *
357  resid = dlansy( '1', 'Upper', p, work, ldv, rwork )
358  result( 4 ) = ( resid / dble( max( 1, p ) ) ) / ulp
359 *
360 * Compute I - Q'*Q
361 *
362  CALL dlaset( 'Full', n, n, zero, one, work, ldq )
363  CALL dsyrk( 'Upper', 'Transpose', n, n, -one, q, ldq, one, work,
364  $ ldq )
365 *
366 * Compute norm( I - Q'*Q ) / ( N * ULP ) .
367 *
368  resid = dlansy( '1', 'Upper', n, work, ldq, rwork )
369  result( 5 ) = ( resid / dble( max( 1, n ) ) ) / ulp
370 *
371 * Check sorting
372 *
373  CALL dcopy( n, alpha, 1, work, 1 )
374  DO 110 i = k + 1, min( k+l, m )
375  j = iwork( i )
376  IF( i.NE.j ) THEN
377  temp = work( i )
378  work( i ) = work( j )
379  work( j ) = temp
380  END IF
381  110 CONTINUE
382 *
383  result( 6 ) = zero
384  DO 120 i = k + 1, min( k+l, m ) - 1
385  IF( work( i ).LT.work( i+1 ) )
386  $ result( 6 ) = ulpinv
387  120 CONTINUE
388 *
389  RETURN
390 *
391 * End of DGSVTS3
392 *
double precision function dlamch(CMACH)
DLAMCH
Definition: dlamch.f:69
subroutine dlacpy(UPLO, M, N, A, LDA, B, LDB)
DLACPY copies all or part of one two-dimensional array to another.
Definition: dlacpy.f:103
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.
Definition: dlaset.f:110
subroutine dcopy(N, DX, INCX, DY, INCY)
DCOPY
Definition: dcopy.f:82
subroutine dsyrk(UPLO, TRANS, N, K, ALPHA, A, LDA, BETA, C, LDC)
DSYRK
Definition: dsyrk.f:169
subroutine dgemm(TRANSA, TRANSB, M, N, K, ALPHA, A, LDA, B, LDB, BETA, C, LDC)
DGEMM
Definition: dgemm.f:187
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 ...
Definition: dlange.f:114
subroutine dggsvd3(JOBU, JOBV, JOBQ, M, N, P, K, L, A, LDA, B, LDB, ALPHA, BETA, U, LDU, V, LDV, Q, LDQ, WORK, LWORK, IWORK, INFO)
DGGSVD3 computes the singular value decomposition (SVD) for OTHER matrices
Definition: dggsvd3.f:349
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,...
Definition: dlansy.f:122
Here is the call graph for this function:
Here is the caller graph for this function: