LAPACK  3.10.0
LAPACK: Linear Algebra PACKage

◆ cgsvts3()

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

CGSVTS3

Purpose:
 CGSVTS3 tests CGGSVD3, 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 COMPLEX array, dimension (LDA,M)
          The M-by-N matrix A.
[out]AF
          AF is COMPLEX array, dimension (LDA,N)
          Details of the GSVD of A and B, as returned by CGGSVD3,
          see CGGSVD3 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 COMPLEX array, dimension (LDB,P)
          On entry, the P-by-N matrix B.
[out]BF
          BF is COMPLEX array, dimension (LDB,N)
          Details of the GSVD of A and B, as returned by CGGSVD3,
          see CGGSVD3 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 COMPLEX array, dimension(LDU,M)
          The M by M unitary matrix U.
[in]LDU
          LDU is INTEGER
          The leading dimension of the array U. LDU >= max(1,M).
[out]V
          V is COMPLEX array, dimension(LDV,M)
          The P by P unitary matrix V.
[in]LDV
          LDV is INTEGER
          The leading dimension of the array V. LDV >= max(1,P).
[out]Q
          Q is COMPLEX array, dimension(LDQ,N)
          The N by N unitary 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 CGGSVD3 for details.
[out]R
          R is COMPLEX 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 COMPLEX 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.

Definition at line 206 of file cgsvts3.f.

209 *
210 * -- LAPACK test routine --
211 * -- LAPACK is a software package provided by Univ. of Tennessee, --
212 * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
213 *
214 * .. Scalar Arguments ..
215  INTEGER LDA, LDB, LDQ, LDR, LDU, LDV, LWORK, M, N, P
216 * ..
217 * .. Array Arguments ..
218  INTEGER IWORK( * )
219  REAL ALPHA( * ), BETA( * ), RESULT( 6 ), RWORK( * )
220  COMPLEX A( LDA, * ), AF( LDA, * ), B( LDB, * ),
221  $ BF( LDB, * ), Q( LDQ, * ), R( LDR, * ),
222  $ U( LDU, * ), V( LDV, * ), WORK( LWORK )
223 * ..
224 *
225 * =====================================================================
226 *
227 * .. Parameters ..
228  REAL ZERO, ONE
229  parameter( zero = 0.0e+0, one = 1.0e+0 )
230  COMPLEX CZERO, CONE
231  parameter( czero = ( 0.0e+0, 0.0e+0 ),
232  $ cone = ( 1.0e+0, 0.0e+0 ) )
233 * ..
234 * .. Local Scalars ..
235  INTEGER I, INFO, J, K, L
236  REAL ANORM, BNORM, RESID, TEMP, ULP, ULPINV, UNFL
237 * ..
238 * .. External Functions ..
239  REAL CLANGE, CLANHE, SLAMCH
240  EXTERNAL clange, clanhe, slamch
241 * ..
242 * .. External Subroutines ..
243  EXTERNAL cgemm, cggsvd3, cherk, clacpy, claset, scopy
244 * ..
245 * .. Intrinsic Functions ..
246  INTRINSIC max, min, real
247 * ..
248 * .. Executable Statements ..
249 *
250  ulp = slamch( 'Precision' )
251  ulpinv = one / ulp
252  unfl = slamch( 'Safe minimum' )
253 *
254 * Copy the matrix A to the array AF.
255 *
256  CALL clacpy( 'Full', m, n, a, lda, af, lda )
257  CALL clacpy( 'Full', p, n, b, ldb, bf, ldb )
258 *
259  anorm = max( clange( '1', m, n, a, lda, rwork ), unfl )
260  bnorm = max( clange( '1', p, n, b, ldb, rwork ), unfl )
261 *
262 * Factorize the matrices A and B in the arrays AF and BF.
263 *
264  CALL cggsvd3( 'U', 'V', 'Q', m, n, p, k, l, af, lda, bf, ldb,
265  $ alpha, beta, u, ldu, v, ldv, q, ldq, work, lwork,
266  $ rwork, iwork, info )
267 *
268 * Copy R
269 *
270  DO 20 i = 1, min( k+l, m )
271  DO 10 j = i, k + l
272  r( i, j ) = af( i, n-k-l+j )
273  10 CONTINUE
274  20 CONTINUE
275 *
276  IF( m-k-l.LT.0 ) THEN
277  DO 40 i = m + 1, k + l
278  DO 30 j = i, k + l
279  r( i, j ) = bf( i-k, n-k-l+j )
280  30 CONTINUE
281  40 CONTINUE
282  END IF
283 *
284 * Compute A:= U'*A*Q - D1*R
285 *
286  CALL cgemm( 'No transpose', 'No transpose', m, n, n, cone, a, lda,
287  $ q, ldq, czero, work, lda )
288 *
289  CALL cgemm( 'Conjugate transpose', 'No transpose', m, n, m, cone,
290  $ u, ldu, work, lda, czero, a, lda )
291 *
292  DO 60 i = 1, k
293  DO 50 j = i, k + l
294  a( i, n-k-l+j ) = a( i, n-k-l+j ) - r( i, j )
295  50 CONTINUE
296  60 CONTINUE
297 *
298  DO 80 i = k + 1, min( k+l, m )
299  DO 70 j = i, k + l
300  a( i, n-k-l+j ) = a( i, n-k-l+j ) - alpha( i )*r( i, j )
301  70 CONTINUE
302  80 CONTINUE
303 *
304 * Compute norm( U'*A*Q - D1*R ) / ( MAX(1,M,N)*norm(A)*ULP ) .
305 *
306  resid = clange( '1', m, n, a, lda, rwork )
307  IF( anorm.GT.zero ) THEN
308  result( 1 ) = ( ( resid / real( 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 cgemm( 'No transpose', 'No transpose', p, n, n, cone, b, ldb,
317  $ q, ldq, czero, work, ldb )
318 *
319  CALL cgemm( 'Conjugate transpose', 'No transpose', p, n, p, cone,
320  $ v, ldv, work, ldb, czero, 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 = clange( '1', p, n, b, ldb, rwork )
331  IF( bnorm.GT.zero ) THEN
332  result( 2 ) = ( ( resid / real( 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 claset( 'Full', m, m, czero, cone, work, ldq )
341  CALL cherk( 'Upper', 'Conjugate transpose', m, m, -one, u, ldu,
342  $ one, work, ldu )
343 *
344 * Compute norm( I - U'*U ) / ( M * ULP ) .
345 *
346  resid = clanhe( '1', 'Upper', m, work, ldu, rwork )
347  result( 3 ) = ( resid / real( max( 1, m ) ) ) / ulp
348 *
349 * Compute I - V'*V
350 *
351  CALL claset( 'Full', p, p, czero, cone, work, ldv )
352  CALL cherk( 'Upper', 'Conjugate transpose', p, p, -one, v, ldv,
353  $ one, work, ldv )
354 *
355 * Compute norm( I - V'*V ) / ( P * ULP ) .
356 *
357  resid = clanhe( '1', 'Upper', p, work, ldv, rwork )
358  result( 4 ) = ( resid / real( max( 1, p ) ) ) / ulp
359 *
360 * Compute I - Q'*Q
361 *
362  CALL claset( 'Full', n, n, czero, cone, work, ldq )
363  CALL cherk( 'Upper', 'Conjugate transpose', n, n, -one, q, ldq,
364  $ one, work, ldq )
365 *
366 * Compute norm( I - Q'*Q ) / ( N * ULP ) .
367 *
368  resid = clanhe( '1', 'Upper', n, work, ldq, rwork )
369  result( 5 ) = ( resid / real( max( 1, n ) ) ) / ulp
370 *
371 * Check sorting
372 *
373  CALL scopy( n, alpha, 1, rwork, 1 )
374  DO 110 i = k + 1, min( k+l, m )
375  j = iwork( i )
376  IF( i.NE.j ) THEN
377  temp = rwork( i )
378  rwork( i ) = rwork( j )
379  rwork( 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( rwork( i ).LT.rwork( i+1 ) )
386  $ result( 6 ) = ulpinv
387  120 CONTINUE
388 *
389  RETURN
390 *
391 * End of CGSVTS3
392 *
subroutine cgemm(TRANSA, TRANSB, M, N, K, ALPHA, A, LDA, B, LDB, BETA, C, LDC)
CGEMM
Definition: cgemm.f:187
subroutine cherk(UPLO, TRANS, N, K, ALPHA, A, LDA, BETA, C, LDC)
CHERK
Definition: cherk.f:173
real function clange(NORM, M, N, A, LDA, WORK)
CLANGE returns the value of the 1-norm, Frobenius norm, infinity-norm, or the largest absolute value ...
Definition: clange.f:115
subroutine cggsvd3(JOBU, JOBV, JOBQ, M, N, P, K, L, A, LDA, B, LDB, ALPHA, BETA, U, LDU, V, LDV, Q, LDQ, WORK, LWORK, RWORK, IWORK, INFO)
CGGSVD3 computes the singular value decomposition (SVD) for OTHER matrices
Definition: cggsvd3.f:354
real function clanhe(NORM, UPLO, N, A, LDA, WORK)
CLANHE returns the value of the 1-norm, or the Frobenius norm, or the infinity norm,...
Definition: clanhe.f:124
subroutine claset(UPLO, M, N, ALPHA, BETA, A, LDA)
CLASET initializes the off-diagonal elements and the diagonal elements of a matrix to given values.
Definition: claset.f:106
subroutine clacpy(UPLO, M, N, A, LDA, B, LDB)
CLACPY copies all or part of one two-dimensional array to another.
Definition: clacpy.f:103
subroutine scopy(N, SX, INCX, SY, INCY)
SCOPY
Definition: scopy.f:82
real function slamch(CMACH)
SLAMCH
Definition: slamch.f:68
Here is the call graph for this function:
Here is the caller graph for this function: