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

## ◆ sgsvts3()

 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.```

Definition at line 207 of file sgsvts3.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 REAL 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 REAL ZERO, ONE
231 parameter( zero = 0.0e+0, one = 1.0e+0 )
232* ..
233* .. Local Scalars ..
234 INTEGER I, INFO, J, K, L
235 REAL ANORM, BNORM, RESID, TEMP, ULP, ULPINV, UNFL
236* ..
237* .. External Functions ..
238 REAL SLAMCH, SLANGE, SLANSY
239 EXTERNAL slamch, slange, slansy
240* ..
241* .. External Subroutines ..
242 EXTERNAL scopy, sgemm, sggsvd3, slacpy, slaset, ssyrk
243* ..
244* .. Intrinsic Functions ..
245 INTRINSIC max, min, real
246* ..
247* .. Executable Statements ..
248*
249 ulp = slamch( 'Precision' )
250 ulpinv = one / ulp
251 unfl = slamch( 'Safe minimum' )
252*
253* Copy the matrix A to the array AF.
254*
255 CALL slacpy( 'Full', m, n, a, lda, af, lda )
256 CALL slacpy( 'Full', p, n, b, ldb, bf, ldb )
257*
258 anorm = max( slange( '1', m, n, a, lda, rwork ), unfl )
259 bnorm = max( slange( '1', p, n, b, ldb, rwork ), unfl )
260*
261* Factorize the matrices A and B in the arrays AF and BF.
262*
263 CALL sggsvd3( '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 sgemm( 'No transpose', 'No transpose', m, n, n, one, a, lda,
286 \$ q, ldq, zero, work, lda )
287*
288 CALL sgemm( '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 = slange( '1', m, n, a, lda, rwork )
306*
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 sgemm( 'No transpose', 'No transpose', p, n, n, one, b, ldb,
317 \$ q, ldq, zero, work, ldb )
318*
319 CALL sgemm( '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 = slange( '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 slaset( 'Full', m, m, zero, one, work, ldq )
341 CALL ssyrk( 'Upper', 'Transpose', m, m, -one, u, ldu, one, work,
342 \$ ldu )
343*
344* Compute norm( I - U'*U ) / ( M * ULP ) .
345*
346 resid = slansy( '1', 'Upper', m, work, ldu, rwork )
347 result( 3 ) = ( resid / real( max( 1, m ) ) ) / ulp
348*
349* Compute I - V'*V
350*
351 CALL slaset( 'Full', p, p, zero, one, work, ldv )
352 CALL ssyrk( 'Upper', 'Transpose', p, p, -one, v, ldv, one, work,
353 \$ ldv )
354*
355* Compute norm( I - V'*V ) / ( P * ULP ) .
356*
357 resid = slansy( '1', 'Upper', p, work, ldv, rwork )
358 result( 4 ) = ( resid / real( max( 1, p ) ) ) / ulp
359*
360* Compute I - Q'*Q
361*
362 CALL slaset( 'Full', n, n, zero, one, work, ldq )
363 CALL ssyrk( 'Upper', 'Transpose', n, n, -one, q, ldq, one, work,
364 \$ ldq )
365*
366* Compute norm( I - Q'*Q ) / ( N * ULP ) .
367*
368 resid = slansy( '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, 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 SGSVTS3
392*
subroutine scopy(n, sx, incx, sy, incy)
SCOPY
Definition scopy.f:82
subroutine sgemm(transa, transb, m, n, k, alpha, a, lda, b, ldb, beta, c, ldc)
SGEMM
Definition sgemm.f:188
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:349
subroutine ssyrk(uplo, trans, n, k, alpha, a, lda, beta, c, ldc)
SSYRK
Definition ssyrk.f:169
subroutine slacpy(uplo, m, n, a, lda, b, ldb)
SLACPY copies all or part of one two-dimensional array to another.
Definition slacpy.f:103
real function slamch(cmach)
SLAMCH
Definition slamch.f:68
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:114
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,...
Definition slansy.f:122
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:110
Here is the call graph for this function:
Here is the caller graph for this function: