LAPACK 3.12.0
LAPACK: Linear Algebra PACKage
Loading...
Searching...
No Matches

◆ 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*
dcopy
subroutine dcopy(n, dx, incx, dy, incy)
DCOPY
Definition dcopy.f:82
dgemm
subroutine dgemm(transa, transb, m, n, k, alpha, a, lda, b, ldb, beta, c, ldc)
DGEMM
Definition dgemm.f:188
dggsvd3
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
dsyrk
subroutine dsyrk(uplo, trans, n, k, alpha, a, lda, beta, c, ldc)
DSYRK
Definition dsyrk.f:169
dlacpy
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
dlamch
double precision function dlamch(cmach)
DLAMCH
Definition dlamch.f:69
dlange
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
dlansy
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
dlaset
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
Here is the call graph for this function:
Here is the caller graph for this function: