Add a BLAS3-based triangular Sylvester equation solver (Reference-LAPACK PR 651)

This commit is contained in:
Martin Kroeker 2022-11-13 23:16:12 +01:00 committed by GitHub
parent 6eb707d941
commit 92174725d9
No known key found for this signature in database
GPG Key ID: 4AEE18F83AFDEB23
14 changed files with 1468 additions and 73 deletions

View File

@ -40,7 +40,7 @@ set(SEIGTST schkee.F
sget54.f sglmts.f sgqrts.f sgrqts.f sgsvts3.f sget54.f sglmts.f sgqrts.f sgrqts.f sgsvts3.f
shst01.f slarfy.f slarhs.f slatm4.f slctes.f slctsx.f slsets.f sort01.f shst01.f slarfy.f slarhs.f slatm4.f slctes.f slctsx.f slsets.f sort01.f
sort03.f ssbt21.f ssgt01.f sslect.f sspt21.f sstt21.f sort03.f ssbt21.f ssgt01.f sslect.f sspt21.f sstt21.f
sstt22.f ssyt21.f ssyt22.f) sstt22.f ssyl01.f ssyt21.f ssyt22.f)
set(CEIGTST cchkee.F set(CEIGTST cchkee.F
cbdt01.f cbdt02.f cbdt03.f cbdt05.f cbdt01.f cbdt02.f cbdt03.f cbdt05.f
@ -56,7 +56,7 @@ set(CEIGTST cchkee.F
cget54.f cglmts.f cgqrts.f cgrqts.f cgsvts3.f cget54.f cglmts.f cgqrts.f cgrqts.f cgsvts3.f
chbt21.f chet21.f chet22.f chpt21.f chst01.f chbt21.f chet21.f chet22.f chpt21.f chst01.f
clarfy.f clarhs.f clatm4.f clctes.f clctsx.f clsets.f csbmv.f clarfy.f clarhs.f clatm4.f clctes.f clctsx.f clsets.f csbmv.f
csgt01.f cslect.f csgt01.f cslect.f csyl01.f
cstt21.f cstt22.f cunt01.f cunt03.f) cstt21.f cstt22.f cunt01.f cunt03.f)
set(DZIGTST dlafts.f dlahd2.f dlasum.f dlatb9.f dstech.f dstect.f set(DZIGTST dlafts.f dlahd2.f dlasum.f dlatb9.f dstech.f dstect.f
@ -77,7 +77,7 @@ set(DEIGTST dchkee.F
dget54.f dglmts.f dgqrts.f dgrqts.f dgsvts3.f dget54.f dglmts.f dgqrts.f dgrqts.f dgsvts3.f
dhst01.f dlarfy.f dlarhs.f dlatm4.f dlctes.f dlctsx.f dlsets.f dort01.f dhst01.f dlarfy.f dlarhs.f dlatm4.f dlctes.f dlctsx.f dlsets.f dort01.f
dort03.f dsbt21.f dsgt01.f dslect.f dspt21.f dstt21.f dort03.f dsbt21.f dsgt01.f dslect.f dspt21.f dstt21.f
dstt22.f dsyt21.f dsyt22.f) dstt22.f dsyl01.f dsyt21.f dsyt22.f)
set(ZEIGTST zchkee.F set(ZEIGTST zchkee.F
zbdt01.f zbdt02.f zbdt03.f zbdt05.f zbdt01.f zbdt02.f zbdt03.f zbdt05.f
@ -93,13 +93,12 @@ set(ZEIGTST zchkee.F
zget54.f zglmts.f zgqrts.f zgrqts.f zgsvts3.f zget54.f zglmts.f zgqrts.f zgrqts.f zgsvts3.f
zhbt21.f zhet21.f zhet22.f zhpt21.f zhst01.f zhbt21.f zhet21.f zhet22.f zhpt21.f zhst01.f
zlarfy.f zlarhs.f zlatm4.f zlctes.f zlctsx.f zlsets.f zsbmv.f zlarfy.f zlarhs.f zlatm4.f zlctes.f zlctsx.f zlsets.f zsbmv.f
zsgt01.f zslect.f zsgt01.f zslect.f zsyl01.f
zstt21.f zstt22.f zunt01.f zunt03.f) zstt21.f zstt22.f zunt01.f zunt03.f)
macro(add_eig_executable name) macro(add_eig_executable name)
add_executable(${name} ${ARGN}) add_executable(${name} ${ARGN})
target_link_libraries(${name} openblas${SUFFIX64_UNDERSCORE}) target_link_libraries(${name} ${TMGLIB} ${LAPACK_LIBRARIES} ${BLAS_LIBRARIES})
#${TMGLIB} ../${LAPACK_LIBRARIES} ${BLAS_LIBRARIES})
endmacro() endmacro()
if(BUILD_SINGLE) if(BUILD_SINGLE)

View File

@ -62,7 +62,7 @@ SEIGTST = schkee.o \
sget54.o sglmts.o sgqrts.o sgrqts.o sgsvts3.o \ sget54.o sglmts.o sgqrts.o sgrqts.o sgsvts3.o \
shst01.o slarfy.o slarhs.o slatm4.o slctes.o slctsx.o slsets.o sort01.o \ shst01.o slarfy.o slarhs.o slatm4.o slctes.o slctsx.o slsets.o sort01.o \
sort03.o ssbt21.o ssgt01.o sslect.o sspt21.o sstt21.o \ sort03.o ssbt21.o ssgt01.o sslect.o sspt21.o sstt21.o \
sstt22.o ssyt21.o ssyt22.o sstt22.o ssyl01.o ssyt21.o ssyt22.o
CEIGTST = cchkee.o \ CEIGTST = cchkee.o \
cbdt01.o cbdt02.o cbdt03.o cbdt05.o \ cbdt01.o cbdt02.o cbdt03.o cbdt05.o \
@ -78,7 +78,7 @@ CEIGTST = cchkee.o \
cget54.o cglmts.o cgqrts.o cgrqts.o cgsvts3.o \ cget54.o cglmts.o cgqrts.o cgrqts.o cgsvts3.o \
chbt21.o chet21.o chet22.o chpt21.o chst01.o \ chbt21.o chet21.o chet22.o chpt21.o chst01.o \
clarfy.o clarhs.o clatm4.o clctes.o clctsx.o clsets.o csbmv.o \ clarfy.o clarhs.o clatm4.o clctes.o clctsx.o clsets.o csbmv.o \
csgt01.o cslect.o \ csgt01.o cslect.o csyl01.o\
cstt21.o cstt22.o cunt01.o cunt03.o cstt21.o cstt22.o cunt01.o cunt03.o
DZIGTST = dlafts.o dlahd2.o dlasum.o dlatb9.o dstech.o dstect.o \ DZIGTST = dlafts.o dlahd2.o dlasum.o dlatb9.o dstech.o dstect.o \
@ -99,7 +99,7 @@ DEIGTST = dchkee.o \
dget54.o dglmts.o dgqrts.o dgrqts.o dgsvts3.o \ dget54.o dglmts.o dgqrts.o dgrqts.o dgsvts3.o \
dhst01.o dlarfy.o dlarhs.o dlatm4.o dlctes.o dlctsx.o dlsets.o dort01.o \ dhst01.o dlarfy.o dlarhs.o dlatm4.o dlctes.o dlctsx.o dlsets.o dort01.o \
dort03.o dsbt21.o dsgt01.o dslect.o dspt21.o dstt21.o \ dort03.o dsbt21.o dsgt01.o dslect.o dspt21.o dstt21.o \
dstt22.o dsyt21.o dsyt22.o dstt22.o dsyl01.o dsyt21.o dsyt22.o
ZEIGTST = zchkee.o \ ZEIGTST = zchkee.o \
zbdt01.o zbdt02.o zbdt03.o zbdt05.o \ zbdt01.o zbdt02.o zbdt03.o zbdt05.o \
@ -115,7 +115,7 @@ ZEIGTST = zchkee.o \
zget54.o zglmts.o zgqrts.o zgrqts.o zgsvts3.o \ zget54.o zglmts.o zgqrts.o zgrqts.o zgsvts3.o \
zhbt21.o zhet21.o zhet22.o zhpt21.o zhst01.o \ zhbt21.o zhet21.o zhet22.o zhpt21.o zhst01.o \
zlarfy.o zlarhs.o zlatm4.o zlctes.o zlctsx.o zlsets.o zsbmv.o \ zlarfy.o zlarhs.o zlatm4.o zlctes.o zlctsx.o zlsets.o zsbmv.o \
zsgt01.o zslect.o \ zsgt01.o zslect.o zsyl01.o\
zstt21.o zstt22.o zunt01.o zunt03.o zstt21.o zstt22.o zunt01.o zunt03.o
.PHONY: all .PHONY: all
@ -127,17 +127,17 @@ complex: xeigtstc
double: xeigtstd double: xeigtstd
complex16: xeigtstz complex16: xeigtstz
xeigtsts: $(SEIGTST) $(SCIGTST) $(AEIGTST) $(TMGLIB) ../$(LAPACKLIB) $(BLASLIB) xeigtsts: $(SEIGTST) $(SCIGTST) $(AEIGTST) $(TMGLIB) $(LAPACKLIB) $(BLASLIB)
$(LOADER) $(FFLAGS) $(LDFLAGS) -o $@ $^ $(FC) $(FFLAGS) $(LDFLAGS) -o $@ $^
xeigtstc: $(CEIGTST) $(SCIGTST) $(AEIGTST) $(TMGLIB) ../$(LAPACKLIB) $(BLASLIB) xeigtstc: $(CEIGTST) $(SCIGTST) $(AEIGTST) $(TMGLIB) $(LAPACKLIB) $(BLASLIB)
$(LOADER) $(FFLAGS) $(LDFLAGS) -o $@ $^ $(FC) $(FFLAGS) $(LDFLAGS) -o $@ $^
xeigtstd: $(DEIGTST) $(DZIGTST) $(AEIGTST) $(TMGLIB) ../$(LAPACKLIB) $(BLASLIB) xeigtstd: $(DEIGTST) $(DZIGTST) $(AEIGTST) $(TMGLIB) $(LAPACKLIB) $(BLASLIB)
$(LOADER) $(FFLAGS) $(LDFLAGS) -o $@ $^ $(FC) $(FFLAGS) $(LDFLAGS) -o $@ $^
xeigtstz: $(ZEIGTST) $(DZIGTST) $(AEIGTST) $(TMGLIB) ../$(LAPACKLIB) $(BLASLIB) xeigtstz: $(ZEIGTST) $(DZIGTST) $(AEIGTST) $(TMGLIB) $(LAPACKLIB) $(BLASLIB)
$(LOADER) $(FFLAGS) $(LDFLAGS) -o $@ $^ $(FC) $(FFLAGS) $(LDFLAGS) -o $@ $^
$(AEIGTST): $(FRC) $(AEIGTST): $(FRC)
$(SCIGTST): $(FRC) $(SCIGTST): $(FRC)

View File

@ -23,7 +23,7 @@
*> \verbatim *> \verbatim
*> *>
*> CCHKEC tests eigen- condition estimation routines *> CCHKEC tests eigen- condition estimation routines
*> CTRSYL, CTREXC, CTRSNA, CTRSEN *> CTRSYL, CTRSYL3, CTREXC, CTRSNA, CTRSEN
*> *>
*> In all cases, the routine runs through a fixed set of numerical *> In all cases, the routine runs through a fixed set of numerical
*> examples, subjects them to various tests, and compares the test *> examples, subjects them to various tests, and compares the test
@ -88,17 +88,17 @@
* .. Local Scalars .. * .. Local Scalars ..
LOGICAL OK LOGICAL OK
CHARACTER*3 PATH CHARACTER*3 PATH
INTEGER KTREXC, KTRSEN, KTRSNA, KTRSYL, LTREXC, LTRSYL, INTEGER KTREXC, KTRSEN, KTRSNA, KTRSYL, KTRSYL3,
$ NTESTS, NTREXC, NTRSYL $ LTREXC, LTRSYL, NTESTS, NTREXC, NTRSYL
REAL EPS, RTREXC, RTRSYL, SFMIN REAL EPS, RTREXC, SFMIN
* .. * ..
* .. Local Arrays .. * .. Local Arrays ..
INTEGER LTRSEN( 3 ), LTRSNA( 3 ), NTRSEN( 3 ), INTEGER FTRSYL( 3 ), ITRSYL( 2 ), LTRSEN( 3 ),
$ NTRSNA( 3 ) $ LTRSNA( 3 ), NTRSEN( 3 ), NTRSNA( 3 )
REAL RTRSEN( 3 ), RTRSNA( 3 ) REAL RTRSEN( 3 ), RTRSNA( 3 ), RTRSYL( 2 )
* .. * ..
* .. External Subroutines .. * .. External Subroutines ..
EXTERNAL CERREC, CGET35, CGET36, CGET37, CGET38 EXTERNAL CERREC, CGET35, CGET36, CGET37, CGET38, CSYL01
* .. * ..
* .. External Functions .. * .. External Functions ..
REAL SLAMCH REAL SLAMCH
@ -120,10 +120,24 @@
$ CALL CERREC( PATH, NOUT ) $ CALL CERREC( PATH, NOUT )
* *
OK = .TRUE. OK = .TRUE.
CALL CGET35( RTRSYL, LTRSYL, NTRSYL, KTRSYL, NIN ) CALL CGET35( RTRSYL( 1 ), LTRSYL, NTRSYL, KTRSYL, NIN )
IF( RTRSYL.GT.THRESH ) THEN IF( RTRSYL( 1 ).GT.THRESH ) THEN
OK = .FALSE. OK = .FALSE.
WRITE( NOUT, FMT = 9999 )RTRSYL, LTRSYL, NTRSYL, KTRSYL WRITE( NOUT, FMT = 9999 )RTRSYL( 1 ), LTRSYL, NTRSYL, KTRSYL
END IF
*
CALL CSYL01( THRESH, FTRSYL, RTRSYL, ITRSYL, KTRSYL3 )
IF( FTRSYL( 1 ).GT.0 ) THEN
OK = .FALSE.
WRITE( NOUT, FMT = 9970 )FTRSYL( 1 ), RTRSYL( 1 ), THRESH
END IF
IF( FTRSYL( 2 ).GT.0 ) THEN
OK = .FALSE.
WRITE( NOUT, FMT = 9971 )FTRSYL( 2 ), RTRSYL( 2 ), THRESH
END IF
IF( FTRSYL( 3 ).GT.0 ) THEN
OK = .FALSE.
WRITE( NOUT, FMT = 9972 )FTRSYL( 3 )
END IF END IF
* *
CALL CGET36( RTREXC, LTREXC, NTREXC, KTREXC, NIN ) CALL CGET36( RTREXC, LTREXC, NTREXC, KTREXC, NIN )
@ -169,6 +183,12 @@
$ / ' Safe minimum (SFMIN) = ', E16.6, / ) $ / ' Safe minimum (SFMIN) = ', E16.6, / )
9992 FORMAT( ' Routines pass computational tests if test ratio is ', 9992 FORMAT( ' Routines pass computational tests if test ratio is ',
$ 'less than', F8.2, / / ) $ 'less than', F8.2, / / )
9972 FORMAT( 'CTRSYL and CTRSYL3 compute an inconsistent scale ',
$ 'factor in ', I8, ' tests.')
9971 FORMAT( 'Error in CTRSYL3: ', I8, ' tests fail the threshold.', /
$ 'Maximum test ratio =', D12.3, ' threshold =', D12.3 )
9970 FORMAT( 'Error in CTRSYL: ', I8, ' tests fail the threshold.', /
$ 'Maximum test ratio =', D12.3, ' threshold =', D12.3 )
RETURN RETURN
* *
* End of CCHKEC * End of CCHKEC

View File

@ -23,7 +23,7 @@
*> *>
*> CERREC tests the error exits for the routines for eigen- condition *> CERREC tests the error exits for the routines for eigen- condition
*> estimation for REAL matrices: *> estimation for REAL matrices:
*> CTRSYL, CTREXC, CTRSNA and CTRSEN. *> CTRSYL, CTRSYL3, CTREXC, CTRSNA and CTRSEN.
*> \endverbatim *> \endverbatim
* *
* Arguments: * Arguments:
@ -77,12 +77,12 @@
* .. * ..
* .. Local Arrays .. * .. Local Arrays ..
LOGICAL SEL( NMAX ) LOGICAL SEL( NMAX )
REAL RW( LW ), S( NMAX ), SEP( NMAX ) REAL RW( LW ), S( NMAX ), SEP( NMAX ), SWORK( NMAX )
COMPLEX A( NMAX, NMAX ), B( NMAX, NMAX ), COMPLEX A( NMAX, NMAX ), B( NMAX, NMAX ),
$ C( NMAX, NMAX ), WORK( LW ), X( NMAX ) $ C( NMAX, NMAX ), WORK( LW ), X( NMAX )
* .. * ..
* .. External Subroutines .. * .. External Subroutines ..
EXTERNAL CHKXER, CTREXC, CTRSEN, CTRSNA, CTRSYL EXTERNAL CHKXER, CTREXC, CTRSEN, CTRSNA, CTRSYL, CTRSYL3
* .. * ..
* .. Scalars in Common .. * .. Scalars in Common ..
LOGICAL LERR, OK LOGICAL LERR, OK
@ -141,6 +141,43 @@
CALL CHKXER( 'CTRSYL', INFOT, NOUT, LERR, OK ) CALL CHKXER( 'CTRSYL', INFOT, NOUT, LERR, OK )
NT = NT + 8 NT = NT + 8
* *
* Test CTRSYL3
*
SRNAMT = 'CTRSYL3'
INFOT = 1
CALL CTRSYL3( 'X', 'N', 1, 0, 0, A, 1, B, 1, C, 1, SCALE,
$ SWORK, NMAX, INFO )
CALL CHKXER( 'CTRSYL3', INFOT, NOUT, LERR, OK )
INFOT = 2
CALL CTRSYL3( 'N', 'X', 1, 0, 0, A, 1, B, 1, C, 1, SCALE,
$ SWORK, NMAX, INFO )
CALL CHKXER( 'CTRSYL3', INFOT, NOUT, LERR, OK )
INFOT = 3
CALL CTRSYL3( 'N', 'N', 0, 0, 0, A, 1, B, 1, C, 1, SCALE,
$ SWORK, NMAX, INFO )
CALL CHKXER( 'CTRSYL3', INFOT, NOUT, LERR, OK )
INFOT = 4
CALL CTRSYL3( 'N', 'N', 1, -1, 0, A, 1, B, 1, C, 1, SCALE,
$ SWORK, NMAX, INFO )
CALL CHKXER( 'CTRSYL3', INFOT, NOUT, LERR, OK )
INFOT = 5
CALL CTRSYL3( 'N', 'N', 1, 0, -1, A, 1, B, 1, C, 1, SCALE,
$ SWORK, NMAX, INFO )
CALL CHKXER( 'CTRSYL3', INFOT, NOUT, LERR, OK )
INFOT = 7
CALL CTRSYL3( 'N', 'N', 1, 2, 0, A, 1, B, 1, C, 2, SCALE,
$ SWORK, NMAX, INFO )
CALL CHKXER( 'CTRSYL3', INFOT, NOUT, LERR, OK )
INFOT = 9
CALL CTRSYL3( 'N', 'N', 1, 0, 2, A, 1, B, 1, C, 1, SCALE,
$ SWORK, NMAX, INFO )
CALL CHKXER( 'CTRSYL3', INFOT, NOUT, LERR, OK )
INFOT = 11
CALL CTRSYL3( 'N', 'N', 1, 2, 0, A, 2, B, 1, C, 1, SCALE,
$ SWORK, NMAX, INFO )
CALL CHKXER( 'CTRSYL3', INFOT, NOUT, LERR, OK )
NT = NT + 8
*
* Test CTREXC * Test CTREXC
* *
SRNAMT = 'CTREXC' SRNAMT = 'CTREXC'

View File

@ -0,0 +1,294 @@
*> \brief \b CSYL01
*
* =========== DOCUMENTATION ===========
*
* Online html documentation available at
* http://www.netlib.org/lapack/explore-html/
*
* Definition:
* ===========
*
* SUBROUTINE CSYL01( THRESH, NFAIL, RMAX, NINFO, KNT )
*
* .. Scalar Arguments ..
* INTEGER KNT
* REAL THRESH
* ..
* .. Array Arguments ..
* INTEGER NFAIL( 3 ), NINFO( 2 )
* REAL RMAX( 2 )
* ..
*
*
*> \par Purpose:
* =============
*>
*> \verbatim
*>
*> CSYL01 tests CTRSYL and CTRSYL3, routines for solving the Sylvester matrix
*> equation
*>
*> op(A)*X + ISGN*X*op(B) = scale*C,
*>
*> where op(A) and op(B) are both upper triangular form, op() represents an
*> optional conjugate transpose, and ISGN can be -1 or +1. Scale is an output
*> less than or equal to 1, chosen to avoid overflow in X.
*>
*> The test code verifies that the following residual does not exceed
*> the provided threshold:
*>
*> norm(op(A)*X + ISGN*X*op(B) - scale*C) /
*> (EPS*max(norm(A),norm(B))*norm(X))
*>
*> This routine complements CGET35 by testing with larger,
*> random matrices, of which some require rescaling of X to avoid overflow.
*>
*> \endverbatim
*
* Arguments:
* ==========
*
*> \param[in] THRESH
*> \verbatim
*> THRESH is REAL
*> A test will count as "failed" if the residual, computed as
*> described above, exceeds THRESH.
*> \endverbatim
*>
*> \param[out] NFAIL
*> \verbatim
*> NFAIL is INTEGER array, dimension (3)
*> NFAIL(1) = No. of times residual CTRSYL exceeds threshold THRESH
*> NFAIL(2) = No. of times residual CTRSYL3 exceeds threshold THRESH
*> NFAIL(3) = No. of times CTRSYL3 and CTRSYL deviate
*> \endverbatim
*>
*> \param[out] RMAX
*> \verbatim
*> RMAX is DOUBLE PRECISION array, dimension (2)
*> RMAX(1) = Value of the largest test ratio of CTRSYL
*> RMAX(2) = Value of the largest test ratio of CTRSYL3
*> \endverbatim
*>
*> \param[out] NINFO
*> \verbatim
*> NINFO is INTEGER array, dimension (2)
*> NINFO(1) = No. of times CTRSYL where INFO is nonzero
*> NINFO(2) = No. of times CTRSYL3 where INFO is nonzero
*> \endverbatim
*>
*> \param[out] KNT
*> \verbatim
*> KNT is INTEGER
*> Total number of examples tested.
*> \endverbatim
*
* -- LAPACK test routine --
SUBROUTINE CSYL01( THRESH, NFAIL, RMAX, NINFO, KNT )
IMPLICIT NONE
*
* -- LAPACK test routine --
* -- LAPACK is a software package provided by Univ. of Tennessee, --
* -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
*
* .. Scalar Arguments ..
INTEGER KNT
REAL THRESH
* ..
* .. Array Arguments ..
INTEGER NFAIL( 3 ), NINFO( 2 )
REAL RMAX( 2 )
* ..
*
* =====================================================================
* ..
* .. Parameters ..
COMPLEX CONE
PARAMETER ( CONE = ( 1.0E+0, 0.0E+0 ) )
REAL ONE, ZERO
PARAMETER ( ZERO = 0.0E+0, ONE = 1.0E+0 )
INTEGER MAXM, MAXN, LDSWORK
PARAMETER ( MAXM = 101, MAXN = 138, LDSWORK = 18 )
* ..
* .. Local Scalars ..
CHARACTER TRANA, TRANB
INTEGER I, INFO, IINFO, ISGN, ITRANA, ITRANB, J, KLA,
$ KUA, KLB, KUB, M, N
REAL ANRM, BNRM, BIGNUM, EPS, RES, RES1,
$ SCALE, SCALE3, SMLNUM, TNRM, XNRM
COMPLEX RMUL
* ..
* .. Local Arrays ..
COMPLEX A( MAXM, MAXM ), B( MAXN, MAXN ),
$ C( MAXM, MAXN ), CC( MAXM, MAXN ),
$ X( MAXM, MAXN ),
$ DUML( MAXM ), DUMR( MAXN ),
$ D( MIN( MAXM, MAXN ) )
REAL SWORK( LDSWORK, 54 ), DUM( MAXN ), VM( 2 )
INTEGER ISEED( 4 ), IWORK( MAXM + MAXN + 2 )
* ..
* .. External Functions ..
LOGICAL SISNAN
REAL SLAMCH, CLANGE
EXTERNAL SISNAN, SLAMCH, CLANGE
* ..
* .. External Subroutines ..
EXTERNAL CLATMR, CLACPY, CGEMM, CTRSYL, CTRSYL3
* ..
* .. Intrinsic Functions ..
INTRINSIC ABS, REAL, MAX
* ..
* .. Executable Statements ..
*
* Get machine parameters
*
EPS = SLAMCH( 'P' )
SMLNUM = SLAMCH( 'S' ) / EPS
BIGNUM = ONE / SMLNUM
*
* Expect INFO = 0
VM( 1 ) = ONE
* Expect INFO = 1
VM( 2 ) = 0.5E+0
*
* Begin test loop
*
NINFO( 1 ) = 0
NINFO( 2 ) = 0
NFAIL( 1 ) = 0
NFAIL( 2 ) = 0
NFAIL( 3 ) = 0
RMAX( 1 ) = ZERO
RMAX( 2 ) = ZERO
KNT = 0
ISEED( 1 ) = 1
ISEED( 2 ) = 1
ISEED( 3 ) = 1
ISEED( 4 ) = 1
SCALE = ONE
SCALE3 = ONE
DO J = 1, 2
DO ISGN = -1, 1, 2
* Reset seed (overwritten by LATMR)
ISEED( 1 ) = 1
ISEED( 2 ) = 1
ISEED( 3 ) = 1
ISEED( 4 ) = 1
DO M = 32, MAXM, 23
KLA = 0
KUA = M - 1
CALL CLATMR( M, M, 'S', ISEED, 'N', D,
$ 6, ONE, CONE, 'T', 'N',
$ DUML, 1, ONE, DUMR, 1, ONE,
$ 'N', IWORK, KLA, KUA, ZERO,
$ ONE, 'NO', A, MAXM, IWORK,
$ IINFO )
DO I = 1, M
A( I, I ) = A( I, I ) * VM( J )
END DO
ANRM = CLANGE( 'M', M, M, A, MAXM, DUM )
DO N = 51, MAXN, 29
KLB = 0
KUB = N - 1
CALL CLATMR( N, N, 'S', ISEED, 'N', D,
$ 6, ONE, CONE, 'T', 'N',
$ DUML, 1, ONE, DUMR, 1, ONE,
$ 'N', IWORK, KLB, KUB, ZERO,
$ ONE, 'NO', B, MAXN, IWORK,
$ IINFO )
DO I = 1, N
B( I, I ) = B( I, I ) * VM ( J )
END DO
BNRM = CLANGE( 'M', N, N, B, MAXN, DUM )
TNRM = MAX( ANRM, BNRM )
CALL CLATMR( M, N, 'S', ISEED, 'N', D,
$ 6, ONE, CONE, 'T', 'N',
$ DUML, 1, ONE, DUMR, 1, ONE,
$ 'N', IWORK, M, N, ZERO, ONE,
$ 'NO', C, MAXM, IWORK, IINFO )
DO ITRANA = 1, 2
IF( ITRANA.EQ.1 )
$ TRANA = 'N'
IF( ITRANA.EQ.2 )
$ TRANA = 'C'
DO ITRANB = 1, 2
IF( ITRANB.EQ.1 )
$ TRANB = 'N'
IF( ITRANB.EQ.2 )
$ TRANB = 'C'
KNT = KNT + 1
*
CALL CLACPY( 'All', M, N, C, MAXM, X, MAXM)
CALL CLACPY( 'All', M, N, C, MAXM, CC, MAXM)
CALL CTRSYL( TRANA, TRANB, ISGN, M, N,
$ A, MAXM, B, MAXN, X, MAXM,
$ SCALE, IINFO )
IF( IINFO.NE.0 )
$ NINFO( 1 ) = NINFO( 1 ) + 1
XNRM = CLANGE( 'M', M, N, X, MAXM, DUM )
RMUL = CONE
IF( XNRM.GT.ONE .AND. TNRM.GT.ONE ) THEN
IF( XNRM.GT.BIGNUM / TNRM ) THEN
RMUL = CONE / MAX( XNRM, TNRM )
END IF
END IF
CALL CGEMM( TRANA, 'N', M, N, M, RMUL,
$ A, MAXM, X, MAXM, -SCALE*RMUL,
$ CC, MAXM )
CALL CGEMM( 'N', TRANB, M, N, N,
$ REAL( ISGN )*RMUL, X, MAXM, B,
$ MAXN, CONE, CC, MAXM )
RES1 = CLANGE( 'M', M, N, CC, MAXM, DUM )
RES = RES1 / MAX( SMLNUM, SMLNUM*XNRM,
$ ( ( ABS( RMUL )*TNRM )*EPS )*XNRM )
IF( RES.GT.THRESH )
$ NFAIL( 1 ) = NFAIL( 1 ) + 1
IF( RES.GT.RMAX( 1 ) )
$ RMAX( 1 ) = RES
*
CALL CLACPY( 'All', M, N, C, MAXM, X, MAXM )
CALL CLACPY( 'All', M, N, C, MAXM, CC, MAXM )
CALL CTRSYL3( TRANA, TRANB, ISGN, M, N,
$ A, MAXM, B, MAXN, X, MAXM,
$ SCALE3, SWORK, LDSWORK, INFO)
IF( INFO.NE.0 )
$ NINFO( 2 ) = NINFO( 2 ) + 1
XNRM = CLANGE( 'M', M, N, X, MAXM, DUM )
RMUL = CONE
IF( XNRM.GT.ONE .AND. TNRM.GT.ONE ) THEN
IF( XNRM.GT.BIGNUM / TNRM ) THEN
RMUL = CONE / MAX( XNRM, TNRM )
END IF
END IF
CALL CGEMM( TRANA, 'N', M, N, M, RMUL,
$ A, MAXM, X, MAXM, -SCALE3*RMUL,
$ CC, MAXM )
CALL CGEMM( 'N', TRANB, M, N, N,
$ REAL( ISGN )*RMUL, X, MAXM, B,
$ MAXN, CONE, CC, MAXM )
RES1 = CLANGE( 'M', M, N, CC, MAXM, DUM )
RES = RES1 / MAX( SMLNUM, SMLNUM*XNRM,
$ ( ( ABS( RMUL )*TNRM )*EPS )*XNRM )
* Verify that TRSYL3 only flushes if TRSYL flushes (but
* there may be cases where TRSYL3 avoid flushing).
IF( SCALE3.EQ.ZERO .AND. SCALE.GT.ZERO .OR.
$ IINFO.NE.INFO ) THEN
NFAIL( 3 ) = NFAIL( 3 ) + 1
END IF
IF( RES.GT.THRESH .OR. SISNAN( RES ) )
$ NFAIL( 2 ) = NFAIL( 2 ) + 1
IF( RES.GT.RMAX( 2 ) )
$ RMAX( 2 ) = RES
END DO
END DO
END DO
END DO
END DO
END DO
*
RETURN
*
* End of CSYL01
*
END

View File

@ -90,21 +90,23 @@
LOGICAL OK LOGICAL OK
CHARACTER*3 PATH CHARACTER*3 PATH
INTEGER KLAEXC, KLALN2, KLANV2, KLAQTR, KLASY2, KTREXC, INTEGER KLAEXC, KLALN2, KLANV2, KLAQTR, KLASY2, KTREXC,
$ KTRSEN, KTRSNA, KTRSYL, LLAEXC, LLALN2, LLANV2, $ KTRSEN, KTRSNA, KTRSYL, KTRSYL3, LLAEXC,
$ LLAQTR, LLASY2, LTREXC, LTRSYL, NLANV2, NLAQTR, $ LLALN2, LLANV2, LLAQTR, LLASY2, LTREXC, LTRSYL,
$ NLASY2, NTESTS, NTRSYL, KTGEXC, NTGEXC, LTGEXC $ NLANV2, NLAQTR, NLASY2, NTESTS, NTRSYL, KTGEXC,
$ LTGEXC
DOUBLE PRECISION EPS, RLAEXC, RLALN2, RLANV2, RLAQTR, RLASY2, DOUBLE PRECISION EPS, RLAEXC, RLALN2, RLANV2, RLAQTR, RLASY2,
$ RTREXC, RTRSYL, SFMIN, RTGEXC $ RTREXC, SFMIN, RTGEXC
* .. * ..
* .. Local Arrays .. * .. Local Arrays ..
INTEGER LTRSEN( 3 ), LTRSNA( 3 ), NLAEXC( 2 ), INTEGER FTRSYL( 3 ), ITRSYL( 2 ), LTRSEN( 3 ),
$ NLALN2( 2 ), NTREXC( 3 ), NTRSEN( 3 ), $ LTRSNA( 3 ), NLAEXC( 2 ), NLALN2( 2 ),
$ NTGEXC( 2 ), NTREXC( 3 ), NTRSEN( 3 ),
$ NTRSNA( 3 ) $ NTRSNA( 3 )
DOUBLE PRECISION RTRSEN( 3 ), RTRSNA( 3 ) DOUBLE PRECISION RTRSEN( 3 ), RTRSNA( 3 ), RTRSYL( 2 )
* .. * ..
* .. External Subroutines .. * .. External Subroutines ..
EXTERNAL DERREC, DGET31, DGET32, DGET33, DGET34, DGET35, EXTERNAL DERREC, DGET31, DGET32, DGET33, DGET34, DGET35,
$ DGET36, DGET37, DGET38, DGET39, DGET40 $ DGET36, DGET37, DGET38, DGET39, DGET40, DSYL01
* .. * ..
* .. External Functions .. * .. External Functions ..
DOUBLE PRECISION DLAMCH DOUBLE PRECISION DLAMCH
@ -153,10 +155,24 @@
WRITE( NOUT, FMT = 9996 )RLAEXC, LLAEXC, NLAEXC, KLAEXC WRITE( NOUT, FMT = 9996 )RLAEXC, LLAEXC, NLAEXC, KLAEXC
END IF END IF
* *
CALL DGET35( RTRSYL, LTRSYL, NTRSYL, KTRSYL ) CALL DGET35( RTRSYL( 1 ), LTRSYL, NTRSYL, KTRSYL )
IF( RTRSYL.GT.THRESH ) THEN IF( RTRSYL( 1 ).GT.THRESH ) THEN
OK = .FALSE. OK = .FALSE.
WRITE( NOUT, FMT = 9995 )RTRSYL, LTRSYL, NTRSYL, KTRSYL WRITE( NOUT, FMT = 9995 )RTRSYL( 1 ), LTRSYL, NTRSYL, KTRSYL
END IF
*
CALL DSYL01( THRESH, FTRSYL, RTRSYL, ITRSYL, KTRSYL3 )
IF( FTRSYL( 1 ).GT.0 ) THEN
OK = .FALSE.
WRITE( NOUT, FMT = 9970 )FTRSYL( 1 ), RTRSYL( 1 ), THRESH
END IF
IF( FTRSYL( 2 ).GT.0 ) THEN
OK = .FALSE.
WRITE( NOUT, FMT = 9971 )FTRSYL( 2 ), RTRSYL( 2 ), THRESH
END IF
IF( FTRSYL( 3 ).GT.0 ) THEN
OK = .FALSE.
WRITE( NOUT, FMT = 9972 )FTRSYL( 3 )
END IF END IF
* *
CALL DGET36( RTREXC, LTREXC, NTREXC, KTREXC, NIN ) CALL DGET36( RTREXC, LTREXC, NTREXC, KTREXC, NIN )
@ -227,7 +243,13 @@
9987 FORMAT( ' Routines pass computational tests if test ratio is les', 9987 FORMAT( ' Routines pass computational tests if test ratio is les',
$ 's than', F8.2, / / ) $ 's than', F8.2, / / )
9986 FORMAT( ' Error in DTGEXC: RMAX =', D12.3, / ' LMAX = ', I8, ' N', 9986 FORMAT( ' Error in DTGEXC: RMAX =', D12.3, / ' LMAX = ', I8, ' N',
$ 'INFO=', I8, ' KNT=', I8 ) $ 'INFO=', 2I8, ' KNT=', I8 )
9972 FORMAT( 'DTRSYL and DTRSYL3 compute an inconsistent result ',
$ 'factor in ', I8, ' tests.')
9971 FORMAT( 'Error in DTRSYL3: ', I8, ' tests fail the threshold.', /
$ 'Maximum test ratio =', D12.3, ' threshold =', D12.3 )
9970 FORMAT( 'Error in DTRSYL: ', I8, ' tests fail the threshold.', /
$ 'Maximum test ratio =', D12.3, ' threshold =', D12.3 )
* *
* End of DCHKEC * End of DCHKEC
* *

View File

@ -23,7 +23,7 @@
*> *>
*> DERREC tests the error exits for the routines for eigen- condition *> DERREC tests the error exits for the routines for eigen- condition
*> estimation for DOUBLE PRECISION matrices: *> estimation for DOUBLE PRECISION matrices:
*> DTRSYL, DTREXC, DTRSNA and DTRSEN. *> DTRSYL, DTRSYL3, DTREXC, DTRSNA and DTRSEN.
*> \endverbatim *> \endverbatim
* *
* Arguments: * Arguments:
@ -82,7 +82,7 @@
$ WI( NMAX ), WORK( NMAX ), WR( NMAX ) $ WI( NMAX ), WORK( NMAX ), WR( NMAX )
* .. * ..
* .. External Subroutines .. * .. External Subroutines ..
EXTERNAL CHKXER, DTREXC, DTRSEN, DTRSNA, DTRSYL EXTERNAL CHKXER, DTREXC, DTRSEN, DTRSNA, DTRSYL, DTRSYL3
* .. * ..
* .. Scalars in Common .. * .. Scalars in Common ..
LOGICAL LERR, OK LOGICAL LERR, OK
@ -141,6 +141,43 @@
CALL CHKXER( 'DTRSYL', INFOT, NOUT, LERR, OK ) CALL CHKXER( 'DTRSYL', INFOT, NOUT, LERR, OK )
NT = NT + 8 NT = NT + 8
* *
* Test DTRSYL3
*
SRNAMT = 'DTRSYL3'
INFOT = 1
CALL DTRSYL3( 'X', 'N', 1, 0, 0, A, 1, B, 1, C, 1, SCALE,
$ IWORK, NMAX, WORK, NMAX, INFO )
CALL CHKXER( 'DTRSYL3', INFOT, NOUT, LERR, OK )
INFOT = 2
CALL DTRSYL3( 'N', 'X', 1, 0, 0, A, 1, B, 1, C, 1, SCALE,
$ IWORK, NMAX, WORK, NMAX, INFO )
CALL CHKXER( 'DTRSYL3', INFOT, NOUT, LERR, OK )
INFOT = 3
CALL DTRSYL3( 'N', 'N', 0, 0, 0, A, 1, B, 1, C, 1, SCALE,
$ IWORK, NMAX, WORK, NMAX, INFO )
CALL CHKXER( 'DTRSYL3', INFOT, NOUT, LERR, OK )
INFOT = 4
CALL DTRSYL3( 'N', 'N', 1, -1, 0, A, 1, B, 1, C, 1, SCALE,
$ IWORK, NMAX, WORK, NMAX, INFO )
CALL CHKXER( 'DTRSYL3', INFOT, NOUT, LERR, OK )
INFOT = 5
CALL DTRSYL3( 'N', 'N', 1, 0, -1, A, 1, B, 1, C, 1, SCALE,
$ IWORK, NMAX, WORK, NMAX, INFO )
CALL CHKXER( 'DTRSYL3', INFOT, NOUT, LERR, OK )
INFOT = 7
CALL DTRSYL3( 'N', 'N', 1, 2, 0, A, 1, B, 1, C, 2, SCALE,
$ IWORK, NMAX, WORK, NMAX, INFO )
CALL CHKXER( 'DTRSYL3', INFOT, NOUT, LERR, OK )
INFOT = 9
CALL DTRSYL3( 'N', 'N', 1, 0, 2, A, 1, B, 1, C, 1, SCALE,
$ IWORK, NMAX, WORK, NMAX, INFO )
CALL CHKXER( 'DTRSYL3', INFOT, NOUT, LERR, OK )
INFOT = 11
CALL DTRSYL3( 'N', 'N', 1, 2, 0, A, 2, B, 1, C, 1, SCALE,
$ IWORK, NMAX, WORK, NMAX, INFO )
CALL CHKXER( 'DTRSYL3', INFOT, NOUT, LERR, OK )
NT = NT + 8
*
* Test DTREXC * Test DTREXC
* *
SRNAMT = 'DTREXC' SRNAMT = 'DTREXC'

View File

@ -0,0 +1,288 @@
*> \brief \b DSYL01
*
* =========== DOCUMENTATION ===========
*
* Online html documentation available at
* http://www.netlib.org/lapack/explore-html/
*
* Definition:
* ===========
*
* SUBROUTINE DSYL01( THRESH, NFAIL, RMAX, NINFO, KNT )
*
* .. Scalar Arguments ..
* INTEGER KNT
* DOUBLE PRECISION THRESH
* ..
* .. Array Arguments ..
* INTEGER NFAIL( 3 ), NINFO( 2 )
* DOUBLE PRECISION RMAX( 2 )
* ..
*
*
*> \par Purpose:
* =============
*>
*> \verbatim
*>
*> DSYL01 tests DTRSYL and DTRSYL3, routines for solving the Sylvester matrix
*> equation
*>
*> op(A)*X + ISGN*X*op(B) = scale*C,
*>
*> A and B are assumed to be in Schur canonical form, op() represents an
*> optional transpose, and ISGN can be -1 or +1. Scale is an output
*> less than or equal to 1, chosen to avoid overflow in X.
*>
*> The test code verifies that the following residual does not exceed
*> the provided threshold:
*>
*> norm(op(A)*X + ISGN*X*op(B) - scale*C) /
*> (EPS*max(norm(A),norm(B))*norm(X))
*>
*> This routine complements DGET35 by testing with larger,
*> random matrices, of which some require rescaling of X to avoid overflow.
*>
*> \endverbatim
*
* Arguments:
* ==========
*
*> \param[in] THRESH
*> \verbatim
*> THRESH is DOUBLE PRECISION
*> A test will count as "failed" if the residual, computed as
*> described above, exceeds THRESH.
*> \endverbatim
*>
*> \param[out] NFAIL
*> \verbatim
*> NFAIL is INTEGER array, dimension (3)
*> NFAIL(1) = No. of times residual DTRSYL exceeds threshold THRESH
*> NFAIL(2) = No. of times residual DTRSYL3 exceeds threshold THRESH
*> NFAIL(3) = No. of times DTRSYL3 and DTRSYL deviate
*> \endverbatim
*>
*> \param[out] RMAX
*> \verbatim
*> RMAX is DOUBLE PRECISION, dimension (2)
*> RMAX(1) = Value of the largest test ratio of DTRSYL
*> RMAX(2) = Value of the largest test ratio of DTRSYL3
*> \endverbatim
*>
*> \param[out] NINFO
*> \verbatim
*> NINFO is INTEGER array, dimension (2)
*> NINFO(1) = No. of times DTRSYL returns an expected INFO
*> NINFO(2) = No. of times DTRSYL3 returns an expected INFO
*> \endverbatim
*>
*> \param[out] KNT
*> \verbatim
*> KNT is INTEGER
*> Total number of examples tested.
*> \endverbatim
*
* -- LAPACK test routine --
SUBROUTINE DSYL01( THRESH, NFAIL, RMAX, NINFO, KNT )
IMPLICIT NONE
*
* -- LAPACK test routine --
* -- LAPACK is a software package provided by Univ. of Tennessee, --
* -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
*
* .. Scalar Arguments ..
INTEGER KNT
DOUBLE PRECISION THRESH
* ..
* .. Array Arguments ..
INTEGER NFAIL( 3 ), NINFO( 2 )
DOUBLE PRECISION RMAX( 2 )
* ..
*
* =====================================================================
* ..
* .. Parameters ..
DOUBLE PRECISION ZERO, ONE
PARAMETER ( ZERO = 0.0D0, ONE = 1.0D0 )
INTEGER MAXM, MAXN, LDSWORK
PARAMETER ( MAXM = 245, MAXN = 192, LDSWORK = 36 )
* ..
* .. Local Scalars ..
CHARACTER TRANA, TRANB
INTEGER I, INFO, IINFO, ISGN, ITRANA, ITRANB, J, KLA,
$ KUA, KLB, KUB, LIWORK, M, N
DOUBLE PRECISION ANRM, BNRM, BIGNUM, EPS, RES, RES1, RMUL,
$ SCALE, SCALE3, SMLNUM, TNRM, XNRM
* ..
* .. Local Arrays ..
DOUBLE PRECISION A( MAXM, MAXM ), B( MAXN, MAXN ),
$ C( MAXM, MAXN ), CC( MAXM, MAXN ),
$ X( MAXM, MAXN ),
$ DUML( MAXM ), DUMR( MAXN ),
$ D( MAX( MAXM, MAXN ) ), DUM( MAXN ),
$ SWORK( LDSWORK, 126 ), VM( 2 )
INTEGER ISEED( 4 ), IWORK( MAXM + MAXN + 2 ), IDUM( 2 )
* ..
* .. External Functions ..
LOGICAL DISNAN
DOUBLE PRECISION DLAMCH, DLANGE
EXTERNAL DLAMCH, DLANGE
* ..
* .. External Subroutines ..
EXTERNAL DLATMR, DLACPY, DGEMM, DTRSYL, DTRSYL3
* ..
* .. Intrinsic Functions ..
INTRINSIC ABS, DBLE, MAX
* ..
* .. Executable Statements ..
*
* Get machine parameters
*
EPS = DLAMCH( 'P' )
SMLNUM = DLAMCH( 'S' ) / EPS
BIGNUM = ONE / SMLNUM
*
VM( 1 ) = ONE
VM( 2 ) = 0.000001D+0
*
* Begin test loop
*
NINFO( 1 ) = 0
NINFO( 2 ) = 0
NFAIL( 1 ) = 0
NFAIL( 2 ) = 0
NFAIL( 3 ) = 0
RMAX( 1 ) = ZERO
RMAX( 2 ) = ZERO
KNT = 0
DO I = 1, 4
ISEED( I ) = 1
END DO
SCALE = ONE
SCALE3 = ONE
LIWORK = MAXM + MAXN + 2
DO J = 1, 2
DO ISGN = -1, 1, 2
* Reset seed (overwritten by LATMR)
DO I = 1, 4
ISEED( I ) = 1
END DO
DO M = 32, MAXM, 71
KLA = 0
KUA = M - 1
CALL DLATMR( M, M, 'S', ISEED, 'N', D,
$ 6, ONE, ONE, 'T', 'N',
$ DUML, 1, ONE, DUMR, 1, ONE,
$ 'N', IWORK, KLA, KUA, ZERO,
$ ONE, 'NO', A, MAXM, IWORK, IINFO )
DO I = 1, M
A( I, I ) = A( I, I ) * VM( J )
END DO
ANRM = DLANGE( 'M', M, M, A, MAXM, DUM )
DO N = 51, MAXN, 47
KLB = 0
KUB = N - 1
CALL DLATMR( N, N, 'S', ISEED, 'N', D,
$ 6, ONE, ONE, 'T', 'N',
$ DUML, 1, ONE, DUMR, 1, ONE,
$ 'N', IWORK, KLB, KUB, ZERO,
$ ONE, 'NO', B, MAXN, IWORK, IINFO )
BNRM = DLANGE( 'M', N, N, B, MAXN, DUM )
TNRM = MAX( ANRM, BNRM )
CALL DLATMR( M, N, 'S', ISEED, 'N', D,
$ 6, ONE, ONE, 'T', 'N',
$ DUML, 1, ONE, DUMR, 1, ONE,
$ 'N', IWORK, M, N, ZERO, ONE,
$ 'NO', C, MAXM, IWORK, IINFO )
DO ITRANA = 1, 2
IF( ITRANA.EQ.1 ) THEN
TRANA = 'N'
END IF
IF( ITRANA.EQ.2 ) THEN
TRANA = 'T'
END IF
DO ITRANB = 1, 2
IF( ITRANB.EQ.1 ) THEN
TRANB = 'N'
END IF
IF( ITRANB.EQ.2 ) THEN
TRANB = 'T'
END IF
KNT = KNT + 1
*
CALL DLACPY( 'All', M, N, C, MAXM, X, MAXM)
CALL DLACPY( 'All', M, N, C, MAXM, CC, MAXM)
CALL DTRSYL( TRANA, TRANB, ISGN, M, N,
$ A, MAXM, B, MAXN, X, MAXM,
$ SCALE, IINFO )
IF( IINFO.NE.0 )
$ NINFO( 1 ) = NINFO( 1 ) + 1
XNRM = DLANGE( 'M', M, N, X, MAXM, DUM )
RMUL = ONE
IF( XNRM.GT.ONE .AND. TNRM.GT.ONE ) THEN
IF( XNRM.GT.BIGNUM / TNRM ) THEN
RMUL = ONE / MAX( XNRM, TNRM )
END IF
END IF
CALL DGEMM( TRANA, 'N', M, N, M, RMUL,
$ A, MAXM, X, MAXM, -SCALE*RMUL,
$ CC, MAXM )
CALL DGEMM( 'N', TRANB, M, N, N,
$ DBLE( ISGN )*RMUL, X, MAXM, B,
$ MAXN, ONE, CC, MAXM )
RES1 = DLANGE( 'M', M, N, CC, MAXM, DUM )
RES = RES1 / MAX( SMLNUM, SMLNUM*XNRM,
$ ( ( RMUL*TNRM )*EPS )*XNRM )
IF( RES.GT.THRESH )
$ NFAIL( 1 ) = NFAIL( 1 ) + 1
IF( RES.GT.RMAX( 1 ) )
$ RMAX( 1 ) = RES
*
CALL DLACPY( 'All', M, N, C, MAXM, X, MAXM )
CALL DLACPY( 'All', M, N, C, MAXM, CC, MAXM )
CALL DTRSYL3( TRANA, TRANB, ISGN, M, N,
$ A, MAXM, B, MAXN, X, MAXM,
$ SCALE3, IWORK, LIWORK,
$ SWORK, LDSWORK, INFO)
IF( INFO.NE.0 )
$ NINFO( 2 ) = NINFO( 2 ) + 1
XNRM = DLANGE( 'M', M, N, X, MAXM, DUM )
RMUL = ONE
IF( XNRM.GT.ONE .AND. TNRM.GT.ONE ) THEN
IF( XNRM.GT.BIGNUM / TNRM ) THEN
RMUL = ONE / MAX( XNRM, TNRM )
END IF
END IF
CALL DGEMM( TRANA, 'N', M, N, M, RMUL,
$ A, MAXM, X, MAXM, -SCALE3*RMUL,
$ CC, MAXM )
CALL DGEMM( 'N', TRANB, M, N, N,
$ DBLE( ISGN )*RMUL, X, MAXM, B,
$ MAXN, ONE, CC, MAXM )
RES1 = DLANGE( 'M', M, N, CC, MAXM, DUM )
RES = RES1 / MAX( SMLNUM, SMLNUM*XNRM,
$ ( ( RMUL*TNRM )*EPS )*XNRM )
* Verify that TRSYL3 only flushes if TRSYL flushes (but
* there may be cases where TRSYL3 avoid flushing).
IF( SCALE3.EQ.ZERO .AND. SCALE.GT.ZERO .OR.
$ IINFO.NE.INFO ) THEN
NFAIL( 3 ) = NFAIL( 3 ) + 1
END IF
IF( RES.GT.THRESH .OR. DISNAN( RES ) )
$ NFAIL( 2 ) = NFAIL( 2 ) + 1
IF( RES.GT.RMAX( 2 ) )
$ RMAX( 2 ) = RES
END DO
END DO
END DO
END DO
END DO
END DO
*
RETURN
*
* End of DSYL01
*
END

View File

@ -90,21 +90,23 @@
LOGICAL OK LOGICAL OK
CHARACTER*3 PATH CHARACTER*3 PATH
INTEGER KLAEXC, KLALN2, KLANV2, KLAQTR, KLASY2, KTREXC, INTEGER KLAEXC, KLALN2, KLANV2, KLAQTR, KLASY2, KTREXC,
$ KTRSEN, KTRSNA, KTRSYL, LLAEXC, LLALN2, LLANV2, $ KTRSEN, KTRSNA, KTRSYL, KTRSYL3, LLAEXC,
$ LLAQTR, LLASY2, LTREXC, LTRSYL, NLANV2, NLAQTR, $ LLALN2, LLANV2, LLAQTR, LLASY2, LTREXC, LTRSYL,
$ NLASY2, NTESTS, NTRSYL, KTGEXC, NTGEXC, LTGEXC $ NLANV2, NLAQTR, NLASY2, NTESTS, NTRSYL, KTGEXC,
$ LTGEXC
REAL EPS, RLAEXC, RLALN2, RLANV2, RLAQTR, RLASY2, REAL EPS, RLAEXC, RLALN2, RLANV2, RLAQTR, RLASY2,
$ RTREXC, RTRSYL, SFMIN, RTGEXC $ RTREXC, SFMIN, RTGEXC
* .. * ..
* .. Local Arrays .. * .. Local Arrays ..
INTEGER LTRSEN( 3 ), LTRSNA( 3 ), NLAEXC( 2 ), INTEGER FTRSYL( 3 ), ITRSYL( 2 ), LTRSEN( 3 ),
$ NLALN2( 2 ), NTREXC( 3 ), NTRSEN( 3 ), $ LTRSNA( 3 ), NLAEXC( 2 ), NLALN2( 2 ),
$ NTGEXC( 2 ), NTREXC( 3 ), NTRSEN( 3 ),
$ NTRSNA( 3 ) $ NTRSNA( 3 )
REAL RTRSEN( 3 ), RTRSNA( 3 ) REAL RTRSEN( 3 ), RTRSNA( 3 ), RTRSYL( 2 )
* .. * ..
* .. External Subroutines .. * .. External Subroutines ..
EXTERNAL SERREC, SGET31, SGET32, SGET33, SGET34, SGET35, EXTERNAL SERREC, SGET31, SGET32, SGET33, SGET34, SGET35,
$ SGET36, SGET37, SGET38, SGET39, SGET40 $ SGET36, SGET37, SGET38, SGET39, SGET40, SSYL01
* .. * ..
* .. External Functions .. * .. External Functions ..
REAL SLAMCH REAL SLAMCH
@ -153,10 +155,24 @@
WRITE( NOUT, FMT = 9996 )RLAEXC, LLAEXC, NLAEXC, KLAEXC WRITE( NOUT, FMT = 9996 )RLAEXC, LLAEXC, NLAEXC, KLAEXC
END IF END IF
* *
CALL SGET35( RTRSYL, LTRSYL, NTRSYL, KTRSYL ) CALL SGET35( RTRSYL( 1 ), LTRSYL, NTRSYL, KTRSYL )
IF( RTRSYL.GT.THRESH ) THEN IF( RTRSYL( 1 ).GT.THRESH ) THEN
OK = .FALSE. OK = .FALSE.
WRITE( NOUT, FMT = 9995 )RTRSYL, LTRSYL, NTRSYL, KTRSYL WRITE( NOUT, FMT = 9995 )RTRSYL( 1 ), LTRSYL, NTRSYL, KTRSYL
END IF
*
CALL SSYL01( THRESH, FTRSYL, RTRSYL, ITRSYL, KTRSYL3 )
IF( FTRSYL( 1 ).GT.0 ) THEN
OK = .FALSE.
WRITE( NOUT, FMT = 9970 )FTRSYL( 1 ), RTRSYL( 1 ), THRESH
END IF
IF( FTRSYL( 2 ).GT.0 ) THEN
OK = .FALSE.
WRITE( NOUT, FMT = 9971 )FTRSYL( 2 ), RTRSYL( 2 ), THRESH
END IF
IF( FTRSYL( 3 ).GT.0 ) THEN
OK = .FALSE.
WRITE( NOUT, FMT = 9972 )FTRSYL( 3 )
END IF END IF
* *
CALL SGET36( RTREXC, LTREXC, NTREXC, KTREXC, NIN ) CALL SGET36( RTREXC, LTREXC, NTREXC, KTREXC, NIN )
@ -227,7 +243,13 @@
9987 FORMAT( ' Routines pass computational tests if test ratio is les', 9987 FORMAT( ' Routines pass computational tests if test ratio is les',
$ 's than', F8.2, / / ) $ 's than', F8.2, / / )
9986 FORMAT( ' Error in STGEXC: RMAX =', E12.3, / ' LMAX = ', I8, ' N', 9986 FORMAT( ' Error in STGEXC: RMAX =', E12.3, / ' LMAX = ', I8, ' N',
$ 'INFO=', I8, ' KNT=', I8 ) $ 'INFO=', 2I8, ' KNT=', I8 )
9972 FORMAT( 'STRSYL and STRSYL3 compute an inconsistent result ',
$ 'factor in ', I8, ' tests.')
9971 FORMAT( 'Error in STRSYL3: ', I8, ' tests fail the threshold.', /
$ 'Maximum test ratio =', D12.3, ' threshold =', D12.3 )
9970 FORMAT( 'Error in STRSYL: ', I8, ' tests fail the threshold.', /
$ 'Maximum test ratio =', D12.3, ' threshold =', D12.3 )
* *
* End of SCHKEC * End of SCHKEC
* *

View File

@ -23,7 +23,7 @@
*> *>
*> SERREC tests the error exits for the routines for eigen- condition *> SERREC tests the error exits for the routines for eigen- condition
*> estimation for REAL matrices: *> estimation for REAL matrices:
*> STRSYL, STREXC, STRSNA and STRSEN. *> STRSYL, STRSYL3, STREXC, STRSNA and STRSEN.
*> \endverbatim *> \endverbatim
* *
* Arguments: * Arguments:
@ -82,7 +82,7 @@
$ WI( NMAX ), WORK( NMAX ), WR( NMAX ) $ WI( NMAX ), WORK( NMAX ), WR( NMAX )
* .. * ..
* .. External Subroutines .. * .. External Subroutines ..
EXTERNAL CHKXER, STREXC, STRSEN, STRSNA, STRSYL EXTERNAL CHKXER, STREXC, STRSEN, STRSNA, STRSYL, STRSYL3
* .. * ..
* .. Scalars in Common .. * .. Scalars in Common ..
LOGICAL LERR, OK LOGICAL LERR, OK
@ -141,6 +141,43 @@
CALL CHKXER( 'STRSYL', INFOT, NOUT, LERR, OK ) CALL CHKXER( 'STRSYL', INFOT, NOUT, LERR, OK )
NT = NT + 8 NT = NT + 8
* *
* Test STRSYL3
*
SRNAMT = 'STRSYL3'
INFOT = 1
CALL STRSYL3( 'X', 'N', 1, 0, 0, A, 1, B, 1, C, 1, SCALE,
$ IWORK, NMAX, WORK, NMAX, INFO )
CALL CHKXER( 'STRSYL3', INFOT, NOUT, LERR, OK )
INFOT = 2
CALL STRSYL3( 'N', 'X', 1, 0, 0, A, 1, B, 1, C, 1, SCALE,
$ IWORK, NMAX, WORK, NMAX, INFO )
CALL CHKXER( 'STRSYL3', INFOT, NOUT, LERR, OK )
INFOT = 3
CALL STRSYL3( 'N', 'N', 0, 0, 0, A, 1, B, 1, C, 1, SCALE,
$ IWORK, NMAX, WORK, NMAX, INFO )
CALL CHKXER( 'STRSYL3', INFOT, NOUT, LERR, OK )
INFOT = 4
CALL STRSYL3( 'N', 'N', 1, -1, 0, A, 1, B, 1, C, 1, SCALE,
$ IWORK, NMAX, WORK, NMAX, INFO )
CALL CHKXER( 'STRSYL3', INFOT, NOUT, LERR, OK )
INFOT = 5
CALL STRSYL3( 'N', 'N', 1, 0, -1, A, 1, B, 1, C, 1, SCALE,
$ IWORK, NMAX, WORK, NMAX, INFO )
CALL CHKXER( 'STRSYL3', INFOT, NOUT, LERR, OK )
INFOT = 7
CALL STRSYL3( 'N', 'N', 1, 2, 0, A, 1, B, 1, C, 2, SCALE,
$ IWORK, NMAX, WORK, NMAX, INFO )
CALL CHKXER( 'STRSYL3', INFOT, NOUT, LERR, OK )
INFOT = 9
CALL STRSYL3( 'N', 'N', 1, 0, 2, A, 1, B, 1, C, 1, SCALE,
$ IWORK, NMAX, WORK, NMAX, INFO )
CALL CHKXER( 'STRSYL3', INFOT, NOUT, LERR, OK )
INFOT = 11
CALL STRSYL3( 'N', 'N', 1, 2, 0, A, 2, B, 1, C, 1, SCALE,
$ IWORK, NMAX, WORK, NMAX, INFO )
CALL CHKXER( 'STRSYL3', INFOT, NOUT, LERR, OK )
NT = NT + 8
*
* Test STREXC * Test STREXC
* *
SRNAMT = 'STREXC' SRNAMT = 'STREXC'

View File

@ -0,0 +1,288 @@
*> \brief \b SSYL01
*
* =========== DOCUMENTATION ===========
*
* Online html documentation available at
* http://www.netlib.org/lapack/explore-html/
*
* Definition:
* ===========
*
* SUBROUTINE SSYL01( THRESH, NFAIL, RMAX, NINFO, KNT )
*
* .. Scalar Arguments ..
* INTEGER KNT
* REAL THRESH
* ..
* .. Array Arguments ..
* INTEGER NFAIL( 3 ), NINFO( 2 )
* REAL RMAX( 2 )
* ..
*
*
*> \par Purpose:
* =============
*>
*> \verbatim
*>
*> SSYL01 tests STRSYL and STRSYL3, routines for solving the Sylvester matrix
*> equation
*>
*> op(A)*X + ISGN*X*op(B) = scale*C,
*>
*> A and B are assumed to be in Schur canonical form, op() represents an
*> optional transpose, and ISGN can be -1 or +1. Scale is an output
*> less than or equal to 1, chosen to avoid overflow in X.
*>
*> The test code verifies that the following residual does not exceed
*> the provided threshold:
*>
*> norm(op(A)*X + ISGN*X*op(B) - scale*C) /
*> (EPS*max(norm(A),norm(B))*norm(X))
*>
*> This routine complements SGET35 by testing with larger,
*> random matrices, of which some require rescaling of X to avoid overflow.
*>
*> \endverbatim
*
* Arguments:
* ==========
*
*> \param[in] THRESH
*> \verbatim
*> THRESH is REAL
*> A test will count as "failed" if the residual, computed as
*> described above, exceeds THRESH.
*> \endverbatim
*>
*> \param[out] NFAIL
*> \verbatim
*> NFAIL is INTEGER array, dimension (3)
*> NFAIL(1) = No. of times residual STRSYL exceeds threshold THRESH
*> NFAIL(2) = No. of times residual STRSYL3 exceeds threshold THRESH
*> NFAIL(3) = No. of times STRSYL3 and STRSYL deviate
*> \endverbatim
*>
*> \param[out] RMAX
*> \verbatim
*> RMAX is REAL, dimension (2)
*> RMAX(1) = Value of the largest test ratio of STRSYL
*> RMAX(2) = Value of the largest test ratio of STRSYL3
*> \endverbatim
*>
*> \param[out] NINFO
*> \verbatim
*> NINFO is INTEGER array, dimension (2)
*> NINFO(1) = No. of times STRSYL returns an expected INFO
*> NINFO(2) = No. of times STRSYL3 returns an expected INFO
*> \endverbatim
*>
*> \param[out] KNT
*> \verbatim
*> KNT is INTEGER
*> Total number of examples tested.
*> \endverbatim
*
* -- LAPACK test routine --
SUBROUTINE SSYL01( THRESH, NFAIL, RMAX, NINFO, KNT )
IMPLICIT NONE
*
* -- LAPACK test routine --
* -- LAPACK is a software package provided by Univ. of Tennessee, --
* -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
*
* .. Scalar Arguments ..
INTEGER KNT
REAL THRESH
* ..
* .. Array Arguments ..
INTEGER NFAIL( 3 ), NINFO( 2 )
REAL RMAX( 2 )
* ..
*
* =====================================================================
* ..
* .. Parameters ..
REAL ZERO, ONE
PARAMETER ( ZERO = 0.0E+0, ONE = 1.0E+0 )
INTEGER MAXM, MAXN, LDSWORK
PARAMETER ( MAXM = 101, MAXN = 138, LDSWORK = 18 )
* ..
* .. Local Scalars ..
CHARACTER TRANA, TRANB
INTEGER I, INFO, IINFO, ISGN, ITRANA, ITRANB, J, KLA,
$ KUA, KLB, KUB, LIWORK, M, N
REAL ANRM, BNRM, BIGNUM, EPS, RES, RES1, RMUL,
$ SCALE, SCALE3, SMLNUM, TNRM, XNRM
* ..
* .. Local Arrays ..
REAL A( MAXM, MAXM ), B( MAXN, MAXN ),
$ C( MAXM, MAXN ), CC( MAXM, MAXN ),
$ X( MAXM, MAXN ),
$ DUML( MAXM ), DUMR( MAXN ),
$ D( MAX( MAXM, MAXN ) ), DUM( MAXN ),
$ SWORK( LDSWORK, 54 ), VM( 2 )
INTEGER ISEED( 4 ), IWORK( MAXM + MAXN + 2 ), IDUM( 2 )
* ..
* .. External Functions ..
LOGICAL SISNAN
REAL SLAMCH, SLANGE
EXTERNAL SISNAN, SLAMCH, SLANGE
* ..
* .. External Subroutines ..
EXTERNAL SLATMR, SLACPY, SGEMM, STRSYL, STRSYL3
* ..
* .. Intrinsic Functions ..
INTRINSIC ABS, REAL, MAX
* ..
* .. Executable Statements ..
*
* Get machine parameters
*
EPS = SLAMCH( 'P' )
SMLNUM = SLAMCH( 'S' ) / EPS
BIGNUM = ONE / SMLNUM
*
VM( 1 ) = ONE
VM( 2 ) = 0.05E+0
*
* Begin test loop
*
NINFO( 1 ) = 0
NINFO( 2 ) = 0
NFAIL( 1 ) = 0
NFAIL( 2 ) = 0
NFAIL( 3 ) = 0
RMAX( 1 ) = ZERO
RMAX( 2 ) = ZERO
KNT = 0
DO I = 1, 4
ISEED( I ) = 1
END DO
SCALE = ONE
SCALE3 = ONE
LIWORK = MAXM + MAXN + 2
DO J = 1, 2
DO ISGN = -1, 1, 2
* Reset seed (overwritten by LATMR)
DO I = 1, 4
ISEED( I ) = 1
END DO
DO M = 32, MAXM, 71
KLA = 0
KUA = M - 1
CALL SLATMR( M, M, 'S', ISEED, 'N', D,
$ 6, ONE, ONE, 'T', 'N',
$ DUML, 1, ONE, DUMR, 1, ONE,
$ 'N', IWORK, KLA, KUA, ZERO,
$ ONE, 'NO', A, MAXM, IWORK, IINFO )
DO I = 1, M
A( I, I ) = A( I, I ) * VM( J )
END DO
ANRM = SLANGE( 'M', M, M, A, MAXM, DUM )
DO N = 51, MAXN, 47
KLB = 0
KUB = N - 1
CALL SLATMR( N, N, 'S', ISEED, 'N', D,
$ 6, ONE, ONE, 'T', 'N',
$ DUML, 1, ONE, DUMR, 1, ONE,
$ 'N', IWORK, KLB, KUB, ZERO,
$ ONE, 'NO', B, MAXN, IWORK, IINFO )
BNRM = SLANGE( 'M', N, N, B, MAXN, DUM )
TNRM = MAX( ANRM, BNRM )
CALL SLATMR( M, N, 'S', ISEED, 'N', D,
$ 6, ONE, ONE, 'T', 'N',
$ DUML, 1, ONE, DUMR, 1, ONE,
$ 'N', IWORK, M, N, ZERO, ONE,
$ 'NO', C, MAXM, IWORK, IINFO )
DO ITRANA = 1, 2
IF( ITRANA.EQ.1 ) THEN
TRANA = 'N'
END IF
IF( ITRANA.EQ.2 ) THEN
TRANA = 'T'
END IF
DO ITRANB = 1, 2
IF( ITRANB.EQ.1 ) THEN
TRANB = 'N'
END IF
IF( ITRANB.EQ.2 ) THEN
TRANB = 'T'
END IF
KNT = KNT + 1
*
CALL SLACPY( 'All', M, N, C, MAXM, X, MAXM)
CALL SLACPY( 'All', M, N, C, MAXM, CC, MAXM)
CALL STRSYL( TRANA, TRANB, ISGN, M, N,
$ A, MAXM, B, MAXN, X, MAXM,
$ SCALE, IINFO )
IF( IINFO.NE.0 )
$ NINFO( 1 ) = NINFO( 1 ) + 1
XNRM = SLANGE( 'M', M, N, X, MAXM, DUM )
RMUL = ONE
IF( XNRM.GT.ONE .AND. TNRM.GT.ONE ) THEN
IF( XNRM.GT.BIGNUM / TNRM ) THEN
RMUL = ONE / MAX( XNRM, TNRM )
END IF
END IF
CALL SGEMM( TRANA, 'N', M, N, M, RMUL,
$ A, MAXM, X, MAXM, -SCALE*RMUL,
$ C, MAXM )
CALL SGEMM( 'N', TRANB, M, N, N,
$ REAL( ISGN )*RMUL, X, MAXM, B,
$ MAXN, ONE, C, MAXM )
RES1 = SLANGE( 'M', M, N, C, MAXM, DUM )
RES = RES1 / MAX( SMLNUM, SMLNUM*XNRM,
$ ( ( RMUL*TNRM )*EPS )*XNRM )
IF( RES.GT.THRESH )
$ NFAIL( 1 ) = NFAIL( 1 ) + 1
IF( RES.GT.RMAX( 1 ) )
$ RMAX( 1 ) = RES
*
CALL SLACPY( 'All', M, N, C, MAXM, X, MAXM )
CALL SLACPY( 'All', M, N, C, MAXM, CC, MAXM )
CALL STRSYL3( TRANA, TRANB, ISGN, M, N,
$ A, MAXM, B, MAXN, X, MAXM,
$ SCALE3, IWORK, LIWORK,
$ SWORK, LDSWORK, INFO)
IF( INFO.NE.0 )
$ NINFO( 2 ) = NINFO( 2 ) + 1
XNRM = SLANGE( 'M', M, N, X, MAXM, DUM )
RMUL = ONE
IF( XNRM.GT.ONE .AND. TNRM.GT.ONE ) THEN
IF( XNRM.GT.BIGNUM / TNRM ) THEN
RMUL = ONE / MAX( XNRM, TNRM )
END IF
END IF
CALL SGEMM( TRANA, 'N', M, N, M, RMUL,
$ A, MAXM, X, MAXM, -SCALE3*RMUL,
$ CC, MAXM )
CALL SGEMM( 'N', TRANB, M, N, N,
$ REAL( ISGN )*RMUL, X, MAXM, B,
$ MAXN, ONE, CC, MAXM )
RES1 = SLANGE( 'M', M, N, CC, MAXM, DUM )
RES = RES1 / MAX( SMLNUM, SMLNUM*XNRM,
$ ( ( RMUL*TNRM )*EPS )*XNRM )
* Verify that TRSYL3 only flushes if TRSYL flushes (but
* there may be cases where TRSYL3 avoid flushing).
IF( SCALE3.EQ.ZERO .AND. SCALE.GT.ZERO .OR.
$ IINFO.NE.INFO ) THEN
NFAIL( 3 ) = NFAIL( 3 ) + 1
END IF
IF( RES.GT.THRESH .OR. SISNAN( RES ) )
$ NFAIL( 2 ) = NFAIL( 2 ) + 1
IF( RES.GT.RMAX( 2 ) )
$ RMAX( 2 ) = RES
END DO
END DO
END DO
END DO
END DO
END DO
*
RETURN
*
* End of SSYL01
*
END

View File

@ -88,17 +88,17 @@
* .. Local Scalars .. * .. Local Scalars ..
LOGICAL OK LOGICAL OK
CHARACTER*3 PATH CHARACTER*3 PATH
INTEGER KTREXC, KTRSEN, KTRSNA, KTRSYL, LTREXC, LTRSYL, INTEGER KTREXC, KTRSEN, KTRSNA, KTRSYL, KTRSYL3,
$ NTESTS, NTREXC, NTRSYL $ LTREXC, LTRSYL, NTESTS, NTREXC, NTRSYL
DOUBLE PRECISION EPS, RTREXC, RTRSYL, SFMIN DOUBLE PRECISION EPS, RTREXC, SFMIN
* .. * ..
* .. Local Arrays .. * .. Local Arrays ..
INTEGER LTRSEN( 3 ), LTRSNA( 3 ), NTRSEN( 3 ), INTEGER FTRSYL( 3 ), ITRSYL( 2 ), LTRSEN( 3 ),
$ NTRSNA( 3 ) $ LTRSNA( 3 ), NTRSEN( 3 ), NTRSNA( 3 )
DOUBLE PRECISION RTRSEN( 3 ), RTRSNA( 3 ) DOUBLE PRECISION RTRSEN( 3 ), RTRSNA( 3 ), RTRSYL( 2 )
* .. * ..
* .. External Subroutines .. * .. External Subroutines ..
EXTERNAL ZERREC, ZGET35, ZGET36, ZGET37, ZGET38 EXTERNAL ZERREC, ZGET35, ZGET36, ZGET37, ZGET38, ZSYL01
* .. * ..
* .. External Functions .. * .. External Functions ..
DOUBLE PRECISION DLAMCH DOUBLE PRECISION DLAMCH
@ -120,10 +120,24 @@
$ CALL ZERREC( PATH, NOUT ) $ CALL ZERREC( PATH, NOUT )
* *
OK = .TRUE. OK = .TRUE.
CALL ZGET35( RTRSYL, LTRSYL, NTRSYL, KTRSYL, NIN ) CALL ZGET35( RTRSYL( 1 ), LTRSYL, NTRSYL, KTRSYL, NIN )
IF( RTRSYL.GT.THRESH ) THEN IF( RTRSYL( 1 ).GT.THRESH ) THEN
OK = .FALSE. OK = .FALSE.
WRITE( NOUT, FMT = 9999 )RTRSYL, LTRSYL, NTRSYL, KTRSYL WRITE( NOUT, FMT = 9999 )RTRSYL( 1 ), LTRSYL, NTRSYL, KTRSYL
END IF
*
CALL ZSYL01( THRESH, FTRSYL, RTRSYL, ITRSYL, KTRSYL3 )
IF( FTRSYL( 1 ).GT.0 ) THEN
OK = .FALSE.
WRITE( NOUT, FMT = 9970 )FTRSYL( 1 ), RTRSYL( 1 ), THRESH
END IF
IF( FTRSYL( 2 ).GT.0 ) THEN
OK = .FALSE.
WRITE( NOUT, FMT = 9971 )FTRSYL( 2 ), RTRSYL( 2 ), THRESH
END IF
IF( FTRSYL( 3 ).GT.0 ) THEN
OK = .FALSE.
WRITE( NOUT, FMT = 9972 )FTRSYL( 3 )
END IF END IF
* *
CALL ZGET36( RTREXC, LTREXC, NTREXC, KTREXC, NIN ) CALL ZGET36( RTREXC, LTREXC, NTREXC, KTREXC, NIN )
@ -148,7 +162,7 @@
WRITE( NOUT, FMT = 9996 )RTRSEN, LTRSEN, NTRSEN, KTRSEN WRITE( NOUT, FMT = 9996 )RTRSEN, LTRSEN, NTRSEN, KTRSEN
END IF END IF
* *
NTESTS = KTRSYL + KTREXC + KTRSNA + KTRSEN NTESTS = KTRSYL + KTRSYL3 + KTREXC + KTRSNA + KTRSEN
IF( OK ) IF( OK )
$ WRITE( NOUT, FMT = 9995 )PATH, NTESTS $ WRITE( NOUT, FMT = 9995 )PATH, NTESTS
* *
@ -169,6 +183,12 @@
$ / ' Safe minimum (SFMIN) = ', D16.6, / ) $ / ' Safe minimum (SFMIN) = ', D16.6, / )
9992 FORMAT( ' Routines pass computational tests if test ratio is ', 9992 FORMAT( ' Routines pass computational tests if test ratio is ',
$ 'less than', F8.2, / / ) $ 'less than', F8.2, / / )
9970 FORMAT( 'Error in ZTRSYL: ', I8, ' tests fail the threshold.', /
$ 'Maximum test ratio =', D12.3, ' threshold =', D12.3 )
9971 FORMAT( 'Error in ZTRSYL3: ', I8, ' tests fail the threshold.', /
$ 'Maximum test ratio =', D12.3, ' threshold =', D12.3 )
9972 FORMAT( 'ZTRSYL and ZTRSYL3 compute an inconsistent scale ',
$ 'factor in ', I8, ' tests.')
RETURN RETURN
* *
* End of ZCHKEC * End of ZCHKEC

View File

@ -23,7 +23,7 @@
*> *>
*> ZERREC tests the error exits for the routines for eigen- condition *> ZERREC tests the error exits for the routines for eigen- condition
*> estimation for DOUBLE PRECISION matrices: *> estimation for DOUBLE PRECISION matrices:
*> ZTRSYL, ZTREXC, ZTRSNA and ZTRSEN. *> ZTRSYL, ZTRSYL3, ZTREXC, ZTRSNA and ZTRSEN.
*> \endverbatim *> \endverbatim
* *
* Arguments: * Arguments:
@ -77,7 +77,7 @@
* .. * ..
* .. Local Arrays .. * .. Local Arrays ..
LOGICAL SEL( NMAX ) LOGICAL SEL( NMAX )
DOUBLE PRECISION RW( LW ), S( NMAX ), SEP( NMAX ) DOUBLE PRECISION RW( LW ), S( NMAX ), SEP( NMAX ), SWORK( NMAX )
COMPLEX*16 A( NMAX, NMAX ), B( NMAX, NMAX ), COMPLEX*16 A( NMAX, NMAX ), B( NMAX, NMAX ),
$ C( NMAX, NMAX ), WORK( LW ), X( NMAX ) $ C( NMAX, NMAX ), WORK( LW ), X( NMAX )
* .. * ..
@ -141,6 +141,43 @@
CALL CHKXER( 'ZTRSYL', INFOT, NOUT, LERR, OK ) CALL CHKXER( 'ZTRSYL', INFOT, NOUT, LERR, OK )
NT = NT + 8 NT = NT + 8
* *
* Test ZTRSYL3
*
SRNAMT = 'ZTRSYL3'
INFOT = 1
CALL ZTRSYL3( 'X', 'N', 1, 0, 0, A, 1, B, 1, C, 1, SCALE,
$ SWORK, NMAX, INFO )
CALL CHKXER( 'ZTRSYL3', INFOT, NOUT, LERR, OK )
INFOT = 2
CALL ZTRSYL3( 'N', 'X', 1, 0, 0, A, 1, B, 1, C, 1, SCALE,
$ SWORK, NMAX, INFO )
CALL CHKXER( 'ZTRSYL3', INFOT, NOUT, LERR, OK )
INFOT = 3
CALL ZTRSYL3( 'N', 'N', 0, 0, 0, A, 1, B, 1, C, 1, SCALE,
$ SWORK, NMAX, INFO )
CALL CHKXER( 'ZTRSYL3', INFOT, NOUT, LERR, OK )
INFOT = 4
CALL ZTRSYL3( 'N', 'N', 1, -1, 0, A, 1, B, 1, C, 1, SCALE,
$ SWORK, NMAX, INFO )
CALL CHKXER( 'ZTRSYL3', INFOT, NOUT, LERR, OK )
INFOT = 5
CALL ZTRSYL3( 'N', 'N', 1, 0, -1, A, 1, B, 1, C, 1, SCALE,
$ SWORK, NMAX, INFO )
CALL CHKXER( 'ZTRSYL3', INFOT, NOUT, LERR, OK )
INFOT = 7
CALL ZTRSYL3( 'N', 'N', 1, 2, 0, A, 1, B, 1, C, 2, SCALE,
$ SWORK, NMAX, INFO )
CALL CHKXER( 'ZTRSYL3', INFOT, NOUT, LERR, OK )
INFOT = 9
CALL ZTRSYL3( 'N', 'N', 1, 0, 2, A, 1, B, 1, C, 1, SCALE,
$ SWORK, NMAX, INFO )
CALL CHKXER( 'ZTRSYL3', INFOT, NOUT, LERR, OK )
INFOT = 11
CALL ZTRSYL3( 'N', 'N', 1, 2, 0, A, 2, B, 1, C, 1, SCALE,
$ SWORK, NMAX, INFO )
CALL CHKXER( 'ZTRSYL3', INFOT, NOUT, LERR, OK )
NT = NT + 8
*
* Test ZTREXC * Test ZTREXC
* *
SRNAMT = 'ZTREXC' SRNAMT = 'ZTREXC'

View File

@ -0,0 +1,294 @@
*> \brief \b ZSYL01
*
* =========== DOCUMENTATION ===========
*
* Online html documentation available at
* http://www.netlib.org/lapack/explore-html/
*
* Definition:
* ===========
*
* SUBROUTINE ZSYL01( THRESH, NFAIL, RMAX, NINFO, KNT )
*
* .. Scalar Arguments ..
* INTEGER KNT
* DOUBLE PRECISION THRESH
* ..
* .. Array Arguments ..
* INTEGER NFAIL( 3 ), NINFO( 2 )
* DOUBLE PRECISION RMAX( 2 )
* ..
*
*
*> \par Purpose:
* =============
*>
*> \verbatim
*>
*> ZSYL01 tests ZTRSYL and ZTRSYL3, routines for solving the Sylvester matrix
*> equation
*>
*> op(A)*X + ISGN*X*op(B) = scale*C,
*>
*> where op(A) and op(B) are both upper triangular form, op() represents an
*> optional conjugate transpose, and ISGN can be -1 or +1. Scale is an output
*> less than or equal to 1, chosen to avoid overflow in X.
*>
*> The test code verifies that the following residual does not exceed
*> the provided threshold:
*>
*> norm(op(A)*X + ISGN*X*op(B) - scale*C) /
*> (EPS*max(norm(A),norm(B))*norm(X))
*>
*> This routine complements ZGET35 by testing with larger,
*> random matrices, of which some require rescaling of X to avoid overflow.
*>
*> \endverbatim
*
* Arguments:
* ==========
*
*> \param[in] THRESH
*> \verbatim
*> THRESH is DOUBLE PRECISION
*> A test will count as "failed" if the residual, computed as
*> described above, exceeds THRESH.
*> \endverbatim
*>
*> \param[out] NFAIL
*> \verbatim
*> NFAIL is INTEGER array, dimension (3)
*> NFAIL(1) = No. of times residual ZTRSYL exceeds threshold THRESH
*> NFAIL(2) = No. of times residual ZTRSYL3 exceeds threshold THRESH
*> NFAIL(3) = No. of times ZTRSYL3 and ZTRSYL deviate
*> \endverbatim
*>
*> \param[out] RMAX
*> \verbatim
*> RMAX is DOUBLE PRECISION array, dimension (2)
*> RMAX(1) = Value of the largest test ratio of ZTRSYL
*> RMAX(2) = Value of the largest test ratio of ZTRSYL3
*> \endverbatim
*>
*> \param[out] NINFO
*> \verbatim
*> NINFO is INTEGER array, dimension (2)
*> NINFO(1) = No. of times ZTRSYL returns an expected INFO
*> NINFO(2) = No. of times ZTRSYL3 returns an expected INFO
*> \endverbatim
*>
*> \param[out] KNT
*> \verbatim
*> KNT is INTEGER
*> Total number of examples tested.
*> \endverbatim
*
* -- LAPACK test routine --
SUBROUTINE ZSYL01( THRESH, NFAIL, RMAX, NINFO, KNT )
IMPLICIT NONE
*
* -- LAPACK test routine --
* -- LAPACK is a software package provided by Univ. of Tennessee, --
* -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
*
* .. Scalar Arguments ..
INTEGER KNT
DOUBLE PRECISION THRESH
* ..
* .. Array Arguments ..
INTEGER NFAIL( 3 ), NINFO( 2 )
DOUBLE PRECISION RMAX( 2 )
* ..
*
* =====================================================================
* ..
* .. Parameters ..
COMPLEX*16 CONE
PARAMETER ( CONE = ( 1.0D0, 0.0D+0 ) )
DOUBLE PRECISION ONE, ZERO
PARAMETER ( ZERO = 0.0D+0, ONE = 1.0D+0 )
INTEGER MAXM, MAXN, LDSWORK
PARAMETER ( MAXM = 185, MAXN = 192, LDSWORK = 36 )
* ..
* .. Local Scalars ..
CHARACTER TRANA, TRANB
INTEGER I, INFO, IINFO, ISGN, ITRANA, ITRANB, J, KLA,
$ KUA, KLB, KUB, M, N
DOUBLE PRECISION ANRM, BNRM, BIGNUM, EPS, RES, RES1,
$ SCALE, SCALE3, SMLNUM, TNRM, XNRM
COMPLEX*16 RMUL
* ..
* .. Local Arrays ..
COMPLEX*16 A( MAXM, MAXM ), B( MAXN, MAXN ),
$ C( MAXM, MAXN ), CC( MAXM, MAXN ),
$ X( MAXM, MAXN ),
$ DUML( MAXM ), DUMR( MAXN ),
$ D( MIN( MAXM, MAXN ) )
DOUBLE PRECISION SWORK( LDSWORK, 103 ), DUM( MAXN ), VM( 2 )
INTEGER ISEED( 4 ), IWORK( MAXM + MAXN + 2 )
* ..
* .. External Functions ..
LOGICAL DISNAN
DOUBLE PRECISION DLAMCH, ZLANGE
EXTERNAL DISNAN, DLAMCH, ZLANGE
* ..
* .. External Subroutines ..
EXTERNAL ZLATMR, ZLACPY, ZGEMM, ZTRSYL, ZTRSYL3
* ..
* .. Intrinsic Functions ..
INTRINSIC ABS, DBLE, MAX, SQRT
* ..
* .. Executable Statements ..
*
* Get machine parameters
*
EPS = DLAMCH( 'P' )
SMLNUM = DLAMCH( 'S' ) / EPS
BIGNUM = ONE / SMLNUM
*
* Expect INFO = 0
VM( 1 ) = ONE
* Expect INFO = 1
VM( 2 ) = 0.05D+0
*
* Begin test loop
*
NINFO( 1 ) = 0
NINFO( 2 ) = 0
NFAIL( 1 ) = 0
NFAIL( 2 ) = 0
NFAIL( 3 ) = 0
RMAX( 1 ) = ZERO
RMAX( 2 ) = ZERO
KNT = 0
ISEED( 1 ) = 1
ISEED( 2 ) = 1
ISEED( 3 ) = 1
ISEED( 4 ) = 1
SCALE = ONE
SCALE3 = ONE
DO J = 1, 2
DO ISGN = -1, 1, 2
* Reset seed (overwritten by LATMR)
ISEED( 1 ) = 1
ISEED( 2 ) = 1
ISEED( 3 ) = 1
ISEED( 4 ) = 1
DO M = 32, MAXM, 51
KLA = 0
KUA = M - 1
CALL ZLATMR( M, M, 'S', ISEED, 'N', D,
$ 6, ONE, CONE, 'T', 'N',
$ DUML, 1, ONE, DUMR, 1, ONE,
$ 'N', IWORK, KLA, KUA, ZERO,
$ ONE, 'NO', A, MAXM, IWORK,
$ IINFO )
DO I = 1, M
A( I, I ) = A( I, I ) * VM( J )
END DO
ANRM = ZLANGE( 'M', M, M, A, MAXM, DUM )
DO N = 51, MAXN, 47
KLB = 0
KUB = N - 1
CALL ZLATMR( N, N, 'S', ISEED, 'N', D,
$ 6, ONE, CONE, 'T', 'N',
$ DUML, 1, ONE, DUMR, 1, ONE,
$ 'N', IWORK, KLB, KUB, ZERO,
$ ONE, 'NO', B, MAXN, IWORK,
$ IINFO )
DO I = 1, N
B( I, I ) = B( I, I ) * VM ( J )
END DO
BNRM = ZLANGE( 'M', N, N, B, MAXN, DUM )
TNRM = MAX( ANRM, BNRM )
CALL ZLATMR( M, N, 'S', ISEED, 'N', D,
$ 6, ONE, CONE, 'T', 'N',
$ DUML, 1, ONE, DUMR, 1, ONE,
$ 'N', IWORK, M, N, ZERO, ONE,
$ 'NO', C, MAXM, IWORK, IINFO )
DO ITRANA = 1, 2
IF( ITRANA.EQ.1 )
$ TRANA = 'N'
IF( ITRANA.EQ.2 )
$ TRANA = 'C'
DO ITRANB = 1, 2
IF( ITRANB.EQ.1 )
$ TRANB = 'N'
IF( ITRANB.EQ.2 )
$ TRANB = 'C'
KNT = KNT + 1
*
CALL ZLACPY( 'All', M, N, C, MAXM, X, MAXM)
CALL ZLACPY( 'All', M, N, C, MAXM, CC, MAXM)
CALL ZTRSYL( TRANA, TRANB, ISGN, M, N,
$ A, MAXM, B, MAXN, X, MAXM,
$ SCALE, IINFO )
IF( IINFO.NE.0 )
$ NINFO( 1 ) = NINFO( 1 ) + 1
XNRM = ZLANGE( 'M', M, N, X, MAXM, DUM )
RMUL = CONE
IF( XNRM.GT.ONE .AND. TNRM.GT.ONE ) THEN
IF( XNRM.GT.BIGNUM / TNRM ) THEN
RMUL = CONE / MAX( XNRM, TNRM )
END IF
END IF
CALL ZGEMM( TRANA, 'N', M, N, M, RMUL,
$ A, MAXM, X, MAXM, -SCALE*RMUL,
$ CC, MAXM )
CALL ZGEMM( 'N', TRANB, M, N, N,
$ DBLE( ISGN )*RMUL, X, MAXM, B,
$ MAXN, CONE, CC, MAXM )
RES1 = ZLANGE( 'M', M, N, CC, MAXM, DUM )
RES = RES1 / MAX( SMLNUM, SMLNUM*XNRM,
$ ( ( ABS( RMUL )*TNRM )*EPS )*XNRM )
IF( RES.GT.THRESH )
$ NFAIL( 1 ) = NFAIL( 1 ) + 1
IF( RES.GT.RMAX( 1 ) )
$ RMAX( 1 ) = RES
*
CALL ZLACPY( 'All', M, N, C, MAXM, X, MAXM )
CALL ZLACPY( 'All', M, N, C, MAXM, CC, MAXM )
CALL ZTRSYL3( TRANA, TRANB, ISGN, M, N,
$ A, MAXM, B, MAXN, X, MAXM,
$ SCALE3, SWORK, LDSWORK, INFO)
IF( INFO.NE.0 )
$ NINFO( 2 ) = NINFO( 2 ) + 1
XNRM = ZLANGE( 'M', M, N, X, MAXM, DUM )
RMUL = CONE
IF( XNRM.GT.ONE .AND. TNRM.GT.ONE ) THEN
IF( XNRM.GT.BIGNUM / TNRM ) THEN
RMUL = CONE / MAX( XNRM, TNRM )
END IF
END IF
CALL ZGEMM( TRANA, 'N', M, N, M, RMUL,
$ A, MAXM, X, MAXM, -SCALE3*RMUL,
$ CC, MAXM )
CALL ZGEMM( 'N', TRANB, M, N, N,
$ DBLE( ISGN )*RMUL, X, MAXM, B,
$ MAXN, CONE, CC, MAXM )
RES1 = ZLANGE( 'M', M, N, CC, MAXM, DUM )
RES = RES1 / MAX( SMLNUM, SMLNUM*XNRM,
$ ( ( ABS( RMUL )*TNRM )*EPS )*XNRM )
* Verify that TRSYL3 only flushes if TRSYL flushes (but
* there may be cases where TRSYL3 avoid flushing).
IF( SCALE3.EQ.ZERO .AND. SCALE.GT.ZERO .OR.
$ IINFO.NE.INFO ) THEN
NFAIL( 3 ) = NFAIL( 3 ) + 1
END IF
IF( RES.GT.THRESH .OR. DISNAN( RES ) )
$ NFAIL( 2 ) = NFAIL( 2 ) + 1
IF( RES.GT.RMAX( 2 ) )
$ RMAX( 2 ) = RES
END DO
END DO
END DO
END DO
END DO
END DO
*
RETURN
*
* End of ZSYL01
*
END