Implement truncated QR with pivoting (Reference-LAPACK PR 891)

This commit is contained in:
Martin Kroeker 2023-11-15 09:53:06 +01:00 committed by GitHub
parent 40109c0392
commit 387830b9d5
No known key found for this signature in database
GPG Key ID: 4AEE18F83AFDEB23
30 changed files with 4132 additions and 238 deletions

View File

@ -9,7 +9,7 @@ set(DZLNTST dlaord.f)
set(SLINTST schkaa.F set(SLINTST schkaa.F
schkeq.f schkgb.f schkge.f schkgt.f schkeq.f schkgb.f schkge.f schkgt.f
schklq.f schkpb.f schkpo.f schkps.f schkpp.f schklq.f schkpb.f schkpo.f schkps.f schkpp.f
schkpt.f schkq3.f schkql.f schkqr.f schkrq.f schkpt.f schkq3.f schkqp3rk.f schkql.f schkqr.f schkrq.f
schksp.f schksy.f schksy_rook.f schksy_rk.f schksp.f schksy.f schksy_rook.f schksy_rk.f
schksy_aa.f schksy_aa_2stage.f schksy_aa.f schksy_aa_2stage.f
schktb.f schktp.f schktr.f schktb.f schktp.f schktr.f
@ -56,7 +56,7 @@ set(CLINTST cchkaa.F
cchkhe.f cchkhe_rook.f cchkhe_rk.f cchkhe.f cchkhe_rook.f cchkhe_rk.f
cchkhe_aa.f cchkhe_aa_2stage.f cchkhe_aa.f cchkhe_aa_2stage.f
cchkhp.f cchklq.f cchkpb.f cchkhp.f cchklq.f cchkpb.f
cchkpo.f cchkps.f cchkpp.f cchkpt.f cchkq3.f cchkql.f cchkpo.f cchkps.f cchkpp.f cchkpt.f cchkq3.f cchkqp3rk.f cchkql.f
cchkqr.f cchkrq.f cchksp.f cchksy.f cchksy_rook.f cchksy_rk.f cchkqr.f cchkrq.f cchksp.f cchksy.f cchksy_rook.f cchksy_rk.f
cchksy_aa.f cchksy_aa_2stage.f cchksy_aa.f cchksy_aa_2stage.f
cchktb.f cchktb.f
@ -110,7 +110,7 @@ endif()
set(DLINTST dchkaa.F set(DLINTST dchkaa.F
dchkeq.f dchkgb.f dchkge.f dchkgt.f dchkeq.f dchkgb.f dchkge.f dchkgt.f
dchklq.f dchkpb.f dchkpo.f dchkps.f dchkpp.f dchklq.f dchkpb.f dchkpo.f dchkps.f dchkpp.f
dchkpt.f dchkq3.f dchkql.f dchkqr.f dchkrq.f dchkpt.f dchkq3.f dchkqp3rk.f dchkql.f dchkqr.f dchkrq.f
dchksp.f dchksy.f dchksy_rook.f dchksy_rk.f dchksp.f dchksy.f dchksy_rook.f dchksy_rk.f
dchksy_aa.f dchksy_aa_2stage.f dchksy_aa.f dchksy_aa_2stage.f
dchktb.f dchktp.f dchktr.f dchktb.f dchktp.f dchktr.f
@ -158,7 +158,7 @@ set(ZLINTST zchkaa.F
zchkhe.f zchkhe_rook.f zchkhe_rk.f zchkhe.f zchkhe_rook.f zchkhe_rk.f
zchkhe_aa.f zchkhe_aa_2stage.f zchkhe_aa.f zchkhe_aa_2stage.f
zchkhp.f zchklq.f zchkpb.f zchkhp.f zchklq.f zchkpb.f
zchkpo.f zchkps.f zchkpp.f zchkpt.f zchkq3.f zchkql.f zchkpo.f zchkps.f zchkpp.f zchkpt.f zchkq3.f zchkqp3rk.f zchkql.f
zchkqr.f zchkrq.f zchksp.f zchksy.f zchksy_rook.f zchksy_rk.f zchkqr.f zchkrq.f zchksp.f zchksy.f zchksy_rook.f zchksy_rk.f
zchksy_aa.f zchksy_aa_2stage.f zchksy_aa.f zchksy_aa_2stage.f
zchktb.f zchktb.f
@ -239,8 +239,7 @@ set(ZLINTSTRFP zchkrfp.f zdrvrfp.f zdrvrf1.f zdrvrf2.f zdrvrf3.f zdrvrf4.f zerrr
macro(add_lin_executable name) macro(add_lin_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

@ -45,7 +45,7 @@ DZLNTST = dlaord.o
SLINTST = schkaa.o \ SLINTST = schkaa.o \
schkeq.o schkgb.o schkge.o schkgt.o \ schkeq.o schkgb.o schkge.o schkgt.o \
schklq.o schkpb.o schkpo.o schkps.o schkpp.o \ schklq.o schkpb.o schkpo.o schkps.o schkpp.o \
schkpt.o schkq3.o schkql.o schkqr.o schkrq.o \ schkpt.o schkq3.o schkqp3rk.o schkql.o schkqr.o schkrq.o \
schksp.o schksy.o schksy_rook.o schksy_rk.o \ schksp.o schksy.o schksy_rook.o schksy_rk.o \
schksy_aa.o schksy_aa_2stage.o schktb.o schktp.o schktr.o \ schksy_aa.o schksy_aa_2stage.o schktb.o schktp.o schktr.o \
schktz.o \ schktz.o \
@ -89,7 +89,7 @@ CLINTST = cchkaa.o \
cchkeq.o cchkgb.o cchkge.o cchkgt.o \ cchkeq.o cchkgb.o cchkge.o cchkgt.o \
cchkhe.o cchkhe_rook.o cchkhe_rk.o \ cchkhe.o cchkhe_rook.o cchkhe_rk.o \
cchkhe_aa.o cchkhe_aa_2stage.o cchkhp.o cchklq.o cchkpb.o \ cchkhe_aa.o cchkhe_aa_2stage.o cchkhp.o cchklq.o cchkpb.o \
cchkpo.o cchkps.o cchkpp.o cchkpt.o cchkq3.o cchkql.o \ cchkpo.o cchkps.o cchkpp.o cchkpt.o cchkq3.o cchkqp3rk.o cchkql.o \
cchkqr.o cchkrq.o cchksp.o cchksy.o cchksy_rook.o cchksy_rk.o \ cchkqr.o cchkrq.o cchksp.o cchksy.o cchksy_rook.o cchksy_rk.o \
cchksy_aa.o cchksy_aa_2stage.o cchktb.o \ cchksy_aa.o cchksy_aa_2stage.o cchktb.o \
cchktp.o cchktr.o cchktz.o \ cchktp.o cchktr.o cchktz.o \
@ -137,7 +137,7 @@ endif
DLINTST = dchkaa.o \ DLINTST = dchkaa.o \
dchkeq.o dchkgb.o dchkge.o dchkgt.o \ dchkeq.o dchkgb.o dchkge.o dchkgt.o \
dchklq.o dchkpb.o dchkpo.o dchkps.o dchkpp.o \ dchklq.o dchkpb.o dchkpo.o dchkps.o dchkpp.o \
dchkpt.o dchkq3.o dchkql.o dchkqr.o dchkrq.o \ dchkpt.o dchkq3.o dchkqp3rk.o dchkql.o dchkqr.o dchkrq.o \
dchksp.o dchksy.o dchksy_rook.o dchksy_rk.o \ dchksp.o dchksy.o dchksy_rook.o dchksy_rk.o \
dchksy_aa.o dchksy_aa_2stage.o dchktb.o dchktp.o dchktr.o \ dchksy_aa.o dchksy_aa_2stage.o dchktb.o dchktp.o dchktr.o \
dchktz.o \ dchktz.o \
@ -182,7 +182,7 @@ ZLINTST = zchkaa.o \
zchkeq.o zchkgb.o zchkge.o zchkgt.o \ zchkeq.o zchkgb.o zchkge.o zchkgt.o \
zchkhe.o zchkhe_rook.o zchkhe_rk.o zchkhe_aa.o zchkhe_aa_2stage.o \ zchkhe.o zchkhe_rook.o zchkhe_rk.o zchkhe_aa.o zchkhe_aa_2stage.o \
zchkhp.o zchklq.o zchkpb.o \ zchkhp.o zchklq.o zchkpb.o \
zchkpo.o zchkps.o zchkpp.o zchkpt.o zchkq3.o zchkql.o \ zchkpo.o zchkps.o zchkpp.o zchkpt.o zchkq3.o zchkqp3rk.o zchkql.o \
zchkqr.o zchkrq.o zchksp.o zchksy.o zchksy_rook.o zchksy_rk.o \ zchkqr.o zchkrq.o zchksp.o zchksy.o zchksy_rook.o zchksy_rk.o \
zchksy_aa.o zchksy_aa_2stage.o zchktb.o \ zchksy_aa.o zchksy_aa_2stage.o zchktb.o \
zchktp.o zchktr.o zchktz.o \ zchktp.o zchktr.o zchktz.o \
@ -269,35 +269,35 @@ proto-double: xlintstds xlintstrfd
proto-complex: xlintstrfc proto-complex: xlintstrfc
proto-complex16: xlintstzc xlintstrfz proto-complex16: xlintstzc xlintstrfz
xlintsts: $(ALINTST) $(SLINTST) $(SCLNTST) $(TMGLIB) $(VARLIB) ../$(LAPACKLIB) $(XBLASLIB) $(BLASLIB) xlintsts: $(ALINTST) $(SLINTST) $(SCLNTST) $(TMGLIB) $(VARLIB) $(LAPACKLIB) $(XBLASLIB) $(BLASLIB)
$(LOADER) $(FFLAGS) $(LDFLAGS) -o $@ $^ $(FC) $(FFLAGS) $(LDFLAGS) -o $@ $^
xlintstc: $(ALINTST) $(CLINTST) $(SCLNTST) $(TMGLIB) $(VARLIB) ../$(LAPACKLIB) $(XBLASLIB) $(BLASLIB) xlintstc: $(ALINTST) $(CLINTST) $(SCLNTST) $(TMGLIB) $(VARLIB) $(LAPACKLIB) $(XBLASLIB) $(BLASLIB)
$(LOADER) $(FFLAGS) $(LDFLAGS) -o $@ $^ $(FC) $(FFLAGS) $(LDFLAGS) -o $@ $^
xlintstd: $(ALINTST) $(DLINTST) $(DZLNTST) $(TMGLIB) $(VARLIB) ../$(LAPACKLIB) $(XBLASLIB) $(BLASLIB) xlintstd: $(ALINTST) $(DLINTST) $(DZLNTST) $(TMGLIB) $(VARLIB) $(LAPACKLIB) $(XBLASLIB) $(BLASLIB)
$(LOADER) $(FFLAGS) $(LDFLAGS) -o $@ $^ $(FC) $(FFLAGS) $(LDFLAGS) -o $@ $^
xlintstz: $(ALINTST) $(ZLINTST) $(DZLNTST) $(TMGLIB) $(VARLIB) ../$(LAPACKLIB) $(XBLASLIB) $(BLASLIB) xlintstz: $(ALINTST) $(ZLINTST) $(DZLNTST) $(TMGLIB) $(VARLIB) $(LAPACKLIB) $(XBLASLIB) $(BLASLIB)
$(LOADER) $(FFLAGS) $(LDFLAGS) -o $@ $^ $(FC) $(FFLAGS) $(LDFLAGS) -o $@ $^
xlintstds: $(DSLINTST) $(TMGLIB) $(VARLIB) ../$(LAPACKLIB) $(BLASLIB) xlintstds: $(DSLINTST) $(TMGLIB) $(VARLIB) $(LAPACKLIB) $(BLASLIB)
$(LOADER) $(FFLAGS) $(LDFLAGS) -o $@ $^ $(FC) $(FFLAGS) $(LDFLAGS) -o $@ $^
xlintstzc: $(ZCLINTST) $(TMGLIB) $(VARLIB) ../$(LAPACKLIB) $(BLASLIB) xlintstzc: $(ZCLINTST) $(TMGLIB) $(VARLIB) $(LAPACKLIB) $(BLASLIB)
$(LOADER) $(FFLAGS) $(LDFLAGS) -o $@ $^ $(FC) $(FFLAGS) $(LDFLAGS) -o $@ $^
xlintstrfs: $(SLINTSTRFP) $(TMGLIB) $(VARLIB) ../$(LAPACKLIB) $(BLASLIB) xlintstrfs: $(SLINTSTRFP) $(TMGLIB) $(VARLIB) $(LAPACKLIB) $(BLASLIB)
$(LOADER) $(FFLAGS) $(LDFLAGS) -o $@ $^ $(FC) $(FFLAGS) $(LDFLAGS) -o $@ $^
xlintstrfd: $(DLINTSTRFP) $(TMGLIB) $(VARLIB) ../$(LAPACKLIB) $(BLASLIB) xlintstrfd: $(DLINTSTRFP) $(TMGLIB) $(VARLIB) $(LAPACKLIB) $(BLASLIB)
$(LOADER) $(FFLAGS) $(LDFLAGS) -o $@ $^ $(FC) $(FFLAGS) $(LDFLAGS) -o $@ $^
xlintstrfc: $(CLINTSTRFP) $(TMGLIB) $(VARLIB) ../$(LAPACKLIB) $(BLASLIB) xlintstrfc: $(CLINTSTRFP) $(TMGLIB) $(VARLIB) $(LAPACKLIB) $(BLASLIB)
$(LOADER) $(FFLAGS) $(LDFLAGS) -o $@ $^ $(FC) $(FFLAGS) $(LDFLAGS) -o $@ $^
xlintstrfz: $(ZLINTSTRFP) $(TMGLIB) $(VARLIB) ../$(LAPACKLIB) $(BLASLIB) xlintstrfz: $(ZLINTSTRFP) $(TMGLIB) $(VARLIB) $(LAPACKLIB) $(BLASLIB)
$(LOADER) $(FFLAGS) $(LDFLAGS) -o $@ $^ $(FC) $(FFLAGS) $(LDFLAGS) -o $@ $^
$(ALINTST): $(FRC) $(ALINTST): $(FRC)
$(SCLNTST): $(FRC) $(SCLNTST): $(FRC)

View File

@ -797,6 +797,18 @@
WRITE( NOUT, FMT = 9978 ) WRITE( NOUT, FMT = 9978 )
$ SUBNAM(1:LEN_TRIM( SUBNAM )), INFO, M, N, IMAT $ SUBNAM(1:LEN_TRIM( SUBNAM )), INFO, M, N, IMAT
END IF END IF
*
ELSE IF( LSAMEN( 2, P2, 'QK' ) ) THEN
*
* xQK: truncated QR factorization with pivoting
*
IF( LSAMEN( 7, SUBNAM( 2: 8 ), 'GEQP3RK' ) ) THEN
WRITE( NOUT, FMT = 9930 )
$ SUBNAM(1:LEN_TRIM( SUBNAM )), INFO, M, N, KL, N5, IMAT
ELSE IF( LSAMEN( 5, SUBNAM( 2: 6 ), 'LATMS' ) ) THEN
WRITE( NOUT, FMT = 9978 )
$ SUBNAM(1:LEN_TRIM( SUBNAM )), INFO, M, N, IMAT
END IF
* *
ELSE IF( LSAMEN( 2, P2, 'LQ' ) ) THEN ELSE IF( LSAMEN( 2, P2, 'LQ' ) ) THEN
* *
@ -1147,6 +1159,11 @@
* What we do next * What we do next
* *
9949 FORMAT( ' ==> Doing only the condition estimate for this case' ) 9949 FORMAT( ' ==> Doing only the condition estimate for this case' )
*
* SUBNAM, INFO, M, N, NB, IMAT
*
9930 FORMAT( ' *** Error code from ', A, '=', I5, / ' ==> M =', I5,
$ ', N =', I5, ', NX =', I5, ', NB =', I4, ', type ', I2 )
* *
RETURN RETURN
* *

View File

@ -584,13 +584,27 @@
* *
* QR decomposition with column pivoting * QR decomposition with column pivoting
* *
WRITE( IOUNIT, FMT = 9986 )PATH WRITE( IOUNIT, FMT = 8006 )PATH
WRITE( IOUNIT, FMT = 9969 ) WRITE( IOUNIT, FMT = 9969 )
WRITE( IOUNIT, FMT = '( '' Test ratios:'' )' ) WRITE( IOUNIT, FMT = '( '' Test ratios:'' )' )
WRITE( IOUNIT, FMT = 9940 )1 WRITE( IOUNIT, FMT = 9940 )1
WRITE( IOUNIT, FMT = 9939 )2 WRITE( IOUNIT, FMT = 9939 )2
WRITE( IOUNIT, FMT = 9938 )3 WRITE( IOUNIT, FMT = 9938 )3
WRITE( IOUNIT, FMT = '( '' Messages:'' )' ) WRITE( IOUNIT, FMT = '( '' Messages:'' )' )
*
ELSE IF( LSAMEN( 2, P2, 'QK' ) ) THEN
*
* truncated QR decomposition with column pivoting
*
WRITE( IOUNIT, FMT = 8006 )PATH
WRITE( IOUNIT, FMT = 9871 )
WRITE( IOUNIT, FMT = '( '' Test ratios:'' )' )
WRITE( IOUNIT, FMT = 8060 )1
WRITE( IOUNIT, FMT = 8061 )2
WRITE( IOUNIT, FMT = 8062 )3
WRITE( IOUNIT, FMT = 8063 )4
WRITE( IOUNIT, FMT = 8064 )5
WRITE( IOUNIT, FMT = '( '' Messages:'' )' )
* *
ELSE IF( LSAMEN( 2, P2, 'TZ' ) ) THEN ELSE IF( LSAMEN( 2, P2, 'TZ' ) ) THEN
* *
@ -779,6 +793,8 @@
$ 'tall-skinny or short-wide matrices' ) $ 'tall-skinny or short-wide matrices' )
8005 FORMAT( / 1X, A3, ': Householder reconstruction from TSQR', 8005 FORMAT( / 1X, A3, ': Householder reconstruction from TSQR',
$ ' factorization output ', /,' for tall-skinny matrices.' ) $ ' factorization output ', /,' for tall-skinny matrices.' )
8006 FORMAT( / 1X, A3, ': truncated QR factorization',
$ ' with column pivoting' )
* *
* GE matrix types * GE matrix types
* *
@ -922,6 +938,36 @@
$ / 4X, '3. Geometric distribution', 10X, $ / 4X, '3. Geometric distribution', 10X,
$ '6. Every second column fixed' ) $ '6. Every second column fixed' )
* *
* QK matrix types
*
9871 FORMAT( 4X, ' 1. Zero matrix', /
$ 4X, ' 2. Random, Diagonal, CNDNUM = 2', /
$ 4X, ' 3. Random, Upper triangular, CNDNUM = 2', /
$ 4X, ' 4. Random, Lower triangular, CNDNUM = 2', /
$ 4X, ' 5. Random, First column is zero, CNDNUM = 2', /
$ 4X, ' 6. Random, Last MINMN column is zero, CNDNUM = 2', /
$ 4X, ' 7. Random, Last N column is zero, CNDNUM = 2', /
$ 4X, ' 8. Random, Middle column in MINMN is zero,',
$ ' CNDNUM = 2', /
$ 4X, ' 9. Random, First half of MINMN columns are zero,',
$ ' CNDNUM = 2', /
$ 4X, '10. Random, Last columns are zero starting from',
$ ' MINMN/2+1, CNDNUM = 2', /
$ 4X, '11. Random, Half MINMN columns in the middle are',
$ ' zero starting from MINMN/2-(MINMN/2)/2+1,'
$ ' CNDNUM = 2', /
$ 4X, '12. Random, Odd columns are ZERO, CNDNUM = 2', /
$ 4X, '13. Random, Even columns are ZERO, CNDNUM = 2', /
$ 4X, '14. Random, CNDNUM = 2', /
$ 4X, '15. Random, CNDNUM = sqrt(0.1/EPS)', /
$ 4X, '16. Random, CNDNUM = 0.1/EPS', /
$ 4X, '17. Random, CNDNUM = 0.1/EPS,',
$ ' one small singular value S(N)=1/CNDNUM', /
$ 4X, '18. Random, CNDNUM = 2, scaled near underflow,',
$ ' NORM = SMALL = SAFMIN', /
$ 4X, '19. Random, CNDNUM = 2, scaled near overflow,',
$ ' NORM = LARGE = 1.0/( 0.25 * ( SAFMIN / EPS ) )' )
*
* TZ matrix types * TZ matrix types
* *
9968 FORMAT( ' Matrix types (2-3 have condition 1/EPS):', / 4X, 9968 FORMAT( ' Matrix types (2-3 have condition 1/EPS):', / 4X,
@ -1030,8 +1076,7 @@
$ ' * norm(C) * EPS )' ) $ ' * norm(C) * EPS )' )
9940 FORMAT( 3X, I2, ': norm(svd(A) - svd(R)) / ', 9940 FORMAT( 3X, I2, ': norm(svd(A) - svd(R)) / ',
$ '( M * norm(svd(R)) * EPS )' ) $ '( M * norm(svd(R)) * EPS )' )
9939 FORMAT( 3X, I2, ': norm( A*P - Q*R ) / ( M * norm(A) * EPS )' 9939 FORMAT( 3X, I2, ': norm( A*P - Q*R ) / ( M * norm(A) * EPS )')
$ )
9938 FORMAT( 3X, I2, ': norm( I - Q''*Q ) / ( M * EPS )' ) 9938 FORMAT( 3X, I2, ': norm( I - Q''*Q ) / ( M * EPS )' )
9937 FORMAT( 3X, I2, ': norm( A - R*Q ) / ( M * norm(A) * EPS )' 9937 FORMAT( 3X, I2, ': norm( A - R*Q ) / ( M * norm(A) * EPS )'
$ ) $ )
@ -1105,6 +1150,15 @@
8054 FORMAT(3X,I2,': norm( C*Q - C*Q ) / ( M * norm(C) * EPS )' ) 8054 FORMAT(3X,I2,': norm( C*Q - C*Q ) / ( M * norm(C) * EPS )' )
8055 FORMAT(3X,I2,': norm( C*Q'' - C*Q'' ) / ( M * norm(C) * EPS )') 8055 FORMAT(3X,I2,': norm( C*Q'' - C*Q'' ) / ( M * norm(C) * EPS )')
8060 FORMAT( 3X, I2, ': 2-norm(svd(A) - svd(R)) / ',
$ '( max(M,N) * 2-norm(svd(R)) * EPS )' )
8061 FORMAT( 3X, I2, ': 1-norm( A*P - Q*R ) / ( max(M,N) * 1-norm(A)',
$ ' * EPS )')
8062 FORMAT( 3X, I2, ': 1-norm( I - Q''*Q ) / ( M * EPS )' )
8063 FORMAT( 3X, I2, ': Returns 1.0D+100, if abs(R(K+1,K+1))',
$ ' > abs(R(K,K)), where K=1:KFACT-1' )
8064 FORMAT( 3X, I2, ': 1-norm(Q**T * B - Q**T * B ) / ( M * EPS )')
* *
RETURN RETURN
* *

View File

@ -28,12 +28,12 @@
*> to evaluate the input line which requested NMATS matrix types for *> to evaluate the input line which requested NMATS matrix types for
*> PATH. The flow of control is as follows: *> PATH. The flow of control is as follows:
*> *>
*> If NMATS = NTYPES then *> IF NMATS = NTYPES THEN
*> DOTYPE(1:NTYPES) = .TRUE. *> DOTYPE(1:NTYPES) = .TRUE.
*> else *> ELSE
*> Read the next input line for NMATS matrix types *> Read the next input line for NMATS matrix types
*> Set DOTYPE(I) = .TRUE. for each valid type I *> Set DOTYPE(I) = .TRUE. for each valid type I
*> endif *> END IF
*> \endverbatim *> \endverbatim
* *
* Arguments: * Arguments:

View File

@ -69,6 +69,7 @@
*> CLQ 8 List types on next line if 0 < NTYPES < 8 *> CLQ 8 List types on next line if 0 < NTYPES < 8
*> CQL 8 List types on next line if 0 < NTYPES < 8 *> CQL 8 List types on next line if 0 < NTYPES < 8
*> CQP 6 List types on next line if 0 < NTYPES < 6 *> CQP 6 List types on next line if 0 < NTYPES < 6
*> ZQK 19 List types on next line if 0 < NTYPES < 19
*> CTZ 3 List types on next line if 0 < NTYPES < 3 *> CTZ 3 List types on next line if 0 < NTYPES < 3
*> CLS 6 List types on next line if 0 < NTYPES < 6 *> CLS 6 List types on next line if 0 < NTYPES < 6
*> CEQ *> CEQ
@ -153,12 +154,11 @@
$ NBVAL( MAXIN ), NBVAL2( MAXIN ), $ NBVAL( MAXIN ), NBVAL2( MAXIN ),
$ NSVAL( MAXIN ), NVAL( MAXIN ), NXVAL( MAXIN ), $ NSVAL( MAXIN ), NVAL( MAXIN ), NXVAL( MAXIN ),
$ RANKVAL( MAXIN ), PIV( NMAX ) $ RANKVAL( MAXIN ), PIV( NMAX )
REAL S( 2*NMAX )
COMPLEX E( NMAX )
* .. * ..
* .. Allocatable Arrays .. * .. Allocatable Arrays ..
INTEGER AllocateStatus INTEGER AllocateStatus
REAL, DIMENSION(:), ALLOCATABLE :: RWORK REAL, DIMENSION(:), ALLOCATABLE :: RWORK, S
COMPLEX, DIMENSION(:), ALLOCATABLE :: E
COMPLEX, DIMENSION(:,:), ALLOCATABLE :: A, B, WORK COMPLEX, DIMENSION(:,:), ALLOCATABLE :: A, B, WORK
* .. * ..
* .. External Functions .. * .. External Functions ..
@ -170,14 +170,14 @@
EXTERNAL ALAREQ, CCHKEQ, CCHKGB, CCHKGE, CCHKGT, CCHKHE, EXTERNAL ALAREQ, CCHKEQ, CCHKGB, CCHKGE, CCHKGT, CCHKHE,
$ CCHKHE_ROOK, CCHKHE_RK, CCHKHE_AA, CCHKHP, $ CCHKHE_ROOK, CCHKHE_RK, CCHKHE_AA, CCHKHP,
$ CCHKLQ, CCHKUNHR_COL, CCHKPB, CCHKPO, CCHKPS, $ CCHKLQ, CCHKUNHR_COL, CCHKPB, CCHKPO, CCHKPS,
$ CCHKPP, CCHKPT, CCHKQ3, CCHKQL, CCHKQR, CCHKRQ, $ CCHKPP, CCHKPT, CCHKQ3, CCHKQP3RK, CCHKQL,
$ CCHKSP, CCHKSY, CCHKSY_ROOK, CCHKSY_RK, $ CCHKQR, CCHKRQ, CCHKSP, CCHKSY, CCHKSY_ROOK,
$ CCHKSY_AA, CCHKTB, CCHKTP, CCHKTR, CCHKTZ, $ CCHKSY_RK, CCHKSY_AA, CCHKTB, CCHKTP, CCHKTR,
$ CDRVGB, CDRVGE, CDRVGT, CDRVHE, CDRVHE_ROOK, $ CCHKTZ, CDRVGB, CDRVGE, CDRVGT, CDRVHE,
$ CDRVHE_RK, CDRVHE_AA, CDRVHP, CDRVLS, CDRVPB, $ CDRVHE_ROOK, CDRVHE_RK, CDRVHE_AA, CDRVHP,
$ CDRVPO, CDRVPP, CDRVPT, CDRVSP, CDRVSY, $ CDRVLS, CDRVPB, CDRVPO, CDRVPP, CDRVPT, CDRVSP,
$ CDRVSY_ROOK, CDRVSY_RK, CDRVSY_AA, ILAVER, $ CDRVSY, CDRVSY_ROOK, CDRVSY_RK, CDRVSY_AA,
$ CCHKQRT, CCHKQRTP $ ILAVER, CCHKQRT, CCHKQRTP
* .. * ..
* .. Scalars in Common .. * .. Scalars in Common ..
LOGICAL LERR, OK LOGICAL LERR, OK
@ -203,6 +203,10 @@
IF (AllocateStatus /= 0) STOP "*** Not enough memory ***" IF (AllocateStatus /= 0) STOP "*** Not enough memory ***"
ALLOCATE ( WORK( NMAX, NMAX+MAXRHS+10 ), STAT = AllocateStatus ) ALLOCATE ( WORK( NMAX, NMAX+MAXRHS+10 ), STAT = AllocateStatus )
IF (AllocateStatus /= 0) STOP "*** Not enough memory ***" IF (AllocateStatus /= 0) STOP "*** Not enough memory ***"
ALLOCATE ( E( NMAX ), STAT = AllocateStatus )
IF (AllocateStatus /= 0) STOP "*** Not enough memory ***"
ALLOCATE ( S( 2*NMAX ), STAT = AllocateStatus)
IF (AllocateStatus /= 0) STOP "*** Not enough memory ***"
ALLOCATE ( RWORK( 150*NMAX+2*MAXRHS ), STAT = AllocateStatus ) ALLOCATE ( RWORK( 150*NMAX+2*MAXRHS ), STAT = AllocateStatus )
IF (AllocateStatus /= 0) STOP "*** Not enough memory ***" IF (AllocateStatus /= 0) STOP "*** Not enough memory ***"
* .. * ..
@ -1109,6 +1113,23 @@
ELSE ELSE
WRITE( NOUT, FMT = 9989 )PATH WRITE( NOUT, FMT = 9989 )PATH
END IF END IF
*
ELSE IF( LSAMEN( 2, C2, 'QK' ) ) THEN
*
* QK: truncated QR factorization with pivoting
*
NTYPES = 19
CALL ALAREQ( PATH, NMATS, DOTYPE, NTYPES, NIN, NOUT )
*
IF( TSTCHK ) THEN
CALL CCHKQP3RK( DOTYPE, NM, MVAL, NN, NVAL, NNS, NSVAL,
$ NNB, NBVAL, NXVAL, THRESH, A( 1, 1 ),
$ A( 1, 2 ), B( 1, 1 ), B( 1, 2 ),
$ S( 1 ), B( 1, 4 ),
$ WORK, RWORK, IWORK, NOUT )
ELSE
WRITE( NOUT, FMT = 9989 )PATH
END IF
* *
ELSE IF( LSAMEN( 2, C2, 'LS' ) ) THEN ELSE IF( LSAMEN( 2, C2, 'LS' ) ) THEN
* *

View File

@ -0,0 +1,836 @@
*> \brief \b CCHKQP3RK
*
* =========== DOCUMENTATION ===========
*
* Online html documentation available at
* http://www.netlib.org/lapack/explore-html/
*
* Definition:
* ===========
*
* SUBROUTINE CCHKQP3RK( DOTYPE, NM, MVAL, NN, NVAL, NNS, NSVAL,
* $ NNB, NBVAL, NXVAL, THRESH, A, COPYA,
* $ B, COPYB, S, TAU,
* $ WORK, RWORK, IWORK, NOUT )
* IMPLICIT NONE
*
* .. Scalar Arguments ..
* INTEGER NM, NN, NNB, NOUT
* REAL THRESH
* ..
* .. Array Arguments ..
* LOGICAL DOTYPE( * )
* INTEGER IWORK( * ), MVAL( * ), NBVAL( * ), NVAL( * ),
* $ NXVAL( * )
* REAL S( * ), RWORK( * )
* COMPLEX A( * ), COPYA( * ), TAU( * ), WORK( * )
* ..
*
*
*> \par Purpose:
* =============
*>
*> \verbatim
*>
*> CCHKQP3RK tests CGEQP3RK.
*> \endverbatim
*
* Arguments:
* ==========
*
*> \param[in] DOTYPE
*> \verbatim
*> DOTYPE is LOGICAL array, dimension (NTYPES)
*> The matrix types to be used for testing. Matrices of type j
*> (for 1 <= j <= NTYPES) are used for testing if DOTYPE(j) =
*> .TRUE.; if DOTYPE(j) = .FALSE., then type j is not used.
*> \endverbatim
*>
*> \param[in] NM
*> \verbatim
*> NM is INTEGER
*> The number of values of M contained in the vector MVAL.
*> \endverbatim
*>
*> \param[in] MVAL
*> \verbatim
*> MVAL is INTEGER array, dimension (NM)
*> The values of the matrix row dimension M.
*> \endverbatim
*>
*> \param[in] NN
*> \verbatim
*> NN is INTEGER
*> The number of values of N contained in the vector NVAL.
*> \endverbatim
*>
*> \param[in] NVAL
*> \verbatim
*> NVAL is INTEGER array, dimension (NN)
*> The values of the matrix column dimension N.
*> \endverbatim
*>
*> \param[in] NNS
*> \verbatim
*> NNS is INTEGER
*> The number of values of NRHS contained in the vector NSVAL.
*> \endverbatim
*>
*> \param[in] NSVAL
*> \verbatim
*> NSVAL is INTEGER array, dimension (NNS)
*> The values of the number of right hand sides NRHS.
*> \endverbatim
*> \param[in] NNB
*> \verbatim
*> NNB is INTEGER
*> The number of values of NB and NX contained in the
*> vectors NBVAL and NXVAL. The blocking parameters are used
*> in pairs (NB,NX).
*> \endverbatim
*>
*> \param[in] NBVAL
*> \verbatim
*> NBVAL is INTEGER array, dimension (NNB)
*> The values of the blocksize NB.
*> \endverbatim
*>
*> \param[in] NXVAL
*> \verbatim
*> NXVAL is INTEGER array, dimension (NNB)
*> The values of the crossover point NX.
*> \endverbatim
*>
*> \param[in] THRESH
*> \verbatim
*> THRESH is REAL
*> The threshold value for the test ratios. A result is
*> included in the output file if RESULT >= THRESH. To have
*> every test ratio printed, use THRESH = 0.
*> \endverbatim
*>
*> \param[out] A
*> \verbatim
*> A is COMPLEX array, dimension (MMAX*NMAX)
*> where MMAX is the maximum value of M in MVAL and NMAX is the
*> maximum value of N in NVAL.
*> \endverbatim
*>
*> \param[out] COPYA
*> \verbatim
*> COPYA is COMPLEX array, dimension (MMAX*NMAX)
*> \endverbatim
*>
*> \param[out] B
*> \verbatim
*> B is COMPLEX array, dimension (MMAX*NSMAX)
*> where MMAX is the maximum value of M in MVAL and NSMAX is the
*> maximum value of NRHS in NSVAL.
*> \endverbatim
*>
*> \param[out] COPYB
*> \verbatim
*> COPYB is COMPLEX array, dimension (MMAX*NSMAX)
*> \endverbatim
*>
*> \param[out] S
*> \verbatim
*> S is REAL array, dimension
*> (min(MMAX,NMAX))
*> \endverbatim
*>
*> \param[out] TAU
*> \verbatim
*> TAU is COMPLEX array, dimension (MMAX)
*> \endverbatim
*>
*> \param[out] WORK
*> \verbatim
*> WORK is COMPLEX array, dimension
*> (max(M*max(M,N) + 4*min(M,N) + max(M,N)))
*> \endverbatim
*>
*> \param[out] RWORK
*> \verbatim
*> RWORK is REAL array, dimension (4*NMAX)
*> \endverbatim
*>
*> \param[out] IWORK
*> \verbatim
*> IWORK is INTEGER array, dimension (2*NMAX)
*> \endverbatim
*>
*> \param[in] NOUT
*> \verbatim
*> NOUT is INTEGER
*> The unit number for output.
*> \endverbatim
*
* Authors:
* ========
*
*> \author Univ. of Tennessee
*> \author Univ. of California Berkeley
*> \author Univ. of Colorado Denver
*> \author NAG Ltd.
*
*> \ingroup complex_lin
*
* =====================================================================
SUBROUTINE CCHKQP3RK( DOTYPE, NM, MVAL, NN, NVAL, NNS, NSVAL,
$ NNB, NBVAL, NXVAL, THRESH, A, COPYA,
$ B, COPYB, S, TAU,
$ WORK, RWORK, IWORK, NOUT )
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 NM, NN, NNB, NNS, NOUT
REAL THRESH
* ..
* .. Array Arguments ..
LOGICAL DOTYPE( * )
INTEGER IWORK( * ), NBVAL( * ), MVAL( * ), NVAL( * ),
$ NSVAL( * ), NXVAL( * )
REAL S( * ), RWORK( * )
COMPLEX A( * ), COPYA( * ), B( * ), COPYB( * ),
$ TAU( * ), WORK( * )
* ..
*
* =====================================================================
*
* .. Parameters ..
INTEGER NTYPES
PARAMETER ( NTYPES = 19 )
INTEGER NTESTS
PARAMETER ( NTESTS = 5 )
REAL ONE, ZERO, BIGNUM
COMPLEX CONE, CZERO
PARAMETER ( ONE = 1.0E+0, ZERO = 0.0E+0,
$ CZERO = ( 0.0E+0, 0.0E+0 ),
$ CONE = ( 1.0E+0, 0.0E+0 ),
$ BIGNUM = 1.0E+38 )
* ..
* .. Local Scalars ..
CHARACTER DIST, TYPE
CHARACTER*3 PATH
INTEGER I, IHIGH, ILOW, IM, IMAT, IN, INC_ZERO,
$ INB, IND_OFFSET_GEN,
$ IND_IN, IND_OUT, INS, INFO,
$ ISTEP, J, J_INC, J_FIRST_NZ, JB_ZERO,
$ KFACT, KL, KMAX, KU, LDA, LW, LWORK,
$ LWORK_MQR, M, MINMN, MINMNB_GEN, MODE, N,
$ NB, NB_ZERO, NERRS, NFAIL, NB_GEN, NRHS,
$ NRUN, NX, T
REAL ANORM, CNDNUM, EPS, ABSTOL, RELTOL,
$ DTEMP, MAXC2NRMK, RELMAXC2NRMK
* ..
* .. Local Arrays ..
INTEGER ISEED( 4 ), ISEEDY( 4 )
REAL RESULT( NTESTS ), RDUMMY( 1 )
* ..
* .. External Functions ..
REAL SLAMCH, CQPT01, CQRT11, CQRT12, CLANGE
EXTERNAL SLAMCH, CQPT01, CQRT11, CQRT12, CLANGE
* ..
* .. External Subroutines ..
EXTERNAL ALAERH, ALAHD, ALASUM, SLAORD, ICOPY, CAXPY,
$ XLAENV, CGEQP3RK, CLACPY, CLASET, CLATB4,
$ CLATMS, CUNMQR, CSWAP
* ..
* .. Intrinsic Functions ..
INTRINSIC ABS, MAX, MIN, MOD, REAL
* ..
* .. Scalars in Common ..
LOGICAL LERR, OK
CHARACTER*32 SRNAMT
INTEGER INFOT, IOUNIT, CUNMQR_LWORK
* ..
* .. Common blocks ..
COMMON / INFOC / INFOT, IOUNIT, OK, LERR
COMMON / SRNAMC / SRNAMT
* ..
* .. Data statements ..
DATA ISEEDY / 1988, 1989, 1990, 1991 /
* ..
* .. Executable Statements ..
*
* Initialize constants and the random number seed.
*
PATH( 1: 1 ) = 'Complex precision'
PATH( 2: 3 ) = 'QK'
NRUN = 0
NFAIL = 0
NERRS = 0
DO I = 1, 4
ISEED( I ) = ISEEDY( I )
END DO
EPS = SLAMCH( 'Epsilon' )
INFOT = 0
*
DO IM = 1, NM
*
* Do for each value of M in MVAL.
*
M = MVAL( IM )
LDA = MAX( 1, M )
*
DO IN = 1, NN
*
* Do for each value of N in NVAL.
*
N = NVAL( IN )
MINMN = MIN( M, N )
LWORK = MAX( 1, M*MAX( M, N )+4*MINMN+MAX( M, N ),
$ M*N + 2*MINMN + 4*N )
*
DO INS = 1, NNS
NRHS = NSVAL( INS )
*
* Set up parameters with CLATB4 and generate
* M-by-NRHS B matrix with CLATMS.
* IMAT = 14:
* Random matrix, CNDNUM = 2, NORM = ONE,
* MODE = 3 (geometric distribution of singular values).
*
CALL CLATB4( PATH, 14, M, NRHS, TYPE, KL, KU, ANORM,
$ MODE, CNDNUM, DIST )
*
SRNAMT = 'CLATMS'
CALL CLATMS( M, NRHS, DIST, ISEED, TYPE, S, MODE,
$ CNDNUM, ANORM, KL, KU, 'No packing',
$ COPYB, LDA, WORK, INFO )
*
* Check error code from CLATMS.
*
IF( INFO.NE.0 ) THEN
CALL ALAERH( PATH, 'CLATMS', INFO, 0, ' ', M,
$ NRHS, -1, -1, -1, 6, NFAIL, NERRS,
$ NOUT )
CYCLE
END IF
*
DO IMAT = 1, NTYPES
*
* Do the tests only if DOTYPE( IMAT ) is true.
*
IF( .NOT.DOTYPE( IMAT ) )
$ CYCLE
*
* The type of distribution used to generate the random
* eigen-/singular values:
* ( 'S' for symmetric distribution ) => UNIFORM( -1, 1 )
*
* Do for each type of NON-SYMMETRIC matrix: CNDNUM NORM MODE
* 1. Zero matrix
* 2. Random, Diagonal, CNDNUM = 2 CNDNUM = 2 ONE 3 ( geometric distribution of singular values )
* 3. Random, Upper triangular, CNDNUM = 2 CNDNUM = 2 ONE 3 ( geometric distribution of singular values )
* 4. Random, Lower triangular, CNDNUM = 2 CNDNUM = 2 ONE 3 ( geometric distribution of singular values )
* 5. Random, First column is zero, CNDNUM = 2 CNDNUM = 2 ONE 3 ( geometric distribution of singular values )
* 6. Random, Last MINMN column is zero, CNDNUM = 2 CNDNUM = 2 ONE 3 ( geometric distribution of singular values )
* 7. Random, Last N column is zero, CNDNUM = 2 CNDNUM = 2 ONE 3 ( geometric distribution of singular values )
* 8. Random, Middle column in MINMN is zero, CNDNUM = 2 CNDNUM = 2 ONE 3 ( geometric distribution of singular values )
* 9. Random, First half of MINMN columns are zero, CNDNUM = 2 CNDNUM = 2 ONE 3 ( geometric distribution of singular values )
* 10. Random, Last columns are zero starting from MINMN/2+1, CNDNUM = 2 CNDNUM = 2 ONE 3 ( geometric distribution of singular values )
* 11. Random, Half MINMN columns in the middle are zero starting
* from MINMN/2-(MINMN/2)/2+1, CNDNUM = 2 CNDNUM = 2 ONE 3 ( geometric distribution of singular values )
* 12. Random, Odd columns are ZERO, CNDNUM = 2 CNDNUM = 2 ONE 3 ( geometric distribution of singular values )
* 13. Random, Even columns are ZERO, CNDNUM = 2 CNDNUM = 2 ONE 3 ( geometric distribution of singular values )
* 14. Random, CNDNUM = 2 CNDNUM = 2 ONE 3 ( geometric distribution of singular values )
* 15. Random, CNDNUM = sqrt(0.1/EPS) CNDNUM = BADC1 = sqrt(0.1/EPS) ONE 3 ( geometric distribution of singular values )
* 16. Random, CNDNUM = 0.1/EPS CNDNUM = BADC2 = 0.1/EPS ONE 3 ( geometric distribution of singular values )
* 17. Random, CNDNUM = 0.1/EPS, CNDNUM = BADC2 = 0.1/EPS ONE 2 ( one small singular value, S(N)=1/CNDNUM )
* one small singular value S(N)=1/CNDNUM
* 18. Random, CNDNUM = 2, scaled near underflow CNDNUM = 2 SMALL = SAFMIN
* 19. Random, CNDNUM = 2, scaled near overflow CNDNUM = 2 LARGE = 1.0/( 0.25 * ( SAFMIN / EPS ) ) 3 ( geometric distribution of singular values )
*
IF( IMAT.EQ.1 ) THEN
*
* Matrix 1: Zero matrix
*
CALL CLASET( 'Full', M, N, CZERO, CZERO, COPYA, LDA )
DO I = 1, MINMN
S( I ) = ZERO
END DO
*
ELSE IF( (IMAT.GE.2 .AND. IMAT.LE.4 )
$ .OR. (IMAT.GE.14 .AND. IMAT.LE.19 ) ) THEN
*
* Matrices 2-5.
*
* Set up parameters with DLATB4 and generate a test
* matrix with CLATMS.
*
CALL CLATB4( PATH, IMAT, M, N, TYPE, KL, KU, ANORM,
$ MODE, CNDNUM, DIST )
*
SRNAMT = 'CLATMS'
CALL CLATMS( M, N, DIST, ISEED, TYPE, S, MODE,
$ CNDNUM, ANORM, KL, KU, 'No packing',
$ COPYA, LDA, WORK, INFO )
*
* Check error code from CLATMS.
*
IF( INFO.NE.0 ) THEN
CALL ALAERH( PATH, 'CLATMS', INFO, 0, ' ', M, N,
$ -1, -1, -1, IMAT, NFAIL, NERRS,
$ NOUT )
CYCLE
END IF
*
CALL SLAORD( 'Decreasing', MINMN, S, 1 )
*
ELSE IF( MINMN.GE.2
$ .AND. IMAT.GE.5 .AND. IMAT.LE.13 ) THEN
*
* Rectangular matrices 5-13 that contain zero columns,
* only for matrices MINMN >=2.
*
* JB_ZERO is the column index of ZERO block.
* NB_ZERO is the column block size of ZERO block.
* NB_GEN is the column blcok size of the
* generated block.
* J_INC in the non_zero column index increment
* for matrix 12 and 13.
* J_FIRS_NZ is the index of the first non-zero
* column.
*
IF( IMAT.EQ.5 ) THEN
*
* First column is zero.
*
JB_ZERO = 1
NB_ZERO = 1
NB_GEN = N - NB_ZERO
*
ELSE IF( IMAT.EQ.6 ) THEN
*
* Last column MINMN is zero.
*
JB_ZERO = MINMN
NB_ZERO = 1
NB_GEN = N - NB_ZERO
*
ELSE IF( IMAT.EQ.7 ) THEN
*
* Last column N is zero.
*
JB_ZERO = N
NB_ZERO = 1
NB_GEN = N - NB_ZERO
*
ELSE IF( IMAT.EQ.8 ) THEN
*
* Middle column in MINMN is zero.
*
JB_ZERO = MINMN / 2 + 1
NB_ZERO = 1
NB_GEN = N - NB_ZERO
*
ELSE IF( IMAT.EQ.9 ) THEN
*
* First half of MINMN columns is zero.
*
JB_ZERO = 1
NB_ZERO = MINMN / 2
NB_GEN = N - NB_ZERO
*
ELSE IF( IMAT.EQ.10 ) THEN
*
* Last columns are zero columns,
* starting from (MINMN / 2 + 1) column.
*
JB_ZERO = MINMN / 2 + 1
NB_ZERO = N - JB_ZERO + 1
NB_GEN = N - NB_ZERO
*
ELSE IF( IMAT.EQ.11 ) THEN
*
* Half of the columns in the middle of MINMN
* columns is zero, starting from
* MINMN/2 - (MINMN/2)/2 + 1 column.
*
JB_ZERO = MINMN / 2 - (MINMN / 2) / 2 + 1
NB_ZERO = MINMN / 2
NB_GEN = N - NB_ZERO
*
ELSE IF( IMAT.EQ.12 ) THEN
*
* Odd-numbered columns are zero,
*
NB_GEN = N / 2
NB_ZERO = N - NB_GEN
J_INC = 2
J_FIRST_NZ = 2
*
ELSE IF( IMAT.EQ.13 ) THEN
*
* Even-numbered columns are zero.
*
NB_ZERO = N / 2
NB_GEN = N - NB_ZERO
J_INC = 2
J_FIRST_NZ = 1
*
END IF
*
*
* 1) Set the first NB_ZERO columns in COPYA(1:M,1:N)
* to zero.
*
CALL CLASET( 'Full', M, NB_ZERO, CZERO, CZERO,
$ COPYA, LDA )
*
* 2) Generate an M-by-(N-NB_ZERO) matrix with the
* chosen singular value distribution
* in COPYA(1:M,NB_ZERO+1:N).
*
CALL CLATB4( PATH, IMAT, M, NB_GEN, TYPE, KL, KU,
$ ANORM, MODE, CNDNUM, DIST )
*
SRNAMT = 'CLATMS'
*
IND_OFFSET_GEN = NB_ZERO * LDA
*
CALL CLATMS( M, NB_GEN, DIST, ISEED, TYPE, S, MODE,
$ CNDNUM, ANORM, KL, KU, 'No packing',
$ COPYA( IND_OFFSET_GEN + 1 ), LDA,
$ WORK, INFO )
*
* Check error code from CLATMS.
*
IF( INFO.NE.0 ) THEN
CALL ALAERH( PATH, 'CLATMS', INFO, 0, ' ', M,
$ NB_GEN, -1, -1, -1, IMAT, NFAIL,
$ NERRS, NOUT )
CYCLE
END IF
*
* 3) Swap the gererated colums from the right side
* NB_GEN-size block in COPYA into correct column
* positions.
*
IF( IMAT.EQ.6
$ .OR. IMAT.EQ.7
$ .OR. IMAT.EQ.8
$ .OR. IMAT.EQ.10
$ .OR. IMAT.EQ.11 ) THEN
*
* Move by swapping the generated columns
* from the right NB_GEN-size block from
* (NB_ZERO+1:NB_ZERO+JB_ZERO)
* into columns (1:JB_ZERO-1).
*
DO J = 1, JB_ZERO-1, 1
CALL CSWAP( M,
$ COPYA( ( NB_ZERO+J-1)*LDA+1), 1,
$ COPYA( (J-1)*LDA + 1 ), 1 )
END DO
*
ELSE IF( IMAT.EQ.12 .OR. IMAT.EQ.13 ) THEN
*
* ( IMAT = 12, Odd-numbered ZERO columns. )
* Swap the generated columns from the right
* NB_GEN-size block into the even zero colums in the
* left NB_ZERO-size block.
*
* ( IMAT = 13, Even-numbered ZERO columns. )
* Swap the generated columns from the right
* NB_GEN-size block into the odd zero colums in the
* left NB_ZERO-size block.
*
DO J = 1, NB_GEN, 1
IND_OUT = ( NB_ZERO+J-1 )*LDA + 1
IND_IN = ( J_INC*(J-1)+(J_FIRST_NZ-1) )*LDA
$ + 1
CALL CSWAP( M,
$ COPYA( IND_OUT ), 1,
$ COPYA( IND_IN), 1 )
END DO
*
END IF
*
* 5) Order the singular values generated by
* DLAMTS in decreasing order and add trailing zeros
* that correspond to zero columns.
* The total number of singular values is MINMN.
*
MINMNB_GEN = MIN( M, NB_GEN )
*
CALL SLAORD( 'Decreasing', MINMNB_GEN, S, 1 )
DO I = MINMNB_GEN+1, MINMN
S( I ) = ZERO
END DO
*
ELSE
*
* IF(MINMN.LT.2) skip this size for this matrix type.
*
CYCLE
END IF
*
* Initialize a copy array for a pivot array for DGEQP3RK.
*
DO I = 1, N
IWORK( I ) = 0
END DO
*
DO INB = 1, NNB
*
* Do for each pair of values (NB,NX) in NBVAL and NXVAL.
*
NB = NBVAL( INB )
CALL XLAENV( 1, NB )
NX = NXVAL( INB )
CALL XLAENV( 3, NX )
*
* We do MIN(M,N)+1 because we need a test for KMAX > N,
* when KMAX is larger than MIN(M,N), KMAX should be
* KMAX = MIN(M,N)
*
DO KMAX = 0, MIN(M,N)+1
*
* Get a working copy of COPYA into A( 1:M,1:N ).
* Get a working copy of COPYB into A( 1:M, (N+1):NRHS ).
* Get a working copy of COPYB into into B( 1:M, 1:NRHS ).
* Get a working copy of IWORK(1:N) awith zeroes into
* which is going to be used as pivot array IWORK( N+1:2N ).
* NOTE: IWORK(2N+1:3N) is going to be used as a WORK array
* for the routine.
*
CALL CLACPY( 'All', M, N, COPYA, LDA, A, LDA )
CALL CLACPY( 'All', M, NRHS, COPYB, LDA,
$ A( LDA*N + 1 ), LDA )
CALL CLACPY( 'All', M, NRHS, COPYB, LDA,
$ B, LDA )
CALL ICOPY( N, IWORK( 1 ), 1, IWORK( N+1 ), 1 )
*
ABSTOL = -1.0
RELTOl = -1.0
*
* Compute the QR factorization with pivoting of A
*
LW = MAX( 1, MAX( 2*N + NB*( N+NRHS+1 ),
$ 3*N + NRHS - 1 ) )
*
* Compute CGEQP3RK factorization of A.
*
SRNAMT = 'CGEQP3RK'
CALL CGEQP3RK( M, N, NRHS, KMAX, ABSTOL, RELTOL,
$ A, LDA, KFACT, MAXC2NRMK,
$ RELMAXC2NRMK, IWORK( N+1 ), TAU,
$ WORK, LW, RWORK, IWORK( 2*N+1 ),
$ INFO )
*
* Check error code from CGEQP3RK.
*
IF( INFO.LT.0 )
$ CALL ALAERH( PATH, 'CGEQP3RK', INFO, 0, ' ',
$ M, N, NX, -1, NB, IMAT,
$ NFAIL, NERRS, NOUT )
*
IF( KFACT.EQ.MINMN ) THEN
*
* Compute test 1:
*
* This test in only for the full rank factorization of
* the matrix A.
*
* Array S(1:min(M,N)) contains svd(A) the sigular values
* of the original matrix A in decreasing absolute value
* order. The test computes svd(R), the vector sigular
* values of the upper trapezoid of A(1:M,1:N) that
* contains the factor R, in decreasing order. The test
* returns the ratio:
*
* 2-norm(svd(R) - svd(A)) / ( max(M,N) * 2-norm(svd(A)) * EPS )
*
RESULT( 1 ) = CQRT12( M, N, A, LDA, S, WORK,
$ LWORK , RWORK )
*
DO T = 1, 1
IF( RESULT( T ).GE.THRESH ) THEN
IF( NFAIL.EQ.0 .AND. NERRS.EQ.0 )
$ CALL ALAHD( NOUT, PATH )
WRITE( NOUT, FMT = 9999 ) 'CGEQP3RK', M, N,
$ NRHS, KMAX, ABSTOL, RELTOL, NB, NX,
$ IMAT, T, RESULT( T )
NFAIL = NFAIL + 1
END IF
END DO
NRUN = NRUN + 1
*
* End test 1
*
END IF
* Compute test 2:
*
* The test returns the ratio:
*
* 1-norm( A*P - Q*R ) / ( max(M,N) * 1-norm(A) * EPS )
*
RESULT( 2 ) = CQPT01( M, N, KFACT, COPYA, A, LDA, TAU,
$ IWORK( N+1 ), WORK, LWORK )
*
* Compute test 3:
*
* The test returns the ratio:
*
* 1-norm( Q**T * Q - I ) / ( M * EPS )
*
RESULT( 3 ) = CQRT11( M, KFACT, A, LDA, TAU, WORK,
$ LWORK )
*
* Print information about the tests that did not pass
* the threshold.
*
DO T = 2, 3
IF( RESULT( T ).GE.THRESH ) THEN
IF( NFAIL.EQ.0 .AND. NERRS.EQ.0 )
$ CALL ALAHD( NOUT, PATH )
WRITE( NOUT, FMT = 9999 ) 'CGEQP3RK', M, N,
$ NRHS, KMAX, ABSTOL, RELTOL,
$ NB, NX, IMAT, T, RESULT( T )
NFAIL = NFAIL + 1
END IF
END DO
NRUN = NRUN + 2
*
* Compute test 4:
*
* This test is only for the factorizations with the
* rank greater than 2.
* The elements on the diagonal of R should be non-
* increasing.
*
* The test returns the ratio:
*
* Returns 1.0D+100 if abs(R(K+1,K+1)) > abs(R(K,K)),
* K=1:KFACT-1
*
IF( MIN(KFACT, MINMN).GE.2 ) THEN
*
DO J = 1, KFACT-1, 1
*
DTEMP = (( ABS( A( (J-1)*M+J ) ) -
$ ABS( A( (J)*M+J+1 ) ) ) /
$ ABS( A(1) ) )
*
IF( DTEMP.LT.ZERO ) THEN
RESULT( 4 ) = BIGNUM
END IF
*
END DO
*
* Print information about the tests that did not
* pass the threshold.
*
DO T = 4, 4
IF( RESULT( T ).GE.THRESH ) THEN
IF( NFAIL.EQ.0 .AND. NERRS.EQ.0 )
$ CALL ALAHD( NOUT, PATH )
WRITE( NOUT, FMT = 9999 ) 'CGEQP3RK',
$ M, N, NRHS, KMAX, ABSTOL, RELTOL,
$ NB, NX, IMAT, T,
$ RESULT( T )
NFAIL = NFAIL + 1
END IF
END DO
NRUN = NRUN + 1
*
* End test 4.
*
END IF
*
* Compute test 5:
*
* This test in only for matrix A with min(M,N) > 0.
*
* The test returns the ratio:
*
* 1-norm(Q**T * B - Q**T * B ) /
* ( M * EPS )
*
* (1) Compute B:=Q**T * B in the matrix B.
*
IF( MINMN.GT.0 ) THEN
*
LWORK_MQR = MAX(1, NRHS)
CALL CUNMQR( 'Left', 'Conjugate transpose',
$ M, NRHS, KFACT, A, LDA, TAU, B, LDA,
$ WORK, LWORK_MQR, INFO )
*
DO I = 1, NRHS
*
* Compare N+J-th column of A and J-column of B.
*
CALL CAXPY( M, -CONE, A( ( N+I-1 )*LDA+1 ), 1,
$ B( ( I-1 )*LDA+1 ), 1 )
END DO
*
RESULT( 5 ) =
$ ABS(
$ CLANGE( 'One-norm', M, NRHS, B, LDA, RDUMMY ) /
$ ( REAL( M )*SLAMCH( 'Epsilon' ) )
$ )
*
* Print information about the tests that did not pass
* the threshold.
*
DO T = 5, 5
IF( RESULT( T ).GE.THRESH ) THEN
IF( NFAIL.EQ.0 .AND. NERRS.EQ.0 )
$ CALL ALAHD( NOUT, PATH )
WRITE( NOUT, FMT = 9999 ) 'CGEQP3RK', M, N,
$ NRHS, KMAX, ABSTOL, RELTOL,
$ NB, NX, IMAT, T, RESULT( T )
NFAIL = NFAIL + 1
END IF
END DO
NRUN = NRUN + 1
*
* End compute test 5.
*
END IF
*
* END DO KMAX = 1, MIN(M,N)+1
*
END DO
*
* END DO for INB = 1, NNB
*
END DO
*
* END DO for IMAT = 1, NTYPES
*
END DO
*
* END DO for INS = 1, NNS
*
END DO
*
* END DO for IN = 1, NN
*
END DO
*
* END DO for IM = 1, NM
*
END DO
*
* Print a summary of the results.
*
CALL ALASUM( PATH, NOUT, NFAIL, NRUN, NERRS )
*
9999 FORMAT( 1X, A, ' M =', I5, ', N =', I5, ', NRHS =', I5,
$ ', KMAX =', I5, ', ABSTOL =', G12.5,
$ ', RELTOL =', G12.5, ', NB =', I4, ', NX =', I4,
$ ', type ', I2, ', test ', I2, ', ratio =', G12.5 )
*
* End of CCHKQP3RK
*
END

View File

@ -154,9 +154,6 @@
* .. Intrinsic Functions .. * .. Intrinsic Functions ..
INTRINSIC ABS, MAX, SQRT INTRINSIC ABS, MAX, SQRT
* .. * ..
* .. External Subroutines ..
EXTERNAL SLABAD
* ..
* .. Save statement .. * .. Save statement ..
SAVE EPS, SMALL, LARGE, BADC1, BADC2, FIRST SAVE EPS, SMALL, LARGE, BADC1, BADC2, FIRST
* .. * ..
@ -174,11 +171,6 @@
BADC1 = SQRT( BADC2 ) BADC1 = SQRT( BADC2 )
SMALL = SLAMCH( 'Safe minimum' ) SMALL = SLAMCH( 'Safe minimum' )
LARGE = ONE / SMALL LARGE = ONE / SMALL
*
* If it looks like we're on a Cray, take the square root of
* SMALL and LARGE to avoid overflow and underflow problems.
*
CALL SLABAD( SMALL, LARGE )
SMALL = SHRINK*( SMALL / EPS ) SMALL = SHRINK*( SMALL / EPS )
LARGE = ONE / SMALL LARGE = ONE / SMALL
END IF END IF
@ -233,6 +225,110 @@
ELSE ELSE
ANORM = ONE ANORM = ONE
END IF END IF
*
ELSE IF( LSAMEN( 2, C2, 'QK' ) ) THEN
*
* xQK: truncated QR with pivoting.
* Set parameters to generate a general
* M x N matrix.
*
* Set TYPE, the type of matrix to be generated. 'N' is nonsymmetric.
*
TYPE = 'N'
*
* Set DIST, the type of distribution for the random
* number generator. 'S' is
*
DIST = 'S'
*
* Set the lower and upper bandwidths.
*
IF( IMAT.EQ.2 ) THEN
*
* 2. Random, Diagonal, CNDNUM = 2
*
KL = 0
KU = 0
CNDNUM = TWO
ANORM = ONE
MODE = 3
ELSE IF( IMAT.EQ.3 ) THEN
*
* 3. Random, Upper triangular, CNDNUM = 2
*
KL = 0
KU = MAX( N-1, 0 )
CNDNUM = TWO
ANORM = ONE
MODE = 3
ELSE IF( IMAT.EQ.4 ) THEN
*
* 4. Random, Lower triangular, CNDNUM = 2
*
KL = MAX( M-1, 0 )
KU = 0
CNDNUM = TWO
ANORM = ONE
MODE = 3
ELSE
*
* 5.-19. Rectangular matrix
*
KL = MAX( M-1, 0 )
KU = MAX( N-1, 0 )
*
IF( IMAT.GE.5 .AND. IMAT.LE.14 ) THEN
*
* 5.-14. Random, CNDNUM = 2.
*
CNDNUM = TWO
ANORM = ONE
MODE = 3
*
ELSE IF( IMAT.EQ.15 ) THEN
*
* 15. Random, CNDNUM = sqrt(0.1/EPS)
*
CNDNUM = BADC1
ANORM = ONE
MODE = 3
*
ELSE IF( IMAT.EQ.16 ) THEN
*
* 16. Random, CNDNUM = 0.1/EPS
*
CNDNUM = BADC2
ANORM = ONE
MODE = 3
*
ELSE IF( IMAT.EQ.17 ) THEN
*
* 17. Random, CNDNUM = 0.1/EPS,
* one small singular value S(N)=1/CNDNUM
*
CNDNUM = BADC2
ANORM = ONE
MODE = 2
*
ELSE IF( IMAT.EQ.18 ) THEN
*
* 18. Random, scaled near underflow
*
CNDNUM = TWO
ANORM = SMALL
MODE = 3
*
ELSE IF( IMAT.EQ.19 ) THEN
*
* 19. Random, scaled near overflow
*
CNDNUM = TWO
ANORM = LARGE
MODE = 3
*
END IF
*
END IF
* *
ELSE IF( LSAMEN( 2, C2, 'GE' ) ) THEN ELSE IF( LSAMEN( 2, C2, 'GE' ) ) THEN
* *
@ -517,17 +613,18 @@
* *
* Set the norm and condition number. * Set the norm and condition number.
* *
IF( IMAT.EQ.2 .OR. IMAT.EQ.8 ) THEN MAT = ABS( IMAT )
IF( MAT.EQ.2 .OR. MAT.EQ.8 ) THEN
CNDNUM = BADC1 CNDNUM = BADC1
ELSE IF( IMAT.EQ.3 .OR. IMAT.EQ.9 ) THEN ELSE IF( MAT.EQ.3 .OR. MAT.EQ.9 ) THEN
CNDNUM = BADC2 CNDNUM = BADC2
ELSE ELSE
CNDNUM = TWO CNDNUM = TWO
END IF END IF
* *
IF( IMAT.EQ.4 ) THEN IF( MAT.EQ.4 ) THEN
ANORM = SMALL ANORM = SMALL
ELSE IF( IMAT.EQ.5 ) THEN ELSE IF( MAT.EQ.5 ) THEN
ANORM = LARGE ANORM = LARGE
ELSE ELSE
ANORM = ONE ANORM = ONE

View File

@ -33,7 +33,8 @@
*> Householder vectors, and the rest of AF contains a partially updated *> Householder vectors, and the rest of AF contains a partially updated
*> matrix. *> matrix.
*> *>
*> This function returns ||A*P - Q*R||/(||norm(A)||*eps*M) *> This function returns ||A*P - Q*R|| / ( ||norm(A)||*eps*max(M,N) )
*> where || . || is matrix one norm.
*> \endverbatim *> \endverbatim
* *
* Arguments: * Arguments:
@ -172,28 +173,28 @@
* *
NORMA = CLANGE( 'One-norm', M, N, A, LDA, RWORK ) NORMA = CLANGE( 'One-norm', M, N, A, LDA, RWORK )
* *
DO 30 J = 1, K DO J = 1, K
DO 10 I = 1, MIN( J, M ) DO I = 1, MIN( J, M )
WORK( ( J-1 )*M+I ) = AF( I, J ) WORK( ( J-1 )*M+I ) = AF( I, J )
10 CONTINUE END DO
DO 20 I = J + 1, M DO I = J + 1, M
WORK( ( J-1 )*M+I ) = ZERO WORK( ( J-1 )*M+I ) = ZERO
20 CONTINUE END DO
30 CONTINUE END DO
DO 40 J = K + 1, N DO J = K + 1, N
CALL CCOPY( M, AF( 1, J ), 1, WORK( ( J-1 )*M+1 ), 1 ) CALL CCOPY( M, AF( 1, J ), 1, WORK( ( J-1 )*M+1 ), 1 )
40 CONTINUE END DO
* *
CALL CUNMQR( 'Left', 'No transpose', M, N, K, AF, LDA, TAU, WORK, CALL CUNMQR( 'Left', 'No transpose', M, N, K, AF, LDA, TAU, WORK,
$ M, WORK( M*N+1 ), LWORK-M*N, INFO ) $ M, WORK( M*N+1 ), LWORK-M*N, INFO )
* *
DO 50 J = 1, N DO J = 1, N
* *
* Compare i-th column of QR and jpvt(i)-th column of A * Compare i-th column of QR and jpvt(i)-th column of A
* *
CALL CAXPY( M, CMPLX( -ONE ), A( 1, JPVT( J ) ), 1, CALL CAXPY( M, CMPLX( -ONE ), A( 1, JPVT( J ) ), 1,
$ WORK( ( J-1 )*M+1 ), 1 ) $ WORK( ( J-1 )*M+1 ), 1 )
50 CONTINUE END DO
* *
CQPT01 = CLANGE( 'One-norm', M, N, WORK, M, RWORK ) / CQPT01 = CLANGE( 'One-norm', M, N, WORK, M, RWORK ) /
$ ( REAL( MAX( M, N ) )*SLAMCH( 'Epsilon' ) ) $ ( REAL( MAX( M, N ) )*SLAMCH( 'Epsilon' ) )

View File

@ -157,9 +157,9 @@
CALL CUNM2R( 'Left', 'Conjugate transpose', M, M, K, A, LDA, TAU, CALL CUNM2R( 'Left', 'Conjugate transpose', M, M, K, A, LDA, TAU,
$ WORK, M, WORK( M*M+1 ), INFO ) $ WORK, M, WORK( M*M+1 ), INFO )
* *
DO 10 J = 1, M DO J = 1, M
WORK( ( J-1 )*M+J ) = WORK( ( J-1 )*M+J ) - ONE WORK( ( J-1 )*M+J ) = WORK( ( J-1 )*M+J ) - ONE
10 CONTINUE END DO
* *
CQRT11 = CLANGE( 'One-norm', M, M, WORK, M, RDUMMY ) / CQRT11 = CLANGE( 'One-norm', M, M, WORK, M, RDUMMY ) /
$ ( REAL( M )*SLAMCH( 'Epsilon' ) ) $ ( REAL( M )*SLAMCH( 'Epsilon' ) )

View File

@ -28,7 +28,7 @@
*> CQRT12 computes the singular values `svlues' of the upper trapezoid *> CQRT12 computes the singular values `svlues' of the upper trapezoid
*> of A(1:M,1:N) and returns the ratio *> of A(1:M,1:N) and returns the ratio
*> *>
*> || s - svlues||/(||svlues||*eps*max(M,N)) *> || svlues -s ||/( ||s||*eps*max(M,N) )
*> \endverbatim *> \endverbatim
* *
* Arguments: * Arguments:
@ -125,8 +125,8 @@
EXTERNAL CLANGE, SASUM, SLAMCH, SNRM2 EXTERNAL CLANGE, SASUM, SLAMCH, SNRM2
* .. * ..
* .. External Subroutines .. * .. External Subroutines ..
EXTERNAL CGEBD2, CLASCL, CLASET, SAXPY, SBDSQR, SLABAD, EXTERNAL CGEBD2, CLASCL, CLASET, SAXPY, SBDSQR, SLASCL,
$ SLASCL, XERBLA $ XERBLA
* .. * ..
* .. Intrinsic Functions .. * .. Intrinsic Functions ..
INTRINSIC CMPLX, MAX, MIN, REAL INTRINSIC CMPLX, MAX, MIN, REAL
@ -153,17 +153,16 @@
* Copy upper triangle of A into work * Copy upper triangle of A into work
* *
CALL CLASET( 'Full', M, N, CMPLX( ZERO ), CMPLX( ZERO ), WORK, M ) CALL CLASET( 'Full', M, N, CMPLX( ZERO ), CMPLX( ZERO ), WORK, M )
DO 20 J = 1, N DO J = 1, N
DO 10 I = 1, MIN( J, M ) DO I = 1, MIN( J, M )
WORK( ( J-1 )*M+I ) = A( I, J ) WORK( ( J-1 )*M+I ) = A( I, J )
10 CONTINUE END DO
20 CONTINUE END DO
* *
* Get machine parameters * Get machine parameters
* *
SMLNUM = SLAMCH( 'S' ) / SLAMCH( 'P' ) SMLNUM = SLAMCH( 'S' ) / SLAMCH( 'P' )
BIGNUM = ONE / SMLNUM BIGNUM = ONE / SMLNUM
CALL SLABAD( SMLNUM, BIGNUM )
* *
* Scale work if max entry outside range [SMLNUM,BIGNUM] * Scale work if max entry outside range [SMLNUM,BIGNUM]
* *
@ -207,9 +206,9 @@
* *
ELSE ELSE
* *
DO 30 I = 1, MN DO I = 1, MN
RWORK( I ) = ZERO RWORK( I ) = ZERO
30 CONTINUE END DO
END IF END IF
* *
* Compare s and singular values of work * Compare s and singular values of work

View File

@ -63,6 +63,7 @@
*> DLQ 8 List types on next line if 0 < NTYPES < 8 *> DLQ 8 List types on next line if 0 < NTYPES < 8
*> DQL 8 List types on next line if 0 < NTYPES < 8 *> DQL 8 List types on next line if 0 < NTYPES < 8
*> DQP 6 List types on next line if 0 < NTYPES < 6 *> DQP 6 List types on next line if 0 < NTYPES < 6
*> DQK 19 List types on next line if 0 < NTYPES < 19
*> DTZ 3 List types on next line if 0 < NTYPES < 3 *> DTZ 3 List types on next line if 0 < NTYPES < 3
*> DLS 6 List types on next line if 0 < NTYPES < 6 *> DLS 6 List types on next line if 0 < NTYPES < 6
*> DEQ *> DEQ
@ -149,11 +150,11 @@
$ NBVAL( MAXIN ), NBVAL2( MAXIN ), $ NBVAL( MAXIN ), NBVAL2( MAXIN ),
$ NSVAL( MAXIN ), NVAL( MAXIN ), NXVAL( MAXIN ), $ NSVAL( MAXIN ), NVAL( MAXIN ), NXVAL( MAXIN ),
$ RANKVAL( MAXIN ), PIV( NMAX ) $ RANKVAL( MAXIN ), PIV( NMAX )
DOUBLE PRECISION E( NMAX ), S( 2*NMAX )
* .. * ..
* .. Allocatable Arrays .. * .. Allocatable Arrays ..
INTEGER AllocateStatus INTEGER AllocateStatus
DOUBLE PRECISION, DIMENSION(:), ALLOCATABLE :: RWORK DOUBLE PRECISION, DIMENSION(:), ALLOCATABLE :: RWORK, S
DOUBLE PRECISION, DIMENSION(:), ALLOCATABLE :: E
DOUBLE PRECISION, DIMENSION(:,:), ALLOCATABLE :: A, B, WORK DOUBLE PRECISION, DIMENSION(:,:), ALLOCATABLE :: A, B, WORK
* .. * ..
* .. External Functions .. * .. External Functions ..
@ -164,13 +165,13 @@
* .. External Subroutines .. * .. External Subroutines ..
EXTERNAL ALAREQ, DCHKEQ, DCHKGB, DCHKGE, DCHKGT, DCHKLQ, EXTERNAL ALAREQ, DCHKEQ, DCHKGB, DCHKGE, DCHKGT, DCHKLQ,
$ DCHKORHR_COL, DCHKPB, DCHKPO, DCHKPS, DCHKPP, $ DCHKORHR_COL, DCHKPB, DCHKPO, DCHKPS, DCHKPP,
$ DCHKPT, DCHKQ3, DCHKQL, DCHKQR, DCHKRQ, DCHKSP, $ DCHKPT, DCHKQ3, DCHKQP3RK, DCHKQL, DCHKQR,
$ DCHKSY, DCHKSY_ROOK, DCHKSY_RK, DCHKSY_AA, $ DCHKRQ, DCHKSP, DCHKSY, DCHKSY_ROOK, DCHKSY_RK,
$ DCHKTB, DCHKTP, DCHKTR, DCHKTZ, DDRVGB, DDRVGE, $ DCHKSY_AA, DCHKTB, DCHKTP, DCHKTR, DCHKTZ,
$ DDRVGT, DDRVLS, DDRVPB, DDRVPO, DDRVPP, DDRVPT, $ DDRVGB, DDRVGE, DDRVGT, DDRVLS, DDRVPB, DDRVPO,
$ DDRVSP, DDRVSY, DDRVSY_ROOK, DDRVSY_RK, $ DDRVPP, DDRVPT, DDRVSP, DDRVSY, DDRVSY_ROOK,
$ DDRVSY_AA, ILAVER, DCHKLQTP, DCHKQRT, DCHKQRTP, $ DDRVSY_RK, DDRVSY_AA, ILAVER, DCHKLQTP, DCHKQRT,
$ DCHKLQT,DCHKTSQR $ DCHKQRTP, DCHKLQT,DCHKTSQR
* .. * ..
* .. Scalars in Common .. * .. Scalars in Common ..
LOGICAL LERR, OK LOGICAL LERR, OK
@ -197,6 +198,10 @@
IF (AllocateStatus /= 0) STOP "*** Not enough memory ***" IF (AllocateStatus /= 0) STOP "*** Not enough memory ***"
ALLOCATE ( WORK( NMAX, 3*NMAX+MAXRHS+30 ), STAT = AllocateStatus ) ALLOCATE ( WORK( NMAX, 3*NMAX+MAXRHS+30 ), STAT = AllocateStatus )
IF (AllocateStatus /= 0) STOP "*** Not enough memory ***" IF (AllocateStatus /= 0) STOP "*** Not enough memory ***"
ALLOCATE ( E( NMAX ), STAT = AllocateStatus )
IF (AllocateStatus /= 0) STOP "*** Not enough memory ***"
ALLOCATE ( S( 2*NMAX ), STAT = AllocateStatus )
IF (AllocateStatus /= 0) STOP "*** Not enough memory ***"
ALLOCATE ( RWORK( 5*NMAX+2*MAXRHS ), STAT = AllocateStatus ) ALLOCATE ( RWORK( 5*NMAX+2*MAXRHS ), STAT = AllocateStatus )
IF (AllocateStatus /= 0) STOP "*** Not enough memory ***" IF (AllocateStatus /= 0) STOP "*** Not enough memory ***"
* *
@ -919,9 +924,26 @@
CALL ALAREQ( PATH, NMATS, DOTYPE, NTYPES, NIN, NOUT ) CALL ALAREQ( PATH, NMATS, DOTYPE, NTYPES, NIN, NOUT )
* *
IF( TSTCHK ) THEN IF( TSTCHK ) THEN
CALL DCHKQ3( DOTYPE, NM, MVAL, NN, NVAL, NNB, NBVAL, NXVAL, CALL DCHKQ3( DOTYPE, NM, MVAL, NN, NVAL, NNB, NBVAL,
$ THRESH, A( 1, 1 ), A( 1, 2 ), B( 1, 1 ), $ NXVAL, THRESH, A( 1, 1 ), A( 1, 2 ),
$ B( 1, 3 ), WORK, IWORK, NOUT ) $ B( 1, 1 ), B( 1, 3 ), WORK, IWORK, NOUT )
ELSE
WRITE( NOUT, FMT = 9989 )PATH
END IF
*
ELSE IF( LSAMEN( 2, C2, 'QK' ) ) THEN
*
* QK: truncated QR factorization with pivoting
*
NTYPES = 19
CALL ALAREQ( PATH, NMATS, DOTYPE, NTYPES, NIN, NOUT )
*
IF( TSTCHK ) THEN
CALL DCHKQP3RK( DOTYPE, NM, MVAL, NN, NVAL, NNS, NSVAL,
$ NNB, NBVAL, NXVAL, THRESH, A( 1, 1 ),
$ A( 1, 2 ), B( 1, 1 ), B( 1, 2 ),
$ B( 1, 3 ), B( 1, 4 ),
$ WORK, IWORK, NOUT )
ELSE ELSE
WRITE( NOUT, FMT = 9989 )PATH WRITE( NOUT, FMT = 9989 )PATH
END IF END IF

View File

@ -0,0 +1,832 @@
*> \brief \b DCHKQP3RK
*
* =========== DOCUMENTATION ===========
*
* Online html documentation available at
* http://www.netlib.org/lapack/explore-html/
*
* Definition:
* ===========
*
* SUBROUTINE DCHKQP3RK( DOTYPE, NM, MVAL, NN, NVAL, NNS, NSVAL,
* $ NNB, NBVAL, NXVAL, THRESH, A, COPYA,
* $ B, COPYB, S, TAU,
* $ WORK, IWORK, NOUT )
* 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 NM, NN, NNS, NNB, NOUT
* DOUBLE PRECISION THRESH
* ..
* .. Array Arguments ..
* LOGICAL DOTYPE( * )
* INTEGER IWORK( * ), MVAL( * ), NBVAL( * ), NSVAL( * ),
* $ NVAL( * ), NXVAL( * )
* DOUBLE PRECISION A( * ), COPYA( * ), B( * ), COPYB( * ),
* $ S( * ), TAU( * ), WORK( * )
* ..
*
*
*> \par Purpose:
* =============
*>
*> \verbatim
*>
*> DCHKQP3RK tests DGEQP3RK.
*> \endverbatim
*
* Arguments:
* ==========
*
*> \param[in] DOTYPE
*> \verbatim
*> DOTYPE is LOGICAL array, dimension (NTYPES)
*> The matrix types to be used for testing. Matrices of type j
*> (for 1 <= j <= NTYPES) are used for testing if DOTYPE(j) =
*> .TRUE.; if DOTYPE(j) = .FALSE., then type j is not used.
*> \endverbatim
*>
*> \param[in] NM
*> \verbatim
*> NM is INTEGER
*> The number of values of M contained in the vector MVAL.
*> \endverbatim
*>
*> \param[in] MVAL
*> \verbatim
*> MVAL is INTEGER array, dimension (NM)
*> The values of the matrix row dimension M.
*> \endverbatim
*>
*> \param[in] NN
*> \verbatim
*> NN is INTEGER
*> The number of values of N contained in the vector NVAL.
*> \endverbatim
*>
*> \param[in] NVAL
*> \verbatim
*> NVAL is INTEGER array, dimension (NN)
*> The values of the matrix column dimension N.
*> \endverbatim
*>
*> \param[in] NNS
*> \verbatim
*> NNS is INTEGER
*> The number of values of NRHS contained in the vector NSVAL.
*> \endverbatim
*>
*> \param[in] NSVAL
*> \verbatim
*> NSVAL is INTEGER array, dimension (NNS)
*> The values of the number of right hand sides NRHS.
*> \endverbatim
*>
*> \param[in] NNB
*> \verbatim
*> NNB is INTEGER
*> The number of values of NB and NX contained in the
*> vectors NBVAL and NXVAL. The blocking parameters are used
*> in pairs (NB,NX).
*> \endverbatim
*>
*> \param[in] NBVAL
*> \verbatim
*> NBVAL is INTEGER array, dimension (NNB)
*> The values of the blocksize NB.
*> \endverbatim
*>
*> \param[in] NXVAL
*> \verbatim
*> NXVAL is INTEGER array, dimension (NNB)
*> The values of the crossover point NX.
*> \endverbatim
*>
*> \param[in] THRESH
*> \verbatim
*> THRESH is DOUBLE PRECISION
*> The threshold value for the test ratios. A result is
*> included in the output file if RESULT >= THRESH. To have
*> every test ratio printed, use THRESH = 0.
*> \endverbatim
*>
*> \param[out] A
*> \verbatim
*> A is DOUBLE PRECISION array, dimension (MMAX*NMAX)
*> where MMAX is the maximum value of M in MVAL and NMAX is the
*> maximum value of N in NVAL.
*> \endverbatim
*>
*> \param[out] COPYA
*> \verbatim
*> COPYA is DOUBLE PRECISION array, dimension (MMAX*NMAX)
*> \endverbatim
*>
*> \param[out] B
*> \verbatim
*> B is DOUBLE PRECISION array, dimension (MMAX*NSMAX)
*> where MMAX is the maximum value of M in MVAL and NSMAX is the
*> maximum value of NRHS in NSVAL.
*> \endverbatim
*>
*> \param[out] COPYB
*> \verbatim
*> COPYB is DOUBLE PRECISION array, dimension (MMAX*NSMAX)
*> \endverbatim
*>
*> \param[out] S
*> \verbatim
*> S is DOUBLE PRECISION array, dimension
*> (min(MMAX,NMAX))
*> \endverbatim
*>
*> \param[out] TAU
*> \verbatim
*> TAU is DOUBLE PRECISION array, dimension (MMAX)
*> \endverbatim
*>
*> \param[out] WORK
*> \verbatim
*> WORK is DOUBLE PRECISION array, dimension
*> (MMAX*NMAX + 4*NMAX + MMAX)
*> \endverbatim
*>
*> \param[out] IWORK
*> \verbatim
*> IWORK is INTEGER array, dimension (2*NMAX)
*> \endverbatim
*>
*> \param[in] NOUT
*> \verbatim
*> NOUT is INTEGER
*> The unit number for output.
*> \endverbatim
*
* Authors:
* ========
*
*> \author Univ. of Tennessee
*> \author Univ. of California Berkeley
*> \author Univ. of Colorado Denver
*> \author NAG Ltd.
*
*> \ingroup double_lin
*
* =====================================================================
SUBROUTINE DCHKQP3RK( DOTYPE, NM, MVAL, NN, NVAL, NNS, NSVAL,
$ NNB, NBVAL, NXVAL, THRESH, A, COPYA,
$ B, COPYB, S, TAU,
$ WORK, IWORK, NOUT )
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 NM, NN, NNB, NNS, NOUT
DOUBLE PRECISION THRESH
* ..
* .. Array Arguments ..
LOGICAL DOTYPE( * )
INTEGER IWORK( * ), NBVAL( * ), MVAL( * ), NVAL( * ),
$ NSVAL( * ), NXVAL( * )
DOUBLE PRECISION A( * ), COPYA( * ), B( * ), COPYB( * ),
$ S( * ), TAU( * ), WORK( * )
* ..
*
* =====================================================================
*
* .. Parameters ..
INTEGER NTYPES
PARAMETER ( NTYPES = 19 )
INTEGER NTESTS
PARAMETER ( NTESTS = 5 )
DOUBLE PRECISION ONE, ZERO, BIGNUM
PARAMETER ( ONE = 1.0D+0, ZERO = 0.0D+0,
$ BIGNUM = 1.0D+38 )
* ..
* .. Local Scalars ..
CHARACTER DIST, TYPE
CHARACTER*3 PATH
INTEGER I, IHIGH, ILOW, IM, IMAT, IN, INC_ZERO,
$ INB, IND_OFFSET_GEN,
$ IND_IN, IND_OUT, INS, INFO,
$ ISTEP, J, J_INC, J_FIRST_NZ, JB_ZERO,
$ KFACT, KL, KMAX, KU, LDA, LW, LWORK,
$ LWORK_MQR, M, MINMN, MINMNB_GEN, MODE, N,
$ NB, NB_ZERO, NERRS, NFAIL, NB_GEN, NRHS,
$ NRUN, NX, T
DOUBLE PRECISION ANORM, CNDNUM, EPS, ABSTOL, RELTOL,
$ DTEMP, MAXC2NRMK, RELMAXC2NRMK
* ..
* .. Local Arrays ..
INTEGER ISEED( 4 ), ISEEDY( 4 )
DOUBLE PRECISION RESULT( NTESTS ), RDUMMY( 1 )
* ..
* .. External Functions ..
DOUBLE PRECISION DLAMCH, DQPT01, DQRT11, DQRT12, DLANGE,
$ DLAPY2
EXTERNAL DLAMCH, DQPT01, DQRT11, DQRT12, DLANGE
* ..
* .. External Subroutines ..
EXTERNAL ALAERH, ALAHD, ALASUM, DAXPY, DGEQP3RK,
$ DLACPY, DLAORD, DLASET, DLATB4, DLATMS,
$ DORMQR, DSWAP, ICOPY, XLAENV
* ..
* .. Intrinsic Functions ..
INTRINSIC ABS, DBLE, MAX, MIN, MOD
* ..
* .. Scalars in Common ..
LOGICAL LERR, OK
CHARACTER*32 SRNAMT
INTEGER INFOT, IOUNIT
* ..
* .. Common blocks ..
COMMON / INFOC / INFOT, IOUNIT, OK, LERR
COMMON / SRNAMC / SRNAMT
* ..
* .. Data statements ..
DATA ISEEDY / 1988, 1989, 1990, 1991 /
* ..
* .. Executable Statements ..
*
* Initialize constants and the random number seed.
*
PATH( 1: 1 ) = 'Double precision'
PATH( 2: 3 ) = 'QK'
NRUN = 0
NFAIL = 0
NERRS = 0
DO I = 1, 4
ISEED( I ) = ISEEDY( I )
END DO
EPS = DLAMCH( 'Epsilon' )
INFOT = 0
*
DO IM = 1, NM
*
* Do for each value of M in MVAL.
*
M = MVAL( IM )
LDA = MAX( 1, M )
*
DO IN = 1, NN
*
* Do for each value of N in NVAL.
*
N = NVAL( IN )
MINMN = MIN( M, N )
LWORK = MAX( 1, M*MAX( M, N )+4*MINMN+MAX( M, N ),
$ M*N + 2*MINMN + 4*N )
*
DO INS = 1, NNS
NRHS = NSVAL( INS )
*
* Set up parameters with DLATB4 and generate
* M-by-NRHS B matrix with DLATMS.
* IMAT = 14:
* Random matrix, CNDNUM = 2, NORM = ONE,
* MODE = 3 (geometric distribution of singular values).
*
CALL DLATB4( PATH, 14, M, NRHS, TYPE, KL, KU, ANORM,
$ MODE, CNDNUM, DIST )
*
SRNAMT = 'DLATMS'
CALL DLATMS( M, NRHS, DIST, ISEED, TYPE, S, MODE,
$ CNDNUM, ANORM, KL, KU, 'No packing',
$ COPYB, LDA, WORK, INFO )
*
* Check error code from DLATMS.
*
IF( INFO.NE.0 ) THEN
CALL ALAERH( PATH, 'DLATMS', INFO, 0, ' ', M,
$ NRHS, -1, -1, -1, 6, NFAIL, NERRS,
$ NOUT )
CYCLE
END IF
*
DO IMAT = 1, NTYPES
*
* Do the tests only if DOTYPE( IMAT ) is true.
*
IF( .NOT.DOTYPE( IMAT ) )
$ CYCLE
*
* The type of distribution used to generate the random
* eigen-/singular values:
* ( 'S' for symmetric distribution ) => UNIFORM( -1, 1 )
*
* Do for each type of NON-SYMMETRIC matrix: CNDNUM NORM MODE
* 1. Zero matrix
* 2. Random, Diagonal, CNDNUM = 2 CNDNUM = 2 ONE 3 ( geometric distribution of singular values )
* 3. Random, Upper triangular, CNDNUM = 2 CNDNUM = 2 ONE 3 ( geometric distribution of singular values )
* 4. Random, Lower triangular, CNDNUM = 2 CNDNUM = 2 ONE 3 ( geometric distribution of singular values )
* 5. Random, First column is zero, CNDNUM = 2 CNDNUM = 2 ONE 3 ( geometric distribution of singular values )
* 6. Random, Last MINMN column is zero, CNDNUM = 2 CNDNUM = 2 ONE 3 ( geometric distribution of singular values )
* 7. Random, Last N column is zero, CNDNUM = 2 CNDNUM = 2 ONE 3 ( geometric distribution of singular values )
* 8. Random, Middle column in MINMN is zero, CNDNUM = 2 CNDNUM = 2 ONE 3 ( geometric distribution of singular values )
* 9. Random, First half of MINMN columns are zero, CNDNUM = 2 CNDNUM = 2 ONE 3 ( geometric distribution of singular values )
* 10. Random, Last columns are zero starting from MINMN/2+1, CNDNUM = 2 CNDNUM = 2 ONE 3 ( geometric distribution of singular values )
* 11. Random, Half MINMN columns in the middle are zero starting
* from MINMN/2-(MINMN/2)/2+1, CNDNUM = 2 CNDNUM = 2 ONE 3 ( geometric distribution of singular values )
* 12. Random, Odd columns are ZERO, CNDNUM = 2 CNDNUM = 2 ONE 3 ( geometric distribution of singular values )
* 13. Random, Even columns are ZERO, CNDNUM = 2 CNDNUM = 2 ONE 3 ( geometric distribution of singular values )
* 14. Random, CNDNUM = 2 CNDNUM = 2 ONE 3 ( geometric distribution of singular values )
* 15. Random, CNDNUM = sqrt(0.1/EPS) CNDNUM = BADC1 = sqrt(0.1/EPS) ONE 3 ( geometric distribution of singular values )
* 16. Random, CNDNUM = 0.1/EPS CNDNUM = BADC2 = 0.1/EPS ONE 3 ( geometric distribution of singular values )
* 17. Random, CNDNUM = 0.1/EPS, CNDNUM = BADC2 = 0.1/EPS ONE 2 ( one small singular value, S(N)=1/CNDNUM )
* one small singular value S(N)=1/CNDNUM
* 18. Random, CNDNUM = 2, scaled near underflow CNDNUM = 2 SMALL = SAFMIN
* 19. Random, CNDNUM = 2, scaled near overflow CNDNUM = 2 LARGE = 1.0/( 0.25 * ( SAFMIN / EPS ) ) 3 ( geometric distribution of singular values )
*
IF( IMAT.EQ.1 ) THEN
*
* Matrix 1: Zero matrix
*
CALL DLASET( 'Full', M, N, ZERO, ZERO, COPYA, LDA )
DO I = 1, MINMN
S( I ) = ZERO
END DO
*
ELSE IF( (IMAT.GE.2 .AND. IMAT.LE.4 )
$ .OR. (IMAT.GE.14 .AND. IMAT.LE.19 ) ) THEN
*
* Matrices 2-5.
*
* Set up parameters with DLATB4 and generate a test
* matrix with DLATMS.
*
CALL DLATB4( PATH, IMAT, M, N, TYPE, KL, KU, ANORM,
$ MODE, CNDNUM, DIST )
*
SRNAMT = 'DLATMS'
CALL DLATMS( M, N, DIST, ISEED, TYPE, S, MODE,
$ CNDNUM, ANORM, KL, KU, 'No packing',
$ COPYA, LDA, WORK, INFO )
*
* Check error code from DLATMS.
*
IF( INFO.NE.0 ) THEN
CALL ALAERH( PATH, 'DLATMS', INFO, 0, ' ', M, N,
$ -1, -1, -1, IMAT, NFAIL, NERRS,
$ NOUT )
CYCLE
END IF
*
CALL DLAORD( 'Decreasing', MINMN, S, 1 )
*
ELSE IF( MINMN.GE.2
$ .AND. IMAT.GE.5 .AND. IMAT.LE.13 ) THEN
*
* Rectangular matrices 5-13 that contain zero columns,
* only for matrices MINMN >=2.
*
* JB_ZERO is the column index of ZERO block.
* NB_ZERO is the column block size of ZERO block.
* NB_GEN is the column blcok size of the
* generated block.
* J_INC in the non_zero column index increment
* for matrix 12 and 13.
* J_FIRS_NZ is the index of the first non-zero
* column.
*
IF( IMAT.EQ.5 ) THEN
*
* First column is zero.
*
JB_ZERO = 1
NB_ZERO = 1
NB_GEN = N - NB_ZERO
*
ELSE IF( IMAT.EQ.6 ) THEN
*
* Last column MINMN is zero.
*
JB_ZERO = MINMN
NB_ZERO = 1
NB_GEN = N - NB_ZERO
*
ELSE IF( IMAT.EQ.7 ) THEN
*
* Last column N is zero.
*
JB_ZERO = N
NB_ZERO = 1
NB_GEN = N - NB_ZERO
*
ELSE IF( IMAT.EQ.8 ) THEN
*
* Middle column in MINMN is zero.
*
JB_ZERO = MINMN / 2 + 1
NB_ZERO = 1
NB_GEN = N - NB_ZERO
*
ELSE IF( IMAT.EQ.9 ) THEN
*
* First half of MINMN columns is zero.
*
JB_ZERO = 1
NB_ZERO = MINMN / 2
NB_GEN = N - NB_ZERO
*
ELSE IF( IMAT.EQ.10 ) THEN
*
* Last columns are zero columns,
* starting from (MINMN / 2 + 1) column.
*
JB_ZERO = MINMN / 2 + 1
NB_ZERO = N - JB_ZERO + 1
NB_GEN = N - NB_ZERO
*
ELSE IF( IMAT.EQ.11 ) THEN
*
* Half of the columns in the middle of MINMN
* columns is zero, starting from
* MINMN/2 - (MINMN/2)/2 + 1 column.
*
JB_ZERO = MINMN / 2 - (MINMN / 2) / 2 + 1
NB_ZERO = MINMN / 2
NB_GEN = N - NB_ZERO
*
ELSE IF( IMAT.EQ.12 ) THEN
*
* Odd-numbered columns are zero,
*
NB_GEN = N / 2
NB_ZERO = N - NB_GEN
J_INC = 2
J_FIRST_NZ = 2
*
ELSE IF( IMAT.EQ.13 ) THEN
*
* Even-numbered columns are zero.
*
NB_ZERO = N / 2
NB_GEN = N - NB_ZERO
J_INC = 2
J_FIRST_NZ = 1
*
END IF
*
*
* 1) Set the first NB_ZERO columns in COPYA(1:M,1:N)
* to zero.
*
CALL DLASET( 'Full', M, NB_ZERO, ZERO, ZERO,
$ COPYA, LDA )
*
* 2) Generate an M-by-(N-NB_ZERO) matrix with the
* chosen singular value distribution
* in COPYA(1:M,NB_ZERO+1:N).
*
CALL DLATB4( PATH, IMAT, M, NB_GEN, TYPE, KL, KU,
$ ANORM, MODE, CNDNUM, DIST )
*
SRNAMT = 'DLATMS'
*
IND_OFFSET_GEN = NB_ZERO * LDA
*
CALL DLATMS( M, NB_GEN, DIST, ISEED, TYPE, S, MODE,
$ CNDNUM, ANORM, KL, KU, 'No packing',
$ COPYA( IND_OFFSET_GEN + 1 ), LDA,
$ WORK, INFO )
*
* Check error code from DLATMS.
*
IF( INFO.NE.0 ) THEN
CALL ALAERH( PATH, 'DLATMS', INFO, 0, ' ', M,
$ NB_GEN, -1, -1, -1, IMAT, NFAIL,
$ NERRS, NOUT )
CYCLE
END IF
*
* 3) Swap the gererated colums from the right side
* NB_GEN-size block in COPYA into correct column
* positions.
*
IF( IMAT.EQ.6
$ .OR. IMAT.EQ.7
$ .OR. IMAT.EQ.8
$ .OR. IMAT.EQ.10
$ .OR. IMAT.EQ.11 ) THEN
*
* Move by swapping the generated columns
* from the right NB_GEN-size block from
* (NB_ZERO+1:NB_ZERO+JB_ZERO)
* into columns (1:JB_ZERO-1).
*
DO J = 1, JB_ZERO-1, 1
CALL DSWAP( M,
$ COPYA( ( NB_ZERO+J-1)*LDA+1), 1,
$ COPYA( (J-1)*LDA + 1 ), 1 )
END DO
*
ELSE IF( IMAT.EQ.12 .OR. IMAT.EQ.13 ) THEN
*
* ( IMAT = 12, Odd-numbered ZERO columns. )
* Swap the generated columns from the right
* NB_GEN-size block into the even zero colums in the
* left NB_ZERO-size block.
*
* ( IMAT = 13, Even-numbered ZERO columns. )
* Swap the generated columns from the right
* NB_GEN-size block into the odd zero colums in the
* left NB_ZERO-size block.
*
DO J = 1, NB_GEN, 1
IND_OUT = ( NB_ZERO+J-1 )*LDA + 1
IND_IN = ( J_INC*(J-1)+(J_FIRST_NZ-1) )*LDA
$ + 1
CALL DSWAP( M,
$ COPYA( IND_OUT ), 1,
$ COPYA( IND_IN), 1 )
END DO
*
END IF
*
* 5) Order the singular values generated by
* DLAMTS in decreasing order and add trailing zeros
* that correspond to zero columns.
* The total number of singular values is MINMN.
*
MINMNB_GEN = MIN( M, NB_GEN )
*
DO I = MINMNB_GEN+1, MINMN
S( I ) = ZERO
END DO
*
ELSE
*
* IF(MINMN.LT.2) skip this size for this matrix type.
*
CYCLE
END IF
*
* Initialize a copy array for a pivot array for DGEQP3RK.
*
DO I = 1, N
IWORK( I ) = 0
END DO
*
DO INB = 1, NNB
*
* Do for each pair of values (NB,NX) in NBVAL and NXVAL.
*
NB = NBVAL( INB )
CALL XLAENV( 1, NB )
NX = NXVAL( INB )
CALL XLAENV( 3, NX )
*
* We do MIN(M,N)+1 because we need a test for KMAX > N,
* when KMAX is larger than MIN(M,N), KMAX should be
* KMAX = MIN(M,N)
*
DO KMAX = 0, MIN(M,N)+1
*
* Get a working copy of COPYA into A( 1:M,1:N ).
* Get a working copy of COPYB into A( 1:M, (N+1):NRHS ).
* Get a working copy of COPYB into into B( 1:M, 1:NRHS ).
* Get a working copy of IWORK(1:N) awith zeroes into
* which is going to be used as pivot array IWORK( N+1:2N ).
* NOTE: IWORK(2N+1:3N) is going to be used as a WORK array
* for the routine.
*
CALL DLACPY( 'All', M, N, COPYA, LDA, A, LDA )
CALL DLACPY( 'All', M, NRHS, COPYB, LDA,
$ A( LDA*N + 1 ), LDA )
CALL DLACPY( 'All', M, NRHS, COPYB, LDA,
$ B, LDA )
CALL ICOPY( N, IWORK( 1 ), 1, IWORK( N+1 ), 1 )
*
ABSTOL = -1.0
RELTOL = -1.0
*
* Compute the QR factorization with pivoting of A
*
LW = MAX( 1, MAX( 2*N + NB*( N+NRHS+1 ),
$ 3*N + NRHS - 1 ) )
*
* Compute DGEQP3RK factorization of A.
*
SRNAMT = 'DGEQP3RK'
CALL DGEQP3RK( M, N, NRHS, KMAX, ABSTOL, RELTOL,
$ A, LDA, KFACT, MAXC2NRMK,
$ RELMAXC2NRMK, IWORK( N+1 ), TAU,
$ WORK, LW, IWORK( 2*N+1 ), INFO )
*
* Check error code from DGEQP3RK.
*
IF( INFO.LT.0 )
$ CALL ALAERH( PATH, 'DGEQP3RK', INFO, 0, ' ',
$ M, N, NX, -1, NB, IMAT,
$ NFAIL, NERRS, NOUT )
*
* Compute test 1:
*
* This test in only for the full rank factorization of
* the matrix A.
*
* Array S(1:min(M,N)) contains svd(A) the sigular values
* of the original matrix A in decreasing absolute value
* order. The test computes svd(R), the vector sigular
* values of the upper trapezoid of A(1:M,1:N) that
* contains the factor R, in decreasing order. The test
* returns the ratio:
*
* 2-norm(svd(R) - svd(A)) / ( max(M,N) * 2-norm(svd(A)) * EPS )
*
IF( KFACT.EQ.MINMN ) THEN
*
RESULT( 1 ) = DQRT12( M, N, A, LDA, S, WORK,
$ LWORK )
*
DO T = 1, 1
IF( RESULT( T ).GE.THRESH ) THEN
IF( NFAIL.EQ.0 .AND. NERRS.EQ.0 )
$ CALL ALAHD( NOUT, PATH )
WRITE( NOUT, FMT = 9999 ) 'DGEQP3RK', M, N,
$ NRHS, KMAX, ABSTOL, RELTOL, NB, NX,
$ IMAT, T, RESULT( T )
NFAIL = NFAIL + 1
END IF
END DO
NRUN = NRUN + 1
*
* End test 1
*
END IF
*
* Compute test 2:
*
* The test returns the ratio:
*
* 1-norm( A*P - Q*R ) / ( max(M,N) * 1-norm(A) * EPS )
*
RESULT( 2 ) = DQPT01( M, N, KFACT, COPYA, A, LDA, TAU,
$ IWORK( N+1 ), WORK, LWORK )
*
* Compute test 3:
*
* The test returns the ratio:
*
* 1-norm( Q**T * Q - I ) / ( M * EPS )
*
RESULT( 3 ) = DQRT11( M, KFACT, A, LDA, TAU, WORK,
$ LWORK )
*
* Print information about the tests that did not pass
* the threshold.
*
DO T = 2, 3
IF( RESULT( T ).GE.THRESH ) THEN
IF( NFAIL.EQ.0 .AND. NERRS.EQ.0 )
$ CALL ALAHD( NOUT, PATH )
WRITE( NOUT, FMT = 9999 ) 'DGEQP3RK', M, N,
$ NRHS, KMAX, ABSTOL, RELTOL,
$ NB, NX, IMAT, T, RESULT( T )
NFAIL = NFAIL + 1
END IF
END DO
NRUN = NRUN + 2
*
* Compute test 4:
*
* This test is only for the factorizations with the
* rank greater than 2.
* The elements on the diagonal of R should be non-
* increasing.
*
* The test returns the ratio:
*
* Returns 1.0D+100 if abs(R(K+1,K+1)) > abs(R(K,K)),
* K=1:KFACT-1
*
IF( MIN(KFACT, MINMN).GE.2 ) THEN
*
DO J = 1, KFACT-1, 1
DTEMP = (( ABS( A( (J-1)*M+J ) ) -
$ ABS( A( (J)*M+J+1 ) ) ) /
$ ABS( A(1) ) )
*
IF( DTEMP.LT.ZERO ) THEN
RESULT( 4 ) = BIGNUM
END IF
*
END DO
*
* Print information about the tests that did not
* pass the threshold.
*
DO T = 4, 4
IF( RESULT( T ).GE.THRESH ) THEN
IF( NFAIL.EQ.0 .AND. NERRS.EQ.0 )
$ CALL ALAHD( NOUT, PATH )
WRITE( NOUT, FMT = 9999 ) 'DGEQP3RK',
$ M, N, NRHS, KMAX, ABSTOL, RELTOL,
$ NB, NX, IMAT, T,
$ RESULT( T )
NFAIL = NFAIL + 1
END IF
END DO
NRUN = NRUN + 1
*
* End test 4.
*
END IF
*
* Compute test 5:
*
* This test in only for matrix A with min(M,N) > 0.
*
* The test returns the ratio:
*
* 1-norm(Q**T * B - Q**T * B ) /
* ( M * EPS )
*
* (1) Compute B:=Q**T * B in the matrix B.
*
IF( MINMN.GT.0 ) THEN
*
LWORK_MQR = MAX(1, NRHS)
CALL DORMQR( 'Left', 'Transpose',
$ M, NRHS, KFACT, A, LDA, TAU, B, LDA,
$ WORK, LWORK_MQR, INFO )
*
DO I = 1, NRHS
*
* Compare N+J-th column of A and J-column of B.
*
CALL DAXPY( M, -ONE, A( ( N+I-1 )*LDA+1 ), 1,
$ B( ( I-1 )*LDA+1 ), 1 )
END DO
*
RESULT( 5 ) =
$ ABS(
$ DLANGE( 'One-norm', M, NRHS, B, LDA, RDUMMY ) /
$ ( DBLE( M )*DLAMCH( 'Epsilon' ) )
$ )
*
* Print information about the tests that did not pass
* the threshold.
*
DO T = 5, 5
IF( RESULT( T ).GE.THRESH ) THEN
IF( NFAIL.EQ.0 .AND. NERRS.EQ.0 )
$ CALL ALAHD( NOUT, PATH )
WRITE( NOUT, FMT = 9999 ) 'DGEQP3RK', M, N,
$ NRHS, KMAX, ABSTOL, RELTOL,
$ NB, NX, IMAT, T, RESULT( T )
NFAIL = NFAIL + 1
END IF
END DO
NRUN = NRUN + 1
*
* End compute test 5.
*
END IF
*
* END DO KMAX = 1, MIN(M,N)+1
*
END DO
*
* END DO for INB = 1, NNB
*
END DO
*
* END DO for IMAT = 1, NTYPES
*
END DO
*
* END DO for INS = 1, NNS
*
END DO
*
* END DO for IN = 1, NN
*
END DO
*
* END DO for IM = 1, NM
*
END DO
*
* Print a summary of the results.
*
CALL ALASUM( PATH, NOUT, NFAIL, NRUN, NERRS )
*
9999 FORMAT( 1X, A, ' M =', I5, ', N =', I5, ', NRHS =', I5,
$ ', KMAX =', I5, ', ABSTOL =', G12.5,
$ ', RELTOL =', G12.5, ', NB =', I4, ', NX =', I4,
$ ', type ', I2, ', test ', I2, ', ratio =', G12.5 )
*
* End of DCHKQP3RK
*
END

View File

@ -133,7 +133,7 @@
* *
* .. Parameters .. * .. Parameters ..
DOUBLE PRECISION SHRINK, TENTH DOUBLE PRECISION SHRINK, TENTH
PARAMETER ( SHRINK = 0.25D0, TENTH = 0.1D+0 ) PARAMETER ( SHRINK = 0.25D+0, TENTH = 0.1D+0 )
DOUBLE PRECISION ONE DOUBLE PRECISION ONE
PARAMETER ( ONE = 1.0D+0 ) PARAMETER ( ONE = 1.0D+0 )
DOUBLE PRECISION TWO DOUBLE PRECISION TWO
@ -153,9 +153,6 @@
* .. Intrinsic Functions .. * .. Intrinsic Functions ..
INTRINSIC ABS, MAX, SQRT INTRINSIC ABS, MAX, SQRT
* .. * ..
* .. External Subroutines ..
EXTERNAL DLABAD
* ..
* .. Save statement .. * .. Save statement ..
SAVE EPS, SMALL, LARGE, BADC1, BADC2, FIRST SAVE EPS, SMALL, LARGE, BADC1, BADC2, FIRST
* .. * ..
@ -173,11 +170,6 @@
BADC1 = SQRT( BADC2 ) BADC1 = SQRT( BADC2 )
SMALL = DLAMCH( 'Safe minimum' ) SMALL = DLAMCH( 'Safe minimum' )
LARGE = ONE / SMALL LARGE = ONE / SMALL
*
* If it looks like we're on a Cray, take the square root of
* SMALL and LARGE to avoid overflow and underflow problems.
*
CALL DLABAD( SMALL, LARGE )
SMALL = SHRINK*( SMALL / EPS ) SMALL = SHRINK*( SMALL / EPS )
LARGE = ONE / SMALL LARGE = ONE / SMALL
END IF END IF
@ -232,6 +224,110 @@
ELSE ELSE
ANORM = ONE ANORM = ONE
END IF END IF
*
ELSE IF( LSAMEN( 2, C2, 'QK' ) ) THEN
*
* xQK: truncated QR with pivoting.
* Set parameters to generate a general
* M x N matrix.
*
* Set TYPE, the type of matrix to be generated. 'N' is nonsymmetric.
*
TYPE = 'N'
*
* Set DIST, the type of distribution for the random
* number generator. 'S' is
*
DIST = 'S'
*
* Set the lower and upper bandwidths.
*
IF( IMAT.EQ.2 ) THEN
*
* 2. Random, Diagonal, CNDNUM = 2
*
KL = 0
KU = 0
CNDNUM = TWO
ANORM = ONE
MODE = 3
ELSE IF( IMAT.EQ.3 ) THEN
*
* 3. Random, Upper triangular, CNDNUM = 2
*
KL = 0
KU = MAX( N-1, 0 )
CNDNUM = TWO
ANORM = ONE
MODE = 3
ELSE IF( IMAT.EQ.4 ) THEN
*
* 4. Random, Lower triangular, CNDNUM = 2
*
KL = MAX( M-1, 0 )
KU = 0
CNDNUM = TWO
ANORM = ONE
MODE = 3
ELSE
*
* 5.-19. Rectangular matrix
*
KL = MAX( M-1, 0 )
KU = MAX( N-1, 0 )
*
IF( IMAT.GE.5 .AND. IMAT.LE.14 ) THEN
*
* 5.-14. Random, CNDNUM = 2.
*
CNDNUM = TWO
ANORM = ONE
MODE = 3
*
ELSE IF( IMAT.EQ.15 ) THEN
*
* 15. Random, CNDNUM = sqrt(0.1/EPS)
*
CNDNUM = BADC1
ANORM = ONE
MODE = 3
*
ELSE IF( IMAT.EQ.16 ) THEN
*
* 16. Random, CNDNUM = 0.1/EPS
*
CNDNUM = BADC2
ANORM = ONE
MODE = 3
*
ELSE IF( IMAT.EQ.17 ) THEN
*
* 17. Random, CNDNUM = 0.1/EPS,
* one small singular value S(N)=1/CNDNUM
*
CNDNUM = BADC2
ANORM = ONE
MODE = 2
*
ELSE IF( IMAT.EQ.18 ) THEN
*
* 18. Random, scaled near underflow
*
CNDNUM = TWO
ANORM = SMALL
MODE = 3
*
ELSE IF( IMAT.EQ.19 ) THEN
*
* 19. Random, scaled near overflow
*
CNDNUM = TWO
ANORM = LARGE
MODE = 3
*
END IF
*
END IF
* *
ELSE IF( LSAMEN( 2, C2, 'GE' ) ) THEN ELSE IF( LSAMEN( 2, C2, 'GE' ) ) THEN
* *
@ -518,17 +614,18 @@
* *
* Set the norm and condition number. * Set the norm and condition number.
* *
IF( IMAT.EQ.2 .OR. IMAT.EQ.8 ) THEN MAT = ABS( IMAT )
IF( MAT.EQ.2 .OR. MAT.EQ.8 ) THEN
CNDNUM = BADC1 CNDNUM = BADC1
ELSE IF( IMAT.EQ.3 .OR. IMAT.EQ.9 ) THEN ELSE IF( MAT.EQ.3 .OR. MAT.EQ.9 ) THEN
CNDNUM = BADC2 CNDNUM = BADC2
ELSE ELSE
CNDNUM = TWO CNDNUM = TWO
END IF END IF
* *
IF( IMAT.EQ.4 ) THEN IF( MAT.EQ.4 ) THEN
ANORM = SMALL ANORM = SMALL
ELSE IF( IMAT.EQ.5 ) THEN ELSE IF( MAT.EQ.5 ) THEN
ANORM = LARGE ANORM = LARGE
ELSE ELSE
ANORM = ONE ANORM = ONE

View File

@ -28,12 +28,13 @@
*> *>
*> DQPT01 tests the QR-factorization with pivoting of a matrix A. The *> DQPT01 tests the QR-factorization with pivoting of a matrix A. The
*> array AF contains the (possibly partial) QR-factorization of A, where *> array AF contains the (possibly partial) QR-factorization of A, where
*> the upper triangle of AF(1:k,1:k) is a partial triangular factor, *> the upper triangle of AF(1:K,1:K) is a partial triangular factor,
*> the entries below the diagonal in the first k columns are the *> the entries below the diagonal in the first K columns are the
*> Householder vectors, and the rest of AF contains a partially updated *> Householder vectors, and the rest of AF contains a partially updated
*> matrix. *> matrix.
*> *>
*> This function returns ||A*P - Q*R||/(||norm(A)||*eps*M) *> This function returns ||A*P - Q*R|| / ( ||norm(A)||*eps*max(M,N) ),
*> where || . || is matrix one norm.
*> \endverbatim *> \endverbatim
* *
* Arguments: * Arguments:
@ -172,28 +173,41 @@
* *
NORMA = DLANGE( 'One-norm', M, N, A, LDA, RWORK ) NORMA = DLANGE( 'One-norm', M, N, A, LDA, RWORK )
* *
DO 30 J = 1, K DO J = 1, K
DO 10 I = 1, MIN( J, M ) *
* Copy the upper triangular part of the factor R stored
* in AF(1:K,1:K) into the work array WORK.
*
DO I = 1, MIN( J, M )
WORK( ( J-1 )*M+I ) = AF( I, J ) WORK( ( J-1 )*M+I ) = AF( I, J )
10 CONTINUE END DO
DO 20 I = J + 1, M *
* Zero out the elements below the diagonal in the work array.
*
DO I = J + 1, M
WORK( ( J-1 )*M+I ) = ZERO WORK( ( J-1 )*M+I ) = ZERO
20 CONTINUE END DO
30 CONTINUE END DO
DO 40 J = K + 1, N *
* Copy columns (K+1,N) from AF into the work array WORK.
* AF(1:K,K+1:N) contains the rectangular block of the upper trapezoidal
* factor R, AF(K+1:M,K+1:N) contains the partially updated residual
* matrix of R.
*
DO J = K + 1, N
CALL DCOPY( M, AF( 1, J ), 1, WORK( ( J-1 )*M+1 ), 1 ) CALL DCOPY( M, AF( 1, J ), 1, WORK( ( J-1 )*M+1 ), 1 )
40 CONTINUE END DO
* *
CALL DORMQR( 'Left', 'No transpose', M, N, K, AF, LDA, TAU, WORK, CALL DORMQR( 'Left', 'No transpose', M, N, K, AF, LDA, TAU, WORK,
$ M, WORK( M*N+1 ), LWORK-M*N, INFO ) $ M, WORK( M*N+1 ), LWORK-M*N, INFO )
* *
DO 50 J = 1, N DO J = 1, N
* *
* Compare i-th column of QR and jpvt(i)-th column of A * Compare J-th column of QR and JPVT(J)-th column of A.
* *
CALL DAXPY( M, -ONE, A( 1, JPVT( J ) ), 1, WORK( ( J-1 )*M+1 ), CALL DAXPY( M, -ONE, A( 1, JPVT( J ) ), 1, WORK( ( J-1 )*M+1 ),
$ 1 ) $ 1 )
50 CONTINUE END DO
* *
DQPT01 = DLANGE( 'One-norm', M, N, WORK, M, RWORK ) / DQPT01 = DLANGE( 'One-norm', M, N, WORK, M, RWORK ) /
$ ( DBLE( MAX( M, N ) )*DLAMCH( 'Epsilon' ) ) $ ( DBLE( MAX( M, N ) )*DLAMCH( 'Epsilon' ) )

View File

@ -157,9 +157,9 @@
CALL DORM2R( 'Left', 'Transpose', M, M, K, A, LDA, TAU, WORK, M, CALL DORM2R( 'Left', 'Transpose', M, M, K, A, LDA, TAU, WORK, M,
$ WORK( M*M+1 ), INFO ) $ WORK( M*M+1 ), INFO )
* *
DO 10 J = 1, M DO J = 1, M
WORK( ( J-1 )*M+J ) = WORK( ( J-1 )*M+J ) - ONE WORK( ( J-1 )*M+J ) = WORK( ( J-1 )*M+J ) - ONE
10 CONTINUE END DO
* *
DQRT11 = DLANGE( 'One-norm', M, M, WORK, M, RDUMMY ) / DQRT11 = DLANGE( 'One-norm', M, M, WORK, M, RDUMMY ) /
$ ( DBLE( M )*DLAMCH( 'Epsilon' ) ) $ ( DBLE( M )*DLAMCH( 'Epsilon' ) )

View File

@ -26,7 +26,7 @@
*> DQRT12 computes the singular values `svlues' of the upper trapezoid *> DQRT12 computes the singular values `svlues' of the upper trapezoid
*> of A(1:M,1:N) and returns the ratio *> of A(1:M,1:N) and returns the ratio
*> *>
*> || s - svlues||/(||svlues||*eps*max(M,N)) *> || svlues - s ||/(||s||*eps*max(M,N))
*> \endverbatim *> \endverbatim
* *
* Arguments: * Arguments:
@ -113,8 +113,7 @@
EXTERNAL DASUM, DLAMCH, DLANGE, DNRM2 EXTERNAL DASUM, DLAMCH, DLANGE, DNRM2
* .. * ..
* .. External Subroutines .. * .. External Subroutines ..
EXTERNAL DAXPY, DBDSQR, DGEBD2, DLABAD, DLASCL, DLASET, EXTERNAL DAXPY, DBDSQR, DGEBD2, DLASCL, DLASET, XERBLA
$ XERBLA
* .. * ..
* .. Intrinsic Functions .. * .. Intrinsic Functions ..
INTRINSIC DBLE, MAX, MIN INTRINSIC DBLE, MAX, MIN
@ -145,17 +144,16 @@
* Copy upper triangle of A into work * Copy upper triangle of A into work
* *
CALL DLASET( 'Full', M, N, ZERO, ZERO, WORK, M ) CALL DLASET( 'Full', M, N, ZERO, ZERO, WORK, M )
DO 20 J = 1, N DO J = 1, N
DO 10 I = 1, MIN( J, M ) DO I = 1, MIN( J, M )
WORK( ( J-1 )*M+I ) = A( I, J ) WORK( ( J-1 )*M+I ) = A( I, J )
10 CONTINUE END DO
20 CONTINUE END DO
* *
* Get machine parameters * Get machine parameters
* *
SMLNUM = DLAMCH( 'S' ) / DLAMCH( 'P' ) SMLNUM = DLAMCH( 'S' ) / DLAMCH( 'P' )
BIGNUM = ONE / SMLNUM BIGNUM = ONE / SMLNUM
CALL DLABAD( SMLNUM, BIGNUM )
* *
* Scale work if max entry outside range [SMLNUM,BIGNUM] * Scale work if max entry outside range [SMLNUM,BIGNUM]
* *
@ -199,16 +197,18 @@
* *
ELSE ELSE
* *
DO 30 I = 1, MN DO I = 1, MN
WORK( M*N+I ) = ZERO WORK( M*N+I ) = ZERO
30 CONTINUE END DO
END IF END IF
* *
* Compare s and singular values of work * Compare s and singular values of work
* *
CALL DAXPY( MN, -ONE, S, 1, WORK( M*N+1 ), 1 ) CALL DAXPY( MN, -ONE, S, 1, WORK( M*N+1 ), 1 )
*
DQRT12 = DASUM( MN, WORK( M*N+1 ), 1 ) / DQRT12 = DASUM( MN, WORK( M*N+1 ), 1 ) /
$ ( DLAMCH( 'Epsilon' )*DBLE( MAX( M, N ) ) ) $ ( DLAMCH('Epsilon') * DBLE( MAX( M, N ) ) )
*
IF( NRMSVL.NE.ZERO ) IF( NRMSVL.NE.ZERO )
$ DQRT12 = DQRT12 / NRMSVL $ DQRT12 = DQRT12 / NRMSVL
* *

View File

@ -63,6 +63,7 @@
*> SLQ 8 List types on next line if 0 < NTYPES < 8 *> SLQ 8 List types on next line if 0 < NTYPES < 8
*> SQL 8 List types on next line if 0 < NTYPES < 8 *> SQL 8 List types on next line if 0 < NTYPES < 8
*> SQP 6 List types on next line if 0 < NTYPES < 6 *> SQP 6 List types on next line if 0 < NTYPES < 6
*> DQK 19 List types on next line if 0 < NTYPES < 19
*> STZ 3 List types on next line if 0 < NTYPES < 3 *> STZ 3 List types on next line if 0 < NTYPES < 3
*> SLS 6 List types on next line if 0 < NTYPES < 6 *> SLS 6 List types on next line if 0 < NTYPES < 6
*> SEQ *> SEQ
@ -147,11 +148,11 @@
$ NBVAL( MAXIN ), NBVAL2( MAXIN ), $ NBVAL( MAXIN ), NBVAL2( MAXIN ),
$ NSVAL( MAXIN ), NVAL( MAXIN ), NXVAL( MAXIN ), $ NSVAL( MAXIN ), NVAL( MAXIN ), NXVAL( MAXIN ),
$ RANKVAL( MAXIN ), PIV( NMAX ) $ RANKVAL( MAXIN ), PIV( NMAX )
REAL E( NMAX ), S( 2*NMAX )
* .. * ..
* .. Allocatable Arrays .. * .. Allocatable Arrays ..
INTEGER AllocateStatus INTEGER AllocateStatus
REAL, DIMENSION(:), ALLOCATABLE :: RWORK REAL, DIMENSION(:), ALLOCATABLE :: RWORK, S
REAL, DIMENSION(:), ALLOCATABLE :: E
REAL, DIMENSION(:,:), ALLOCATABLE :: A, B, WORK REAL, DIMENSION(:,:), ALLOCATABLE :: A, B, WORK
* .. * ..
* .. External Functions .. * .. External Functions ..
@ -162,13 +163,13 @@
* .. External Subroutines .. * .. External Subroutines ..
EXTERNAL ALAREQ, SCHKEQ, SCHKGB, SCHKGE, SCHKGT, SCHKLQ, EXTERNAL ALAREQ, SCHKEQ, SCHKGB, SCHKGE, SCHKGT, SCHKLQ,
$ SCHKORHR_COL, SCHKPB, SCHKPO, SCHKPS, SCHKPP, $ SCHKORHR_COL, SCHKPB, SCHKPO, SCHKPS, SCHKPP,
$ SCHKPT, SCHKQ3, SCHKQL, SCHKQR, SCHKRQ, SCHKSP, $ SCHKPT, SCHKQ3, SCHKQP3RK, SCHKQL, SCHKQR,
$ SCHKSY, SCHKSY_ROOK, SCHKSY_RK, SCHKSY_AA, $ SCHKRQ, SCHKSP, SCHKSY, SCHKSY_ROOK, SCHKSY_RK,
$ SCHKTB, SCHKTP, SCHKTR, SCHKTZ, SDRVGB, SDRVGE, $ SCHKSY_AA, SCHKTB, SCHKTP, SCHKTR, SCHKTZ,
$ SDRVGT, SDRVLS, SDRVPB, SDRVPO, SDRVPP, SDRVPT, $ SDRVGB, SDRVGE, SDRVGT, SDRVLS, SDRVPB, SDRVPO,
$ SDRVSP, SDRVSY, SDRVSY_ROOK, SDRVSY_RK, $ SDRVPP, SDRVPT, SDRVSP, SDRVSY, SDRVSY_ROOK,
$ SDRVSY_AA, ILAVER, SCHKLQTP, SCHKQRT, SCHKQRTP, $ SDRVSY_RK, SDRVSY_AA, ILAVER, SCHKLQTP, SCHKQRT,
$ SCHKLQT, SCHKTSQR $ SCHKQRTP, SCHKLQT, SCHKTSQR
* .. * ..
* .. Scalars in Common .. * .. Scalars in Common ..
LOGICAL LERR, OK LOGICAL LERR, OK
@ -188,13 +189,17 @@
* .. * ..
* .. Allocate memory dynamically .. * .. Allocate memory dynamically ..
* *
ALLOCATE (A( ( KDMAX+1 )*NMAX, 7 ), STAT = AllocateStatus ) ALLOCATE ( A( ( KDMAX+1 )*NMAX, 7 ), STAT = AllocateStatus )
IF (AllocateStatus /= 0) STOP "*** Not enough memory ***" IF (AllocateStatus /= 0) STOP "*** Not enough memory ***"
ALLOCATE (B( NMAX*MAXRHS, 4 ), STAT = AllocateStatus ) ALLOCATE ( B( NMAX*MAXRHS, 4 ), STAT = AllocateStatus )
IF (AllocateStatus /= 0) STOP "*** Not enough memory ***" IF (AllocateStatus /= 0) STOP "*** Not enough memory ***"
ALLOCATE (WORK( NMAX, NMAX+MAXRHS+30 ) , STAT = AllocateStatus ) ALLOCATE ( WORK( NMAX, 3*NMAX+MAXRHS+30 ), STAT = AllocateStatus )
IF (AllocateStatus /= 0) STOP "*** Not enough memory ***" IF (AllocateStatus /= 0) STOP "*** Not enough memory ***"
ALLOCATE (RWORK( 5*NMAX+2*MAXRHS ), STAT = AllocateStatus ) ALLOCATE ( E( NMAX ), STAT = AllocateStatus )
IF (AllocateStatus /= 0) STOP "*** Not enough memory ***"
ALLOCATE ( S( 2*NMAX ), STAT = AllocateStatus )
IF (AllocateStatus /= 0) STOP "*** Not enough memory ***"
ALLOCATE ( RWORK( 5*NMAX+2*MAXRHS ), STAT = AllocateStatus )
IF (AllocateStatus /= 0) STOP "*** Not enough memory ***" IF (AllocateStatus /= 0) STOP "*** Not enough memory ***"
* .. * ..
* .. Executable Statements .. * .. Executable Statements ..
@ -920,6 +925,23 @@
ELSE ELSE
WRITE( NOUT, FMT = 9989 )PATH WRITE( NOUT, FMT = 9989 )PATH
END IF END IF
*
ELSE IF( LSAMEN( 2, C2, 'QK' ) ) THEN
*
* QK: truncated QR factorization with pivoting
*
NTYPES = 19
CALL ALAREQ( PATH, NMATS, DOTYPE, NTYPES, NIN, NOUT )
*
IF( TSTCHK ) THEN
CALL SCHKQP3RK( DOTYPE, NM, MVAL, NN, NVAL, NNS, NSVAL,
$ NNB, NBVAL, NXVAL, THRESH, A( 1, 1 ),
$ A( 1, 2 ), B( 1, 1 ), B( 1, 2 ),
$ B( 1, 3 ), B( 1, 4 ),
$ WORK, IWORK, NOUT )
ELSE
WRITE( NOUT, FMT = 9989 )PATH
END IF
* *
ELSE IF( LSAMEN( 2, C2, 'TZ' ) ) THEN ELSE IF( LSAMEN( 2, C2, 'TZ' ) ) THEN
* *

View File

@ -0,0 +1,831 @@
*> \brief \b SCHKQP3RK
*
* =========== DOCUMENTATION ===========
*
* Online html documentation available at
* http://www.netlib.org/lapack/explore-html/
*
* Definition:
* ===========
*
* SUBROUTINE SCHKQP3RK( DOTYPE, NM, MVAL, NN, NVAL, NNS, NSVAL,
* $ NNB, NBVAL, NXVAL, THRESH, A, COPYA,
* $ B, COPYB, S, TAU,
* $ WORK, IWORK, NOUT )
* 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 NM, NN, NNS, NNB, NOUT
* REAL THRESH
* ..
* .. Array Arguments ..
* LOGICAL DOTYPE( * )
* INTEGER IWORK( * ), MVAL( * ), NBVAL( * ), NSVAL( * ),
* $ NVAL( * ), NXVAL( * )
* REAL A( * ), COPYA( * ), B( * ), COPYB( * ),
* $ S( * ), TAU( * ), WORK( * )
* ..
*
*
*> \par Purpose:
* =============
*>
*> \verbatim
*>
*> SCHKQP3RK tests SGEQP3RK.
*> \endverbatim
*
* Arguments:
* ==========
*
*> \param[in] DOTYPE
*> \verbatim
*> DOTYPE is LOGICAL array, dimension (NTYPES)
*> The matrix types to be used for testing. Matrices of type j
*> (for 1 <= j <= NTYPES) are used for testing if DOTYPE(j) =
*> .TRUE.; if DOTYPE(j) = .FALSE., then type j is not used.
*> \endverbatim
*>
*> \param[in] NM
*> \verbatim
*> NM is INTEGER
*> The number of values of M contained in the vector MVAL.
*> \endverbatim
*>
*> \param[in] MVAL
*> \verbatim
*> MVAL is INTEGER array, dimension (NM)
*> The values of the matrix row dimension M.
*> \endverbatim
*>
*> \param[in] NN
*> \verbatim
*> NN is INTEGER
*> The number of values of N contained in the vector NVAL.
*> \endverbatim
*>
*> \param[in] NVAL
*> \verbatim
*> NVAL is INTEGER array, dimension (NN)
*> The values of the matrix column dimension N.
*> \endverbatim
*>
*> \param[in] NNS
*> \verbatim
*> NNS is INTEGER
*> The number of values of NRHS contained in the vector NSVAL.
*> \endverbatim
*>
*> \param[in] NSVAL
*> \verbatim
*> NSVAL is INTEGER array, dimension (NNS)
*> The values of the number of right hand sides NRHS.
*> \endverbatim
*>
*> \param[in] NNB
*> \verbatim
*> NNB is INTEGER
*> The number of values of NB and NX contained in the
*> vectors NBVAL and NXVAL. The blocking parameters are used
*> in pairs (NB,NX).
*> \endverbatim
*>
*> \param[in] NBVAL
*> \verbatim
*> NBVAL is INTEGER array, dimension (NNB)
*> The values of the blocksize NB.
*> \endverbatim
*>
*> \param[in] NXVAL
*> \verbatim
*> NXVAL is INTEGER array, dimension (NNB)
*> The values of the crossover point NX.
*> \endverbatim
*>
*> \param[in] THRESH
*> \verbatim
*> THRESH is REAL
*> The threshold value for the test ratios. A result is
*> included in the output file if RESULT >= THRESH. To have
*> every test ratio printed, use THRESH = 0.
*> \endverbatim
*>
*> \param[out] A
*> \verbatim
*> A is REAL array, dimension (MMAX*NMAX)
*> where MMAX is the maximum value of M in MVAL and NMAX is the
*> maximum value of N in NVAL.
*> \endverbatim
*>
*> \param[out] COPYA
*> \verbatim
*> COPYA is REAL array, dimension (MMAX*NMAX)
*> \endverbatim
*>
*> \param[out] B
*> \verbatim
*> B is REAL array, dimension (MMAX*NSMAX)
*> where MMAX is the maximum value of M in MVAL and NSMAX is the
*> maximum value of NRHS in NSVAL.
*> \endverbatim
*>
*> \param[out] COPYB
*> \verbatim
*> COPYB is REAL array, dimension (MMAX*NSMAX)
*> \endverbatim
*>
*> \param[out] S
*> \verbatim
*> S is REAL array, dimension
*> (min(MMAX,NMAX))
*> \endverbatim
*>
*> \param[out] TAU
*> \verbatim
*> TAU is REAL array, dimension (MMAX)
*> \endverbatim
*>
*> \param[out] WORK
*> \verbatim
*> WORK is REAL array, dimension
*> (MMAX*NMAX + 4*NMAX + MMAX)
*> \endverbatim
*>
*> \param[out] IWORK
*> \verbatim
*> IWORK is INTEGER array, dimension (2*NMAX)
*> \endverbatim
*>
*> \param[in] NOUT
*> \verbatim
*> NOUT is INTEGER
*> The unit number for output.
*> \endverbatim
*
* Authors:
* ========
*
*> \author Univ. of Tennessee
*> \author Univ. of California Berkeley
*> \author Univ. of Colorado Denver
*> \author NAG Ltd.
*
*> \ingroup single_lin
*
* =====================================================================
SUBROUTINE SCHKQP3RK( DOTYPE, NM, MVAL, NN, NVAL, NNS, NSVAL,
$ NNB, NBVAL, NXVAL, THRESH, A, COPYA,
$ B, COPYB, S, TAU,
$ WORK, IWORK, NOUT )
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 NM, NN, NNB, NNS, NOUT
REAL THRESH
* ..
* .. Array Arguments ..
LOGICAL DOTYPE( * )
INTEGER IWORK( * ), NBVAL( * ), MVAL( * ), NVAL( * ),
$ NSVAL( * ), NXVAL( * )
REAL A( * ), COPYA( * ), B( * ), COPYB( * ),
$ S( * ), TAU( * ), WORK( * )
* ..
*
* =====================================================================
*
* .. Parameters ..
INTEGER NTYPES
PARAMETER ( NTYPES = 19 )
INTEGER NTESTS
PARAMETER ( NTESTS = 5 )
REAL ONE, ZERO, BIGNUM
PARAMETER ( ONE = 1.0E+0, ZERO = 0.0E+0,
$ BIGNUM = 1.0E+38 )
* ..
* .. Local Scalars ..
CHARACTER DIST, TYPE
CHARACTER*3 PATH
INTEGER I, IHIGH, ILOW, IM, IMAT, IN, INC_ZERO,
$ INB, IND_OFFSET_GEN,
$ IND_IN, IND_OUT, INS, INFO,
$ ISTEP, J, J_INC, J_FIRST_NZ, JB_ZERO,
$ KFACT, KL, KMAX, KU, LDA, LW, LWORK,
$ LWORK_MQR, M, MINMN, MINMNB_GEN, MODE, N,
$ NB, NB_ZERO, NERRS, NFAIL, NB_GEN, NRHS,
$ NRUN, NX, T
REAL ANORM, CNDNUM, EPS, ABSTOL, RELTOL,
$ DTEMP, MAXC2NRMK, RELMAXC2NRMK
* ..
* .. Local Arrays ..
INTEGER ISEED( 4 ), ISEEDY( 4 )
REAL RESULT( NTESTS ), RDUMMY( 1 )
* ..
* .. External Functions ..
REAL SLAMCH, SQPT01, SQRT11, SQRT12, SLANGE
EXTERNAL SLAMCH, SQPT01, SQRT11, SQRT12, SLANGE
* ..
* .. External Subroutines ..
EXTERNAL ALAERH, ALAHD, ALASUM, SAXPY, SGEQP3RK,
$ SLACPY, SLAORD, SLASET, SLATB4, SLATMS,
$ SORMQR, SSWAP, ICOPY, XLAENV
* ..
* .. Intrinsic Functions ..
INTRINSIC ABS, MAX, MIN, MOD, REAL
* ..
* .. Scalars in Common ..
LOGICAL LERR, OK
CHARACTER*32 SRNAMT
INTEGER INFOT, IOUNIT
* ..
* .. Common blocks ..
COMMON / INFOC / INFOT, IOUNIT, OK, LERR
COMMON / SRNAMC / SRNAMT
* ..
* .. Data statements ..
DATA ISEEDY / 1988, 1989, 1990, 1991 /
* ..
* .. Executable Statements ..
*
* Initialize constants and the random number seed.
*
PATH( 1: 1 ) = 'Single precision'
PATH( 2: 3 ) = 'QK'
NRUN = 0
NFAIL = 0
NERRS = 0
DO I = 1, 4
ISEED( I ) = ISEEDY( I )
END DO
EPS = SLAMCH( 'Epsilon' )
INFOT = 0
*
DO IM = 1, NM
*
* Do for each value of M in MVAL.
*
M = MVAL( IM )
LDA = MAX( 1, M )
*
DO IN = 1, NN
*
* Do for each value of N in NVAL.
*
N = NVAL( IN )
MINMN = MIN( M, N )
LWORK = MAX( 1, M*MAX( M, N )+4*MINMN+MAX( M, N ),
$ M*N + 2*MINMN + 4*N )
*
DO INS = 1, NNS
NRHS = NSVAL( INS )
*
* Set up parameters with SLATB4 and generate
* M-by-NRHS B matrix with SLATMS.
* IMAT = 14:
* Random matrix, CNDNUM = 2, NORM = ONE,
* MODE = 3 (geometric distribution of singular values).
*
CALL SLATB4( PATH, 14, M, NRHS, TYPE, KL, KU, ANORM,
$ MODE, CNDNUM, DIST )
*
SRNAMT = 'SLATMS'
CALL SLATMS( M, NRHS, DIST, ISEED, TYPE, S, MODE,
$ CNDNUM, ANORM, KL, KU, 'No packing',
$ COPYB, LDA, WORK, INFO )
*
* Check error code from SLATMS.
*
IF( INFO.NE.0 ) THEN
CALL ALAERH( PATH, 'SLATMS', INFO, 0, ' ', M,
$ NRHS, -1, -1, -1, 6, NFAIL, NERRS,
$ NOUT )
CYCLE
END IF
*
DO IMAT = 1, NTYPES
*
* Do the tests only if DOTYPE( IMAT ) is true.
*
IF( .NOT.DOTYPE( IMAT ) )
$ CYCLE
*
* The type of distribution used to generate the random
* eigen-/singular values:
* ( 'S' for symmetric distribution ) => UNIFORM( -1, 1 )
*
* Do for each type of NON-SYMMETRIC matrix: CNDNUM NORM MODE
* 1. Zero matrix
* 2. Random, Diagonal, CNDNUM = 2 CNDNUM = 2 ONE 3 ( geometric distribution of singular values )
* 3. Random, Upper triangular, CNDNUM = 2 CNDNUM = 2 ONE 3 ( geometric distribution of singular values )
* 4. Random, Lower triangular, CNDNUM = 2 CNDNUM = 2 ONE 3 ( geometric distribution of singular values )
* 5. Random, First column is zero, CNDNUM = 2 CNDNUM = 2 ONE 3 ( geometric distribution of singular values )
* 6. Random, Last MINMN column is zero, CNDNUM = 2 CNDNUM = 2 ONE 3 ( geometric distribution of singular values )
* 7. Random, Last N column is zero, CNDNUM = 2 CNDNUM = 2 ONE 3 ( geometric distribution of singular values )
* 8. Random, Middle column in MINMN is zero, CNDNUM = 2 CNDNUM = 2 ONE 3 ( geometric distribution of singular values )
* 9. Random, First half of MINMN columns are zero, CNDNUM = 2 CNDNUM = 2 ONE 3 ( geometric distribution of singular values )
* 10. Random, Last columns are zero starting from MINMN/2+1, CNDNUM = 2 CNDNUM = 2 ONE 3 ( geometric distribution of singular values )
* 11. Random, Half MINMN columns in the middle are zero starting
* from MINMN/2-(MINMN/2)/2+1, CNDNUM = 2 CNDNUM = 2 ONE 3 ( geometric distribution of singular values )
* 12. Random, Odd columns are ZERO, CNDNUM = 2 CNDNUM = 2 ONE 3 ( geometric distribution of singular values )
* 13. Random, Even columns are ZERO, CNDNUM = 2 CNDNUM = 2 ONE 3 ( geometric distribution of singular values )
* 14. Random, CNDNUM = 2 CNDNUM = 2 ONE 3 ( geometric distribution of singular values )
* 15. Random, CNDNUM = sqrt(0.1/EPS) CNDNUM = BADC1 = sqrt(0.1/EPS) ONE 3 ( geometric distribution of singular values )
* 16. Random, CNDNUM = 0.1/EPS CNDNUM = BADC2 = 0.1/EPS ONE 3 ( geometric distribution of singular values )
* 17. Random, CNDNUM = 0.1/EPS, CNDNUM = BADC2 = 0.1/EPS ONE 2 ( one small singular value, S(N)=1/CNDNUM )
* one small singular value S(N)=1/CNDNUM
* 18. Random, CNDNUM = 2, scaled near underflow CNDNUM = 2 SMALL = SAFMIN
* 19. Random, CNDNUM = 2, scaled near overflow CNDNUM = 2 LARGE = 1.0/( 0.25 * ( SAFMIN / EPS ) ) 3 ( geometric distribution of singular values )
*
IF( IMAT.EQ.1 ) THEN
*
* Matrix 1: Zero matrix
*
CALL SLASET( 'Full', M, N, ZERO, ZERO, COPYA, LDA )
DO I = 1, MINMN
S( I ) = ZERO
END DO
*
ELSE IF( (IMAT.GE.2 .AND. IMAT.LE.4 )
$ .OR. (IMAT.GE.14 .AND. IMAT.LE.19 ) ) THEN
*
* Matrices 2-5.
*
* Set up parameters with SLATB4 and generate a test
* matrix with SLATMS.
*
CALL SLATB4( PATH, IMAT, M, N, TYPE, KL, KU, ANORM,
$ MODE, CNDNUM, DIST )
*
SRNAMT = 'SLATMS'
CALL SLATMS( M, N, DIST, ISEED, TYPE, S, MODE,
$ CNDNUM, ANORM, KL, KU, 'No packing',
$ COPYA, LDA, WORK, INFO )
*
* Check error code from SLATMS.
*
IF( INFO.NE.0 ) THEN
CALL ALAERH( PATH, 'SLATMS', INFO, 0, ' ', M, N,
$ -1, -1, -1, IMAT, NFAIL, NERRS,
$ NOUT )
CYCLE
END IF
*
CALL SLAORD( 'Decreasing', MINMN, S, 1 )
*
ELSE IF( MINMN.GE.2
$ .AND. IMAT.GE.5 .AND. IMAT.LE.13 ) THEN
*
* Rectangular matrices 5-13 that contain zero columns,
* only for matrices MINMN >=2.
*
* JB_ZERO is the column index of ZERO block.
* NB_ZERO is the column block size of ZERO block.
* NB_GEN is the column blcok size of the
* generated block.
* J_INC in the non_zero column index increment
* for matrix 12 and 13.
* J_FIRS_NZ is the index of the first non-zero
* column.
*
IF( IMAT.EQ.5 ) THEN
*
* First column is zero.
*
JB_ZERO = 1
NB_ZERO = 1
NB_GEN = N - NB_ZERO
*
ELSE IF( IMAT.EQ.6 ) THEN
*
* Last column MINMN is zero.
*
JB_ZERO = MINMN
NB_ZERO = 1
NB_GEN = N - NB_ZERO
*
ELSE IF( IMAT.EQ.7 ) THEN
*
* Last column N is zero.
*
JB_ZERO = N
NB_ZERO = 1
NB_GEN = N - NB_ZERO
*
ELSE IF( IMAT.EQ.8 ) THEN
*
* Middle column in MINMN is zero.
*
JB_ZERO = MINMN / 2 + 1
NB_ZERO = 1
NB_GEN = N - NB_ZERO
*
ELSE IF( IMAT.EQ.9 ) THEN
*
* First half of MINMN columns is zero.
*
JB_ZERO = 1
NB_ZERO = MINMN / 2
NB_GEN = N - NB_ZERO
*
ELSE IF( IMAT.EQ.10 ) THEN
*
* Last columns are zero columns,
* starting from (MINMN / 2 + 1) column.
*
JB_ZERO = MINMN / 2 + 1
NB_ZERO = N - JB_ZERO + 1
NB_GEN = N - NB_ZERO
*
ELSE IF( IMAT.EQ.11 ) THEN
*
* Half of the columns in the middle of MINMN
* columns is zero, starting from
* MINMN/2 - (MINMN/2)/2 + 1 column.
*
JB_ZERO = MINMN / 2 - (MINMN / 2) / 2 + 1
NB_ZERO = MINMN / 2
NB_GEN = N - NB_ZERO
*
ELSE IF( IMAT.EQ.12 ) THEN
*
* Odd-numbered columns are zero,
*
NB_GEN = N / 2
NB_ZERO = N - NB_GEN
J_INC = 2
J_FIRST_NZ = 2
*
ELSE IF( IMAT.EQ.13 ) THEN
*
* Even-numbered columns are zero.
*
NB_ZERO = N / 2
NB_GEN = N - NB_ZERO
J_INC = 2
J_FIRST_NZ = 1
*
END IF
*
*
* 1) Set the first NB_ZERO columns in COPYA(1:M,1:N)
* to zero.
*
CALL SLASET( 'Full', M, NB_ZERO, ZERO, ZERO,
$ COPYA, LDA )
*
* 2) Generate an M-by-(N-NB_ZERO) matrix with the
* chosen singular value distribution
* in COPYA(1:M,NB_ZERO+1:N).
*
CALL SLATB4( PATH, IMAT, M, NB_GEN, TYPE, KL, KU,
$ ANORM, MODE, CNDNUM, DIST )
*
SRNAMT = 'SLATMS'
*
IND_OFFSET_GEN = NB_ZERO * LDA
*
CALL SLATMS( M, NB_GEN, DIST, ISEED, TYPE, S, MODE,
$ CNDNUM, ANORM, KL, KU, 'No packing',
$ COPYA( IND_OFFSET_GEN + 1 ), LDA,
$ WORK, INFO )
*
* Check error code from SLATMS.
*
IF( INFO.NE.0 ) THEN
CALL ALAERH( PATH, 'SLATMS', INFO, 0, ' ', M,
$ NB_GEN, -1, -1, -1, IMAT, NFAIL,
$ NERRS, NOUT )
CYCLE
END IF
*
* 3) Swap the gererated colums from the right side
* NB_GEN-size block in COPYA into correct column
* positions.
*
IF( IMAT.EQ.6
$ .OR. IMAT.EQ.7
$ .OR. IMAT.EQ.8
$ .OR. IMAT.EQ.10
$ .OR. IMAT.EQ.11 ) THEN
*
* Move by swapping the generated columns
* from the right NB_GEN-size block from
* (NB_ZERO+1:NB_ZERO+JB_ZERO)
* into columns (1:JB_ZERO-1).
*
DO J = 1, JB_ZERO-1, 1
CALL SSWAP( M,
$ COPYA( ( NB_ZERO+J-1)*LDA+1), 1,
$ COPYA( (J-1)*LDA + 1 ), 1 )
END DO
*
ELSE IF( IMAT.EQ.12 .OR. IMAT.EQ.13 ) THEN
*
* ( IMAT = 12, Odd-numbered ZERO columns. )
* Swap the generated columns from the right
* NB_GEN-size block into the even zero colums in the
* left NB_ZERO-size block.
*
* ( IMAT = 13, Even-numbered ZERO columns. )
* Swap the generated columns from the right
* NB_GEN-size block into the odd zero colums in the
* left NB_ZERO-size block.
*
DO J = 1, NB_GEN, 1
IND_OUT = ( NB_ZERO+J-1 )*LDA + 1
IND_IN = ( J_INC*(J-1)+(J_FIRST_NZ-1) )*LDA
$ + 1
CALL SSWAP( M,
$ COPYA( IND_OUT ), 1,
$ COPYA( IND_IN), 1 )
END DO
*
END IF
*
* 5) Order the singular values generated by
* DLAMTS in decreasing order and add trailing zeros
* that correspond to zero columns.
* The total number of singular values is MINMN.
*
MINMNB_GEN = MIN( M, NB_GEN )
*
DO I = MINMNB_GEN+1, MINMN
S( I ) = ZERO
END DO
*
ELSE
*
* IF(MINMN.LT.2) skip this size for this matrix type.
*
CYCLE
END IF
*
* Initialize a copy array for a pivot array for SGEQP3RK.
*
DO I = 1, N
IWORK( I ) = 0
END DO
*
DO INB = 1, NNB
*
* Do for each pair of values (NB,NX) in NBVAL and NXVAL.
*
NB = NBVAL( INB )
CALL XLAENV( 1, NB )
NX = NXVAL( INB )
CALL XLAENV( 3, NX )
*
* We do MIN(M,N)+1 because we need a test for KMAX > N,
* when KMAX is larger than MIN(M,N), KMAX should be
* KMAX = MIN(M,N)
*
DO KMAX = 0, MIN(M,N)+1
*
* Get a working copy of COPYA into A( 1:M,1:N ).
* Get a working copy of COPYB into A( 1:M, (N+1):NRHS ).
* Get a working copy of COPYB into into B( 1:M, 1:NRHS ).
* Get a working copy of IWORK(1:N) awith zeroes into
* which is going to be used as pivot array IWORK( N+1:2N ).
* NOTE: IWORK(2N+1:3N) is going to be used as a WORK array
* for the routine.
*
CALL SLACPY( 'All', M, N, COPYA, LDA, A, LDA )
CALL SLACPY( 'All', M, NRHS, COPYB, LDA,
$ A( LDA*N + 1 ), LDA )
CALL SLACPY( 'All', M, NRHS, COPYB, LDA,
$ B, LDA )
CALL ICOPY( N, IWORK( 1 ), 1, IWORK( N+1 ), 1 )
*
ABSTOL = -1.0
RELTOL = -1.0
*
* Compute the QR factorization with pivoting of A
*
LW = MAX( 1, MAX( 2*N + NB*( N+NRHS+1 ),
$ 3*N + NRHS - 1 ) )
*
* Compute SGEQP3RK factorization of A.
*
SRNAMT = 'SGEQP3RK'
CALL SGEQP3RK( M, N, NRHS, KMAX, ABSTOL, RELTOL,
$ A, LDA, KFACT, MAXC2NRMK,
$ RELMAXC2NRMK, IWORK( N+1 ), TAU,
$ WORK, LW, IWORK( 2*N+1 ), INFO )
*
* Check error code from SGEQP3RK.
*
IF( INFO.LT.0 )
$ CALL ALAERH( PATH, 'SGEQP3RK', INFO, 0, ' ',
$ M, N, NX, -1, NB, IMAT,
$ NFAIL, NERRS, NOUT )
*
* Compute test 1:
*
* This test in only for the full rank factorization of
* the matrix A.
*
* Array S(1:min(M,N)) contains svd(A) the sigular values
* of the original matrix A in decreasing absolute value
* order. The test computes svd(R), the vector sigular
* values of the upper trapezoid of A(1:M,1:N) that
* contains the factor R, in decreasing order. The test
* returns the ratio:
*
* 2-norm(svd(R) - svd(A)) / ( max(M,N) * 2-norm(svd(A)) * EPS )
*
IF( KFACT.EQ.MINMN ) THEN
*
RESULT( 1 ) = SQRT12( M, N, A, LDA, S, WORK,
$ LWORK )
*
DO T = 1, 1
IF( RESULT( T ).GE.THRESH ) THEN
IF( NFAIL.EQ.0 .AND. NERRS.EQ.0 )
$ CALL ALAHD( NOUT, PATH )
WRITE( NOUT, FMT = 9999 ) 'SGEQP3RK', M, N,
$ NRHS, KMAX, ABSTOL, RELTOL, NB, NX,
$ IMAT, T, RESULT( T )
NFAIL = NFAIL + 1
END IF
END DO
NRUN = NRUN + 1
*
* End test 1
*
END IF
*
* Compute test 2:
*
* The test returns the ratio:
*
* 1-norm( A*P - Q*R ) / ( max(M,N) * 1-norm(A) * EPS )
*
RESULT( 2 ) = SQPT01( M, N, KFACT, COPYA, A, LDA, TAU,
$ IWORK( N+1 ), WORK, LWORK )
*
* Compute test 3:
*
* The test returns the ratio:
*
* 1-norm( Q**T * Q - I ) / ( M * EPS )
*
RESULT( 3 ) = SQRT11( M, KFACT, A, LDA, TAU, WORK,
$ LWORK )
*
* Print information about the tests that did not pass
* the threshold.
*
DO T = 2, 3
IF( RESULT( T ).GE.THRESH ) THEN
IF( NFAIL.EQ.0 .AND. NERRS.EQ.0 )
$ CALL ALAHD( NOUT, PATH )
WRITE( NOUT, FMT = 9999 ) 'SGEQP3RK', M, N,
$ NRHS, KMAX, ABSTOL, RELTOL,
$ NB, NX, IMAT, T, RESULT( T )
NFAIL = NFAIL + 1
END IF
END DO
NRUN = NRUN + 2
*
* Compute test 4:
*
* This test is only for the factorizations with the
* rank greater than 2.
* The elements on the diagonal of R should be non-
* increasing.
*
* The test returns the ratio:
*
* Returns 1.0D+100 if abs(R(K+1,K+1)) > abs(R(K,K)),
* K=1:KFACT-1
*
IF( MIN(KFACT, MINMN).GE.2 ) THEN
*
DO J = 1, KFACT-1, 1
DTEMP = (( ABS( A( (J-1)*M+J ) ) -
$ ABS( A( (J)*M+J+1 ) ) ) /
$ ABS( A(1) ) )
*
IF( DTEMP.LT.ZERO ) THEN
RESULT( 4 ) = BIGNUM
END IF
*
END DO
*
* Print information about the tests that did not
* pass the threshold.
*
DO T = 4, 4
IF( RESULT( T ).GE.THRESH ) THEN
IF( NFAIL.EQ.0 .AND. NERRS.EQ.0 )
$ CALL ALAHD( NOUT, PATH )
WRITE( NOUT, FMT = 9999 ) 'SGEQP3RK',
$ M, N, NRHS, KMAX, ABSTOL, RELTOL,
$ NB, NX, IMAT, T,
$ RESULT( T )
NFAIL = NFAIL + 1
END IF
END DO
NRUN = NRUN + 1
*
* End test 4.
*
END IF
*
* Compute test 5:
*
* This test in only for matrix A with min(M,N) > 0.
*
* The test returns the ratio:
*
* 1-norm(Q**T * B - Q**T * B ) /
* ( M * EPS )
*
* (1) Compute B:=Q**T * B in the matrix B.
*
IF( MINMN.GT.0 ) THEN
*
LWORK_MQR = MAX(1, NRHS)
CALL SORMQR( 'Left', 'Transpose',
$ M, NRHS, KFACT, A, LDA, TAU, B, LDA,
$ WORK, LWORK_MQR, INFO )
*
DO I = 1, NRHS
*
* Compare N+J-th column of A and J-column of B.
*
CALL SAXPY( M, -ONE, A( ( N+I-1 )*LDA+1 ), 1,
$ B( ( I-1 )*LDA+1 ), 1 )
END DO
*
RESULT( 5 ) =
$ ABS(
$ SLANGE( 'One-norm', M, NRHS, B, LDA, RDUMMY ) /
$ ( REAL( M )*SLAMCH( 'Epsilon' ) )
$ )
*
* Print information about the tests that did not pass
* the threshold.
*
DO T = 5, 5
IF( RESULT( T ).GE.THRESH ) THEN
IF( NFAIL.EQ.0 .AND. NERRS.EQ.0 )
$ CALL ALAHD( NOUT, PATH )
WRITE( NOUT, FMT = 9999 ) 'SGEQP3RK', M, N,
$ NRHS, KMAX, ABSTOL, RELTOL,
$ NB, NX, IMAT, T, RESULT( T )
NFAIL = NFAIL + 1
END IF
END DO
NRUN = NRUN + 1
*
* End compute test 5.
*
END IF
*
* END DO KMAX = 1, MIN(M,N)+1
*
END DO
*
* END DO for INB = 1, NNB
*
END DO
*
* END DO for IMAT = 1, NTYPES
*
END DO
*
* END DO for INS = 1, NNS
*
END DO
*
* END DO for IN = 1, NN
*
END DO
*
* END DO for IM = 1, NM
*
END DO
*
* Print a summary of the results.
*
CALL ALASUM( PATH, NOUT, NFAIL, NRUN, NERRS )
*
9999 FORMAT( 1X, A, ' M =', I5, ', N =', I5, ', NRHS =', I5,
$ ', KMAX =', I5, ', ABSTOL =', G12.5,
$ ', RELTOL =', G12.5, ', NB =', I4, ', NX =', I4,
$ ', type ', I2, ', test ', I2, ', ratio =', G12.5 )
*
* End of SCHKQP3RK
*
END

View File

@ -153,9 +153,6 @@
* .. Intrinsic Functions .. * .. Intrinsic Functions ..
INTRINSIC ABS, MAX, SQRT INTRINSIC ABS, MAX, SQRT
* .. * ..
* .. External Subroutines ..
EXTERNAL SLABAD
* ..
* .. Save statement .. * .. Save statement ..
SAVE EPS, SMALL, LARGE, BADC1, BADC2, FIRST SAVE EPS, SMALL, LARGE, BADC1, BADC2, FIRST
* .. * ..
@ -173,11 +170,6 @@
BADC1 = SQRT( BADC2 ) BADC1 = SQRT( BADC2 )
SMALL = SLAMCH( 'Safe minimum' ) SMALL = SLAMCH( 'Safe minimum' )
LARGE = ONE / SMALL LARGE = ONE / SMALL
*
* If it looks like we're on a Cray, take the square root of
* SMALL and LARGE to avoid overflow and underflow problems.
*
CALL SLABAD( SMALL, LARGE )
SMALL = SHRINK*( SMALL / EPS ) SMALL = SHRINK*( SMALL / EPS )
LARGE = ONE / SMALL LARGE = ONE / SMALL
END IF END IF
@ -232,6 +224,110 @@
ELSE ELSE
ANORM = ONE ANORM = ONE
END IF END IF
*
ELSE IF( LSAMEN( 2, C2, 'QK' ) ) THEN
*
* xQK: truncated QR with pivoting.
* Set parameters to generate a general
* M x N matrix.
*
* Set TYPE, the type of matrix to be generated. 'N' is nonsymmetric.
*
TYPE = 'N'
*
* Set DIST, the type of distribution for the random
* number generator. 'S' is
*
DIST = 'S'
*
* Set the lower and upper bandwidths.
*
IF( IMAT.EQ.2 ) THEN
*
* 2. Random, Diagonal, CNDNUM = 2
*
KL = 0
KU = 0
CNDNUM = TWO
ANORM = ONE
MODE = 3
ELSE IF( IMAT.EQ.3 ) THEN
*
* 3. Random, Upper triangular, CNDNUM = 2
*
KL = 0
KU = MAX( N-1, 0 )
CNDNUM = TWO
ANORM = ONE
MODE = 3
ELSE IF( IMAT.EQ.4 ) THEN
*
* 4. Random, Lower triangular, CNDNUM = 2
*
KL = MAX( M-1, 0 )
KU = 0
CNDNUM = TWO
ANORM = ONE
MODE = 3
ELSE
*
* 5.-19. Rectangular matrix
*
KL = MAX( M-1, 0 )
KU = MAX( N-1, 0 )
*
IF( IMAT.GE.5 .AND. IMAT.LE.14 ) THEN
*
* 5.-14. Random, CNDNUM = 2.
*
CNDNUM = TWO
ANORM = ONE
MODE = 3
*
ELSE IF( IMAT.EQ.15 ) THEN
*
* 15. Random, CNDNUM = sqrt(0.1/EPS)
*
CNDNUM = BADC1
ANORM = ONE
MODE = 3
*
ELSE IF( IMAT.EQ.16 ) THEN
*
* 16. Random, CNDNUM = 0.1/EPS
*
CNDNUM = BADC2
ANORM = ONE
MODE = 3
*
ELSE IF( IMAT.EQ.17 ) THEN
*
* 17. Random, CNDNUM = 0.1/EPS,
* one small singular value S(N)=1/CNDNUM
*
CNDNUM = BADC2
ANORM = ONE
MODE = 2
*
ELSE IF( IMAT.EQ.18 ) THEN
*
* 18. Random, scaled near underflow
*
CNDNUM = TWO
ANORM = SMALL
MODE = 3
*
ELSE IF( IMAT.EQ.19 ) THEN
*
* 19. Random, scaled near overflow
*
CNDNUM = TWO
ANORM = LARGE
MODE = 3
*
END IF
*
END IF
* *
ELSE IF( LSAMEN( 2, C2, 'GE' ) ) THEN ELSE IF( LSAMEN( 2, C2, 'GE' ) ) THEN
* *
@ -518,17 +614,18 @@
* *
* Set the norm and condition number. * Set the norm and condition number.
* *
IF( IMAT.EQ.2 .OR. IMAT.EQ.8 ) THEN MAT = ABS( IMAT )
IF( MAT.EQ.2 .OR. MAT.EQ.8 ) THEN
CNDNUM = BADC1 CNDNUM = BADC1
ELSE IF( IMAT.EQ.3 .OR. IMAT.EQ.9 ) THEN ELSE IF( MAT.EQ.3 .OR. MAT.EQ.9 ) THEN
CNDNUM = BADC2 CNDNUM = BADC2
ELSE ELSE
CNDNUM = TWO CNDNUM = TWO
END IF END IF
* *
IF( IMAT.EQ.4 ) THEN IF( MAT.EQ.4 ) THEN
ANORM = SMALL ANORM = SMALL
ELSE IF( IMAT.EQ.5 ) THEN ELSE IF( MAT.EQ.5 ) THEN
ANORM = LARGE ANORM = LARGE
ELSE ELSE
ANORM = ONE ANORM = ONE

View File

@ -33,7 +33,8 @@
*> Householder vectors, and the rest of AF contains a partially updated *> Householder vectors, and the rest of AF contains a partially updated
*> matrix. *> matrix.
*> *>
*> This function returns ||A*P - Q*R||/(||norm(A)||*eps*M) *> This function returns ||A*P - Q*R|| / ( ||norm(A)||*eps*max(M,N) )
*> where || . || is matrix one norm.
*> \endverbatim *> \endverbatim
* *
* Arguments: * Arguments:
@ -172,28 +173,28 @@
* *
NORMA = SLANGE( 'One-norm', M, N, A, LDA, RWORK ) NORMA = SLANGE( 'One-norm', M, N, A, LDA, RWORK )
* *
DO 30 J = 1, K DO J = 1, K
DO 10 I = 1, MIN( J, M ) DO I = 1, MIN( J, M )
WORK( ( J-1 )*M+I ) = AF( I, J ) WORK( ( J-1 )*M+I ) = AF( I, J )
10 CONTINUE END DO
DO 20 I = J + 1, M DO I = J + 1, M
WORK( ( J-1 )*M+I ) = ZERO WORK( ( J-1 )*M+I ) = ZERO
20 CONTINUE END DO
30 CONTINUE END DO
DO 40 J = K + 1, N DO J = K + 1, N
CALL SCOPY( M, AF( 1, J ), 1, WORK( ( J-1 )*M+1 ), 1 ) CALL SCOPY( M, AF( 1, J ), 1, WORK( ( J-1 )*M+1 ), 1 )
40 CONTINUE END DO
* *
CALL SORMQR( 'Left', 'No transpose', M, N, K, AF, LDA, TAU, WORK, CALL SORMQR( 'Left', 'No transpose', M, N, K, AF, LDA, TAU, WORK,
$ M, WORK( M*N+1 ), LWORK-M*N, INFO ) $ M, WORK( M*N+1 ), LWORK-M*N, INFO )
* *
DO 50 J = 1, N DO J = 1, N
* *
* Compare i-th column of QR and jpvt(i)-th column of A * Compare i-th column of QR and jpvt(i)-th column of A
* *
CALL SAXPY( M, -ONE, A( 1, JPVT( J ) ), 1, WORK( ( J-1 )*M+1 ), CALL SAXPY( M, -ONE, A( 1, JPVT( J ) ), 1, WORK( ( J-1 )*M+1 ),
$ 1 ) $ 1 )
50 CONTINUE END DO
* *
SQPT01 = SLANGE( 'One-norm', M, N, WORK, M, RWORK ) / SQPT01 = SLANGE( 'One-norm', M, N, WORK, M, RWORK ) /
$ ( REAL( MAX( M, N ) )*SLAMCH( 'Epsilon' ) ) $ ( REAL( MAX( M, N ) )*SLAMCH( 'Epsilon' ) )

View File

@ -157,9 +157,9 @@
CALL SORM2R( 'Left', 'Transpose', M, M, K, A, LDA, TAU, WORK, M, CALL SORM2R( 'Left', 'Transpose', M, M, K, A, LDA, TAU, WORK, M,
$ WORK( M*M+1 ), INFO ) $ WORK( M*M+1 ), INFO )
* *
DO 10 J = 1, M DO J = 1, M
WORK( ( J-1 )*M+J ) = WORK( ( J-1 )*M+J ) - ONE WORK( ( J-1 )*M+J ) = WORK( ( J-1 )*M+J ) - ONE
10 CONTINUE END DO
* *
SQRT11 = SLANGE( 'One-norm', M, M, WORK, M, RDUMMY ) / SQRT11 = SLANGE( 'One-norm', M, M, WORK, M, RDUMMY ) /
$ ( REAL( M )*SLAMCH( 'Epsilon' ) ) $ ( REAL( M )*SLAMCH( 'Epsilon' ) )

View File

@ -26,7 +26,7 @@
*> SQRT12 computes the singular values `svlues' of the upper trapezoid *> SQRT12 computes the singular values `svlues' of the upper trapezoid
*> of A(1:M,1:N) and returns the ratio *> of A(1:M,1:N) and returns the ratio
*> *>
*> || s - svlues||/(||svlues||*eps*max(M,N)) *> || svlues - s ||/(||s||*eps*max(M,N))
*> \endverbatim *> \endverbatim
* *
* Arguments: * Arguments:
@ -113,8 +113,7 @@
EXTERNAL SASUM, SLAMCH, SLANGE, SNRM2 EXTERNAL SASUM, SLAMCH, SLANGE, SNRM2
* .. * ..
* .. External Subroutines .. * .. External Subroutines ..
EXTERNAL SAXPY, SBDSQR, SGEBD2, SLABAD, SLASCL, SLASET, EXTERNAL SAXPY, SBDSQR, SGEBD2, SLASCL, SLASET, XERBLA
$ XERBLA
* .. * ..
* .. Intrinsic Functions .. * .. Intrinsic Functions ..
INTRINSIC MAX, MIN, REAL INTRINSIC MAX, MIN, REAL
@ -145,17 +144,16 @@
* Copy upper triangle of A into work * Copy upper triangle of A into work
* *
CALL SLASET( 'Full', M, N, ZERO, ZERO, WORK, M ) CALL SLASET( 'Full', M, N, ZERO, ZERO, WORK, M )
DO 20 J = 1, N DO J = 1, N
DO 10 I = 1, MIN( J, M ) DO I = 1, MIN( J, M )
WORK( ( J-1 )*M+I ) = A( I, J ) WORK( ( J-1 )*M+I ) = A( I, J )
10 CONTINUE END DO
20 CONTINUE END DO
* *
* Get machine parameters * Get machine parameters
* *
SMLNUM = SLAMCH( 'S' ) / SLAMCH( 'P' ) SMLNUM = SLAMCH( 'S' ) / SLAMCH( 'P' )
BIGNUM = ONE / SMLNUM BIGNUM = ONE / SMLNUM
CALL SLABAD( SMLNUM, BIGNUM )
* *
* Scale work if max entry outside range [SMLNUM,BIGNUM] * Scale work if max entry outside range [SMLNUM,BIGNUM]
* *
@ -199,9 +197,9 @@
* *
ELSE ELSE
* *
DO 30 I = 1, MN DO I = 1, MN
WORK( M*N+I ) = ZERO WORK( M*N+I ) = ZERO
30 CONTINUE END DO
END IF END IF
* *
* Compare s and singular values of work * Compare s and singular values of work

View File

@ -69,6 +69,7 @@
*> ZLQ 8 List types on next line if 0 < NTYPES < 8 *> ZLQ 8 List types on next line if 0 < NTYPES < 8
*> ZQL 8 List types on next line if 0 < NTYPES < 8 *> ZQL 8 List types on next line if 0 < NTYPES < 8
*> ZQP 6 List types on next line if 0 < NTYPES < 6 *> ZQP 6 List types on next line if 0 < NTYPES < 6
*> ZQK 19 List types on next line if 0 < NTYPES < 19
*> ZTZ 3 List types on next line if 0 < NTYPES < 3 *> ZTZ 3 List types on next line if 0 < NTYPES < 3
*> ZLS 6 List types on next line if 0 < NTYPES < 6 *> ZLS 6 List types on next line if 0 < NTYPES < 6
*> ZEQ *> ZEQ
@ -153,12 +154,11 @@
$ NBVAL( MAXIN ), NBVAL2( MAXIN ), $ NBVAL( MAXIN ), NBVAL2( MAXIN ),
$ NSVAL( MAXIN ), NVAL( MAXIN ), NXVAL( MAXIN ), $ NSVAL( MAXIN ), NVAL( MAXIN ), NXVAL( MAXIN ),
$ RANKVAL( MAXIN ), PIV( NMAX ) $ RANKVAL( MAXIN ), PIV( NMAX )
DOUBLE PRECISION S( 2*NMAX ) * ..
COMPLEX*16 E( NMAX )
*
* .. Allocatable Arrays .. * .. Allocatable Arrays ..
INTEGER AllocateStatus INTEGER AllocateStatus
DOUBLE PRECISION, DIMENSION(:), ALLOCATABLE:: RWORK DOUBLE PRECISION, DIMENSION(:), ALLOCATABLE:: RWORK, S
COMPLEX*16, DIMENSION(:), ALLOCATABLE :: E
COMPLEX*16, DIMENSION(:,:), ALLOCATABLE:: A, B, WORK COMPLEX*16, DIMENSION(:,:), ALLOCATABLE:: A, B, WORK
* .. * ..
* .. External Functions .. * .. External Functions ..
@ -170,15 +170,16 @@
EXTERNAL ALAREQ, ZCHKEQ, ZCHKGB, ZCHKGE, ZCHKGT, ZCHKHE, EXTERNAL ALAREQ, ZCHKEQ, ZCHKGB, ZCHKGE, ZCHKGT, ZCHKHE,
$ ZCHKHE_ROOK, ZCHKHE_RK, ZCHKHE_AA, ZCHKHP, $ ZCHKHE_ROOK, ZCHKHE_RK, ZCHKHE_AA, ZCHKHP,
$ ZCHKLQ, ZCHKUNHR_COL, ZCHKPB, ZCHKPO, ZCHKPS, $ ZCHKLQ, ZCHKUNHR_COL, ZCHKPB, ZCHKPO, ZCHKPS,
$ ZCHKPP, ZCHKPT, ZCHKQ3, ZCHKQL, ZCHKQR, ZCHKRQ, $ ZCHKPP, ZCHKPT, ZCHKQ3, ZCHKQP3RK, ZCHKQL,
$ ZCHKSP, ZCHKSY, ZCHKSY_ROOK, ZCHKSY_RK, $ ZCHKQR, ZCHKRQ, ZCHKSP, ZCHKSY, ZCHKSY_ROOK,
$ ZCHKSY_AA, ZCHKTB, ZCHKTP, ZCHKTR, ZCHKTZ, $ ZCHKSY_RK, ZCHKSY_AA, ZCHKTB, ZCHKTP, ZCHKTR,
$ ZDRVGB, ZDRVGE, ZDRVGT, ZDRVHE, ZDRVHE_ROOK, $ ZCHKTZ, ZDRVGB, ZDRVGE, ZDRVGT, ZDRVHE,
$ ZDRVHE_RK, ZDRVHE_AA, ZDRVHE_AA_2STAGE, ZDRVHP, $ ZDRVHE_ROOK, ZDRVHE_RK, ZDRVHE_AA,
$ ZDRVLS, ZDRVPB, ZDRVPO, ZDRVPP, ZDRVPT, $ ZDRVHE_AA_2STAGE, ZDRVHP, ZDRVLS, ZDRVPB,
$ ZDRVSP, ZDRVSY, ZDRVSY_ROOK, ZDRVSY_RK, $ ZDRVPO, ZDRVPP, ZDRVPT, ZDRVSP, ZDRVSY,
$ ZDRVSY_AA, ZDRVSY_AA_2STAGE, ILAVER, ZCHKQRT, $ ZDRVSY_ROOK, ZDRVSY_RK, ZDRVSY_AA,
$ ZCHKQRTP, ZCHKLQT, ZCHKLQTP, ZCHKTSQR $ ZDRVSY_AA_2STAGE, ILAVER, ZCHKQRT, ZCHKQRTP,
$ ZCHKLQT, ZCHKLQTP, ZCHKTSQR
* .. * ..
* .. Scalars in Common .. * .. Scalars in Common ..
LOGICAL LERR, OK LOGICAL LERR, OK
@ -197,13 +198,18 @@
DATA THREQ / 2.0D0 / , INTSTR / '0123456789' / DATA THREQ / 2.0D0 / , INTSTR / '0123456789' /
* *
* .. Allocate memory dynamically .. * .. Allocate memory dynamically ..
ALLOCATE (RWORK( 150*NMAX+2*MAXRHS ), STAT = AllocateStatus) *
ALLOCATE ( A ( (KDMAX+1) * NMAX, 7 ), STAT = AllocateStatus)
IF (AllocateStatus /= 0) STOP "*** Not enough memory ***" IF (AllocateStatus /= 0) STOP "*** Not enough memory ***"
ALLOCATE (A ((KDMAX+1) * NMAX, 7), STAT = AllocateStatus) ALLOCATE ( B ( NMAX * MAXRHS, 4 ), STAT = AllocateStatus)
IF (AllocateStatus /= 0 ) STOP "*** Not enough memory ***"
ALLOCATE ( WORK ( NMAX, NMAX+MAXRHS+10 ), STAT = AllocateStatus)
IF (AllocateStatus /= 0) STOP "*** Not enough memory ***" IF (AllocateStatus /= 0) STOP "*** Not enough memory ***"
ALLOCATE (B (NMAX * MAXRHS, 4), STAT = AllocateStatus) ALLOCATE ( E( NMAX ), STAT = AllocateStatus )
IF (AllocateStatus /= 0) STOP "*** Not enough memory ***" IF (AllocateStatus /= 0) STOP "*** Not enough memory ***"
ALLOCATE (WORK (NMAX, NMAX+MAXRHS+10), STAT = AllocateStatus) ALLOCATE ( S( 2*NMAX ), STAT = AllocateStatus)
IF (AllocateStatus /= 0) STOP "*** Not enough memory ***"
ALLOCATE ( RWORK( 150*NMAX+2*MAXRHS ), STAT = AllocateStatus)
IF (AllocateStatus /= 0) STOP "*** Not enough memory ***" IF (AllocateStatus /= 0) STOP "*** Not enough memory ***"
* .. * ..
* .. Executable Statements .. * .. Executable Statements ..
@ -1109,6 +1115,23 @@
ELSE ELSE
WRITE( NOUT, FMT = 9989 )PATH WRITE( NOUT, FMT = 9989 )PATH
END IF END IF
*
ELSE IF( LSAMEN( 2, C2, 'QK' ) ) THEN
*
* QK: truncated QR factorization with pivoting
*
NTYPES = 19
CALL ALAREQ( PATH, NMATS, DOTYPE, NTYPES, NIN, NOUT )
*
IF( TSTCHK ) THEN
CALL ZCHKQP3RK( DOTYPE, NM, MVAL, NN, NVAL, NNS, NSVAL,
$ NNB, NBVAL, NXVAL, THRESH, A( 1, 1 ),
$ A( 1, 2 ), B( 1, 1 ), B( 1, 2 ),
$ S( 1 ), B( 1, 4 ),
$ WORK, RWORK, IWORK, NOUT )
ELSE
WRITE( NOUT, FMT = 9989 )PATH
END IF
* *
ELSE IF( LSAMEN( 2, C2, 'LS' ) ) THEN ELSE IF( LSAMEN( 2, C2, 'LS' ) ) THEN
* *

View File

@ -0,0 +1,836 @@
*> \brief \b ZCHKQP3RK
*
* =========== DOCUMENTATION ===========
*
* Online html documentation available at
* http://www.netlib.org/lapack/explore-html/
*
* Definition:
* ===========
*
* SUBROUTINE ZCHKQP3RK( DOTYPE, NM, MVAL, NN, NVAL, NNS, NSVAL,
* $ NNB, NBVAL, NXVAL, THRESH, A, COPYA,
* $ B, COPYB, S, TAU,
* $ WORK, RWORK, IWORK, NOUT )
* IMPLICIT NONE
*
* .. Scalar Arguments ..
* INTEGER NM, NN, NNB, NOUT
* DOUBLE PRECISION THRESH
* ..
* .. Array Arguments ..
* LOGICAL DOTYPE( * )
* INTEGER IWORK( * ), MVAL( * ), NBVAL( * ), NVAL( * ),
* $ NXVAL( * )
* DOUBLE PRECISION S( * ), RWORK( * )
* COMPLEX*16 A( * ), COPYA( * ), TAU( * ), WORK( * )
* ..
*
*
*> \par Purpose:
* =============
*>
*> \verbatim
*>
*> ZCHKQP3RK tests ZGEQP3RK.
*> \endverbatim
*
* Arguments:
* ==========
*
*> \param[in] DOTYPE
*> \verbatim
*> DOTYPE is LOGICAL array, dimension (NTYPES)
*> The matrix types to be used for testing. Matrices of type j
*> (for 1 <= j <= NTYPES) are used for testing if DOTYPE(j) =
*> .TRUE.; if DOTYPE(j) = .FALSE., then type j is not used.
*> \endverbatim
*>
*> \param[in] NM
*> \verbatim
*> NM is INTEGER
*> The number of values of M contained in the vector MVAL.
*> \endverbatim
*>
*> \param[in] MVAL
*> \verbatim
*> MVAL is INTEGER array, dimension (NM)
*> The values of the matrix row dimension M.
*> \endverbatim
*>
*> \param[in] NN
*> \verbatim
*> NN is INTEGER
*> The number of values of N contained in the vector NVAL.
*> \endverbatim
*>
*> \param[in] NVAL
*> \verbatim
*> NVAL is INTEGER array, dimension (NN)
*> The values of the matrix column dimension N.
*> \endverbatim
*>
*> \param[in] NNS
*> \verbatim
*> NNS is INTEGER
*> The number of values of NRHS contained in the vector NSVAL.
*> \endverbatim
*>
*> \param[in] NSVAL
*> \verbatim
*> NSVAL is INTEGER array, dimension (NNS)
*> The values of the number of right hand sides NRHS.
*> \endverbatim
*> \param[in] NNB
*> \verbatim
*> NNB is INTEGER
*> The number of values of NB and NX contained in the
*> vectors NBVAL and NXVAL. The blocking parameters are used
*> in pairs (NB,NX).
*> \endverbatim
*>
*> \param[in] NBVAL
*> \verbatim
*> NBVAL is INTEGER array, dimension (NNB)
*> The values of the blocksize NB.
*> \endverbatim
*>
*> \param[in] NXVAL
*> \verbatim
*> NXVAL is INTEGER array, dimension (NNB)
*> The values of the crossover point NX.
*> \endverbatim
*>
*> \param[in] THRESH
*> \verbatim
*> THRESH is DOUBLE PRECISION
*> The threshold value for the test ratios. A result is
*> included in the output file if RESULT >= THRESH. To have
*> every test ratio printed, use THRESH = 0.
*> \endverbatim
*>
*> \param[out] A
*> \verbatim
*> A is COMPLEX*16 array, dimension (MMAX*NMAX)
*> where MMAX is the maximum value of M in MVAL and NMAX is the
*> maximum value of N in NVAL.
*> \endverbatim
*>
*> \param[out] COPYA
*> \verbatim
*> COPYA is COMPLEX*16 array, dimension (MMAX*NMAX)
*> \endverbatim
*>
*> \param[out] B
*> \verbatim
*> B is COMPLEX*16 array, dimension (MMAX*NSMAX)
*> where MMAX is the maximum value of M in MVAL and NSMAX is the
*> maximum value of NRHS in NSVAL.
*> \endverbatim
*>
*> \param[out] COPYB
*> \verbatim
*> COPYB is COMPLEX*16 array, dimension (MMAX*NSMAX)
*> \endverbatim
*>
*> \param[out] S
*> \verbatim
*> S is DOUBLE PRECISION array, dimension
*> (min(MMAX,NMAX))
*> \endverbatim
*>
*> \param[out] TAU
*> \verbatim
*> TAU is COMPLEX*16 array, dimension (MMAX)
*> \endverbatim
*>
*> \param[out] WORK
*> \verbatim
*> WORK is COMPLEX*16 array, dimension
*> (max(M*max(M,N) + 4*min(M,N) + max(M,N)))
*> \endverbatim
*>
*> \param[out] RWORK
*> \verbatim
*> RWORK is DOUBLE PRECISION array, dimension (4*NMAX)
*> \endverbatim
*>
*> \param[out] IWORK
*> \verbatim
*> IWORK is INTEGER array, dimension (2*NMAX)
*> \endverbatim
*>
*> \param[in] NOUT
*> \verbatim
*> NOUT is INTEGER
*> The unit number for output.
*> \endverbatim
*
* Authors:
* ========
*
*> \author Univ. of Tennessee
*> \author Univ. of California Berkeley
*> \author Univ. of Colorado Denver
*> \author NAG Ltd.
*
*> \ingroup complex16_lin
*
* =====================================================================
SUBROUTINE ZCHKQP3RK( DOTYPE, NM, MVAL, NN, NVAL, NNS, NSVAL,
$ NNB, NBVAL, NXVAL, THRESH, A, COPYA,
$ B, COPYB, S, TAU,
$ WORK, RWORK, IWORK, NOUT )
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 NM, NN, NNB, NNS, NOUT
DOUBLE PRECISION THRESH
* ..
* .. Array Arguments ..
LOGICAL DOTYPE( * )
INTEGER IWORK( * ), NBVAL( * ), MVAL( * ), NVAL( * ),
$ NSVAL( * ), NXVAL( * )
DOUBLE PRECISION S( * ), RWORK( * )
COMPLEX*16 A( * ), COPYA( * ), B( * ), COPYB( * ),
$ TAU( * ), WORK( * )
* ..
*
* =====================================================================
*
* .. Parameters ..
INTEGER NTYPES
PARAMETER ( NTYPES = 19 )
INTEGER NTESTS
PARAMETER ( NTESTS = 5 )
DOUBLE PRECISION ONE, ZERO, BIGNUM
COMPLEX*16 CONE, CZERO
PARAMETER ( ONE = 1.0D+0, ZERO = 0.0D+0,
$ CZERO = ( 0.0D+0, 0.0D+0 ),
$ CONE = ( 1.0D+0, 0.0D+0 ),
$ BIGNUM = 1.0D+38 )
* ..
* .. Local Scalars ..
CHARACTER DIST, TYPE
CHARACTER*3 PATH
INTEGER I, IHIGH, ILOW, IM, IMAT, IN, INC_ZERO,
$ INB, IND_OFFSET_GEN,
$ IND_IN, IND_OUT, INS, INFO,
$ ISTEP, J, J_INC, J_FIRST_NZ, JB_ZERO,
$ KFACT, KL, KMAX, KU, LDA, LW, LWORK,
$ LWORK_MQR, M, MINMN, MINMNB_GEN, MODE, N,
$ NB, NB_ZERO, NERRS, NFAIL, NB_GEN, NRHS,
$ NRUN, NX, T
DOUBLE PRECISION ANORM, CNDNUM, EPS, ABSTOL, RELTOL,
$ DTEMP, MAXC2NRMK, RELMAXC2NRMK
* ..
* .. Local Arrays ..
INTEGER ISEED( 4 ), ISEEDY( 4 )
DOUBLE PRECISION RESULT( NTESTS ), RDUMMY( 1 )
* ..
* .. External Functions ..
DOUBLE PRECISION DLAMCH, ZQPT01, ZQRT11, ZQRT12, ZLANGE
EXTERNAL DLAMCH, ZQPT01, ZQRT11, ZQRT12, ZLANGE
* ..
* .. External Subroutines ..
EXTERNAL ALAERH, ALAHD, ALASUM, DLAORD, ICOPY, ZAXPY,
$ XLAENV, ZGEQP3RK, ZLACPY, ZLASET, ZLATB4,
$ ZLATMS, ZUNMQR, ZSWAP
* ..
* .. Intrinsic Functions ..
INTRINSIC ABS, DBLE, MAX, MIN, MOD
* ..
* .. Scalars in Common ..
LOGICAL LERR, OK
CHARACTER*32 SRNAMT
INTEGER INFOT, IOUNIT, ZUNMQR_LWORK
* ..
* .. Common blocks ..
COMMON / INFOC / INFOT, IOUNIT, OK, LERR
COMMON / SRNAMC / SRNAMT
* ..
* .. Data statements ..
DATA ISEEDY / 1988, 1989, 1990, 1991 /
* ..
* .. Executable Statements ..
*
* Initialize constants and the random number seed.
*
PATH( 1: 1 ) = 'Zomplex precision'
PATH( 2: 3 ) = 'QK'
NRUN = 0
NFAIL = 0
NERRS = 0
DO I = 1, 4
ISEED( I ) = ISEEDY( I )
END DO
EPS = DLAMCH( 'Epsilon' )
INFOT = 0
*
DO IM = 1, NM
*
* Do for each value of M in MVAL.
*
M = MVAL( IM )
LDA = MAX( 1, M )
*
DO IN = 1, NN
*
* Do for each value of N in NVAL.
*
N = NVAL( IN )
MINMN = MIN( M, N )
LWORK = MAX( 1, M*MAX( M, N )+4*MINMN+MAX( M, N ),
$ M*N + 2*MINMN + 4*N )
*
DO INS = 1, NNS
NRHS = NSVAL( INS )
*
* Set up parameters with ZLATB4 and generate
* M-by-NRHS B matrix with ZLATMS.
* IMAT = 14:
* Random matrix, CNDNUM = 2, NORM = ONE,
* MODE = 3 (geometric distribution of singular values).
*
CALL ZLATB4( PATH, 14, M, NRHS, TYPE, KL, KU, ANORM,
$ MODE, CNDNUM, DIST )
*
SRNAMT = 'ZLATMS'
CALL ZLATMS( M, NRHS, DIST, ISEED, TYPE, S, MODE,
$ CNDNUM, ANORM, KL, KU, 'No packing',
$ COPYB, LDA, WORK, INFO )
*
* Check error code from ZLATMS.
*
IF( INFO.NE.0 ) THEN
CALL ALAERH( PATH, 'ZLATMS', INFO, 0, ' ', M,
$ NRHS, -1, -1, -1, 6, NFAIL, NERRS,
$ NOUT )
CYCLE
END IF
*
DO IMAT = 1, NTYPES
*
* Do the tests only if DOTYPE( IMAT ) is true.
*
IF( .NOT.DOTYPE( IMAT ) )
$ CYCLE
*
* The type of distribution used to generate the random
* eigen-/singular values:
* ( 'S' for symmetric distribution ) => UNIFORM( -1, 1 )
*
* Do for each type of NON-SYMMETRIC matrix: CNDNUM NORM MODE
* 1. Zero matrix
* 2. Random, Diagonal, CNDNUM = 2 CNDNUM = 2 ONE 3 ( geometric distribution of singular values )
* 3. Random, Upper triangular, CNDNUM = 2 CNDNUM = 2 ONE 3 ( geometric distribution of singular values )
* 4. Random, Lower triangular, CNDNUM = 2 CNDNUM = 2 ONE 3 ( geometric distribution of singular values )
* 5. Random, First column is zero, CNDNUM = 2 CNDNUM = 2 ONE 3 ( geometric distribution of singular values )
* 6. Random, Last MINMN column is zero, CNDNUM = 2 CNDNUM = 2 ONE 3 ( geometric distribution of singular values )
* 7. Random, Last N column is zero, CNDNUM = 2 CNDNUM = 2 ONE 3 ( geometric distribution of singular values )
* 8. Random, Middle column in MINMN is zero, CNDNUM = 2 CNDNUM = 2 ONE 3 ( geometric distribution of singular values )
* 9. Random, First half of MINMN columns are zero, CNDNUM = 2 CNDNUM = 2 ONE 3 ( geometric distribution of singular values )
* 10. Random, Last columns are zero starting from MINMN/2+1, CNDNUM = 2 CNDNUM = 2 ONE 3 ( geometric distribution of singular values )
* 11. Random, Half MINMN columns in the middle are zero starting
* from MINMN/2-(MINMN/2)/2+1, CNDNUM = 2 CNDNUM = 2 ONE 3 ( geometric distribution of singular values )
* 12. Random, Odd columns are ZERO, CNDNUM = 2 CNDNUM = 2 ONE 3 ( geometric distribution of singular values )
* 13. Random, Even columns are ZERO, CNDNUM = 2 CNDNUM = 2 ONE 3 ( geometric distribution of singular values )
* 14. Random, CNDNUM = 2 CNDNUM = 2 ONE 3 ( geometric distribution of singular values )
* 15. Random, CNDNUM = sqrt(0.1/EPS) CNDNUM = BADC1 = sqrt(0.1/EPS) ONE 3 ( geometric distribution of singular values )
* 16. Random, CNDNUM = 0.1/EPS CNDNUM = BADC2 = 0.1/EPS ONE 3 ( geometric distribution of singular values )
* 17. Random, CNDNUM = 0.1/EPS, CNDNUM = BADC2 = 0.1/EPS ONE 2 ( one small singular value, S(N)=1/CNDNUM )
* one small singular value S(N)=1/CNDNUM
* 18. Random, CNDNUM = 2, scaled near underflow CNDNUM = 2 SMALL = SAFMIN
* 19. Random, CNDNUM = 2, scaled near overflow CNDNUM = 2 LARGE = 1.0/( 0.25 * ( SAFMIN / EPS ) ) 3 ( geometric distribution of singular values )
*
IF( IMAT.EQ.1 ) THEN
*
* Matrix 1: Zero matrix
*
CALL ZLASET( 'Full', M, N, CZERO, CZERO, COPYA, LDA )
DO I = 1, MINMN
S( I ) = ZERO
END DO
*
ELSE IF( (IMAT.GE.2 .AND. IMAT.LE.4 )
$ .OR. (IMAT.GE.14 .AND. IMAT.LE.19 ) ) THEN
*
* Matrices 2-5.
*
* Set up parameters with DLATB4 and generate a test
* matrix with ZLATMS.
*
CALL ZLATB4( PATH, IMAT, M, N, TYPE, KL, KU, ANORM,
$ MODE, CNDNUM, DIST )
*
SRNAMT = 'ZLATMS'
CALL ZLATMS( M, N, DIST, ISEED, TYPE, S, MODE,
$ CNDNUM, ANORM, KL, KU, 'No packing',
$ COPYA, LDA, WORK, INFO )
*
* Check error code from ZLATMS.
*
IF( INFO.NE.0 ) THEN
CALL ALAERH( PATH, 'ZLATMS', INFO, 0, ' ', M, N,
$ -1, -1, -1, IMAT, NFAIL, NERRS,
$ NOUT )
CYCLE
END IF
*
CALL DLAORD( 'Decreasing', MINMN, S, 1 )
*
ELSE IF( MINMN.GE.2
$ .AND. IMAT.GE.5 .AND. IMAT.LE.13 ) THEN
*
* Rectangular matrices 5-13 that contain zero columns,
* only for matrices MINMN >=2.
*
* JB_ZERO is the column index of ZERO block.
* NB_ZERO is the column block size of ZERO block.
* NB_GEN is the column blcok size of the
* generated block.
* J_INC in the non_zero column index increment
* for matrix 12 and 13.
* J_FIRS_NZ is the index of the first non-zero
* column.
*
IF( IMAT.EQ.5 ) THEN
*
* First column is zero.
*
JB_ZERO = 1
NB_ZERO = 1
NB_GEN = N - NB_ZERO
*
ELSE IF( IMAT.EQ.6 ) THEN
*
* Last column MINMN is zero.
*
JB_ZERO = MINMN
NB_ZERO = 1
NB_GEN = N - NB_ZERO
*
ELSE IF( IMAT.EQ.7 ) THEN
*
* Last column N is zero.
*
JB_ZERO = N
NB_ZERO = 1
NB_GEN = N - NB_ZERO
*
ELSE IF( IMAT.EQ.8 ) THEN
*
* Middle column in MINMN is zero.
*
JB_ZERO = MINMN / 2 + 1
NB_ZERO = 1
NB_GEN = N - NB_ZERO
*
ELSE IF( IMAT.EQ.9 ) THEN
*
* First half of MINMN columns is zero.
*
JB_ZERO = 1
NB_ZERO = MINMN / 2
NB_GEN = N - NB_ZERO
*
ELSE IF( IMAT.EQ.10 ) THEN
*
* Last columns are zero columns,
* starting from (MINMN / 2 + 1) column.
*
JB_ZERO = MINMN / 2 + 1
NB_ZERO = N - JB_ZERO + 1
NB_GEN = N - NB_ZERO
*
ELSE IF( IMAT.EQ.11 ) THEN
*
* Half of the columns in the middle of MINMN
* columns is zero, starting from
* MINMN/2 - (MINMN/2)/2 + 1 column.
*
JB_ZERO = MINMN / 2 - (MINMN / 2) / 2 + 1
NB_ZERO = MINMN / 2
NB_GEN = N - NB_ZERO
*
ELSE IF( IMAT.EQ.12 ) THEN
*
* Odd-numbered columns are zero,
*
NB_GEN = N / 2
NB_ZERO = N - NB_GEN
J_INC = 2
J_FIRST_NZ = 2
*
ELSE IF( IMAT.EQ.13 ) THEN
*
* Even-numbered columns are zero.
*
NB_ZERO = N / 2
NB_GEN = N - NB_ZERO
J_INC = 2
J_FIRST_NZ = 1
*
END IF
*
*
* 1) Set the first NB_ZERO columns in COPYA(1:M,1:N)
* to zero.
*
CALL ZLASET( 'Full', M, NB_ZERO, CZERO, CZERO,
$ COPYA, LDA )
*
* 2) Generate an M-by-(N-NB_ZERO) matrix with the
* chosen singular value distribution
* in COPYA(1:M,NB_ZERO+1:N).
*
CALL ZLATB4( PATH, IMAT, M, NB_GEN, TYPE, KL, KU,
$ ANORM, MODE, CNDNUM, DIST )
*
SRNAMT = 'ZLATMS'
*
IND_OFFSET_GEN = NB_ZERO * LDA
*
CALL ZLATMS( M, NB_GEN, DIST, ISEED, TYPE, S, MODE,
$ CNDNUM, ANORM, KL, KU, 'No packing',
$ COPYA( IND_OFFSET_GEN + 1 ), LDA,
$ WORK, INFO )
*
* Check error code from ZLATMS.
*
IF( INFO.NE.0 ) THEN
CALL ALAERH( PATH, 'ZLATMS', INFO, 0, ' ', M,
$ NB_GEN, -1, -1, -1, IMAT, NFAIL,
$ NERRS, NOUT )
CYCLE
END IF
*
* 3) Swap the gererated colums from the right side
* NB_GEN-size block in COPYA into correct column
* positions.
*
IF( IMAT.EQ.6
$ .OR. IMAT.EQ.7
$ .OR. IMAT.EQ.8
$ .OR. IMAT.EQ.10
$ .OR. IMAT.EQ.11 ) THEN
*
* Move by swapping the generated columns
* from the right NB_GEN-size block from
* (NB_ZERO+1:NB_ZERO+JB_ZERO)
* into columns (1:JB_ZERO-1).
*
DO J = 1, JB_ZERO-1, 1
CALL ZSWAP( M,
$ COPYA( ( NB_ZERO+J-1)*LDA+1), 1,
$ COPYA( (J-1)*LDA + 1 ), 1 )
END DO
*
ELSE IF( IMAT.EQ.12 .OR. IMAT.EQ.13 ) THEN
*
* ( IMAT = 12, Odd-numbered ZERO columns. )
* Swap the generated columns from the right
* NB_GEN-size block into the even zero colums in the
* left NB_ZERO-size block.
*
* ( IMAT = 13, Even-numbered ZERO columns. )
* Swap the generated columns from the right
* NB_GEN-size block into the odd zero colums in the
* left NB_ZERO-size block.
*
DO J = 1, NB_GEN, 1
IND_OUT = ( NB_ZERO+J-1 )*LDA + 1
IND_IN = ( J_INC*(J-1)+(J_FIRST_NZ-1) )*LDA
$ + 1
CALL ZSWAP( M,
$ COPYA( IND_OUT ), 1,
$ COPYA( IND_IN), 1 )
END DO
*
END IF
*
* 5) Order the singular values generated by
* DLAMTS in decreasing order and add trailing zeros
* that correspond to zero columns.
* The total number of singular values is MINMN.
*
MINMNB_GEN = MIN( M, NB_GEN )
*
CALL DLAORD( 'Decreasing', MINMNB_GEN, S, 1 )
DO I = MINMNB_GEN+1, MINMN
S( I ) = ZERO
END DO
*
ELSE
*
* IF(MINMN.LT.2) skip this size for this matrix type.
*
CYCLE
END IF
*
* Initialize a copy array for a pivot array for DGEQP3RK.
*
DO I = 1, N
IWORK( I ) = 0
END DO
*
DO INB = 1, NNB
*
* Do for each pair of values (NB,NX) in NBVAL and NXVAL.
*
NB = NBVAL( INB )
CALL XLAENV( 1, NB )
NX = NXVAL( INB )
CALL XLAENV( 3, NX )
*
* We do MIN(M,N)+1 because we need a test for KMAX > N,
* when KMAX is larger than MIN(M,N), KMAX should be
* KMAX = MIN(M,N)
*
DO KMAX = 0, MIN(M,N)+1
*
* Get a working copy of COPYA into A( 1:M,1:N ).
* Get a working copy of COPYB into A( 1:M, (N+1):NRHS ).
* Get a working copy of COPYB into into B( 1:M, 1:NRHS ).
* Get a working copy of IWORK(1:N) awith zeroes into
* which is going to be used as pivot array IWORK( N+1:2N ).
* NOTE: IWORK(2N+1:3N) is going to be used as a WORK array
* for the routine.
*
CALL ZLACPY( 'All', M, N, COPYA, LDA, A, LDA )
CALL ZLACPY( 'All', M, NRHS, COPYB, LDA,
$ A( LDA*N + 1 ), LDA )
CALL ZLACPY( 'All', M, NRHS, COPYB, LDA,
$ B, LDA )
CALL ICOPY( N, IWORK( 1 ), 1, IWORK( N+1 ), 1 )
*
ABSTOL = -1.0
RELTOl = -1.0
*
* Compute the QR factorization with pivoting of A
*
LW = MAX( 1, MAX( 2*N + NB*( N+NRHS+1 ),
$ 3*N + NRHS - 1 ) )
*
* Compute ZGEQP3RK factorization of A.
*
SRNAMT = 'ZGEQP3RK'
CALL ZGEQP3RK( M, N, NRHS, KMAX, ABSTOL, RELTOL,
$ A, LDA, KFACT, MAXC2NRMK,
$ RELMAXC2NRMK, IWORK( N+1 ), TAU,
$ WORK, LW, RWORK, IWORK( 2*N+1 ),
$ INFO )
*
* Check error code from ZGEQP3RK.
*
IF( INFO.LT.0 )
$ CALL ALAERH( PATH, 'ZGEQP3RK', INFO, 0, ' ',
$ M, N, NX, -1, NB, IMAT,
$ NFAIL, NERRS, NOUT )
*
IF( KFACT.EQ.MINMN ) THEN
*
* Compute test 1:
*
* This test in only for the full rank factorization of
* the matrix A.
*
* Array S(1:min(M,N)) contains svd(A) the sigular values
* of the original matrix A in decreasing absolute value
* order. The test computes svd(R), the vector sigular
* values of the upper trapezoid of A(1:M,1:N) that
* contains the factor R, in decreasing order. The test
* returns the ratio:
*
* 2-norm(svd(R) - svd(A)) / ( max(M,N) * 2-norm(svd(A)) * EPS )
*
RESULT( 1 ) = ZQRT12( M, N, A, LDA, S, WORK,
$ LWORK , RWORK )
*
DO T = 1, 1
IF( RESULT( T ).GE.THRESH ) THEN
IF( NFAIL.EQ.0 .AND. NERRS.EQ.0 )
$ CALL ALAHD( NOUT, PATH )
WRITE( NOUT, FMT = 9999 ) 'ZGEQP3RK', M, N,
$ NRHS, KMAX, ABSTOL, RELTOL, NB, NX,
$ IMAT, T, RESULT( T )
NFAIL = NFAIL + 1
END IF
END DO
NRUN = NRUN + 1
*
* End test 1
*
END IF
* Compute test 2:
*
* The test returns the ratio:
*
* 1-norm( A*P - Q*R ) / ( max(M,N) * 1-norm(A) * EPS )
*
RESULT( 2 ) = ZQPT01( M, N, KFACT, COPYA, A, LDA, TAU,
$ IWORK( N+1 ), WORK, LWORK )
*
* Compute test 3:
*
* The test returns the ratio:
*
* 1-norm( Q**T * Q - I ) / ( M * EPS )
*
RESULT( 3 ) = ZQRT11( M, KFACT, A, LDA, TAU, WORK,
$ LWORK )
*
* Print information about the tests that did not pass
* the threshold.
*
DO T = 2, 3
IF( RESULT( T ).GE.THRESH ) THEN
IF( NFAIL.EQ.0 .AND. NERRS.EQ.0 )
$ CALL ALAHD( NOUT, PATH )
WRITE( NOUT, FMT = 9999 ) 'ZGEQP3RK', M, N,
$ NRHS, KMAX, ABSTOL, RELTOL,
$ NB, NX, IMAT, T, RESULT( T )
NFAIL = NFAIL + 1
END IF
END DO
NRUN = NRUN + 2
*
* Compute test 4:
*
* This test is only for the factorizations with the
* rank greater than 2.
* The elements on the diagonal of R should be non-
* increasing.
*
* The test returns the ratio:
*
* Returns 1.0D+100 if abs(R(K+1,K+1)) > abs(R(K,K)),
* K=1:KFACT-1
*
IF( MIN(KFACT, MINMN).GE.2 ) THEN
*
DO J = 1, KFACT-1, 1
*
DTEMP = (( ABS( A( (J-1)*M+J ) ) -
$ ABS( A( (J)*M+J+1 ) ) ) /
$ ABS( A(1) ) )
*
IF( DTEMP.LT.ZERO ) THEN
RESULT( 4 ) = BIGNUM
END IF
*
END DO
*
* Print information about the tests that did not
* pass the threshold.
*
DO T = 4, 4
IF( RESULT( T ).GE.THRESH ) THEN
IF( NFAIL.EQ.0 .AND. NERRS.EQ.0 )
$ CALL ALAHD( NOUT, PATH )
WRITE( NOUT, FMT = 9999 ) 'ZGEQP3RK',
$ M, N, NRHS, KMAX, ABSTOL, RELTOL,
$ NB, NX, IMAT, T,
$ RESULT( T )
NFAIL = NFAIL + 1
END IF
END DO
NRUN = NRUN + 1
*
* End test 4.
*
END IF
*
* Compute test 5:
*
* This test in only for matrix A with min(M,N) > 0.
*
* The test returns the ratio:
*
* 1-norm(Q**T * B - Q**T * B ) /
* ( M * EPS )
*
* (1) Compute B:=Q**T * B in the matrix B.
*
IF( MINMN.GT.0 ) THEN
*
LWORK_MQR = MAX(1, NRHS)
CALL ZUNMQR( 'Left', 'Conjugate transpose',
$ M, NRHS, KFACT, A, LDA, TAU, B, LDA,
$ WORK, LWORK_MQR, INFO )
*
DO I = 1, NRHS
*
* Compare N+J-th column of A and J-column of B.
*
CALL ZAXPY( M, -CONE, A( ( N+I-1 )*LDA+1 ), 1,
$ B( ( I-1 )*LDA+1 ), 1 )
END DO
*
RESULT( 5 ) =
$ ABS(
$ ZLANGE( 'One-norm', M, NRHS, B, LDA, RDUMMY ) /
$ ( DBLE( M )*DLAMCH( 'Epsilon' ) )
$ )
*
* Print information about the tests that did not pass
* the threshold.
*
DO T = 5, 5
IF( RESULT( T ).GE.THRESH ) THEN
IF( NFAIL.EQ.0 .AND. NERRS.EQ.0 )
$ CALL ALAHD( NOUT, PATH )
WRITE( NOUT, FMT = 9999 ) 'ZGEQP3RK', M, N,
$ NRHS, KMAX, ABSTOL, RELTOL,
$ NB, NX, IMAT, T, RESULT( T )
NFAIL = NFAIL + 1
END IF
END DO
NRUN = NRUN + 1
*
* End compute test 5.
*
END IF
*
* END DO KMAX = 1, MIN(M,N)+1
*
END DO
*
* END DO for INB = 1, NNB
*
END DO
*
* END DO for IMAT = 1, NTYPES
*
END DO
*
* END DO for INS = 1, NNS
*
END DO
*
* END DO for IN = 1, NN
*
END DO
*
* END DO for IM = 1, NM
*
END DO
*
* Print a summary of the results.
*
CALL ALASUM( PATH, NOUT, NFAIL, NRUN, NERRS )
*
9999 FORMAT( 1X, A, ' M =', I5, ', N =', I5, ', NRHS =', I5,
$ ', KMAX =', I5, ', ABSTOL =', G12.5,
$ ', RELTOL =', G12.5, ', NB =', I4, ', NX =', I4,
$ ', type ', I2, ', test ', I2, ', ratio =', G12.5 )
*
* End of ZCHKQP3RK
*
END

View File

@ -154,9 +154,6 @@
* .. Intrinsic Functions .. * .. Intrinsic Functions ..
INTRINSIC ABS, MAX, SQRT INTRINSIC ABS, MAX, SQRT
* .. * ..
* .. External Subroutines ..
EXTERNAL DLABAD
* ..
* .. Save statement .. * .. Save statement ..
SAVE EPS, SMALL, LARGE, BADC1, BADC2, FIRST SAVE EPS, SMALL, LARGE, BADC1, BADC2, FIRST
* .. * ..
@ -174,11 +171,6 @@
BADC1 = SQRT( BADC2 ) BADC1 = SQRT( BADC2 )
SMALL = DLAMCH( 'Safe minimum' ) SMALL = DLAMCH( 'Safe minimum' )
LARGE = ONE / SMALL LARGE = ONE / SMALL
*
* If it looks like we're on a Cray, take the square root of
* SMALL and LARGE to avoid overflow and underflow problems.
*
CALL DLABAD( SMALL, LARGE )
SMALL = SHRINK*( SMALL / EPS ) SMALL = SHRINK*( SMALL / EPS )
LARGE = ONE / SMALL LARGE = ONE / SMALL
END IF END IF
@ -233,6 +225,110 @@
ELSE ELSE
ANORM = ONE ANORM = ONE
END IF END IF
*
ELSE IF( LSAMEN( 2, C2, 'QK' ) ) THEN
*
* xQK: truncated QR with pivoting.
* Set parameters to generate a general
* M x N matrix.
*
* Set TYPE, the type of matrix to be generated. 'N' is nonsymmetric.
*
TYPE = 'N'
*
* Set DIST, the type of distribution for the random
* number generator. 'S' is
*
DIST = 'S'
*
* Set the lower and upper bandwidths.
*
IF( IMAT.EQ.2 ) THEN
*
* 2. Random, Diagonal, CNDNUM = 2
*
KL = 0
KU = 0
CNDNUM = TWO
ANORM = ONE
MODE = 3
ELSE IF( IMAT.EQ.3 ) THEN
*
* 3. Random, Upper triangular, CNDNUM = 2
*
KL = 0
KU = MAX( N-1, 0 )
CNDNUM = TWO
ANORM = ONE
MODE = 3
ELSE IF( IMAT.EQ.4 ) THEN
*
* 4. Random, Lower triangular, CNDNUM = 2
*
KL = MAX( M-1, 0 )
KU = 0
CNDNUM = TWO
ANORM = ONE
MODE = 3
ELSE
*
* 5.-19. Rectangular matrix
*
KL = MAX( M-1, 0 )
KU = MAX( N-1, 0 )
*
IF( IMAT.GE.5 .AND. IMAT.LE.14 ) THEN
*
* 5.-14. Random, CNDNUM = 2.
*
CNDNUM = TWO
ANORM = ONE
MODE = 3
*
ELSE IF( IMAT.EQ.15 ) THEN
*
* 15. Random, CNDNUM = sqrt(0.1/EPS)
*
CNDNUM = BADC1
ANORM = ONE
MODE = 3
*
ELSE IF( IMAT.EQ.16 ) THEN
*
* 16. Random, CNDNUM = 0.1/EPS
*
CNDNUM = BADC2
ANORM = ONE
MODE = 3
*
ELSE IF( IMAT.EQ.17 ) THEN
*
* 17. Random, CNDNUM = 0.1/EPS,
* one small singular value S(N)=1/CNDNUM
*
CNDNUM = BADC2
ANORM = ONE
MODE = 2
*
ELSE IF( IMAT.EQ.18 ) THEN
*
* 18. Random, scaled near underflow
*
CNDNUM = TWO
ANORM = SMALL
MODE = 3
*
ELSE IF( IMAT.EQ.19 ) THEN
*
* 19. Random, scaled near overflow
*
CNDNUM = TWO
ANORM = LARGE
MODE = 3
*
END IF
*
END IF
* *
ELSE IF( LSAMEN( 2, C2, 'GE' ) ) THEN ELSE IF( LSAMEN( 2, C2, 'GE' ) ) THEN
* *
@ -517,17 +613,18 @@
* *
* Set the norm and condition number. * Set the norm and condition number.
* *
IF( IMAT.EQ.2 .OR. IMAT.EQ.8 ) THEN MAT = ABS( IMAT )
IF( MAT.EQ.2 .OR. MAT.EQ.8 ) THEN
CNDNUM = BADC1 CNDNUM = BADC1
ELSE IF( IMAT.EQ.3 .OR. IMAT.EQ.9 ) THEN ELSE IF( MAT.EQ.3 .OR. MAT.EQ.9 ) THEN
CNDNUM = BADC2 CNDNUM = BADC2
ELSE ELSE
CNDNUM = TWO CNDNUM = TWO
END IF END IF
* *
IF( IMAT.EQ.4 ) THEN IF( MAT.EQ.4 ) THEN
ANORM = SMALL ANORM = SMALL
ELSE IF( IMAT.EQ.5 ) THEN ELSE IF( MAT.EQ.5 ) THEN
ANORM = LARGE ANORM = LARGE
ELSE ELSE
ANORM = ONE ANORM = ONE

View File

@ -33,7 +33,7 @@
*> Householder vectors, and the rest of AF contains a partially updated *> Householder vectors, and the rest of AF contains a partially updated
*> matrix. *> matrix.
*> *>
*> This function returns ||A*P - Q*R||/(||norm(A)||*eps*M) *> This function returns ||A*P - Q*R|| / ( ||norm(A)||*eps*max(M,N) )
*> \endverbatim *> \endverbatim
* *
* Arguments: * Arguments:
@ -172,28 +172,28 @@
* *
NORMA = ZLANGE( 'One-norm', M, N, A, LDA, RWORK ) NORMA = ZLANGE( 'One-norm', M, N, A, LDA, RWORK )
* *
DO 30 J = 1, K DO J = 1, K
DO 10 I = 1, MIN( J, M ) DO I = 1, MIN( J, M )
WORK( ( J-1 )*M+I ) = AF( I, J ) WORK( ( J-1 )*M+I ) = AF( I, J )
10 CONTINUE END DO
DO 20 I = J + 1, M DO I = J + 1, M
WORK( ( J-1 )*M+I ) = ZERO WORK( ( J-1 )*M+I ) = ZERO
20 CONTINUE END DO
30 CONTINUE END DO
DO 40 J = K + 1, N DO J = K + 1, N
CALL ZCOPY( M, AF( 1, J ), 1, WORK( ( J-1 )*M+1 ), 1 ) CALL ZCOPY( M, AF( 1, J ), 1, WORK( ( J-1 )*M+1 ), 1 )
40 CONTINUE END DO
* *
CALL ZUNMQR( 'Left', 'No transpose', M, N, K, AF, LDA, TAU, WORK, CALL ZUNMQR( 'Left', 'No transpose', M, N, K, AF, LDA, TAU, WORK,
$ M, WORK( M*N+1 ), LWORK-M*N, INFO ) $ M, WORK( M*N+1 ), LWORK-M*N, INFO )
* *
DO 50 J = 1, N DO J = 1, N
* *
* Compare i-th column of QR and jpvt(i)-th column of A * Compare i-th column of QR and jpvt(i)-th column of A
* *
CALL ZAXPY( M, DCMPLX( -ONE ), A( 1, JPVT( J ) ), 1, CALL ZAXPY( M, DCMPLX( -ONE ), A( 1, JPVT( J ) ), 1,
$ WORK( ( J-1 )*M+1 ), 1 ) $ WORK( ( J-1 )*M+1 ), 1 )
50 CONTINUE END DO
* *
ZQPT01 = ZLANGE( 'One-norm', M, N, WORK, M, RWORK ) / ZQPT01 = ZLANGE( 'One-norm', M, N, WORK, M, RWORK ) /
$ ( DBLE( MAX( M, N ) )*DLAMCH( 'Epsilon' ) ) $ ( DBLE( MAX( M, N ) )*DLAMCH( 'Epsilon' ) )

View File

@ -158,9 +158,9 @@
CALL ZUNM2R( 'Left', 'Conjugate transpose', M, M, K, A, LDA, TAU, CALL ZUNM2R( 'Left', 'Conjugate transpose', M, M, K, A, LDA, TAU,
$ WORK, M, WORK( M*M+1 ), INFO ) $ WORK, M, WORK( M*M+1 ), INFO )
* *
DO 10 J = 1, M DO J = 1, M
WORK( ( J-1 )*M+J ) = WORK( ( J-1 )*M+J ) - ONE WORK( ( J-1 )*M+J ) = WORK( ( J-1 )*M+J ) - ONE
10 CONTINUE END DO
* *
ZQRT11 = ZLANGE( 'One-norm', M, M, WORK, M, RDUMMY ) / ZQRT11 = ZLANGE( 'One-norm', M, M, WORK, M, RDUMMY ) /
$ ( DBLE( M )*DLAMCH( 'Epsilon' ) ) $ ( DBLE( M )*DLAMCH( 'Epsilon' ) )

View File

@ -28,7 +28,7 @@
*> ZQRT12 computes the singular values `svlues' of the upper trapezoid *> ZQRT12 computes the singular values `svlues' of the upper trapezoid
*> of A(1:M,1:N) and returns the ratio *> of A(1:M,1:N) and returns the ratio
*> *>
*> || s - svlues||/(||svlues||*eps*max(M,N)) *> || svlues - s||/(||s||*eps*max(M,N))
*> \endverbatim *> \endverbatim
* *
* Arguments: * Arguments:
@ -125,8 +125,8 @@
EXTERNAL DASUM, DLAMCH, DNRM2, ZLANGE EXTERNAL DASUM, DLAMCH, DNRM2, ZLANGE
* .. * ..
* .. External Subroutines .. * .. External Subroutines ..
EXTERNAL DAXPY, DBDSQR, DLABAD, DLASCL, XERBLA, ZGEBD2, EXTERNAL DAXPY, DBDSQR, DLASCL, XERBLA, ZGEBD2, ZLASCL,
$ ZLASCL, ZLASET $ ZLASET
* .. * ..
* .. Intrinsic Functions .. * .. Intrinsic Functions ..
INTRINSIC DBLE, DCMPLX, MAX, MIN INTRINSIC DBLE, DCMPLX, MAX, MIN
@ -154,17 +154,16 @@
* *
CALL ZLASET( 'Full', M, N, DCMPLX( ZERO ), DCMPLX( ZERO ), WORK, CALL ZLASET( 'Full', M, N, DCMPLX( ZERO ), DCMPLX( ZERO ), WORK,
$ M ) $ M )
DO 20 J = 1, N DO J = 1, N
DO 10 I = 1, MIN( J, M ) DO I = 1, MIN( J, M )
WORK( ( J-1 )*M+I ) = A( I, J ) WORK( ( J-1 )*M+I ) = A( I, J )
10 CONTINUE END DO
20 CONTINUE END DO
* *
* Get machine parameters * Get machine parameters
* *
SMLNUM = DLAMCH( 'S' ) / DLAMCH( 'P' ) SMLNUM = DLAMCH( 'S' ) / DLAMCH( 'P' )
BIGNUM = ONE / SMLNUM BIGNUM = ONE / SMLNUM
CALL DLABAD( SMLNUM, BIGNUM )
* *
* Scale work if max entry outside range [SMLNUM,BIGNUM] * Scale work if max entry outside range [SMLNUM,BIGNUM]
* *
@ -208,9 +207,9 @@
* *
ELSE ELSE
* *
DO 30 I = 1, MN DO I = 1, MN
RWORK( I ) = ZERO RWORK( I ) = ZERO
30 CONTINUE END DO
END IF END IF
* *
* Compare s and singular values of work * Compare s and singular values of work
@ -218,6 +217,7 @@
CALL DAXPY( MN, -ONE, S, 1, RWORK( 1 ), 1 ) CALL DAXPY( MN, -ONE, S, 1, RWORK( 1 ), 1 )
ZQRT12 = DASUM( MN, RWORK( 1 ), 1 ) / ZQRT12 = DASUM( MN, RWORK( 1 ), 1 ) /
$ ( DLAMCH( 'Epsilon' )*DBLE( MAX( M, N ) ) ) $ ( DLAMCH( 'Epsilon' )*DBLE( MAX( M, N ) ) )
*
IF( NRMSVL.NE.ZERO ) IF( NRMSVL.NE.ZERO )
$ ZQRT12 = ZQRT12 / NRMSVL $ ZQRT12 = ZQRT12 / NRMSVL
* *