LAPACK  3.6.1
LAPACK: Linear Algebra PACKage
subroutine sgsvts3 ( integer  M,
integer  P,
integer  N,
real, dimension( lda, * )  A,
real, dimension( lda, * )  AF,
integer  LDA,
real, dimension( ldb, * )  B,
real, dimension( ldb, * )  BF,
integer  LDB,
real, dimension( ldu, * )  U,
integer  LDU,
real, dimension( ldv, * )  V,
integer  LDV,
real, dimension( ldq, * )  Q,
integer  LDQ,
real, dimension( * )  ALPHA,
real, dimension( * )  BETA,
real, dimension( ldr, * )  R,
integer  LDR,
integer, dimension( * )  IWORK,
real, dimension( lwork )  WORK,
integer  LWORK,
real, dimension( * )  RWORK,
real, dimension( 6 )  RESULT 
)

SGSVTS3

Purpose: 
 SGSVTS3 tests SGGSVD3, 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 REAL array, dimension (LDA,M)
          The M-by-N matrix A.
[out]AF
          AF is REAL array, dimension (LDA,N)
          Details of the GSVD of A and B, as returned by SGGSVD3,
          see SGGSVD3 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 REAL array, dimension (LDB,P)
          On entry, the P-by-N matrix B.
[out]BF
          BF is REAL array, dimension (LDB,N)
          Details of the GSVD of A and B, as returned by SGGSVD3,
          see SGGSVD3 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 REAL 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 REAL 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 REAL 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 REAL array, dimension (N)
[out]BETA
          BETA is REAL 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 SGGSVD3 for details.
[out]R
          R is REAL 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 REAL 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 REAL array, dimension (max(M,P,N))
[out]RESULT
          RESULT is REAL 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.
Date
August 2015

Definition at line 212 of file sgsvts3.f.

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

Here is the call graph for this function:

Here is the caller graph for this function: