Update LAPACK to 3.9.0

This commit is contained in:
Martin Kroeker 2019-12-30 16:41:31 +01:00 committed by GitHub
parent ab74361a0c
commit abbc7941ff
No known key found for this signature in database
GPG Key ID: 4AEE18F83AFDEB23
48 changed files with 1038 additions and 455 deletions

View File

@ -1,5 +1,3 @@
include ../../make.inc
########################################################################
# This is the makefile for the eigenvalue test program from LAPACK.
# The test files are organized as follows:
@ -33,6 +31,9 @@ include ../../make.inc
#
########################################################################
TOPSRCDIR = ../..
include $(TOPSRCDIR)/make.inc
AEIGTST = \
alahdg.o \
alasum.o \
@ -117,24 +118,26 @@ ZEIGTST = zchkee.o \
zsgt01.o zslect.o \
zstt21.o zstt22.o zunt01.o zunt03.o
.PHONY: all
all: single complex double complex16
.PHONY: single complex double complex16
single: xeigtsts
complex: xeigtstc
double: xeigtstd
complex16: xeigtstz
xeigtsts: $(SEIGTST) $(SCIGTST) $(AEIGTST) ../../$(TMGLIB) ../../$(LAPACKLIB) $(BLASLIB)
$(LOADER) $(LOADOPTS) -o $@ $^
xeigtsts: $(SEIGTST) $(SCIGTST) $(AEIGTST) $(TMGLIB) $(LAPACKLIB) $(BLASLIB)
$(FC) $(FFLAGS) $(LDFLAGS) -o $@ $^
xeigtstc: $(CEIGTST) $(SCIGTST) $(AEIGTST) ../../$(TMGLIB) ../../$(LAPACKLIB) $(BLASLIB)
$(LOADER) $(LOADOPTS) -o $@ $^
xeigtstc: $(CEIGTST) $(SCIGTST) $(AEIGTST) $(TMGLIB) $(LAPACKLIB) $(BLASLIB)
$(FC) $(FFLAGS) $(LDFLAGS) -o $@ $^
xeigtstd: $(DEIGTST) $(DZIGTST) $(AEIGTST) ../../$(TMGLIB) ../../$(LAPACKLIB) $(BLASLIB)
$(LOADER) $(LOADOPTS) -o $@ $^
xeigtstd: $(DEIGTST) $(DZIGTST) $(AEIGTST) $(TMGLIB) $(LAPACKLIB) $(BLASLIB)
$(FC) $(FFLAGS) $(LDFLAGS) -o $@ $^
xeigtstz: $(ZEIGTST) $(DZIGTST) $(AEIGTST) ../../$(TMGLIB) ../../$(LAPACKLIB) $(BLASLIB)
$(LOADER) $(LOADOPTS) -o $@ $^
xeigtstz: $(ZEIGTST) $(DZIGTST) $(AEIGTST) $(TMGLIB) $(LAPACKLIB) $(BLASLIB)
$(FC) $(FFLAGS) $(LDFLAGS) -o $@ $^
$(AEIGTST): $(FRC)
$(SCIGTST): $(FRC)
@ -147,6 +150,7 @@ $(ZEIGTST): $(FRC)
FRC:
@FRC=$(FRC)
.PHONY: clean cleanobj cleanexe
clean: cleanobj cleanexe
cleanobj:
rm -f *.o
@ -154,13 +158,10 @@ cleanexe:
rm -f xeigtst*
schkee.o: schkee.f
$(FORTRAN) $(DRVOPTS) -c -o $@ $<
$(FC) $(FFLAGS_DRV) -c -o $@ $<
dchkee.o: dchkee.f
$(FORTRAN) $(DRVOPTS) -c -o $@ $<
$(FC) $(FFLAGS_DRV) -c -o $@ $<
cchkee.o: cchkee.f
$(FORTRAN) $(DRVOPTS) -c -o $@ $<
$(FC) $(FFLAGS_DRV) -c -o $@ $<
zchkee.o: zchkee.f
$(FORTRAN) $(DRVOPTS) -c -o $@ $<
.f.o:
$(FORTRAN) $(OPTS) -c -o $@ $<
$(FC) $(FFLAGS_DRV) -c -o $@ $<

View File

@ -52,6 +52,7 @@
*> \verbatim
*> A is COMPLEX array, dimension (LDA,N)
*> The m by n matrix A.
*> \endverbatim
*>
*> \param[in] LDA
*> \verbatim

View File

@ -167,7 +167,7 @@
*> CSTEMR('V', 'I')
*>
*> Tests 29 through 34 are disable at present because CSTEMR
*> does not handle partial specturm requests.
*> does not handle partial spectrum requests.
*>
*> (29) | S - Z D Z* | / ( |S| n ulp ) CSTEMR('V', 'I')
*>

View File

@ -188,7 +188,7 @@
*> CSTEMR('V', 'I')
*>
*> Tests 29 through 34 are disable at present because CSTEMR
*> does not handle partial specturm requests.
*> does not handle partial spectrum requests.
*>
*> (29) | S - Z D Z* | / ( |S| n ulp ) CSTEMR('V', 'I')
*>

View File

@ -737,7 +737,7 @@
CALL CLACPY( 'Full', MPLUSN, MPLUSN, AI, LDA, A, LDA )
CALL CLACPY( 'Full', MPLUSN, MPLUSN, BI, LDA, B, LDA )
*
* Compute the Schur factorization while swaping the
* Compute the Schur factorization while swapping the
* m-by-m (1,1)-blocks with n-by-n (2,2)-blocks.
*
CALL CGGESX( 'V', 'V', 'S', CLCTSX, 'B', MPLUSN, AI, LDA, BI, LDA,

View File

@ -33,8 +33,9 @@
*>
*> \verbatim
*>
*> CDRVBD checks the singular value decomposition (SVD) driver CGESVD
*> and CGESDD.
*> CDRVBD checks the singular value decomposition (SVD) driver CGESVD,
*> CGESDD, CGESVJ, CGEJSV, CGESVDX, and CGESVDQ.
*>
*> CGESVD and CGESDD factors A = U diag(S) VT, where U and VT are
*> unitary and diag(S) is diagonal with the entries of the array S on
*> its diagonal. The entries of S are the singular values, nonnegative
@ -73,81 +74,92 @@
*>
*> Test for CGESDD:
*>
*> (1) | A - U diag(S) VT | / ( |A| max(M,N) ulp )
*>
*> (2) | I - U'U | / ( M ulp )
*>
*> (3) | I - VT VT' | / ( N ulp )
*>
*> (4) S contains MNMIN nonnegative values in decreasing order.
*> (Return 0 if true, 1/ULP if false.)
*>
*> (5) | U - Upartial | / ( M ulp ) where Upartial is a partially
*> computed U.
*>
*> (6) | VT - VTpartial | / ( N ulp ) where VTpartial is a partially
*> computed VT.
*>
*> (7) | S - Spartial | / ( MNMIN ulp |S| ) where Spartial is the
*> vector of singular values from the partial SVD
*>
*> Test for CGESVJ:
*>
*> (1) | A - U diag(S) VT | / ( |A| max(M,N) ulp )
*>
*> (2) | I - U'U | / ( M ulp )
*>
*> (3) | I - VT VT' | / ( N ulp )
*>
*> (4) S contains MNMIN nonnegative values in decreasing order.
*> (Return 0 if true, 1/ULP if false.)
*>
*> Test for CGEJSV:
*>
*> (1) | A - U diag(S) VT | / ( |A| max(M,N) ulp )
*>
*> (2) | I - U'U | / ( M ulp )
*>
*> (3) | I - VT VT' | / ( N ulp )
*>
*> (4) S contains MNMIN nonnegative values in decreasing order.
*> (Return 0 if true, 1/ULP if false.)
*>
*> Test for CGESVDX( 'V', 'V', 'A' )/CGESVDX( 'N', 'N', 'A' )
*>
*> (1) | A - U diag(S) VT | / ( |A| max(M,N) ulp )
*>
*> (2) | I - U'U | / ( M ulp )
*>
*> (3) | I - VT VT' | / ( N ulp )
*>
*> (4) S contains MNMIN nonnegative values in decreasing order.
*> (Return 0 if true, 1/ULP if false.)
*>
*> (5) | U - Upartial | / ( M ulp ) where Upartial is a partially
*> computed U.
*>
*> (6) | VT - VTpartial | / ( N ulp ) where VTpartial is a partially
*> computed VT.
*>
*> (7) | S - Spartial | / ( MNMIN ulp |S| ) where Spartial is the
*> vector of singular values from the partial SVD
*>
*> Test for CGESVDX( 'V', 'V', 'I' )
*>
*> (8) | U' A VT''' - diag(S) | / ( |A| max(M,N) ulp )
*> (8) | A - U diag(S) VT | / ( |A| max(M,N) ulp )
*>
*> (9) | I - U'U | / ( M ulp )
*>
*> (10) | I - VT VT' | / ( N ulp )
*>
*> (11) S contains MNMIN nonnegative values in decreasing order.
*> (Return 0 if true, 1/ULP if false.)
*>
*> (12) | U - Upartial | / ( M ulp ) where Upartial is a partially
*> computed U.
*>
*> (13) | VT - VTpartial | / ( N ulp ) where VTpartial is a partially
*> computed VT.
*>
*> (14) | S - Spartial | / ( MNMIN ulp |S| ) where Spartial is the
*> vector of singular values from the partial SVD
*>
*> Test for CGESVDQ:
*>
*> (36) | A - U diag(S) VT | / ( |A| max(M,N) ulp )
*>
*> (37) | I - U'U | / ( M ulp )
*>
*> (38) | I - VT VT' | / ( N ulp )
*>
*> (39) S contains MNMIN nonnegative values in decreasing order.
*> (Return 0 if true, 1/ULP if false.)
*>
*> Test for CGESVJ:
*>
*> (15) | A - U diag(S) VT | / ( |A| max(M,N) ulp )
*>
*> (16) | I - U'U | / ( M ulp )
*>
*> (17) | I - VT VT' | / ( N ulp )
*>
*> (18) S contains MNMIN nonnegative values in decreasing order.
*> (Return 0 if true, 1/ULP if false.)
*>
*> Test for CGEJSV:
*>
*> (19) | A - U diag(S) VT | / ( |A| max(M,N) ulp )
*>
*> (20) | I - U'U | / ( M ulp )
*>
*> (21) | I - VT VT' | / ( N ulp )
*>
*> (22) S contains MNMIN nonnegative values in decreasing order.
*> (Return 0 if true, 1/ULP if false.)
*>
*> Test for CGESVDX( 'V', 'V', 'A' )/CGESVDX( 'N', 'N', 'A' )
*>
*> (23) | A - U diag(S) VT | / ( |A| max(M,N) ulp )
*>
*> (24) | I - U'U | / ( M ulp )
*>
*> (25) | I - VT VT' | / ( N ulp )
*>
*> (26) S contains MNMIN nonnegative values in decreasing order.
*> (Return 0 if true, 1/ULP if false.)
*>
*> (27) | U - Upartial | / ( M ulp ) where Upartial is a partially
*> computed U.
*>
*> (28) | VT - VTpartial | / ( N ulp ) where VTpartial is a partially
*> computed VT.
*>
*> (29) | S - Spartial | / ( MNMIN ulp |S| ) where Spartial is the
*> vector of singular values from the partial SVD
*>
*> Test for CGESVDX( 'V', 'V', 'I' )
*>
*> (30) | U' A VT''' - diag(S) | / ( |A| max(M,N) ulp )
*>
*> (31) | I - U'U | / ( M ulp )
*>
*> (32) | I - VT VT' | / ( N ulp )
*>
*> Test for CGESVDX( 'V', 'V', 'V' )
*>
*> (11) | U' A VT''' - diag(S) | / ( |A| max(M,N) ulp )
*> (33) | U' A VT''' - diag(S) | / ( |A| max(M,N) ulp )
*>
*> (12) | I - U'U | / ( M ulp )
*> (34) | I - U'U | / ( M ulp )
*>
*> (13) | I - VT VT' | / ( N ulp )
*> (35) | I - VT VT' | / ( N ulp )
*>
*> The "sizes" are specified by the arrays MM(1:NSIZES) and
*> NN(1:NSIZES); the value of each element pair (MM(j),NN(j))
@ -393,6 +405,8 @@
* -- LAPACK is a software package provided by Univ. of Tennessee, --
* -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
* June 2016
*
IMPLICIT NONE
*
* .. Scalar Arguments ..
INTEGER INFO, LDA, LDU, LDVT, LWORK, NOUNIT, NSIZES,
@ -411,7 +425,7 @@
* =====================================================================
*
* .. Parameters ..
REAL ZERO, ONE, TWO, HALF
REAL ZERO, ONE, TWO, HALF
PARAMETER ( ZERO = 0.0E0, ONE = 1.0E0, TWO = 2.0E0,
$ HALF = 0.5E0 )
COMPLEX CZERO, CONE
@ -431,10 +445,13 @@
REAL ANORM, DIF, DIV, OVFL, RTUNFL, ULP, ULPINV,
$ UNFL, VL, VU
* ..
* .. Local Scalars for CGESVDQ ..
INTEGER LIWORK, NUMRANK
* ..
* .. Local Arrays ..
CHARACTER CJOB( 4 ), CJOBR( 3 ), CJOBV( 2 )
INTEGER IOLDSD( 4 ), ISEED2( 4 )
REAL RESULT( 35 )
REAL RESULT( 39 )
* ..
* .. External Functions ..
REAL SLAMCH, SLARND
@ -442,8 +459,8 @@
* ..
* .. External Subroutines ..
EXTERNAL ALASVM, XERBLA, CBDT01, CBDT05, CGESDD,
$ CGESVD, CGESVJ, CGEJSV, CGESVDX, CLACPY,
$ CLASET, CLATMS, CUNT01, CUNT03
$ CGESVD, CGESVDQ, CGESVJ, CGEJSV, CGESVDX,
$ CLACPY, CLASET, CLATMS, CUNT01, CUNT03
* ..
* .. Intrinsic Functions ..
INTRINSIC ABS, REAL, MAX, MIN
@ -838,8 +855,64 @@
130 CONTINUE
*
* Test CGESVJ: Factorize A
* Note: CGESVJ does not work for M < N
* Test CGESVDQ
* Note: CGESVDQ only works for M >= N
*
RESULT( 36 ) = ZERO
RESULT( 37 ) = ZERO
RESULT( 38 ) = ZERO
RESULT( 39 ) = ZERO
*
IF( M.GE.N ) THEN
IWTMP = 2*MNMIN*MNMIN + 2*MNMIN + MAX( M, N )
LSWORK = IWTMP + ( IWSPC-1 )*( LWORK-IWTMP ) / 3
LSWORK = MIN( LSWORK, LWORK )
LSWORK = MAX( LSWORK, 1 )
IF( IWSPC.EQ.4 )
$ LSWORK = LWORK
*
CALL CLACPY( 'F', M, N, ASAV, LDA, A, LDA )
SRNAMT = 'CGESVDQ'
*
LRWORK = MAX(2, M, 5*N)
LIWORK = MAX( N, 1 )
CALL CGESVDQ( 'H', 'N', 'N', 'A', 'A',
$ M, N, A, LDA, SSAV, USAV, LDU,
$ VTSAV, LDVT, NUMRANK, IWORK, LIWORK,
$ WORK, LWORK, RWORK, LRWORK, IINFO )
*
IF( IINFO.NE.0 ) THEN
WRITE( NOUNIT, FMT = 9995 )'CGESVDQ', IINFO, M, N,
$ JTYPE, LSWORK, IOLDSD
INFO = ABS( IINFO )
RETURN
END IF
*
* Do tests 36--39
*
CALL CBDT01( M, N, 0, ASAV, LDA, USAV, LDU, SSAV, E,
$ VTSAV, LDVT, WORK, RWORK, RESULT( 36 ) )
IF( M.NE.0 .AND. N.NE.0 ) THEN
CALL CUNT01( 'Columns', M, M, USAV, LDU, WORK,
$ LWORK, RWORK, RESULT( 37 ) )
CALL CUNT01( 'Rows', N, N, VTSAV, LDVT, WORK,
$ LWORK, RWORK, RESULT( 38 ) )
END IF
RESULT( 39 ) = ZERO
DO 199 I = 1, MNMIN - 1
IF( SSAV( I ).LT.SSAV( I+1 ) )
$ RESULT( 39 ) = ULPINV
IF( SSAV( I ).LT.ZERO )
$ RESULT( 39 ) = ULPINV
199 CONTINUE
IF( MNMIN.GE.1 ) THEN
IF( SSAV( MNMIN ).LT.ZERO )
$ RESULT( 39 ) = ULPINV
END IF
END IF
*
* Test CGESVJ
* Note: CGESVJ only works for M >= N
*
RESULT( 15 ) = ZERO
RESULT( 16 ) = ZERO
@ -847,13 +920,13 @@
RESULT( 18 ) = ZERO
*
IF( M.GE.N ) THEN
IWTMP = 2*MNMIN*MNMIN + 2*MNMIN + MAX( M, N )
LSWORK = IWTMP + ( IWSPC-1 )*( LWORK-IWTMP ) / 3
LSWORK = MIN( LSWORK, LWORK )
LSWORK = MAX( LSWORK, 1 )
LRWORK = MAX(6,N)
IF( IWSPC.EQ.4 )
$ LSWORK = LWORK
IWTMP = 2*MNMIN*MNMIN + 2*MNMIN + MAX( M, N )
LSWORK = IWTMP + ( IWSPC-1 )*( LWORK-IWTMP ) / 3
LSWORK = MIN( LSWORK, LWORK )
LSWORK = MAX( LSWORK, 1 )
LRWORK = MAX(6,N)
IF( IWSPC.EQ.4 )
$ LSWORK = LWORK
*
CALL CLACPY( 'F', M, N, ASAV, LDA, USAV, LDA )
SRNAMT = 'CGESVJ'
@ -861,8 +934,7 @@
& 0, A, LDVT, WORK, LWORK, RWORK,
& LRWORK, IINFO )
*
* CGESVJ retuns V not VT, so we transpose to use the same
* test suite.
* CGESVJ returns V not VH
*
DO J=1,N
DO I=1,N
@ -900,31 +972,30 @@
END IF
END IF
*
* Test CGEJSV: Factorize A
* Note: CGEJSV does not work for M < N
* Test CGEJSV
* Note: CGEJSV only works for M >= N
*
RESULT( 19 ) = ZERO
RESULT( 20 ) = ZERO
RESULT( 21 ) = ZERO
RESULT( 22 ) = ZERO
IF( M.GE.N ) THEN
IWTMP = 2*MNMIN*MNMIN + 2*MNMIN + MAX( M, N )
LSWORK = IWTMP + ( IWSPC-1 )*( LWORK-IWTMP ) / 3
LSWORK = MIN( LSWORK, LWORK )
LSWORK = MAX( LSWORK, 1 )
IF( IWSPC.EQ.4 )
$ LSWORK = LWORK
LRWORK = MAX( 7, N + 2*M)
IWTMP = 2*MNMIN*MNMIN + 2*MNMIN + MAX( M, N )
LSWORK = IWTMP + ( IWSPC-1 )*( LWORK-IWTMP ) / 3
LSWORK = MIN( LSWORK, LWORK )
LSWORK = MAX( LSWORK, 1 )
IF( IWSPC.EQ.4 )
$ LSWORK = LWORK
LRWORK = MAX( 7, N + 2*M)
*
CALL CLACPY( 'F', M, N, ASAV, LDA, VTSAV, LDA )
CALL CLACPY( 'F', M, N, ASAV, LDA, VTSAV, LDA )
SRNAMT = 'CGEJSV'
CALL CGEJSV( 'G', 'U', 'V', 'R', 'N', 'N',
& M, N, VTSAV, LDA, SSAV, USAV, LDU, A, LDVT,
& WORK, LWORK, RWORK,
& LRWORK, IWORK, IINFO )
*
* CGEJSV retuns V not VT, so we transpose to use the same
* test suite.
* CGEJSV returns V not VH
*
DO 133 J=1,N
DO 132 I=1,N
@ -933,7 +1004,7 @@
133 END DO
*
IF( IINFO.NE.0 ) THEN
WRITE( NOUNIT, FMT = 9995 )'GESVJ', IINFO, M, N,
WRITE( NOUNIT, FMT = 9995 )'GEJSV', IINFO, M, N,
$ JTYPE, LSWORK, IOLDSD
INFO = ABS( IINFO )
RETURN
@ -1160,7 +1231,7 @@
*
NTEST = 0
NFAIL = 0
DO 190 J = 1, 35
DO 190 J = 1, 39
IF( RESULT( J ).GE.ZERO )
$ NTEST = NTEST + 1
IF( RESULT( J ).GE.THRESH )
@ -1175,7 +1246,7 @@
NTESTF = 2
END IF
*
DO 200 J = 1, 35
DO 200 J = 1, 39
IF( RESULT( J ).GE.THRESH ) THEN
WRITE( NOUNIT, FMT = 9997 )M, N, JTYPE, IWSPC,
$ IOLDSD, J, RESULT( J )
@ -1251,6 +1322,12 @@
$ / '33 = | U**T A VT**T - diag(S) | / ( |A| max(M,N) ulp )',
$ / '34 = | I - U**T U | / ( M ulp ) ',
$ / '35 = | I - VT VT**T | / ( N ulp ) ',
$ ' CGESVDQ(H,N,N,A,A',
$ / '36 = | A - U diag(S) VT | / ( |A| max(M,N) ulp ) ',
$ / '37 = | I - U**T U | / ( M ulp ) ',
$ / '38 = | I - VT VT**T | / ( N ulp ) ',
$ / '39 = 0 if S contains min(M,N) nonnegative values in',
$ ' decreasing order, else 1/ulp',
$ / / )
9997 FORMAT( ' M=', I5, ', N=', I5, ', type ', I1, ', IWS=', I1,
$ ', seed=', 4( I4, ',' ), ' test(', I2, ')=', G11.4 )

View File

@ -36,6 +36,8 @@
*> CGEJSV compute SVD of an M-by-N matrix A where M >= N
*> CGESVDX compute SVD of an M-by-N matrix A(by bisection
*> and inverse iteration)
*> CGESVDQ compute SVD of an M-by-N matrix A(with a
*> QR-Preconditioned )
*> \endverbatim
*
* Arguments:
@ -101,7 +103,7 @@
* ..
* .. External Subroutines ..
EXTERNAL CHKXER, CGEES, CGEESX, CGEEV, CGEEVX, CGEJSV,
$ CGESDD, CGESVD
$ CGESDD, CGESVD, CGESVDX, CGESVDQ
* ..
* .. External Functions ..
LOGICAL LSAMEN, CSLECT
@ -495,6 +497,61 @@
ELSE
WRITE( NOUT, FMT = 9998 )
END IF
*
* Test CGESVDQ
*
SRNAMT = 'CGESVDQ'
INFOT = 1
CALL CGESVDQ( 'X', 'P', 'T', 'A', 'A', 0, 0, A, 1, S, U,
$ 0, VT, 0, NS, IW, 1, W, 1, RW, 1, INFO )
CALL CHKXER( 'CGESVDQ', INFOT, NOUT, LERR, OK )
INFOT = 2
CALL CGESVDQ( 'A', 'X', 'T', 'A', 'A', 0, 0, A, 1, S, U,
$ 0, VT, 0, NS, IW, 1, W, 1, RW, 1, INFO )
CALL CHKXER( 'CGESVDQ', INFOT, NOUT, LERR, OK )
INFOT = 3
CALL CGESVDQ( 'A', 'P', 'X', 'A', 'A', 0, 0, A, 1, S, U,
$ 0, VT, 0, NS, IW, 1, W, 1, RW, 1, INFO )
CALL CHKXER( 'CGESVDQ', INFOT, NOUT, LERR, OK )
INFOT = 4
CALL CGESVDQ( 'A', 'P', 'T', 'X', 'A', 0, 0, A, 1, S, U,
$ 0, VT, 0, NS, IW, 1, W, 1, RW, 1, INFO )
CALL CHKXER( 'CGESVDQ', INFOT, NOUT, LERR, OK )
INFOT = 5
CALL CGESVDQ( 'A', 'P', 'T', 'A', 'X', 0, 0, A, 1, S, U,
$ 0, VT, 0, NS, IW, 1, W, 1, RW, 1, INFO )
CALL CHKXER( 'CGESVDQ', INFOT, NOUT, LERR, OK )
INFOT = 6
CALL CGESVDQ( 'A', 'P', 'T', 'A', 'A', -1, 0, A, 1, S, U,
$ 0, VT, 0, NS, IW, 1, W, 1, RW, 1, INFO )
CALL CHKXER( 'CGESVDQ', INFOT, NOUT, LERR, OK )
INFOT = 7
CALL CGESVDQ( 'A', 'P', 'T', 'A', 'A', 0, 1, A, 1, S, U,
$ 0, VT, 0, NS, IW, 1, W, 1, RW, 1, INFO )
CALL CHKXER( 'CGESVDQ', INFOT, NOUT, LERR, OK )
INFOT = 9
CALL CGESVDQ( 'A', 'P', 'T', 'A', 'A', 1, 1, A, 0, S, U,
$ 0, VT, 0, NS, IW, 1, W, 1, RW, 1, INFO )
CALL CHKXER( 'CGESVDQ', INFOT, NOUT, LERR, OK )
INFOT = 12
CALL CGESVDQ( 'A', 'P', 'T', 'A', 'A', 1, 1, A, 1, S, U,
$ -1, VT, 0, NS, IW, 1, W, 1, RW, 1, INFO )
CALL CHKXER( 'CGESVDQ', INFOT, NOUT, LERR, OK )
INFOT = 14
CALL CGESVDQ( 'A', 'P', 'T', 'A', 'A', 1, 1, A, 1, S, U,
$ 1, VT, -1, NS, IW, 1, W, 1, RW, 1, INFO )
CALL CHKXER( 'CGESVDQ', INFOT, NOUT, LERR, OK )
INFOT = 17
CALL CGESVDQ( 'A', 'P', 'T', 'A', 'A', 1, 1, A, 1, S, U,
$ 1, VT, 1, NS, IW, -5, W, 1, RW, 1, INFO )
CALL CHKXER( 'CGESVDQ', INFOT, NOUT, LERR, OK )
NT = 11
IF( OK ) THEN
WRITE( NOUT, FMT = 9999 )SRNAMT( 1:LEN_TRIM( SRNAMT ) ),
$ NT
ELSE
WRITE( NOUT, FMT = 9998 )
END IF
END IF
*
* Print a summary line.

View File

@ -29,12 +29,13 @@
*>
*> CGET51 generally checks a decomposition of the form
*>
*> A = U B VC>
*> where * means conjugate transpose and U and V are unitary.
*> A = U B V**H
*>
*> where **H means conjugate transpose and U and V are unitary.
*>
*> Specifically, if ITYPE=1
*>
*> RESULT = | A - U B V* | / ( |A| n ulp )
*> RESULT = | A - U B V**H | / ( |A| n ulp )
*>
*> If ITYPE=2, then:
*>
@ -42,7 +43,7 @@
*>
*> If ITYPE=3, then:
*>
*> RESULT = | I - UU* | / ( n ulp )
*> RESULT = | I - U U**H | / ( n ulp )
*> \endverbatim
*
* Arguments:
@ -52,9 +53,9 @@
*> \verbatim
*> ITYPE is INTEGER
*> Specifies the type of tests to be performed.
*> =1: RESULT = | A - U B V* | / ( |A| n ulp )
*> =1: RESULT = | A - U B V**H | / ( |A| n ulp )
*> =2: RESULT = | A - B | / ( |A| n ulp )
*> =3: RESULT = | I - UU* | / ( n ulp )
*> =3: RESULT = | I - U U**H | / ( n ulp )
*> \endverbatim
*>
*> \param[in] N
@ -218,7 +219,7 @@
*
IF( ITYPE.EQ.1 ) THEN
*
* ITYPE=1: Compute W = A - UBV'
* ITYPE=1: Compute W = A - U B V**H
*
CALL CLACPY( ' ', N, N, A, LDA, WORK, N )
CALL CGEMM( 'N', 'N', N, N, N, CONE, U, LDU, B, LDB, CZERO,
@ -259,7 +260,7 @@
*
* Tests not scaled by norm(A)
*
* ITYPE=3: Compute UU' - I
* ITYPE=3: Compute U U**H - I
*
CALL CGEMM( 'N', 'C', N, N, N, CONE, U, LDU, U, LDU, CZERO,
$ WORK, N )

View File

@ -28,14 +28,16 @@
*>
*> CHBT21 generally checks a decomposition of the form
*>
*> A = U S UC>
*> where * means conjugate transpose, A is hermitian banded, U is
*> A = U S U**H
*>
*> where **H means conjugate transpose, A is hermitian banded, U is
*> unitary, and S is diagonal (if KS=0) or symmetric
*> tridiagonal (if KS=1).
*>
*> Specifically:
*>
*> RESULT(1) = | A - U S U* | / ( |A| n ulp ) *andC> RESULT(2) = | I - UU* | / ( n ulp )
*> RESULT(1) = | A - U S U**H | / ( |A| n ulp ) and
*> RESULT(2) = | I - U U**H | / ( n ulp )
*> \endverbatim
*
* Arguments:
@ -220,7 +222,7 @@
*
ANORM = MAX( CLANHB( '1', CUPLO, N, IKA, A, LDA, RWORK ), UNFL )
*
* Compute error matrix: Error = A - U S U*
* Compute error matrix: Error = A - U S U**H
*
* Copy A from SB to SP storage format.
*
@ -271,7 +273,7 @@
*
* Do Test 2
*
* Compute UU* - I
* Compute U U**H - I
*
CALL CGEMM( 'N', 'C', N, N, N, CONE, U, LDU, U, LDU, CZERO, WORK,
$ N )

View File

@ -29,8 +29,9 @@
*>
*> CHET21 generally checks a decomposition of the form
*>
*> A = U S UC>
*> where * means conjugate transpose, A is hermitian, U is unitary, and
*> A = U S U**H
*>
*> where **H means conjugate transpose, A is hermitian, U is unitary, and
*> S is diagonal (if KBAND=0) or (real) symmetric tridiagonal (if
*> KBAND=1).
*>
@ -42,18 +43,19 @@
*>
*> Specifically, if ITYPE=1, then:
*>
*> RESULT(1) = | A - U S U* | / ( |A| n ulp ) *andC> RESULT(2) = | I - UU* | / ( n ulp )
*> RESULT(1) = | A - U S U**H | / ( |A| n ulp ) and
*> RESULT(2) = | I - U U**H | / ( n ulp )
*>
*> If ITYPE=2, then:
*>
*> RESULT(1) = | A - V S V* | / ( |A| n ulp )
*> RESULT(1) = | A - V S V**H | / ( |A| n ulp )
*>
*> If ITYPE=3, then:
*>
*> RESULT(1) = | I - UV* | / ( n ulp )
*> RESULT(1) = | I - U V**H | / ( n ulp )
*>
*> For ITYPE > 1, the transformation U is expressed as a product
*> V = H(1)...H(n-2), where H(j) = I - tau(j) v(j) v(j)C> and each
*> V = H(1)...H(n-2), where H(j) = I - tau(j) v(j) v(j)**H and each
*> vector v(j) has its first j elements 0 and the remaining n-j elements
*> stored in V(j+1:n,j).
*> \endverbatim
@ -66,14 +68,15 @@
*> ITYPE is INTEGER
*> Specifies the type of tests to be performed.
*> 1: U expressed as a dense unitary matrix:
*> RESULT(1) = | A - U S U* | / ( |A| n ulp ) *andC> RESULT(2) = | I - UU* | / ( n ulp )
*> RESULT(1) = | A - U S U**H | / ( |A| n ulp ) and
*> RESULT(2) = | I - U U**H | / ( n ulp )
*>
*> 2: U expressed as a product V of Housholder transformations:
*> RESULT(1) = | A - V S V* | / ( |A| n ulp )
*> RESULT(1) = | A - V S V**H | / ( |A| n ulp )
*>
*> 3: U expressed both as a dense unitary matrix and
*> as a product of Housholder transformations:
*> RESULT(1) = | I - UV* | / ( n ulp )
*> RESULT(1) = | I - U V**H | / ( n ulp )
*> \endverbatim
*>
*> \param[in] UPLO
@ -171,7 +174,7 @@
*> \verbatim
*> TAU is COMPLEX array, dimension (N)
*> If ITYPE >= 2, then TAU(j) is the scalar factor of
*> v(j) v(j)* in the Householder transformation H(j) of
*> v(j) v(j)**H in the Householder transformation H(j) of
*> the product U = H(1)...H(n-2)
*> If ITYPE < 2, then TAU is not referenced.
*> \endverbatim
@ -294,7 +297,7 @@
*
IF( ITYPE.EQ.1 ) THEN
*
* ITYPE=1: error = A - U S U*
* ITYPE=1: error = A - U S U**H
*
CALL CLASET( 'Full', N, N, CZERO, CZERO, WORK, N )
CALL CLACPY( CUPLO, N, N, A, LDA, WORK, N )
@ -304,8 +307,7 @@
10 CONTINUE
*
IF( N.GT.1 .AND. KBAND.EQ.1 ) THEN
CMK DO 20 J = 1, N - 1
DO 20 J = 2, N - 1
DO 20 J = 1, N - 1
CALL CHER2( CUPLO, N, -CMPLX( E( J ) ), U( 1, J ), 1,
$ U( 1, J-1 ), 1, WORK, N )
20 CONTINUE
@ -314,7 +316,7 @@ CMK DO 20 J = 1, N - 1
*
ELSE IF( ITYPE.EQ.2 ) THEN
*
* ITYPE=2: error = V S V* - A
* ITYPE=2: error = V S V**H - A
*
CALL CLASET( 'Full', N, N, CZERO, CZERO, WORK, N )
*
@ -371,7 +373,7 @@ CMK DO 20 J = 1, N - 1
*
ELSE IF( ITYPE.EQ.3 ) THEN
*
* ITYPE=3: error = U V* - I
* ITYPE=3: error = U V**H - I
*
IF( N.LT.2 )
$ RETURN
@ -407,7 +409,7 @@ CMK DO 20 J = 1, N - 1
*
* Do Test 2
*
* Compute UU* - I
* Compute U U**H - I
*
IF( ITYPE.EQ.1 ) THEN
CALL CGEMM( 'N', 'C', N, N, N, CONE, U, LDU, U, LDU, CZERO,

View File

@ -42,7 +42,8 @@
*>
*> Specifically, if ITYPE=1, then:
*>
*> RESULT(1) = | U' A U - S | / ( |A| m ulp ) *andC> RESULT(2) = | I - U'U | / ( m ulp )
*> RESULT(1) = | U**H A U - S | / ( |A| m ulp ) and
*> RESULT(2) = | I - U**H U | / ( m ulp )
*> \endverbatim
*
* Arguments:
@ -52,7 +53,8 @@
*> ITYPE INTEGER
*> Specifies the type of tests to be performed.
*> 1: U expressed as a dense orthogonal matrix:
*> RESULT(1) = | A - U S U' | / ( |A| n ulp ) *andC> RESULT(2) = | I - UU' | / ( n ulp )
*> RESULT(1) = | A - U S U**H | / ( |A| n ulp ) and
*> RESULT(2) = | I - U U**H | / ( n ulp )
*>
*> UPLO CHARACTER
*> If UPLO='U', the upper triangle of A will be used and the
@ -122,7 +124,7 @@
*>
*> TAU COMPLEX array, dimension (N)
*> If ITYPE >= 2, then TAU(j) is the scalar factor of
*> v(j) v(j)' in the Householder transformation H(j) of
*> v(j) v(j)**H in the Householder transformation H(j) of
*> the product U = H(1)...H(n-2)
*> If ITYPE < 2, then TAU is not referenced.
*> Not modified.
@ -215,7 +217,7 @@
*
* Compute error matrix:
*
* ITYPE=1: error = U' A U - S
* ITYPE=1: error = U**H A U - S
*
CALL CHEMM( 'L', UPLO, N, M, CONE, A, LDA, U, LDU, CZERO, WORK,
$ N )
@ -249,7 +251,7 @@
*
* Do Test 2
*
* Compute U'U - I
* Compute U**H U - I
*
IF( ITYPE.EQ.1 )
$ CALL CUNT01( 'Columns', N, M, U, LDU, WORK, 2*N*N, RWORK,

View File

@ -29,8 +29,9 @@
*>
*> CHPT21 generally checks a decomposition of the form
*>
*> A = U S UC>
*> where * means conjugate transpose, A is hermitian, U is
*> A = U S U**H
*>
*> where **H means conjugate transpose, A is hermitian, U is
*> unitary, and S is diagonal (if KBAND=0) or (real) symmetric
*> tridiagonal (if KBAND=1). If ITYPE=1, then U is represented as
*> a dense matrix, otherwise the U is expressed as a product of
@ -41,15 +42,16 @@
*>
*> Specifically, if ITYPE=1, then:
*>
*> RESULT(1) = | A - U S U* | / ( |A| n ulp ) *andC> RESULT(2) = | I - UU* | / ( n ulp )
*> RESULT(1) = | A - U S U**H | / ( |A| n ulp ) and
*> RESULT(2) = | I - U U**H | / ( n ulp )
*>
*> If ITYPE=2, then:
*>
*> RESULT(1) = | A - V S V* | / ( |A| n ulp )
*> RESULT(1) = | A - V S V**H | / ( |A| n ulp )
*>
*> If ITYPE=3, then:
*>
*> RESULT(1) = | I - UV* | / ( n ulp )
*> RESULT(1) = | I - U V**H | / ( n ulp )
*>
*> Packed storage means that, for example, if UPLO='U', then the columns
*> of the upper triangle of A are stored one after another, so that
@ -70,14 +72,16 @@
*>
*> If UPLO='U', then V = H(n-1)...H(1), where
*>
*> H(j) = I - tau(j) v(j) v(j)C>
*> H(j) = I - tau(j) v(j) v(j)**H
*>
*> and the first j-1 elements of v(j) are stored in V(1:j-1,j+1),
*> (i.e., VP( j*(j+1)/2 + 1 : j*(j+1)/2 + j-1 ) ),
*> the j-th element is 1, and the last n-j elements are 0.
*>
*> If UPLO='L', then V = H(1)...H(n-1), where
*>
*> H(j) = I - tau(j) v(j) v(j)C>
*> H(j) = I - tau(j) v(j) v(j)**H
*>
*> and the first j elements of v(j) are 0, the (j+1)-st is 1, and the
*> (j+2)-nd through n-th elements are stored in V(j+2:n,j) (i.e.,
*> in VP( (2*n-j)*(j-1)/2 + j+2 : (2*n-j)*(j-1)/2 + n ) .)
@ -91,14 +95,15 @@
*> ITYPE is INTEGER
*> Specifies the type of tests to be performed.
*> 1: U expressed as a dense unitary matrix:
*> RESULT(1) = | A - U S U* | / ( |A| n ulp ) *andC> RESULT(2) = | I - UU* | / ( n ulp )
*> RESULT(1) = | A - U S U**H | / ( |A| n ulp ) and
*> RESULT(2) = | I - U U**H | / ( n ulp )
*>
*> 2: U expressed as a product V of Housholder transformations:
*> RESULT(1) = | A - V S V* | / ( |A| n ulp )
*> RESULT(1) = | A - V S V**H | / ( |A| n ulp )
*>
*> 3: U expressed both as a dense unitary matrix and
*> as a product of Housholder transformations:
*> RESULT(1) = | I - UV* | / ( n ulp )
*> RESULT(1) = | I - U V**H | / ( n ulp )
*> \endverbatim
*>
*> \param[in] UPLO
@ -181,7 +186,7 @@
*> \verbatim
*> TAU is COMPLEX array, dimension (N)
*> If ITYPE >= 2, then TAU(j) is the scalar factor of
*> v(j) v(j)* in the Householder transformation H(j) of
*> v(j) v(j)**H in the Householder transformation H(j) of
*> the product U = H(1)...H(n-2)
*> If ITYPE < 2, then TAU is not referenced.
*> \endverbatim
@ -313,7 +318,7 @@
*
IF( ITYPE.EQ.1 ) THEN
*
* ITYPE=1: error = A - U S U*
* ITYPE=1: error = A - U S U**H
*
CALL CLASET( 'Full', N, N, CZERO, CZERO, WORK, N )
CALL CCOPY( LAP, AP, 1, WORK, 1 )
@ -323,7 +328,7 @@
10 CONTINUE
*
IF( N.GT.1 .AND. KBAND.EQ.1 ) THEN
DO 20 J = 2, N - 1
DO 20 J = 1, N - 1
CALL CHPR2( CUPLO, N, -CMPLX( E( J ) ), U( 1, J ), 1,
$ U( 1, J-1 ), 1, WORK )
20 CONTINUE
@ -332,7 +337,7 @@
*
ELSE IF( ITYPE.EQ.2 ) THEN
*
* ITYPE=2: error = V S V* - A
* ITYPE=2: error = V S V**H - A
*
CALL CLASET( 'Full', N, N, CZERO, CZERO, WORK, N )
*
@ -400,7 +405,7 @@
*
ELSE IF( ITYPE.EQ.3 ) THEN
*
* ITYPE=3: error = U V* - I
* ITYPE=3: error = U V**H - I
*
IF( N.LT.2 )
$ RETURN
@ -431,7 +436,7 @@
*
* Do Test 2
*
* Compute UU* - I
* Compute U U**H - I
*
IF( ITYPE.EQ.1 ) THEN
CALL CGEMM( 'N', 'C', N, N, N, CONE, U, LDU, U, LDU, CZERO,

View File

@ -28,14 +28,15 @@
*>
*> CSTT21 checks a decomposition of the form
*>
*> A = U S UC>
*> where * means conjugate transpose, A is real symmetric tridiagonal,
*> A = U S U**H
*>
*> where **H means conjugate transpose, A is real symmetric tridiagonal,
*> U is unitary, and S is real and diagonal (if KBAND=0) or symmetric
*> tridiagonal (if KBAND=1). Two tests are performed:
*>
*> RESULT(1) = | A - U S U* | / ( |A| n ulp )
*> RESULT(1) = | A - U S U**H | / ( |A| n ulp )
*>
*> RESULT(2) = | I - UU* | / ( n ulp )
*> RESULT(2) = | I - U U**H | / ( n ulp )
*> \endverbatim
*
* Arguments:
@ -201,7 +202,7 @@
WORK( N**2 ) = AD( N )
ANORM = MAX( ANORM, ABS( AD( N ) )+TEMP1, UNFL )
*
* Norm of A - USU*
* Norm of A - U S U**H
*
DO 20 J = 1, N
CALL CHER( 'L', N, -SD( J ), U( 1, J ), 1, WORK, N )
@ -228,7 +229,7 @@
*
* Do Test 2
*
* Compute UU* - I
* Compute U U**H - I
*
CALL CGEMM( 'N', 'C', N, N, N, CONE, U, LDU, U, LDU, CZERO, WORK,
$ N )

View File

@ -52,6 +52,7 @@
*> \verbatim
*> A is DOUBLE PRECISION array, dimension (LDA,N)
*> The m by n matrix A.
*> \endverbatim
*>
*> \param[in] LDA
*> \verbatim

View File

@ -166,7 +166,7 @@
*> DSTEMR('V', 'I')
*>
*> Tests 29 through 34 are disable at present because DSTEMR
*> does not handle partial specturm requests.
*> does not handle partial spectrum requests.
*>
*> (29) | S - Z D Z' | / ( |S| n ulp ) DSTEMR('V', 'I')
*>

View File

@ -187,7 +187,7 @@
*> DSTEMR('V', 'I')
*>
*> Tests 29 through 34 are disable at present because DSTEMR
*> does not handle partial specturm requests.
*> does not handle partial spectrum requests.
*>
*> (29) | S - Z D Z' | / ( |S| n ulp ) DSTEMR('V', 'I')
*>

View File

@ -769,7 +769,7 @@
CALL DLACPY( 'Full', MPLUSN, MPLUSN, AI, LDA, A, LDA )
CALL DLACPY( 'Full', MPLUSN, MPLUSN, BI, LDA, B, LDA )
*
* Compute the Schur factorization while swaping the
* Compute the Schur factorization while swapping the
* m-by-m (1,1)-blocks with n-by-n (2,2)-blocks.
*
CALL DGGESX( 'V', 'V', 'S', DLCTSX, 'B', MPLUSN, AI, LDA, BI, LDA,

View File

@ -32,7 +32,7 @@
*> \verbatim
*>
*> DDRVBD checks the singular value decomposition (SVD) drivers
*> DGESVD, DGESDD, DGESVJ, and DGEJSV.
*> DGESVD, DGESDD, DGESVDQ, DGESVJ, DGEJSV, and DGESVDX.
*>
*> Both DGESVD and DGESDD factor A = U diag(S) VT, where U and VT are
*> orthogonal and diag(S) is diagonal with the entries of the array S
@ -90,6 +90,17 @@
*> (14) | S - Spartial | / ( MNMIN ulp |S| ) where Spartial is the
*> vector of singular values from the partial SVD
*>
*> Test for DGESVDQ:
*>
*> (36) | A - U diag(S) VT | / ( |A| max(M,N) ulp )
*>
*> (37) | I - U'U | / ( M ulp )
*>
*> (38) | I - VT VT' | / ( N ulp )
*>
*> (39) S contains MNMIN nonnegative values in decreasing order.
*> (Return 0 if true, 1/ULP if false.)
*>
*> Test for DGESVJ:
*>
*> (15) | A - U diag(S) VT | / ( |A| max(M,N) ulp )
@ -354,6 +365,8 @@
SUBROUTINE DDRVBD( NSIZES, MM, NN, NTYPES, DOTYPE, ISEED, THRESH,
$ A, LDA, U, LDU, VT, LDVT, ASAV, USAV, VTSAV, S,
$ SSAV, E, WORK, LWORK, IWORK, NOUT, INFO )
*
IMPLICIT NONE
*
* -- LAPACK test routine (version 3.7.0) --
* -- LAPACK is a software package provided by Univ. of Tennessee, --
@ -390,13 +403,19 @@
$ ITEMP, J, JSIZE, JTYPE, LSWORK, M, MINWRK,
$ MMAX, MNMAX, MNMIN, MTYPES, N, NFAIL,
$ NMAX, NS, NSI, NSV, NTEST
DOUBLE PRECISION ANORM, DIF, DIV, OVFL, RTUNFL, ULP,
$ ULPINV, UNFL, VL, VU
DOUBLE PRECISION ANORM, DIF, DIV, OVFL, RTUNFL, ULP,
$ ULPINV, UNFL, VL, VU
* ..
* .. Local Scalars for DGESVDQ ..
INTEGER LIWORK, LRWORK, NUMRANK
* ..
* .. Local Arrays for DGESVDQ ..
DOUBLE PRECISION RWORK( 2 )
* ..
* .. Local Arrays ..
CHARACTER CJOB( 4 ), CJOBR( 3 ), CJOBV( 2 )
INTEGER IOLDSD( 4 ), ISEED2( 4 )
DOUBLE PRECISION RESULT( 40 )
DOUBLE PRECISION RESULT( 39 )
* ..
* .. External Functions ..
DOUBLE PRECISION DLAMCH, DLARND
@ -404,8 +423,8 @@
* ..
* .. External Subroutines ..
EXTERNAL ALASVM, DBDT01, DGEJSV, DGESDD, DGESVD,
$ DGESVDX, DGESVJ, DLABAD, DLACPY, DLASET,
$ DLATMS, DORT01, DORT03, XERBLA
$ DGESVDQ, DGESVDX, DGESVJ, DLABAD, DLACPY,
$ DLASET, DLATMS, DORT01, DORT03, XERBLA
* ..
* .. Intrinsic Functions ..
INTRINSIC ABS, DBLE, INT, MAX, MIN
@ -781,8 +800,64 @@
RESULT( 14 ) = MAX( RESULT( 14 ), DIF )
110 CONTINUE
*
* Test DGESVJ: Factorize A
* Note: DGESVJ does not work for M < N
* Test DGESVDQ
* Note: DGESVDQ only works for M >= N
*
RESULT( 36 ) = ZERO
RESULT( 37 ) = ZERO
RESULT( 38 ) = ZERO
RESULT( 39 ) = ZERO
*
IF( M.GE.N ) THEN
IWTMP = 5*MNMIN*MNMIN + 9*MNMIN + MAX( M, N )
LSWORK = IWTMP + ( IWS-1 )*( LWORK-IWTMP ) / 3
LSWORK = MIN( LSWORK, LWORK )
LSWORK = MAX( LSWORK, 1 )
IF( IWS.EQ.4 )
$ LSWORK = LWORK
*
CALL DLACPY( 'F', M, N, ASAV, LDA, A, LDA )
SRNAMT = 'DGESVDQ'
*
LRWORK = 2
LIWORK = MAX( N, 1 )
CALL DGESVDQ( 'H', 'N', 'N', 'A', 'A',
$ M, N, A, LDA, SSAV, USAV, LDU,
$ VTSAV, LDVT, NUMRANK, IWORK, LIWORK,
$ WORK, LWORK, RWORK, LRWORK, IINFO )
*
IF( IINFO.NE.0 ) THEN
WRITE( NOUT, FMT = 9995 )'DGESVDQ', IINFO, M, N,
$ JTYPE, LSWORK, IOLDSD
INFO = ABS( IINFO )
RETURN
END IF
*
* Do tests 36--39
*
CALL DBDT01( M, N, 0, ASAV, LDA, USAV, LDU, SSAV, E,
$ VTSAV, LDVT, WORK, RESULT( 36 ) )
IF( M.NE.0 .AND. N.NE.0 ) THEN
CALL DORT01( 'Columns', M, M, USAV, LDU, WORK,
$ LWORK, RESULT( 37 ) )
CALL DORT01( 'Rows', N, N, VTSAV, LDVT, WORK,
$ LWORK, RESULT( 38 ) )
END IF
RESULT( 39 ) = ZERO
DO 199 I = 1, MNMIN - 1
IF( SSAV( I ).LT.SSAV( I+1 ) )
$ RESULT( 39 ) = ULPINV
IF( SSAV( I ).LT.ZERO )
$ RESULT( 39 ) = ULPINV
199 CONTINUE
IF( MNMIN.GE.1 ) THEN
IF( SSAV( MNMIN ).LT.ZERO )
$ RESULT( 39 ) = ULPINV
END IF
END IF
*
* Test DGESVJ
* Note: DGESVJ only works for M >= N
*
RESULT( 15 ) = ZERO
RESULT( 16 ) = ZERO
@ -802,8 +877,7 @@
CALL DGESVJ( 'G', 'U', 'V', M, N, USAV, LDA, SSAV,
& 0, A, LDVT, WORK, LWORK, INFO )
*
* DGESVJ retuns V not VT, so we transpose to use the same
* test suite.
* DGESVJ returns V not VT
*
DO J=1,N
DO I=1,N
@ -841,8 +915,8 @@
END IF
END IF
*
* Test DGEJSV: Factorize A
* Note: DGEJSV does not work for M < N
* Test DGEJSV
* Note: DGEJSV only works for M >= N
*
RESULT( 19 ) = ZERO
RESULT( 20 ) = ZERO
@ -862,8 +936,7 @@
& M, N, VTSAV, LDA, SSAV, USAV, LDU, A, LDVT,
& WORK, LWORK, IWORK, INFO )
*
* DGEJSV retuns V not VT, so we transpose to use the same
* test suite.
* DGEJSV returns V not VT
*
DO 140 J=1,N
DO 130 I=1,N
@ -872,7 +945,7 @@
140 END DO
*
IF( IINFO.NE.0 ) THEN
WRITE( NOUT, FMT = 9995 )'GESVJ', IINFO, M, N,
WRITE( NOUT, FMT = 9995 )'GEJSV', IINFO, M, N,
$ JTYPE, LSWORK, IOLDSD
INFO = ABS( IINFO )
RETURN
@ -1086,7 +1159,7 @@
*
* End of Loop -- Check for RESULT(j) > THRESH
*
DO 210 J = 1, 35
DO 210 J = 1, 39
IF( RESULT( J ).GE.THRESH ) THEN
IF( NFAIL.EQ.0 ) THEN
WRITE( NOUT, FMT = 9999 )
@ -1097,7 +1170,7 @@
NFAIL = NFAIL + 1
END IF
210 CONTINUE
NTEST = NTEST + 35
NTEST = NTEST + 39
220 CONTINUE
230 CONTINUE
240 CONTINUE
@ -1158,6 +1231,12 @@
$ ' DGESVDX(V,V,V) ',
$ / '34 = | I - U**T U | / ( M ulp ) ',
$ / '35 = | I - VT VT**T | / ( N ulp ) ',
$ ' DGESVDQ(H,N,N,A,A',
$ / '36 = | A - U diag(S) VT | / ( |A| max(M,N) ulp ) ',
$ / '37 = | I - U**T U | / ( M ulp ) ',
$ / '38 = | I - VT VT**T | / ( N ulp ) ',
$ / '39 = 0 if S contains min(M,N) nonnegative values in',
$ ' decreasing order, else 1/ulp',
$ / / )
9997 FORMAT( ' M=', I5, ', N=', I5, ', type ', I1, ', IWS=', I1,
$ ', seed=', 4( I4, ',' ), ' test(', I2, ')=', G11.4 )

View File

@ -36,6 +36,8 @@
*> DGEJSV compute SVD of an M-by-N matrix A where M >= N
*> DGESVDX compute SVD of an M-by-N matrix A(by bisection
*> and inverse iteration)
*> DGESVDQ compute SVD of an M-by-N matrix A(with a
*> QR-Preconditioned )
*> \endverbatim
*
* Arguments:
@ -100,7 +102,7 @@
* ..
* .. External Subroutines ..
EXTERNAL CHKXER, DGEES, DGEESX, DGEEV, DGEEVX, DGEJSV,
$ DGESDD, DGESVD
$ DGESDD, DGESVD, DGESVDX, DGESVQ
* ..
* .. External Functions ..
LOGICAL DSLECT, LSAMEN
@ -486,6 +488,61 @@
ELSE
WRITE( NOUT, FMT = 9998 )
END IF
*
* Test DGESVDQ
*
SRNAMT = 'DGESVDQ'
INFOT = 1
CALL DGESVDQ( 'X', 'P', 'T', 'A', 'A', 0, 0, A, 1, S, U,
$ 0, VT, 0, NS, IW, 1, W, 1, W, 1, INFO )
CALL CHKXER( 'DGESVDQ', INFOT, NOUT, LERR, OK )
INFOT = 2
CALL DGESVDQ( 'A', 'X', 'T', 'A', 'A', 0, 0, A, 1, S, U,
$ 0, VT, 0, NS, IW, 1, W, 1, W, 1, INFO )
CALL CHKXER( 'DGESVDQ', INFOT, NOUT, LERR, OK )
INFOT = 3
CALL DGESVDQ( 'A', 'P', 'X', 'A', 'A', 0, 0, A, 1, S, U,
$ 0, VT, 0, NS, IW, 1, W, 1, W, 1, INFO )
CALL CHKXER( 'DGESVDQ', INFOT, NOUT, LERR, OK )
INFOT = 4
CALL DGESVDQ( 'A', 'P', 'T', 'X', 'A', 0, 0, A, 1, S, U,
$ 0, VT, 0, NS, IW, 1, W, 1, W, 1, INFO )
CALL CHKXER( 'DGESVDQ', INFOT, NOUT, LERR, OK )
INFOT = 5
CALL DGESVDQ( 'A', 'P', 'T', 'A', 'X', 0, 0, A, 1, S, U,
$ 0, VT, 0, NS, IW, 1, W, 1, W, 1, INFO )
CALL CHKXER( 'DGESVDQ', INFOT, NOUT, LERR, OK )
INFOT = 6
CALL DGESVDQ( 'A', 'P', 'T', 'A', 'A', -1, 0, A, 1, S, U,
$ 0, VT, 0, NS, IW, 1, W, 1, W, 1, INFO )
CALL CHKXER( 'DGESVDQ', INFOT, NOUT, LERR, OK )
INFOT = 7
CALL DGESVDQ( 'A', 'P', 'T', 'A', 'A', 0, 1, A, 1, S, U,
$ 0, VT, 0, NS, IW, 1, W, 1, W, 1, INFO )
CALL CHKXER( 'DGESVDQ', INFOT, NOUT, LERR, OK )
INFOT = 9
CALL DGESVDQ( 'A', 'P', 'T', 'A', 'A', 1, 1, A, 0, S, U,
$ 0, VT, 0, NS, IW, 1, W, 1, W, 1, INFO )
CALL CHKXER( 'DGESVDQ', INFOT, NOUT, LERR, OK )
INFOT = 12
CALL DGESVDQ( 'A', 'P', 'T', 'A', 'A', 1, 1, A, 1, S, U,
$ -1, VT, 0, NS, IW, 1, W, 1, W, 1, INFO )
CALL CHKXER( 'DGESVDQ', INFOT, NOUT, LERR, OK )
INFOT = 14
CALL DGESVDQ( 'A', 'P', 'T', 'A', 'A', 1, 1, A, 1, S, U,
$ 1, VT, -1, NS, IW, 1, W, 1, W, 1, INFO )
CALL CHKXER( 'DGESVDQ', INFOT, NOUT, LERR, OK )
INFOT = 17
CALL DGESVDQ( 'A', 'P', 'T', 'A', 'A', 1, 1, A, 1, S, U,
$ 1, VT, 1, NS, IW, -5, W, 1, W, 1, INFO )
CALL CHKXER( 'DGESVDQ', INFOT, NOUT, LERR, OK )
NT = 11
IF( OK ) THEN
WRITE( NOUT, FMT = 9999 )SRNAMT( 1:LEN_TRIM( SRNAMT ) ),
$ NT
ELSE
WRITE( NOUT, FMT = 9998 )
END IF
END IF
*
* Print a summary line.

View File

@ -194,7 +194,7 @@
VM5( 2 ) = EPS
VM5( 3 ) = SQRT( SMLNUM )
*
* Initalization
* Initialization
*
KNT = 0
RMAX = ZERO

View File

@ -28,15 +28,16 @@
*>
*> DSBT21 generally checks a decomposition of the form
*>
*> A = U S U'
*> A = U S U**T
*>
*> where ' means transpose, A is symmetric banded, U is
*> where **T means transpose, A is symmetric banded, U is
*> orthogonal, and S is diagonal (if KS=0) or symmetric
*> tridiagonal (if KS=1).
*>
*> Specifically:
*>
*> RESULT(1) = | A - U S U' | / ( |A| n ulp ) *andC> RESULT(2) = | I - UU' | / ( n ulp )
*> RESULT(1) = | A - U S U**T | / ( |A| n ulp ) and
*> RESULT(2) = | I - U U**T | / ( n ulp )
*> \endverbatim
*
* Arguments:
@ -214,7 +215,7 @@
*
ANORM = MAX( DLANSB( '1', CUPLO, N, IKA, A, LDA, WORK ), UNFL )
*
* Compute error matrix: Error = A - U S U'
* Compute error matrix: Error = A - U S U**T
*
* Copy A from SB to SP storage format.
*
@ -265,7 +266,7 @@
*
* Do Test 2
*
* Compute UU' - I
* Compute U U**T - I
*
CALL DGEMM( 'N', 'C', N, N, N, ONE, U, LDU, U, LDU, ZERO, WORK,
$ N )

View File

@ -28,9 +28,9 @@
*>
*> DSPT21 generally checks a decomposition of the form
*>
*> A = U S U'
*> A = U S U**T
*>
*> where ' means transpose, A is symmetric (stored in packed format), U
*> where **T means transpose, A is symmetric (stored in packed format), U
*> is orthogonal, and S is diagonal (if KBAND=0) or symmetric
*> tridiagonal (if KBAND=1). If ITYPE=1, then U is represented as a
*> dense matrix, otherwise the U is expressed as a product of
@ -41,15 +41,16 @@
*>
*> Specifically, if ITYPE=1, then:
*>
*> RESULT(1) = | A - U S U' | / ( |A| n ulp ) *andC> RESULT(2) = | I - UU' | / ( n ulp )
*> RESULT(1) = | A - U S U**T | / ( |A| n ulp ) and
*> RESULT(2) = | I - U U**T | / ( n ulp )
*>
*> If ITYPE=2, then:
*>
*> RESULT(1) = | A - V S V' | / ( |A| n ulp )
*> RESULT(1) = | A - V S V**T | / ( |A| n ulp )
*>
*> If ITYPE=3, then:
*>
*> RESULT(1) = | I - VU' | / ( n ulp )
*> RESULT(1) = | I - V U**T | / ( n ulp )
*>
*> Packed storage means that, for example, if UPLO='U', then the columns
*> of the upper triangle of A are stored one after another, so that
@ -70,7 +71,7 @@
*>
*> If UPLO='U', then V = H(n-1)...H(1), where
*>
*> H(j) = I - tau(j) v(j) v(j)'
*> H(j) = I - tau(j) v(j) v(j)**T
*>
*> and the first j-1 elements of v(j) are stored in V(1:j-1,j+1),
*> (i.e., VP( j*(j+1)/2 + 1 : j*(j+1)/2 + j-1 ) ),
@ -78,7 +79,7 @@
*>
*> If UPLO='L', then V = H(1)...H(n-1), where
*>
*> H(j) = I - tau(j) v(j) v(j)'
*> H(j) = I - tau(j) v(j) v(j)**T
*>
*> and the first j elements of v(j) are 0, the (j+1)-st is 1, and the
*> (j+2)-nd through n-th elements are stored in V(j+2:n,j) (i.e.,
@ -93,14 +94,15 @@
*> ITYPE is INTEGER
*> Specifies the type of tests to be performed.
*> 1: U expressed as a dense orthogonal matrix:
*> RESULT(1) = | A - U S U' | / ( |A| n ulp ) *andC> RESULT(2) = | I - UU' | / ( n ulp )
*> RESULT(1) = | A - U S U**T | / ( |A| n ulp ) and
*> RESULT(2) = | I - U U**T | / ( n ulp )
*>
*> 2: U expressed as a product V of Housholder transformations:
*> RESULT(1) = | A - V S V' | / ( |A| n ulp )
*> RESULT(1) = | A - V S V**T | / ( |A| n ulp )
*>
*> 3: U expressed both as a dense orthogonal matrix and
*> as a product of Housholder transformations:
*> RESULT(1) = | I - VU' | / ( n ulp )
*> RESULT(1) = | I - V U**T | / ( n ulp )
*> \endverbatim
*>
*> \param[in] UPLO
@ -183,7 +185,7 @@
*> \verbatim
*> TAU is DOUBLE PRECISION array, dimension (N)
*> If ITYPE >= 2, then TAU(j) is the scalar factor of
*> v(j) v(j)' in the Householder transformation H(j) of
*> v(j) v(j)**T in the Householder transformation H(j) of
*> the product U = H(1)...H(n-2)
*> If ITYPE < 2, then TAU is not referenced.
*> \endverbatim
@ -303,7 +305,7 @@
*
IF( ITYPE.EQ.1 ) THEN
*
* ITYPE=1: error = A - U S U'
* ITYPE=1: error = A - U S U**T
*
CALL DLASET( 'Full', N, N, ZERO, ZERO, WORK, N )
CALL DCOPY( LAP, AP, 1, WORK, 1 )
@ -322,7 +324,7 @@
*
ELSE IF( ITYPE.EQ.2 ) THEN
*
* ITYPE=2: error = V S V' - A
* ITYPE=2: error = V S V**T - A
*
CALL DLASET( 'Full', N, N, ZERO, ZERO, WORK, N )
*
@ -389,7 +391,7 @@
*
ELSE IF( ITYPE.EQ.3 ) THEN
*
* ITYPE=3: error = U V' - I
* ITYPE=3: error = U V**T - I
*
IF( N.LT.2 )
$ RETURN
@ -420,7 +422,7 @@
*
* Do Test 2
*
* Compute UU' - I
* Compute U U**T - I
*
IF( ITYPE.EQ.1 ) THEN
CALL DGEMM( 'N', 'C', N, N, N, ONE, U, LDU, U, LDU, ZERO, WORK,

View File

@ -28,9 +28,9 @@
*>
*> DSYT21 generally checks a decomposition of the form
*>
*> A = U S U'
*> A = U S U**T
*>
*> where ' means transpose, A is symmetric, U is orthogonal, and S is
*> where **T means transpose, A is symmetric, U is orthogonal, and S is
*> diagonal (if KBAND=0) or symmetric tridiagonal (if KBAND=1).
*>
*> If ITYPE=1, then U is represented as a dense matrix; otherwise U is
@ -41,18 +41,19 @@
*>
*> Specifically, if ITYPE=1, then:
*>
*> RESULT(1) = | A - U S U' | / ( |A| n ulp ) *andC> RESULT(2) = | I - UU' | / ( n ulp )
*> RESULT(1) = | A - U S U**T | / ( |A| n ulp ) and
*> RESULT(2) = | I - U U**T | / ( n ulp )
*>
*> If ITYPE=2, then:
*>
*> RESULT(1) = | A - V S V' | / ( |A| n ulp )
*> RESULT(1) = | A - V S V**T | / ( |A| n ulp )
*>
*> If ITYPE=3, then:
*>
*> RESULT(1) = | I - VU' | / ( n ulp )
*> RESULT(1) = | I - V U**T | / ( n ulp )
*>
*> For ITYPE > 1, the transformation U is expressed as a product
*> V = H(1)...H(n-2), where H(j) = I - tau(j) v(j) v(j)' and each
*> V = H(1)...H(n-2), where H(j) = I - tau(j) v(j) v(j)**T and each
*> vector v(j) has its first j elements 0 and the remaining n-j elements
*> stored in V(j+1:n,j).
*> \endverbatim
@ -65,14 +66,15 @@
*> ITYPE is INTEGER
*> Specifies the type of tests to be performed.
*> 1: U expressed as a dense orthogonal matrix:
*> RESULT(1) = | A - U S U' | / ( |A| n ulp ) *andC> RESULT(2) = | I - UU' | / ( n ulp )
*> RESULT(1) = | A - U S U**T | / ( |A| n ulp ) and
*> RESULT(2) = | I - U U**T | / ( n ulp )
*>
*> 2: U expressed as a product V of Housholder transformations:
*> RESULT(1) = | A - V S V' | / ( |A| n ulp )
*> RESULT(1) = | A - V S V**T | / ( |A| n ulp )
*>
*> 3: U expressed both as a dense orthogonal matrix and
*> as a product of Housholder transformations:
*> RESULT(1) = | I - VU' | / ( n ulp )
*> RESULT(1) = | I - V U**T | / ( n ulp )
*> \endverbatim
*>
*> \param[in] UPLO
@ -170,7 +172,7 @@
*> \verbatim
*> TAU is DOUBLE PRECISION array, dimension (N)
*> If ITYPE >= 2, then TAU(j) is the scalar factor of
*> v(j) v(j)' in the Householder transformation H(j) of
*> v(j) v(j)**T in the Householder transformation H(j) of
*> the product U = H(1)...H(n-2)
*> If ITYPE < 2, then TAU is not referenced.
*> \endverbatim
@ -283,7 +285,7 @@
*
IF( ITYPE.EQ.1 ) THEN
*
* ITYPE=1: error = A - U S U'
* ITYPE=1: error = A - U S U**T
*
CALL DLASET( 'Full', N, N, ZERO, ZERO, WORK, N )
CALL DLACPY( CUPLO, N, N, A, LDA, WORK, N )
@ -302,7 +304,7 @@
*
ELSE IF( ITYPE.EQ.2 ) THEN
*
* ITYPE=2: error = V S V' - A
* ITYPE=2: error = V S V**T - A
*
CALL DLASET( 'Full', N, N, ZERO, ZERO, WORK, N )
*
@ -359,7 +361,7 @@
*
ELSE IF( ITYPE.EQ.3 ) THEN
*
* ITYPE=3: error = U V' - I
* ITYPE=3: error = U V**T - I
*
IF( N.LT.2 )
$ RETURN
@ -395,7 +397,7 @@
*
* Do Test 2
*
* Compute UU' - I
* Compute U U**T - I
*
IF( ITYPE.EQ.1 ) THEN
CALL DGEMM( 'N', 'C', N, N, N, ONE, U, LDU, U, LDU, ZERO, WORK,

View File

@ -41,7 +41,8 @@
*>
*> Specifically, if ITYPE=1, then:
*>
*> RESULT(1) = | U' A U - S | / ( |A| m ulp ) *andC> RESULT(2) = | I - U'U | / ( m ulp )
*> RESULT(1) = | U**T A U - S | / ( |A| m ulp ) and
*> RESULT(2) = | I - U**T U | / ( m ulp )
*> \endverbatim
*
* Arguments:
@ -51,7 +52,8 @@
*> ITYPE INTEGER
*> Specifies the type of tests to be performed.
*> 1: U expressed as a dense orthogonal matrix:
*> RESULT(1) = | A - U S U' | / ( |A| n ulp ) *andC> RESULT(2) = | I - UU' | / ( n ulp )
*> RESULT(1) = | A - U S U**T | / ( |A| n ulp ) and
*> RESULT(2) = | I - U U**T | / ( n ulp )
*>
*> UPLO CHARACTER
*> If UPLO='U', the upper triangle of A will be used and the
@ -122,7 +124,7 @@
*>
*> TAU DOUBLE PRECISION array, dimension (N)
*> If ITYPE >= 2, then TAU(j) is the scalar factor of
*> v(j) v(j)' in the Householder transformation H(j) of
*> v(j) v(j)**T in the Householder transformation H(j) of
*> the product U = H(1)...H(n-2)
*> If ITYPE < 2, then TAU is not referenced.
*> Not modified.
@ -207,7 +209,7 @@
*
* Compute error matrix:
*
* ITYPE=1: error = U' A U - S
* ITYPE=1: error = U**T A U - S
*
CALL DSYMM( 'L', UPLO, N, M, ONE, A, LDA, U, LDU, ZERO, WORK, N )
NN = N*N
@ -240,7 +242,7 @@
*
* Do Test 2
*
* Compute U'U - I
* Compute U**T U - I
*
IF( ITYPE.EQ.1 )
$ CALL DORT01( 'Columns', N, M, U, LDU, WORK, 2*N*N,

View File

@ -52,6 +52,7 @@
*> \verbatim
*> A is REAL array, dimension (LDA,N)
*> The m by n matrix A.
*> \endverbatim
*>
*> \param[in] LDA
*> \verbatim

View File

@ -166,7 +166,7 @@
*> SSTEMR('V', 'I')
*>
*> Tests 29 through 34 are disable at present because SSTEMR
*> does not handle partial specturm requests.
*> does not handle partial spectrum requests.
*>
*> (29) | S - Z D Z' | / ( |S| n ulp ) SSTEMR('V', 'I')
*>

View File

@ -187,7 +187,7 @@
*> SSTEMR('V', 'I')
*>
*> Tests 29 through 34 are disable at present because SSTEMR
*> does not handle partial specturm requests.
*> does not handle partial spectrum requests.
*>
*> (29) | S - Z D Z' | / ( |S| n ulp ) SSTEMR('V', 'I')
*>

View File

@ -770,7 +770,7 @@ c MINWRK = MAX( 10*( NSIZE+1 ), 5*NSIZE*NSIZE / 2-2 )
CALL SLACPY( 'Full', MPLUSN, MPLUSN, AI, LDA, A, LDA )
CALL SLACPY( 'Full', MPLUSN, MPLUSN, BI, LDA, B, LDA )
*
* Compute the Schur factorization while swaping the
* Compute the Schur factorization while swapping the
* m-by-m (1,1)-blocks with n-by-n (2,2)-blocks.
*
CALL SGGESX( 'V', 'V', 'S', SLCTSX, 'B', MPLUSN, AI, LDA, BI, LDA,

View File

@ -32,7 +32,7 @@
*> \verbatim
*>
*> SDRVBD checks the singular value decomposition (SVD) drivers
*> SGESVD, SGESDD, SGESVJ, and SGEJSV.
*> SGESVD, SGESDD, SGESVDQ, SGESVJ, SGEJSV, and DGESVDX.
*>
*> Both SGESVD and SGESDD factor A = U diag(S) VT, where U and VT are
*> orthogonal and diag(S) is diagonal with the entries of the array S
@ -90,6 +90,17 @@
*> (14) | S - Spartial | / ( MNMIN ulp |S| ) where Spartial is the
*> vector of singular values from the partial SVD
*>
*> Test for SGESVDQ:
*>
*> (36) | A - U diag(S) VT | / ( |A| max(M,N) ulp )
*>
*> (37) | I - U'U | / ( M ulp )
*>
*> (38) | I - VT VT' | / ( N ulp )
*>
*> (39) S contains MNMIN nonnegative values in decreasing order.
*> (Return 0 if true, 1/ULP if false.)
*>
*> Test for SGESVJ:
*>
*> (15) | A - U diag(S) VT | / ( |A| max(M,N) ulp )
@ -359,6 +370,8 @@
* -- LAPACK is a software package provided by Univ. of Tennessee, --
* -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
* June 2016
*
IMPLICIT NONE
*
* .. Scalar Arguments ..
INTEGER INFO, LDA, LDU, LDVT, LWORK, NOUT, NSIZES,
@ -391,12 +404,18 @@
$ MMAX, MNMAX, MNMIN, MTYPES, N, NFAIL,
$ NMAX, NS, NSI, NSV, NTEST
REAL ANORM, DIF, DIV, OVFL, RTUNFL, ULP,
$ ULPINV, UNFL, VL, VU
$ ULPINV, UNFL, VL, VU
* ..
* .. Local Scalars for DGESVDQ ..
INTEGER LIWORK, LRWORK, NUMRANK
* ..
* .. Local Arrays for DGESVDQ ..
REAL RWORK( 2 )
* ..
* .. Local Arrays ..
CHARACTER CJOB( 4 ), CJOBR( 3 ), CJOBV( 2 )
INTEGER IOLDSD( 4 ), ISEED2( 4 )
REAL RESULT( 40 )
REAL RESULT( 39 )
* ..
* .. External Functions ..
REAL SLAMCH, SLARND
@ -404,8 +423,8 @@
* ..
* .. External Subroutines ..
EXTERNAL ALASVM, SBDT01, SGEJSV, SGESDD, SGESVD,
$ SGESVDX, SGESVJ, SLABAD, SLACPY, SLASET,
$ SLATMS, SORT01, SORT03, XERBLA
$ SGESVDQ, SGESVDX, SGESVJ, SLABAD, SLACPY,
$ SLASET, SLATMS, SORT01, SORT03, XERBLA
* ..
* .. Intrinsic Functions ..
INTRINSIC ABS, REAL, INT, MAX, MIN
@ -781,8 +800,64 @@
RESULT( 14 ) = MAX( RESULT( 14 ), DIF )
110 CONTINUE
*
* Test SGESVJ: Factorize A
* Note: SGESVJ does not work for M < N
* Test SGESVDQ
* Note: SGESVDQ only works for M >= N
*
RESULT( 36 ) = ZERO
RESULT( 37 ) = ZERO
RESULT( 38 ) = ZERO
RESULT( 39 ) = ZERO
*
IF( M.GE.N ) THEN
IWTMP = 5*MNMIN*MNMIN + 9*MNMIN + MAX( M, N )
LSWORK = IWTMP + ( IWS-1 )*( LWORK-IWTMP ) / 3
LSWORK = MIN( LSWORK, LWORK )
LSWORK = MAX( LSWORK, 1 )
IF( IWS.EQ.4 )
$ LSWORK = LWORK
*
CALL SLACPY( 'F', M, N, ASAV, LDA, A, LDA )
SRNAMT = 'SGESVDQ'
*
LRWORK = 2
LIWORK = MAX( N, 1 )
CALL SGESVDQ( 'H', 'N', 'N', 'A', 'A',
$ M, N, A, LDA, SSAV, USAV, LDU,
$ VTSAV, LDVT, NUMRANK, IWORK, LIWORK,
$ WORK, LWORK, RWORK, LRWORK, IINFO )
*
IF( IINFO.NE.0 ) THEN
WRITE( NOUT, FMT = 9995 )'SGESVDQ', IINFO, M, N,
$ JTYPE, LSWORK, IOLDSD
INFO = ABS( IINFO )
RETURN
END IF
*
* Do tests 36--39
*
CALL SBDT01( M, N, 0, ASAV, LDA, USAV, LDU, SSAV, E,
$ VTSAV, LDVT, WORK, RESULT( 36 ) )
IF( M.NE.0 .AND. N.NE.0 ) THEN
CALL SORT01( 'Columns', M, M, USAV, LDU, WORK,
$ LWORK, RESULT( 37 ) )
CALL SORT01( 'Rows', N, N, VTSAV, LDVT, WORK,
$ LWORK, RESULT( 38 ) )
END IF
RESULT( 39 ) = ZERO
DO 199 I = 1, MNMIN - 1
IF( SSAV( I ).LT.SSAV( I+1 ) )
$ RESULT( 39 ) = ULPINV
IF( SSAV( I ).LT.ZERO )
$ RESULT( 39 ) = ULPINV
199 CONTINUE
IF( MNMIN.GE.1 ) THEN
IF( SSAV( MNMIN ).LT.ZERO )
$ RESULT( 39 ) = ULPINV
END IF
END IF
*
* Test SGESVJ
* Note: SGESVJ only works for M >= N
*
RESULT( 15 ) = ZERO
RESULT( 16 ) = ZERO
@ -802,8 +877,7 @@
CALL SGESVJ( 'G', 'U', 'V', M, N, USAV, LDA, SSAV,
& 0, A, LDVT, WORK, LWORK, INFO )
*
* SGESVJ retuns V not VT, so we transpose to use the same
* test suite.
* SGESVJ returns V not VT
*
DO J=1,N
DO I=1,N
@ -841,8 +915,8 @@
END IF
END IF
*
* Test SGEJSV: Factorize A
* Note: SGEJSV does not work for M < N
* Test SGEJSV
* Note: SGEJSV only works for M >= N
*
RESULT( 19 ) = ZERO
RESULT( 20 ) = ZERO
@ -862,8 +936,7 @@
& M, N, VTSAV, LDA, SSAV, USAV, LDU, A, LDVT,
& WORK, LWORK, IWORK, INFO )
*
* SGEJSV retuns V not VT, so we transpose to use the same
* test suite.
* SGEJSV returns V not VT
*
DO 140 J=1,N
DO 130 I=1,N
@ -872,7 +945,7 @@
140 END DO
*
IF( IINFO.NE.0 ) THEN
WRITE( NOUT, FMT = 9995 )'GESVJ', IINFO, M, N,
WRITE( NOUT, FMT = 9995 )'GEJSV', IINFO, M, N,
$ JTYPE, LSWORK, IOLDSD
INFO = ABS( IINFO )
RETURN
@ -1086,7 +1159,7 @@
*
* End of Loop -- Check for RESULT(j) > THRESH
*
DO 210 J = 1, 35
DO 210 J = 1, 39
IF( RESULT( J ).GE.THRESH ) THEN
IF( NFAIL.EQ.0 ) THEN
WRITE( NOUT, FMT = 9999 )
@ -1097,7 +1170,7 @@
NFAIL = NFAIL + 1
END IF
210 CONTINUE
NTEST = NTEST + 35
NTEST = NTEST + 39
220 CONTINUE
230 CONTINUE
240 CONTINUE
@ -1158,6 +1231,12 @@
$ ' SGESVDX(V,V,V) ',
$ / '34 = | I - U**T U | / ( M ulp ) ',
$ / '35 = | I - VT VT**T | / ( N ulp ) ',
$ ' SGESVDQ(H,N,N,A,A',
$ / '36 = | A - U diag(S) VT | / ( |A| max(M,N) ulp ) ',
$ / '37 = | I - U**T U | / ( M ulp ) ',
$ / '38 = | I - VT VT**T | / ( N ulp ) ',
$ / '39 = 0 if S contains min(M,N) nonnegative values in',
$ ' decreasing order, else 1/ulp',
$ / / )
9997 FORMAT( ' M=', I5, ', N=', I5, ', type ', I1, ', IWS=', I1,
$ ', seed=', 4( I4, ',' ), ' test(', I2, ')=', G11.4 )

View File

@ -36,6 +36,8 @@
*> SGEJSV compute SVD of an M-by-N matrix A where M >= N
*> SGESVDX compute SVD of an M-by-N matrix A(by bisection
*> and inverse iteration)
*> SGESVDQ compute SVD of an M-by-N matrix A(with a
*> QR-Preconditioned )
*> \endverbatim
*
* Arguments:
@ -100,7 +102,7 @@
* ..
* .. External Subroutines ..
EXTERNAL CHKXER, SGEES, SGEESX, SGEEV, SGEEVX, SGEJSV,
$ SGESDD, SGESVD
$ SGESDD, SGESVD, SGESVDX, SGESVDQ
* ..
* .. External Functions ..
LOGICAL SSLECT, LSAMEN
@ -486,6 +488,61 @@
ELSE
WRITE( NOUT, FMT = 9998 )
END IF
*
* Test SGESVDQ
*
SRNAMT = 'SGESVDQ'
INFOT = 1
CALL SGESVDQ( 'X', 'P', 'T', 'A', 'A', 0, 0, A, 1, S, U,
$ 0, VT, 0, NS, IW, 1, W, 1, W, 1, INFO )
CALL CHKXER( 'SGESVDQ', INFOT, NOUT, LERR, OK )
INFOT = 2
CALL SGESVDQ( 'A', 'X', 'T', 'A', 'A', 0, 0, A, 1, S, U,
$ 0, VT, 0, NS, IW, 1, W, 1, W, 1, INFO )
CALL CHKXER( 'SGESVDQ', INFOT, NOUT, LERR, OK )
INFOT = 3
CALL SGESVDQ( 'A', 'P', 'X', 'A', 'A', 0, 0, A, 1, S, U,
$ 0, VT, 0, NS, IW, 1, W, 1, W, 1, INFO )
CALL CHKXER( 'SGESVDQ', INFOT, NOUT, LERR, OK )
INFOT = 4
CALL SGESVDQ( 'A', 'P', 'T', 'X', 'A', 0, 0, A, 1, S, U,
$ 0, VT, 0, NS, IW, 1, W, 1, W, 1, INFO )
CALL CHKXER( 'SGESVDQ', INFOT, NOUT, LERR, OK )
INFOT = 5
CALL SGESVDQ( 'A', 'P', 'T', 'A', 'X', 0, 0, A, 1, S, U,
$ 0, VT, 0, NS, IW, 1, W, 1, W, 1, INFO )
CALL CHKXER( 'SGESVDQ', INFOT, NOUT, LERR, OK )
INFOT = 6
CALL SGESVDQ( 'A', 'P', 'T', 'A', 'A', -1, 0, A, 1, S, U,
$ 0, VT, 0, NS, IW, 1, W, 1, W, 1, INFO )
CALL CHKXER( 'SGESVDQ', INFOT, NOUT, LERR, OK )
INFOT = 7
CALL SGESVDQ( 'A', 'P', 'T', 'A', 'A', 0, 1, A, 1, S, U,
$ 0, VT, 0, NS, IW, 1, W, 1, W, 1, INFO )
CALL CHKXER( 'SGESVDQ', INFOT, NOUT, LERR, OK )
INFOT = 9
CALL SGESVDQ( 'A', 'P', 'T', 'A', 'A', 1, 1, A, 0, S, U,
$ 0, VT, 0, NS, IW, 1, W, 1, W, 1, INFO )
CALL CHKXER( 'SGESVDQ', INFOT, NOUT, LERR, OK )
INFOT = 12
CALL SGESVDQ( 'A', 'P', 'T', 'A', 'A', 1, 1, A, 1, S, U,
$ -1, VT, 0, NS, IW, 1, W, 1, W, 1, INFO )
CALL CHKXER( 'SGESVDQ', INFOT, NOUT, LERR, OK )
INFOT = 14
CALL SGESVDQ( 'A', 'P', 'T', 'A', 'A', 1, 1, A, 1, S, U,
$ 1, VT, -1, NS, IW, 1, W, 1, W, 1, INFO )
CALL CHKXER( 'SGESVDQ', INFOT, NOUT, LERR, OK )
INFOT = 17
CALL SGESVDQ( 'A', 'P', 'T', 'A', 'A', 1, 1, A, 1, S, U,
$ 1, VT, 1, NS, IW, -5, W, 1, W, 1, INFO )
CALL CHKXER( 'SGESVDQ', INFOT, NOUT, LERR, OK )
NT = 11
IF( OK ) THEN
WRITE( NOUT, FMT = 9999 )SRNAMT( 1:LEN_TRIM( SRNAMT ) ),
$ NT
ELSE
WRITE( NOUT, FMT = 9998 )
END IF
END IF
*
* Print a summary line.

View File

@ -194,7 +194,7 @@
VM5( 2 ) = EPS
VM5( 3 ) = SQRT( SMLNUM )
*
* Initalization
* Initialization
*
KNT = 0
RMAX = ZERO

View File

@ -28,15 +28,16 @@
*>
*> SSBT21 generally checks a decomposition of the form
*>
*> A = U S U'
*> A = U S U**T
*>
*> where ' means transpose, A is symmetric banded, U is
*> where **T means transpose, A is symmetric banded, U is
*> orthogonal, and S is diagonal (if KS=0) or symmetric
*> tridiagonal (if KS=1).
*>
*> Specifically:
*>
*> RESULT(1) = | A - U S U' | / ( |A| n ulp ) *andC> RESULT(2) = | I - UU' | / ( n ulp )
*> RESULT(1) = | A - U S U**T | / ( |A| n ulp ) and
*> RESULT(2) = | I - U U**T | / ( n ulp )
*> \endverbatim
*
* Arguments:
@ -214,7 +215,7 @@
*
ANORM = MAX( SLANSB( '1', CUPLO, N, IKA, A, LDA, WORK ), UNFL )
*
* Compute error matrix: Error = A - U S U'
* Compute error matrix: Error = A - U S U**T
*
* Copy A from SB to SP storage format.
*
@ -265,7 +266,7 @@
*
* Do Test 2
*
* Compute UU' - I
* Compute U U**T - I
*
CALL SGEMM( 'N', 'C', N, N, N, ONE, U, LDU, U, LDU, ZERO, WORK,
$ N )

View File

@ -28,9 +28,9 @@
*>
*> SSPT21 generally checks a decomposition of the form
*>
*> A = U S U'
*> A = U S U**T
*>
*> where ' means transpose, A is symmetric (stored in packed format), U
*> where **T means transpose, A is symmetric (stored in packed format), U
*> is orthogonal, and S is diagonal (if KBAND=0) or symmetric
*> tridiagonal (if KBAND=1). If ITYPE=1, then U is represented as a
*> dense matrix, otherwise the U is expressed as a product of
@ -41,15 +41,16 @@
*>
*> Specifically, if ITYPE=1, then:
*>
*> RESULT(1) = | A - U S U' | / ( |A| n ulp ) *andC> RESULT(2) = | I - UU' | / ( n ulp )
*> RESULT(1) = | A - U S U**T | / ( |A| n ulp ) and
*> RESULT(2) = | I - U U**T | / ( n ulp )
*>
*> If ITYPE=2, then:
*>
*> RESULT(1) = | A - V S V' | / ( |A| n ulp )
*> RESULT(1) = | A - V S V**T | / ( |A| n ulp )
*>
*> If ITYPE=3, then:
*>
*> RESULT(1) = | I - VU' | / ( n ulp )
*> RESULT(1) = | I - V U**T | / ( n ulp )
*>
*> Packed storage means that, for example, if UPLO='U', then the columns
*> of the upper triangle of A are stored one after another, so that
@ -70,7 +71,7 @@
*>
*> If UPLO='U', then V = H(n-1)...H(1), where
*>
*> H(j) = I - tau(j) v(j) v(j)'
*> H(j) = I - tau(j) v(j) v(j)**T
*>
*> and the first j-1 elements of v(j) are stored in V(1:j-1,j+1),
*> (i.e., VP( j*(j+1)/2 + 1 : j*(j+1)/2 + j-1 ) ),
@ -78,7 +79,7 @@
*>
*> If UPLO='L', then V = H(1)...H(n-1), where
*>
*> H(j) = I - tau(j) v(j) v(j)'
*> H(j) = I - tau(j) v(j) v(j)**T
*>
*> and the first j elements of v(j) are 0, the (j+1)-st is 1, and the
*> (j+2)-nd through n-th elements are stored in V(j+2:n,j) (i.e.,
@ -93,14 +94,15 @@
*> ITYPE is INTEGER
*> Specifies the type of tests to be performed.
*> 1: U expressed as a dense orthogonal matrix:
*> RESULT(1) = | A - U S U' | / ( |A| n ulp ) *andC> RESULT(2) = | I - UU' | / ( n ulp )
*> RESULT(1) = | A - U S U**T | / ( |A| n ulp ) and
*> RESULT(2) = | I - U U**T | / ( n ulp )
*>
*> 2: U expressed as a product V of Housholder transformations:
*> RESULT(1) = | A - V S V' | / ( |A| n ulp )
*> RESULT(1) = | A - V S V**T | / ( |A| n ulp )
*>
*> 3: U expressed both as a dense orthogonal matrix and
*> as a product of Housholder transformations:
*> RESULT(1) = | I - VU' | / ( n ulp )
*> RESULT(1) = | I - V U**T | / ( n ulp )
*> \endverbatim
*>
*> \param[in] UPLO
@ -183,7 +185,7 @@
*> \verbatim
*> TAU is REAL array, dimension (N)
*> If ITYPE >= 2, then TAU(j) is the scalar factor of
*> v(j) v(j)' in the Householder transformation H(j) of
*> v(j) v(j)**T in the Householder transformation H(j) of
*> the product U = H(1)...H(n-2)
*> If ITYPE < 2, then TAU is not referenced.
*> \endverbatim
@ -303,7 +305,7 @@
*
IF( ITYPE.EQ.1 ) THEN
*
* ITYPE=1: error = A - U S U'
* ITYPE=1: error = A - U S U**T
*
CALL SLASET( 'Full', N, N, ZERO, ZERO, WORK, N )
CALL SCOPY( LAP, AP, 1, WORK, 1 )
@ -322,7 +324,7 @@
*
ELSE IF( ITYPE.EQ.2 ) THEN
*
* ITYPE=2: error = V S V' - A
* ITYPE=2: error = V S V**T - A
*
CALL SLASET( 'Full', N, N, ZERO, ZERO, WORK, N )
*
@ -389,7 +391,7 @@
*
ELSE IF( ITYPE.EQ.3 ) THEN
*
* ITYPE=3: error = U V' - I
* ITYPE=3: error = U V**T - I
*
IF( N.LT.2 )
$ RETURN
@ -420,7 +422,7 @@
*
* Do Test 2
*
* Compute UU' - I
* Compute U U**T - I
*
IF( ITYPE.EQ.1 ) THEN
CALL SGEMM( 'N', 'C', N, N, N, ONE, U, LDU, U, LDU, ZERO, WORK,

View File

@ -28,9 +28,9 @@
*>
*> SSYT21 generally checks a decomposition of the form
*>
*> A = U S U'
*> A = U S U**T
*>
*> where ' means transpose, A is symmetric, U is orthogonal, and S is
*> where **T means transpose, A is symmetric, U is orthogonal, and S is
*> diagonal (if KBAND=0) or symmetric tridiagonal (if KBAND=1).
*>
*> If ITYPE=1, then U is represented as a dense matrix; otherwise U is
@ -41,18 +41,19 @@
*>
*> Specifically, if ITYPE=1, then:
*>
*> RESULT(1) = | A - U S U' | / ( |A| n ulp ) *andC> RESULT(2) = | I - UU' | / ( n ulp )
*> RESULT(1) = | A - U S U**T | / ( |A| n ulp ) and
*> RESULT(2) = | I - U U**T | / ( n ulp )
*>
*> If ITYPE=2, then:
*>
*> RESULT(1) = | A - V S V' | / ( |A| n ulp )
*> RESULT(1) = | A - V S V**T | / ( |A| n ulp )
*>
*> If ITYPE=3, then:
*>
*> RESULT(1) = | I - VU' | / ( n ulp )
*> RESULT(1) = | I - V U**T | / ( n ulp )
*>
*> For ITYPE > 1, the transformation U is expressed as a product
*> V = H(1)...H(n-2), where H(j) = I - tau(j) v(j) v(j)' and each
*> V = H(1)...H(n-2), where H(j) = I - tau(j) v(j) v(j)**T and each
*> vector v(j) has its first j elements 0 and the remaining n-j elements
*> stored in V(j+1:n,j).
*> \endverbatim
@ -65,14 +66,15 @@
*> ITYPE is INTEGER
*> Specifies the type of tests to be performed.
*> 1: U expressed as a dense orthogonal matrix:
*> RESULT(1) = | A - U S U' | / ( |A| n ulp ) *andC> RESULT(2) = | I - UU' | / ( n ulp )
*> RESULT(1) = | A - U S U**T | / ( |A| n ulp ) and
*> RESULT(2) = | I - U U**T | / ( n ulp )
*>
*> 2: U expressed as a product V of Housholder transformations:
*> RESULT(1) = | A - V S V' | / ( |A| n ulp )
*> RESULT(1) = | A - V S V**T | / ( |A| n ulp )
*>
*> 3: U expressed both as a dense orthogonal matrix and
*> as a product of Housholder transformations:
*> RESULT(1) = | I - VU' | / ( n ulp )
*> RESULT(1) = | I - V U**T | / ( n ulp )
*> \endverbatim
*>
*> \param[in] UPLO
@ -170,7 +172,7 @@
*> \verbatim
*> TAU is REAL array, dimension (N)
*> If ITYPE >= 2, then TAU(j) is the scalar factor of
*> v(j) v(j)' in the Householder transformation H(j) of
*> v(j) v(j)**T in the Householder transformation H(j) of
*> the product U = H(1)...H(n-2)
*> If ITYPE < 2, then TAU is not referenced.
*> \endverbatim
@ -283,7 +285,7 @@
*
IF( ITYPE.EQ.1 ) THEN
*
* ITYPE=1: error = A - U S U'
* ITYPE=1: error = A - U S U**T
*
CALL SLASET( 'Full', N, N, ZERO, ZERO, WORK, N )
CALL SLACPY( CUPLO, N, N, A, LDA, WORK, N )
@ -302,7 +304,7 @@
*
ELSE IF( ITYPE.EQ.2 ) THEN
*
* ITYPE=2: error = V S V' - A
* ITYPE=2: error = V S V**T - A
*
CALL SLASET( 'Full', N, N, ZERO, ZERO, WORK, N )
*
@ -359,7 +361,7 @@
*
ELSE IF( ITYPE.EQ.3 ) THEN
*
* ITYPE=3: error = U V' - I
* ITYPE=3: error = U V**T - I
*
IF( N.LT.2 )
$ RETURN
@ -395,7 +397,7 @@
*
* Do Test 2
*
* Compute UU' - I
* Compute U U**T - I
*
IF( ITYPE.EQ.1 ) THEN
CALL SGEMM( 'N', 'C', N, N, N, ONE, U, LDU, U, LDU, ZERO, WORK,

View File

@ -41,7 +41,8 @@
*>
*> Specifically, if ITYPE=1, then:
*>
*> RESULT(1) = | U' A U - S | / ( |A| m ulp ) *andC> RESULT(2) = | I - U'U | / ( m ulp )
*> RESULT(1) = | U**T A U - S | / ( |A| m ulp ) and
*> RESULT(2) = | I - U**T U | / ( m ulp )
*> \endverbatim
*
* Arguments:
@ -51,7 +52,8 @@
*> ITYPE INTEGER
*> Specifies the type of tests to be performed.
*> 1: U expressed as a dense orthogonal matrix:
*> RESULT(1) = | A - U S U' | / ( |A| n ulp ) *andC> RESULT(2) = | I - UU' | / ( n ulp )
*> RESULT(1) = | A - U S U**T | / ( |A| n ulp ) and
*> RESULT(2) = | I - U U**T | / ( n ulp )
*>
*> UPLO CHARACTER
*> If UPLO='U', the upper triangle of A will be used and the
@ -122,7 +124,7 @@
*>
*> TAU REAL array, dimension (N)
*> If ITYPE >= 2, then TAU(j) is the scalar factor of
*> v(j) v(j)' in the Householder transformation H(j) of
*> v(j) v(j)**T in the Householder transformation H(j) of
*> the product U = H(1)...H(n-2)
*> If ITYPE < 2, then TAU is not referenced.
*> Not modified.
@ -207,7 +209,7 @@
*
* Compute error matrix:
*
* ITYPE=1: error = U' A U - S
* ITYPE=1: error = U**T A U - S
*
CALL SSYMM( 'L', UPLO, N, M, ONE, A, LDA, U, LDU, ZERO, WORK, N )
NN = N*N
@ -240,7 +242,7 @@
*
* Do Test 2
*
* Compute U'U - I
* Compute U**T U - I
*
IF( ITYPE.EQ.1 )
$ CALL SORT01( 'Columns', N, M, U, LDU, WORK, 2*N*N,

View File

@ -52,6 +52,7 @@
*> \verbatim
*> A is COMPLEX*16 array, dimension (LDA,N)
*> The m by n matrix A.
*> \endverbatim
*>
*> \param[in] LDA
*> \verbatim

View File

@ -167,7 +167,7 @@
*> ZSTEMR('V', 'I')
*>
*> Tests 29 through 34 are disable at present because ZSTEMR
*> does not handle partial specturm requests.
*> does not handle partial spectrum requests.
*>
*> (29) | S - Z D Z* | / ( |S| n ulp ) ZSTEMR('V', 'I')
*>

View File

@ -188,7 +188,7 @@
*> ZSTEMR('V', 'I')
*>
*> Tests 29 through 34 are disable at present because ZSTEMR
*> does not handle partial specturm requests.
*> does not handle partial spectrum requests.
*>
*> (29) | S - Z D Z* | / ( |S| n ulp ) ZSTEMR('V', 'I')
*>

View File

@ -389,7 +389,7 @@
*> \author Univ. of Colorado Denver
*> \author NAG Ltd.
*
*> \date Febuary 2015
*> \date February 2015
*
*> \ingroup complex16_eig
*

View File

@ -738,7 +738,7 @@
CALL ZLACPY( 'Full', MPLUSN, MPLUSN, AI, LDA, A, LDA )
CALL ZLACPY( 'Full', MPLUSN, MPLUSN, BI, LDA, B, LDA )
*
* Compute the Schur factorization while swaping the
* Compute the Schur factorization while swapping the
* m-by-m (1,1)-blocks with n-by-n (2,2)-blocks.
*
CALL ZGGESX( 'V', 'V', 'S', ZLCTSX, 'B', MPLUSN, AI, LDA, BI, LDA,

View File

@ -33,8 +33,9 @@
*>
*> \verbatim
*>
*> ZDRVBD checks the singular value decomposition (SVD) driver ZGESVD
*> and ZGESDD.
*> ZDRVBD checks the singular value decomposition (SVD) driver ZGESVD,
*> ZGESDD, ZGESVJ, ZGEJSV, ZGESVDX, and ZGESVDQ.
*>
*> ZGESVD and ZGESDD factors A = U diag(S) VT, where U and VT are
*> unitary and diag(S) is diagonal with the entries of the array S on
*> its diagonal. The entries of S are the singular values, nonnegative
@ -73,81 +74,92 @@
*>
*> Test for ZGESDD:
*>
*> (1) | A - U diag(S) VT | / ( |A| max(M,N) ulp )
*>
*> (2) | I - U'U | / ( M ulp )
*>
*> (3) | I - VT VT' | / ( N ulp )
*>
*> (4) S contains MNMIN nonnegative values in decreasing order.
*> (Return 0 if true, 1/ULP if false.)
*>
*> (5) | U - Upartial | / ( M ulp ) where Upartial is a partially
*> computed U.
*>
*> (6) | VT - VTpartial | / ( N ulp ) where VTpartial is a partially
*> computed VT.
*>
*> (7) | S - Spartial | / ( MNMIN ulp |S| ) where Spartial is the
*> vector of singular values from the partial SVD
*>
*> Test for ZGESVJ:
*>
*> (1) | A - U diag(S) VT | / ( |A| max(M,N) ulp )
*>
*> (2) | I - U'U | / ( M ulp )
*>
*> (3) | I - VT VT' | / ( N ulp )
*>
*> (4) S contains MNMIN nonnegative values in decreasing order.
*> (Return 0 if true, 1/ULP if false.)
*>
*> Test for ZGEJSV:
*>
*> (1) | A - U diag(S) VT | / ( |A| max(M,N) ulp )
*>
*> (2) | I - U'U | / ( M ulp )
*>
*> (3) | I - VT VT' | / ( N ulp )
*>
*> (4) S contains MNMIN nonnegative values in decreasing order.
*> (Return 0 if true, 1/ULP if false.)
*>
*> Test for ZGESVDX( 'V', 'V', 'A' )/ZGESVDX( 'N', 'N', 'A' )
*>
*> (1) | A - U diag(S) VT | / ( |A| max(M,N) ulp )
*>
*> (2) | I - U'U | / ( M ulp )
*>
*> (3) | I - VT VT' | / ( N ulp )
*>
*> (4) S contains MNMIN nonnegative values in decreasing order.
*> (Return 0 if true, 1/ULP if false.)
*>
*> (5) | U - Upartial | / ( M ulp ) where Upartial is a partially
*> computed U.
*>
*> (6) | VT - VTpartial | / ( N ulp ) where VTpartial is a partially
*> computed VT.
*>
*> (7) | S - Spartial | / ( MNMIN ulp |S| ) where Spartial is the
*> vector of singular values from the partial SVD
*>
*> Test for ZGESVDX( 'V', 'V', 'I' )
*>
*> (8) | U' A VT''' - diag(S) | / ( |A| max(M,N) ulp )
*> (8) | A - U diag(S) VT | / ( |A| max(M,N) ulp )
*>
*> (9) | I - U'U | / ( M ulp )
*>
*> (10) | I - VT VT' | / ( N ulp )
*>
*> (11) S contains MNMIN nonnegative values in decreasing order.
*> (Return 0 if true, 1/ULP if false.)
*>
*> (12) | U - Upartial | / ( M ulp ) where Upartial is a partially
*> computed U.
*>
*> (13) | VT - VTpartial | / ( N ulp ) where VTpartial is a partially
*> computed VT.
*>
*> (14) | S - Spartial | / ( MNMIN ulp |S| ) where Spartial is the
*> vector of singular values from the partial SVD
*>
*> Test for ZGESVDQ:
*>
*> (36) | A - U diag(S) VT | / ( |A| max(M,N) ulp )
*>
*> (37) | I - U'U | / ( M ulp )
*>
*> (38) | I - VT VT' | / ( N ulp )
*>
*> (39) S contains MNMIN nonnegative values in decreasing order.
*> (Return 0 if true, 1/ULP if false.)
*>
*> Test for ZGESVJ:
*>
*> (15) | A - U diag(S) VT | / ( |A| max(M,N) ulp )
*>
*> (16) | I - U'U | / ( M ulp )
*>
*> (17) | I - VT VT' | / ( N ulp )
*>
*> (18) S contains MNMIN nonnegative values in decreasing order.
*> (Return 0 if true, 1/ULP if false.)
*>
*> Test for ZGEJSV:
*>
*> (19) | A - U diag(S) VT | / ( |A| max(M,N) ulp )
*>
*> (20) | I - U'U | / ( M ulp )
*>
*> (21) | I - VT VT' | / ( N ulp )
*>
*> (22) S contains MNMIN nonnegative values in decreasing order.
*> (Return 0 if true, 1/ULP if false.)
*>
*> Test for ZGESVDX( 'V', 'V', 'A' )/ZGESVDX( 'N', 'N', 'A' )
*>
*> (23) | A - U diag(S) VT | / ( |A| max(M,N) ulp )
*>
*> (24) | I - U'U | / ( M ulp )
*>
*> (25) | I - VT VT' | / ( N ulp )
*>
*> (26) S contains MNMIN nonnegative values in decreasing order.
*> (Return 0 if true, 1/ULP if false.)
*>
*> (27) | U - Upartial | / ( M ulp ) where Upartial is a partially
*> computed U.
*>
*> (28) | VT - VTpartial | / ( N ulp ) where VTpartial is a partially
*> computed VT.
*>
*> (29) | S - Spartial | / ( MNMIN ulp |S| ) where Spartial is the
*> vector of singular values from the partial SVD
*>
*> Test for ZGESVDX( 'V', 'V', 'I' )
*>
*> (30) | U' A VT''' - diag(S) | / ( |A| max(M,N) ulp )
*>
*> (31) | I - U'U | / ( M ulp )
*>
*> (32) | I - VT VT' | / ( N ulp )
*>
*> Test for ZGESVDX( 'V', 'V', 'V' )
*>
*> (11) | U' A VT''' - diag(S) | / ( |A| max(M,N) ulp )
*> (33) | U' A VT''' - diag(S) | / ( |A| max(M,N) ulp )
*>
*> (12) | I - U'U | / ( M ulp )
*> (34) | I - U'U | / ( M ulp )
*>
*> (13) | I - VT VT' | / ( N ulp )
*> (35) | I - VT VT' | / ( N ulp )
*>
*> The "sizes" are specified by the arrays MM(1:NSIZES) and
*> NN(1:NSIZES); the value of each element pair (MM(j),NN(j))
@ -393,6 +405,8 @@
* -- LAPACK is a software package provided by Univ. of Tennessee, --
* -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
* June 2016
*
IMPLICIT NONE
*
* .. Scalar Arguments ..
INTEGER INFO, LDA, LDU, LDVT, LWORK, NOUNIT, NSIZES,
@ -411,7 +425,7 @@
* =====================================================================
*
* .. Parameters ..
DOUBLE PRECISION ZERO, ONE, TWO, HALF
DOUBLE PRECISION ZERO, ONE, TWO, HALF
PARAMETER ( ZERO = 0.0D0, ONE = 1.0D0, TWO = 2.0D0,
$ HALF = 0.5D0 )
COMPLEX*16 CZERO, CONE
@ -431,10 +445,13 @@
DOUBLE PRECISION ANORM, DIF, DIV, OVFL, RTUNFL, ULP, ULPINV,
$ UNFL, VL, VU
* ..
* .. Local Scalars for ZGESVDQ ..
INTEGER LIWORK, NUMRANK
* ..
* .. Local Arrays ..
CHARACTER CJOB( 4 ), CJOBR( 3 ), CJOBV( 2 )
INTEGER IOLDSD( 4 ), ISEED2( 4 )
DOUBLE PRECISION RESULT( 35 )
DOUBLE PRECISION RESULT( 39 )
* ..
* .. External Functions ..
DOUBLE PRECISION DLAMCH, DLARND
@ -442,8 +459,8 @@
* ..
* .. External Subroutines ..
EXTERNAL ALASVM, XERBLA, ZBDT01, ZBDT05, ZGESDD,
$ ZGESVD, ZGESVJ, ZGEJSV, ZGESVDX, ZLACPY,
$ ZLASET, ZLATMS, ZUNT01, ZUNT03
$ ZGESVD, ZGESVDQ, ZGESVJ, ZGEJSV, ZGESVDX,
$ ZLACPY, ZLASET, ZLATMS, ZUNT01, ZUNT03
* ..
* .. Intrinsic Functions ..
INTRINSIC ABS, DBLE, MAX, MIN
@ -836,10 +853,65 @@
120 CONTINUE
RESULT( 14 ) = MAX( RESULT( 14 ), DIF )
130 CONTINUE
*
* Test ZGESVJ: Factorize A
* Note: ZGESVJ does not work for M < N
* Test ZGESVDQ
* Note: ZGESVDQ only works for M >= N
*
RESULT( 36 ) = ZERO
RESULT( 37 ) = ZERO
RESULT( 38 ) = ZERO
RESULT( 39 ) = ZERO
*
IF( M.GE.N ) THEN
IWTMP = 2*MNMIN*MNMIN + 2*MNMIN + MAX( M, N )
LSWORK = IWTMP + ( IWSPC-1 )*( LWORK-IWTMP ) / 3
LSWORK = MIN( LSWORK, LWORK )
LSWORK = MAX( LSWORK, 1 )
IF( IWSPC.EQ.4 )
$ LSWORK = LWORK
*
CALL ZLACPY( 'F', M, N, ASAV, LDA, A, LDA )
SRNAMT = 'ZGESVDQ'
*
LRWORK = MAX(2, M, 5*N)
LIWORK = MAX( N, 1 )
CALL ZGESVDQ( 'H', 'N', 'N', 'A', 'A',
$ M, N, A, LDA, SSAV, USAV, LDU,
$ VTSAV, LDVT, NUMRANK, IWORK, LIWORK,
$ WORK, LWORK, RWORK, LRWORK, IINFO )
*
IF( IINFO.NE.0 ) THEN
WRITE( NOUNIT, FMT = 9995 )'ZGESVDQ', IINFO, M, N,
$ JTYPE, LSWORK, IOLDSD
INFO = ABS( IINFO )
RETURN
END IF
*
* Do tests 36--39
*
CALL ZBDT01( M, N, 0, ASAV, LDA, USAV, LDU, SSAV, E,
$ VTSAV, LDVT, WORK, RWORK, RESULT( 36 ) )
IF( M.NE.0 .AND. N.NE.0 ) THEN
CALL ZUNT01( 'Columns', M, M, USAV, LDU, WORK,
$ LWORK, RWORK, RESULT( 37 ) )
CALL ZUNT01( 'Rows', N, N, VTSAV, LDVT, WORK,
$ LWORK, RWORK, RESULT( 38 ) )
END IF
RESULT( 39 ) = ZERO
DO 199 I = 1, MNMIN - 1
IF( SSAV( I ).LT.SSAV( I+1 ) )
$ RESULT( 39 ) = ULPINV
IF( SSAV( I ).LT.ZERO )
$ RESULT( 39 ) = ULPINV
199 CONTINUE
IF( MNMIN.GE.1 ) THEN
IF( SSAV( MNMIN ).LT.ZERO )
$ RESULT( 39 ) = ULPINV
END IF
END IF
*
* Test ZGESVJ
* Note: ZGESVJ only works for M >= N
*
RESULT( 15 ) = ZERO
RESULT( 16 ) = ZERO
@ -847,13 +919,13 @@
RESULT( 18 ) = ZERO
*
IF( M.GE.N ) THEN
IWTMP = 2*MNMIN*MNMIN + 2*MNMIN + MAX( M, N )
LSWORK = IWTMP + ( IWSPC-1 )*( LWORK-IWTMP ) / 3
LSWORK = MIN( LSWORK, LWORK )
LSWORK = MAX( LSWORK, 1 )
LRWORK = MAX(6,N)
IF( IWSPC.EQ.4 )
$ LSWORK = LWORK
IWTMP = 2*MNMIN*MNMIN + 2*MNMIN + MAX( M, N )
LSWORK = IWTMP + ( IWSPC-1 )*( LWORK-IWTMP ) / 3
LSWORK = MIN( LSWORK, LWORK )
LSWORK = MAX( LSWORK, 1 )
LRWORK = MAX(6,N)
IF( IWSPC.EQ.4 )
$ LSWORK = LWORK
*
CALL ZLACPY( 'F', M, N, ASAV, LDA, USAV, LDA )
SRNAMT = 'ZGESVJ'
@ -861,8 +933,7 @@
& 0, A, LDVT, WORK, LWORK, RWORK,
& LRWORK, IINFO )
*
* ZGESVJ retuns V not VT, so we transpose to use the same
* test suite.
* ZGESVJ returns V not VH
*
DO J=1,N
DO I=1,N
@ -900,21 +971,21 @@
END IF
END IF
*
* Test ZGEJSV: Factorize A
* Note: ZGEJSV does not work for M < N
* Test ZGEJSV
* Note: ZGEJSV only works for M >= N
*
RESULT( 19 ) = ZERO
RESULT( 20 ) = ZERO
RESULT( 21 ) = ZERO
RESULT( 22 ) = ZERO
IF( M.GE.N ) THEN
IWTMP = 2*MNMIN*MNMIN + 2*MNMIN + MAX( M, N )
LSWORK = IWTMP + ( IWSPC-1 )*( LWORK-IWTMP ) / 3
LSWORK = MIN( LSWORK, LWORK )
LSWORK = MAX( LSWORK, 1 )
IF( IWSPC.EQ.4 )
$ LSWORK = LWORK
LRWORK = MAX( 7, N + 2*M)
IWTMP = 2*MNMIN*MNMIN + 2*MNMIN + MAX( M, N )
LSWORK = IWTMP + ( IWSPC-1 )*( LWORK-IWTMP ) / 3
LSWORK = MIN( LSWORK, LWORK )
LSWORK = MAX( LSWORK, 1 )
IF( IWSPC.EQ.4 )
$ LSWORK = LWORK
LRWORK = MAX( 7, N + 2*M)
*
CALL ZLACPY( 'F', M, N, ASAV, LDA, VTSAV, LDA )
SRNAMT = 'ZGEJSV'
@ -923,8 +994,7 @@
& WORK, LWORK, RWORK,
& LRWORK, IWORK, IINFO )
*
* ZGEJSV retuns V not VT, so we transpose to use the same
* test suite.
* ZGEJSV returns V not VH
*
DO 133 J=1,N
DO 132 I=1,N
@ -933,7 +1003,7 @@
133 END DO
*
IF( IINFO.NE.0 ) THEN
WRITE( NOUNIT, FMT = 9995 )'GESVJ', IINFO, M, N,
WRITE( NOUNIT, FMT = 9995 )'GEJSV', IINFO, M, N,
$ JTYPE, LSWORK, IOLDSD
INFO = ABS( IINFO )
RETURN
@ -1160,7 +1230,7 @@
*
NTEST = 0
NFAIL = 0
DO 190 J = 1, 35
DO 190 J = 1, 39
IF( RESULT( J ).GE.ZERO )
$ NTEST = NTEST + 1
IF( RESULT( J ).GE.THRESH )
@ -1175,7 +1245,7 @@
NTESTF = 2
END IF
*
DO 200 J = 1, 35
DO 200 J = 1, 39
IF( RESULT( J ).GE.THRESH ) THEN
WRITE( NOUNIT, FMT = 9997 )M, N, JTYPE, IWSPC,
$ IOLDSD, J, RESULT( J )
@ -1251,6 +1321,12 @@
$ / '33 = | U**T A VT**T - diag(S) | / ( |A| max(M,N) ulp )',
$ / '34 = | I - U**T U | / ( M ulp ) ',
$ / '35 = | I - VT VT**T | / ( N ulp ) ',
$ ' ZGESVDQ(H,N,N,A,A',
$ / '36 = | A - U diag(S) VT | / ( |A| max(M,N) ulp ) ',
$ / '37 = | I - U**T U | / ( M ulp ) ',
$ / '38 = | I - VT VT**T | / ( N ulp ) ',
$ / '39 = 0 if S contains min(M,N) nonnegative values in',
$ ' decreasing order, else 1/ulp',
$ / / )
9997 FORMAT( ' M=', I5, ', N=', I5, ', type ', I1, ', IWS=', I1,
$ ', seed=', 4( I4, ',' ), ' test(', I2, ')=', G11.4 )

View File

@ -36,6 +36,8 @@
*> ZGEJSV compute SVD of an M-by-N matrix A where M >= N
*> ZGESVDX compute SVD of an M-by-N matrix A(by bisection
*> and inverse iteration)
*> ZGESVDQ compute SVD of an M-by-N matrix A(with a
*> QR-Preconditioned )
*> \endverbatim
*
* Arguments:
@ -101,7 +103,7 @@
* ..
* .. External Subroutines ..
EXTERNAL CHKXER, ZGEES, ZGEESX, ZGEEV, ZGEEVX, ZGESVJ,
$ ZGESDD, ZGESVD
$ ZGESDD, ZGESVD, ZGESVDX, ZGESVQ
* ..
* .. External Functions ..
LOGICAL LSAMEN, ZSLECT
@ -495,6 +497,61 @@
ELSE
WRITE( NOUT, FMT = 9998 )
END IF
*
* Test ZGESVDQ
*
SRNAMT = 'ZGESVDQ'
INFOT = 1
CALL ZGESVDQ( 'X', 'P', 'T', 'A', 'A', 0, 0, A, 1, S, U,
$ 0, VT, 0, NS, IW, 1, W, 1, RW, 1, INFO )
CALL CHKXER( 'ZGESVDQ', INFOT, NOUT, LERR, OK )
INFOT = 2
CALL ZGESVDQ( 'A', 'X', 'T', 'A', 'A', 0, 0, A, 1, S, U,
$ 0, VT, 0, NS, IW, 1, W, 1, RW, 1, INFO )
CALL CHKXER( 'ZGESVDQ', INFOT, NOUT, LERR, OK )
INFOT = 3
CALL ZGESVDQ( 'A', 'P', 'X', 'A', 'A', 0, 0, A, 1, S, U,
$ 0, VT, 0, NS, IW, 1, W, 1, RW, 1, INFO )
CALL CHKXER( 'ZGESVDQ', INFOT, NOUT, LERR, OK )
INFOT = 4
CALL ZGESVDQ( 'A', 'P', 'T', 'X', 'A', 0, 0, A, 1, S, U,
$ 0, VT, 0, NS, IW, 1, W, 1, RW, 1, INFO )
CALL CHKXER( 'ZGESVDQ', INFOT, NOUT, LERR, OK )
INFOT = 5
CALL ZGESVDQ( 'A', 'P', 'T', 'A', 'X', 0, 0, A, 1, S, U,
$ 0, VT, 0, NS, IW, 1, W, 1, RW, 1, INFO )
CALL CHKXER( 'ZGESVDQ', INFOT, NOUT, LERR, OK )
INFOT = 6
CALL ZGESVDQ( 'A', 'P', 'T', 'A', 'A', -1, 0, A, 1, S, U,
$ 0, VT, 0, NS, IW, 1, W, 1, RW, 1, INFO )
CALL CHKXER( 'ZGESVDQ', INFOT, NOUT, LERR, OK )
INFOT = 7
CALL ZGESVDQ( 'A', 'P', 'T', 'A', 'A', 0, 1, A, 1, S, U,
$ 0, VT, 0, NS, IW, 1, W, 1, RW, 1, INFO )
CALL CHKXER( 'ZGESVDQ', INFOT, NOUT, LERR, OK )
INFOT = 9
CALL ZGESVDQ( 'A', 'P', 'T', 'A', 'A', 1, 1, A, 0, S, U,
$ 0, VT, 0, NS, IW, 1, W, 1, RW, 1, INFO )
CALL CHKXER( 'ZGESVDQ', INFOT, NOUT, LERR, OK )
INFOT = 12
CALL ZGESVDQ( 'A', 'P', 'T', 'A', 'A', 1, 1, A, 1, S, U,
$ -1, VT, 0, NS, IW, 1, W, 1, RW, 1, INFO )
CALL CHKXER( 'ZGESVDQ', INFOT, NOUT, LERR, OK )
INFOT = 14
CALL ZGESVDQ( 'A', 'P', 'T', 'A', 'A', 1, 1, A, 1, S, U,
$ 1, VT, -1, NS, IW, 1, W, 1, RW, 1, INFO )
CALL CHKXER( 'ZGESVDQ', INFOT, NOUT, LERR, OK )
INFOT = 17
CALL ZGESVDQ( 'A', 'P', 'T', 'A', 'A', 1, 1, A, 1, S, U,
$ 1, VT, 1, NS, IW, -5, W, 1, RW, 1, INFO )
CALL CHKXER( 'ZGESVDQ', INFOT, NOUT, LERR, OK )
NT = 11
IF( OK ) THEN
WRITE( NOUT, FMT = 9999 )SRNAMT( 1:LEN_TRIM( SRNAMT ) ),
$ NT
ELSE
WRITE( NOUT, FMT = 9998 )
END IF
END IF
*
* Print a summary line.

View File

@ -29,12 +29,13 @@
*>
*> ZGET51 generally checks a decomposition of the form
*>
*> A = U B VC>
*> where * means conjugate transpose and U and V are unitary.
*> A = U B V**H
*>
*> where **H means conjugate transpose and U and V are unitary.
*>
*> Specifically, if ITYPE=1
*>
*> RESULT = | A - U B V* | / ( |A| n ulp )
*> RESULT = | A - U B V**H | / ( |A| n ulp )
*>
*> If ITYPE=2, then:
*>
@ -42,7 +43,7 @@
*>
*> If ITYPE=3, then:
*>
*> RESULT = | I - UU* | / ( n ulp )
*> RESULT = | I - U U**H | / ( n ulp )
*> \endverbatim
*
* Arguments:
@ -52,9 +53,9 @@
*> \verbatim
*> ITYPE is INTEGER
*> Specifies the type of tests to be performed.
*> =1: RESULT = | A - U B V* | / ( |A| n ulp )
*> =1: RESULT = | A - U B V**H | / ( |A| n ulp )
*> =2: RESULT = | A - B | / ( |A| n ulp )
*> =3: RESULT = | I - UU* | / ( n ulp )
*> =3: RESULT = | I - U U**H | / ( n ulp )
*> \endverbatim
*>
*> \param[in] N
@ -218,7 +219,7 @@
*
IF( ITYPE.EQ.1 ) THEN
*
* ITYPE=1: Compute W = A - UBV'
* ITYPE=1: Compute W = A - U B V**H
*
CALL ZLACPY( ' ', N, N, A, LDA, WORK, N )
CALL ZGEMM( 'N', 'N', N, N, N, CONE, U, LDU, B, LDB, CZERO,
@ -259,7 +260,7 @@
*
* Tests not scaled by norm(A)
*
* ITYPE=3: Compute UU' - I
* ITYPE=3: Compute U U**H - I
*
CALL ZGEMM( 'N', 'C', N, N, N, CONE, U, LDU, U, LDU, CZERO,
$ WORK, N )

View File

@ -28,14 +28,16 @@
*>
*> ZHBT21 generally checks a decomposition of the form
*>
*> A = U S UC>
*> where * means conjugate transpose, A is hermitian banded, U is
*> A = U S U**H
*>
*> where **H means conjugate transpose, A is hermitian banded, U is
*> unitary, and S is diagonal (if KS=0) or symmetric
*> tridiagonal (if KS=1).
*>
*> Specifically:
*>
*> RESULT(1) = | A - U S U* | / ( |A| n ulp ) *andC> RESULT(2) = | I - UU* | / ( n ulp )
*> RESULT(1) = | A - U S U**H | / ( |A| n ulp ) and
*> RESULT(2) = | I - U U**H | / ( n ulp )
*> \endverbatim
*
* Arguments:
@ -220,7 +222,7 @@
*
ANORM = MAX( ZLANHB( '1', CUPLO, N, IKA, A, LDA, RWORK ), UNFL )
*
* Compute error matrix: Error = A - U S U*
* Compute error matrix: Error = A - U S U**H
*
* Copy A from SB to SP storage format.
*
@ -271,7 +273,7 @@
*
* Do Test 2
*
* Compute UU* - I
* Compute U U**H - I
*
CALL ZGEMM( 'N', 'C', N, N, N, CONE, U, LDU, U, LDU, CZERO, WORK,
$ N )

View File

@ -29,8 +29,9 @@
*>
*> ZHET21 generally checks a decomposition of the form
*>
*> A = U S UC>
*> where * means conjugate transpose, A is hermitian, U is unitary, and
*> A = U S U**H
*>
*> where **H means conjugate transpose, A is hermitian, U is unitary, and
*> S is diagonal (if KBAND=0) or (real) symmetric tridiagonal (if
*> KBAND=1).
*>
@ -42,18 +43,19 @@
*>
*> Specifically, if ITYPE=1, then:
*>
*> RESULT(1) = | A - U S U* | / ( |A| n ulp ) *andC> RESULT(2) = | I - UU* | / ( n ulp )
*> RESULT(1) = | A - U S U**H | / ( |A| n ulp ) and
*> RESULT(2) = | I - U U**H | / ( n ulp )
*>
*> If ITYPE=2, then:
*>
*> RESULT(1) = | A - V S V* | / ( |A| n ulp )
*> RESULT(1) = | A - V S V**H | / ( |A| n ulp )
*>
*> If ITYPE=3, then:
*>
*> RESULT(1) = | I - UV* | / ( n ulp )
*> RESULT(1) = | I - U V**H | / ( n ulp )
*>
*> For ITYPE > 1, the transformation U is expressed as a product
*> V = H(1)...H(n-2), where H(j) = I - tau(j) v(j) v(j)C> and each
*> V = H(1)...H(n-2), where H(j) = I - tau(j) v(j) v(j)**H and each
*> vector v(j) has its first j elements 0 and the remaining n-j elements
*> stored in V(j+1:n,j).
*> \endverbatim
@ -66,14 +68,15 @@
*> ITYPE is INTEGER
*> Specifies the type of tests to be performed.
*> 1: U expressed as a dense unitary matrix:
*> RESULT(1) = | A - U S U* | / ( |A| n ulp ) *andC> RESULT(2) = | I - UU* | / ( n ulp )
*> RESULT(1) = | A - U S U**H | / ( |A| n ulp ) and
*> RESULT(2) = | I - U U**H | / ( n ulp )
*>
*> 2: U expressed as a product V of Housholder transformations:
*> RESULT(1) = | A - V S V* | / ( |A| n ulp )
*> RESULT(1) = | A - V S V**H | / ( |A| n ulp )
*>
*> 3: U expressed both as a dense unitary matrix and
*> as a product of Housholder transformations:
*> RESULT(1) = | I - UV* | / ( n ulp )
*> RESULT(1) = | I - U V**H | / ( n ulp )
*> \endverbatim
*>
*> \param[in] UPLO
@ -171,7 +174,7 @@
*> \verbatim
*> TAU is COMPLEX*16 array, dimension (N)
*> If ITYPE >= 2, then TAU(j) is the scalar factor of
*> v(j) v(j)* in the Householder transformation H(j) of
*> v(j) v(j)**H in the Householder transformation H(j) of
*> the product U = H(1)...H(n-2)
*> If ITYPE < 2, then TAU is not referenced.
*> \endverbatim
@ -294,7 +297,7 @@
*
IF( ITYPE.EQ.1 ) THEN
*
* ITYPE=1: error = A - U S U*
* ITYPE=1: error = A - U S U**H
*
CALL ZLASET( 'Full', N, N, CZERO, CZERO, WORK, N )
CALL ZLACPY( CUPLO, N, N, A, LDA, WORK, N )
@ -304,8 +307,7 @@
10 CONTINUE
*
IF( N.GT.1 .AND. KBAND.EQ.1 ) THEN
CMK DO 20 J = 1, N - 1
DO 20 J = 2, N - 1
DO 20 J = 1, N - 1
CALL ZHER2( CUPLO, N, -DCMPLX( E( J ) ), U( 1, J ), 1,
$ U( 1, J-1 ), 1, WORK, N )
20 CONTINUE
@ -314,7 +316,7 @@ CMK DO 20 J = 1, N - 1
*
ELSE IF( ITYPE.EQ.2 ) THEN
*
* ITYPE=2: error = V S V* - A
* ITYPE=2: error = V S V**H - A
*
CALL ZLASET( 'Full', N, N, CZERO, CZERO, WORK, N )
*
@ -371,7 +373,7 @@ CMK DO 20 J = 1, N - 1
*
ELSE IF( ITYPE.EQ.3 ) THEN
*
* ITYPE=3: error = U V* - I
* ITYPE=3: error = U V**H - I
*
IF( N.LT.2 )
$ RETURN
@ -407,7 +409,7 @@ CMK DO 20 J = 1, N - 1
*
* Do Test 2
*
* Compute UU* - I
* Compute U U**H - I
*
IF( ITYPE.EQ.1 ) THEN
CALL ZGEMM( 'N', 'C', N, N, N, CONE, U, LDU, U, LDU, CZERO,

View File

@ -42,7 +42,8 @@
*>
*> Specifically, if ITYPE=1, then:
*>
*> RESULT(1) = | U' A U - S | / ( |A| m ulp ) *andC> RESULT(2) = | I - U'U | / ( m ulp )
*> RESULT(1) = | U**H A U - S | / ( |A| m ulp ) and
*> RESULT(2) = | I - U**H U | / ( m ulp )
*> \endverbatim
*
* Arguments:
@ -52,7 +53,8 @@
*> ITYPE INTEGER
*> Specifies the type of tests to be performed.
*> 1: U expressed as a dense orthogonal matrix:
*> RESULT(1) = | A - U S U' | / ( |A| n ulp ) *andC> RESULT(2) = | I - UU' | / ( n ulp )
*> RESULT(1) = | A - U S U**H | / ( |A| n ulp ) *and
*> RESULT(2) = | I - U U**H | / ( n ulp )
*>
*> UPLO CHARACTER
*> If UPLO='U', the upper triangle of A will be used and the
@ -122,7 +124,7 @@
*>
*> TAU COMPLEX*16 array, dimension (N)
*> If ITYPE >= 2, then TAU(j) is the scalar factor of
*> v(j) v(j)' in the Householder transformation H(j) of
*> v(j) v(j)**H in the Householder transformation H(j) of
*> the product U = H(1)...H(n-2)
*> If ITYPE < 2, then TAU is not referenced.
*> Not modified.
@ -215,7 +217,7 @@
*
* Compute error matrix:
*
* ITYPE=1: error = U' A U - S
* ITYPE=1: error = U**H A U - S
*
CALL ZHEMM( 'L', UPLO, N, M, CONE, A, LDA, U, LDU, CZERO, WORK,
$ N )
@ -249,7 +251,7 @@
*
* Do Test 2
*
* Compute U'U - I
* Compute U**H U - I
*
IF( ITYPE.EQ.1 )
$ CALL ZUNT01( 'Columns', N, M, U, LDU, WORK, 2*N*N, RWORK,

View File

@ -29,8 +29,9 @@
*>
*> ZHPT21 generally checks a decomposition of the form
*>
*> A = U S UC>
*> where * means conjugate transpose, A is hermitian, U is
*> A = U S U**H
*>
*> where **H means conjugate transpose, A is hermitian, U is
*> unitary, and S is diagonal (if KBAND=0) or (real) symmetric
*> tridiagonal (if KBAND=1). If ITYPE=1, then U is represented as
*> a dense matrix, otherwise the U is expressed as a product of
@ -41,15 +42,16 @@
*>
*> Specifically, if ITYPE=1, then:
*>
*> RESULT(1) = | A - U S U* | / ( |A| n ulp ) *andC> RESULT(2) = | I - UU* | / ( n ulp )
*> RESULT(1) = | A - U S U**H | / ( |A| n ulp ) and
*> RESULT(2) = | I - U U**H | / ( n ulp )
*>
*> If ITYPE=2, then:
*>
*> RESULT(1) = | A - V S V* | / ( |A| n ulp )
*> RESULT(1) = | A - V S V**H | / ( |A| n ulp )
*>
*> If ITYPE=3, then:
*>
*> RESULT(1) = | I - UV* | / ( n ulp )
*> RESULT(1) = | I - U V**H | / ( n ulp )
*>
*> Packed storage means that, for example, if UPLO='U', then the columns
*> of the upper triangle of A are stored one after another, so that
@ -70,14 +72,16 @@
*>
*> If UPLO='U', then V = H(n-1)...H(1), where
*>
*> H(j) = I - tau(j) v(j) v(j)C>
*> H(j) = I - tau(j) v(j) v(j)**H
*>
*> and the first j-1 elements of v(j) are stored in V(1:j-1,j+1),
*> (i.e., VP( j*(j+1)/2 + 1 : j*(j+1)/2 + j-1 ) ),
*> the j-th element is 1, and the last n-j elements are 0.
*>
*> If UPLO='L', then V = H(1)...H(n-1), where
*>
*> H(j) = I - tau(j) v(j) v(j)C>
*> H(j) = I - tau(j) v(j) v(j)**H
*>
*> and the first j elements of v(j) are 0, the (j+1)-st is 1, and the
*> (j+2)-nd through n-th elements are stored in V(j+2:n,j) (i.e.,
*> in VP( (2*n-j)*(j-1)/2 + j+2 : (2*n-j)*(j-1)/2 + n ) .)
@ -91,14 +95,15 @@
*> ITYPE is INTEGER
*> Specifies the type of tests to be performed.
*> 1: U expressed as a dense unitary matrix:
*> RESULT(1) = | A - U S U* | / ( |A| n ulp ) *andC> RESULT(2) = | I - UU* | / ( n ulp )
*> RESULT(1) = | A - U S U**H | / ( |A| n ulp ) and
*> RESULT(2) = | I - U U**H | / ( n ulp )
*>
*> 2: U expressed as a product V of Housholder transformations:
*> RESULT(1) = | A - V S V* | / ( |A| n ulp )
*> RESULT(1) = | A - V S V**H | / ( |A| n ulp )
*>
*> 3: U expressed both as a dense unitary matrix and
*> as a product of Housholder transformations:
*> RESULT(1) = | I - UV* | / ( n ulp )
*> RESULT(1) = | I - U V**H | / ( n ulp )
*> \endverbatim
*>
*> \param[in] UPLO
@ -181,7 +186,7 @@
*> \verbatim
*> TAU is COMPLEX*16 array, dimension (N)
*> If ITYPE >= 2, then TAU(j) is the scalar factor of
*> v(j) v(j)* in the Householder transformation H(j) of
*> v(j) v(j)**H in the Householder transformation H(j) of
*> the product U = H(1)...H(n-2)
*> If ITYPE < 2, then TAU is not referenced.
*> \endverbatim
@ -313,7 +318,7 @@
*
IF( ITYPE.EQ.1 ) THEN
*
* ITYPE=1: error = A - U S U*
* ITYPE=1: error = A - U S U**H
*
CALL ZLASET( 'Full', N, N, CZERO, CZERO, WORK, N )
CALL ZCOPY( LAP, AP, 1, WORK, 1 )
@ -323,8 +328,7 @@
10 CONTINUE
*
IF( N.GT.1 .AND. KBAND.EQ.1 ) THEN
CMK DO 20 J = 1, N - 1
DO 20 J = 2, N - 1
DO 20 J = 1, N - 1
CALL ZHPR2( CUPLO, N, -DCMPLX( E( J ) ), U( 1, J ), 1,
$ U( 1, J-1 ), 1, WORK )
20 CONTINUE
@ -333,7 +337,7 @@ CMK DO 20 J = 1, N - 1
*
ELSE IF( ITYPE.EQ.2 ) THEN
*
* ITYPE=2: error = V S V* - A
* ITYPE=2: error = V S V**H - A
*
CALL ZLASET( 'Full', N, N, CZERO, CZERO, WORK, N )
*
@ -401,7 +405,7 @@ CMK DO 20 J = 1, N - 1
*
ELSE IF( ITYPE.EQ.3 ) THEN
*
* ITYPE=3: error = U V* - I
* ITYPE=3: error = U V**H - I
*
IF( N.LT.2 )
$ RETURN
@ -432,7 +436,7 @@ CMK DO 20 J = 1, N - 1
*
* Do Test 2
*
* Compute UU* - I
* Compute U U**H - I
*
IF( ITYPE.EQ.1 ) THEN
CALL ZGEMM( 'N', 'C', N, N, N, CONE, U, LDU, U, LDU, CZERO,

View File

@ -28,14 +28,15 @@
*>
*> ZSTT21 checks a decomposition of the form
*>
*> A = U S UC>
*> where * means conjugate transpose, A is real symmetric tridiagonal,
*> A = U S U**H
*>
*> where **H means conjugate transpose, A is real symmetric tridiagonal,
*> U is unitary, and S is real and diagonal (if KBAND=0) or symmetric
*> tridiagonal (if KBAND=1). Two tests are performed:
*>
*> RESULT(1) = | A - U S U* | / ( |A| n ulp )
*> RESULT(1) = | A - U S U**H | / ( |A| n ulp )
*>
*> RESULT(2) = | I - UU* | / ( n ulp )
*> RESULT(2) = | I - U U**H | / ( n ulp )
*> \endverbatim
*
* Arguments:
@ -228,7 +229,7 @@
*
* Do Test 2
*
* Compute UU* - I
* Compute U U**H - I
*
CALL ZGEMM( 'N', 'C', N, N, N, CONE, U, LDU, U, LDU, CZERO, WORK,
$ N )