diff --git a/lapack-netlib/TESTING/LIN/CMakeLists.txt b/lapack-netlib/TESTING/LIN/CMakeLists.txt index 676857a80..e28818c76 100644 --- a/lapack-netlib/TESTING/LIN/CMakeLists.txt +++ b/lapack-netlib/TESTING/LIN/CMakeLists.txt @@ -9,7 +9,7 @@ set(DZLNTST dlaord.f) set(SLINTST schkaa.F schkeq.f schkgb.f schkge.f schkgt.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 schksy_aa.f schksy_aa_2stage.f schktb.f schktp.f schktr.f @@ -56,7 +56,7 @@ set(CLINTST cchkaa.F cchkhe.f cchkhe_rook.f cchkhe_rk.f cchkhe_aa.f cchkhe_aa_2stage.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 cchksy_aa.f cchksy_aa_2stage.f cchktb.f @@ -110,7 +110,7 @@ endif() set(DLINTST dchkaa.F dchkeq.f dchkgb.f dchkge.f dchkgt.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 dchksy_aa.f dchksy_aa_2stage.f dchktb.f dchktp.f dchktr.f @@ -158,7 +158,7 @@ set(ZLINTST zchkaa.F zchkhe.f zchkhe_rook.f zchkhe_rk.f zchkhe_aa.f zchkhe_aa_2stage.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 zchksy_aa.f zchksy_aa_2stage.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) add_executable(${name} ${ARGN}) - target_link_libraries(${name} openblas${SUFFIX64_UNDERSCORE}) -#${TMGLIB} ${LAPACK_LIBRARIES} ${BLAS_LIBRARIES}) + target_link_libraries(${name} ${TMGLIB} ${LAPACK_LIBRARIES} ${BLAS_LIBRARIES}) endmacro() if(BUILD_SINGLE) diff --git a/lapack-netlib/TESTING/LIN/Makefile b/lapack-netlib/TESTING/LIN/Makefile index 64abc4dba..46e096c2f 100644 --- a/lapack-netlib/TESTING/LIN/Makefile +++ b/lapack-netlib/TESTING/LIN/Makefile @@ -45,7 +45,7 @@ DZLNTST = dlaord.o SLINTST = schkaa.o \ schkeq.o schkgb.o schkge.o schkgt.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 \ schksy_aa.o schksy_aa_2stage.o schktb.o schktp.o schktr.o \ schktz.o \ @@ -89,7 +89,7 @@ CLINTST = cchkaa.o \ cchkeq.o cchkgb.o cchkge.o cchkgt.o \ cchkhe.o cchkhe_rook.o cchkhe_rk.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 \ cchksy_aa.o cchksy_aa_2stage.o cchktb.o \ cchktp.o cchktr.o cchktz.o \ @@ -137,7 +137,7 @@ endif DLINTST = dchkaa.o \ dchkeq.o dchkgb.o dchkge.o dchkgt.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 \ dchksy_aa.o dchksy_aa_2stage.o dchktb.o dchktp.o dchktr.o \ dchktz.o \ @@ -182,7 +182,7 @@ ZLINTST = zchkaa.o \ zchkeq.o zchkgb.o zchkge.o zchkgt.o \ zchkhe.o zchkhe_rook.o zchkhe_rk.o zchkhe_aa.o zchkhe_aa_2stage.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 \ zchksy_aa.o zchksy_aa_2stage.o zchktb.o \ zchktp.o zchktr.o zchktz.o \ @@ -269,35 +269,35 @@ proto-double: xlintstds xlintstrfd proto-complex: xlintstrfc proto-complex16: xlintstzc xlintstrfz -xlintsts: $(ALINTST) $(SLINTST) $(SCLNTST) $(TMGLIB) $(VARLIB) ../$(LAPACKLIB) $(XBLASLIB) $(BLASLIB) - $(LOADER) $(FFLAGS) $(LDFLAGS) -o $@ $^ +xlintsts: $(ALINTST) $(SLINTST) $(SCLNTST) $(TMGLIB) $(VARLIB) $(LAPACKLIB) $(XBLASLIB) $(BLASLIB) + $(FC) $(FFLAGS) $(LDFLAGS) -o $@ $^ -xlintstc: $(ALINTST) $(CLINTST) $(SCLNTST) $(TMGLIB) $(VARLIB) ../$(LAPACKLIB) $(XBLASLIB) $(BLASLIB) - $(LOADER) $(FFLAGS) $(LDFLAGS) -o $@ $^ +xlintstc: $(ALINTST) $(CLINTST) $(SCLNTST) $(TMGLIB) $(VARLIB) $(LAPACKLIB) $(XBLASLIB) $(BLASLIB) + $(FC) $(FFLAGS) $(LDFLAGS) -o $@ $^ -xlintstd: $(ALINTST) $(DLINTST) $(DZLNTST) $(TMGLIB) $(VARLIB) ../$(LAPACKLIB) $(XBLASLIB) $(BLASLIB) - $(LOADER) $(FFLAGS) $(LDFLAGS) -o $@ $^ +xlintstd: $(ALINTST) $(DLINTST) $(DZLNTST) $(TMGLIB) $(VARLIB) $(LAPACKLIB) $(XBLASLIB) $(BLASLIB) + $(FC) $(FFLAGS) $(LDFLAGS) -o $@ $^ -xlintstz: $(ALINTST) $(ZLINTST) $(DZLNTST) $(TMGLIB) $(VARLIB) ../$(LAPACKLIB) $(XBLASLIB) $(BLASLIB) - $(LOADER) $(FFLAGS) $(LDFLAGS) -o $@ $^ +xlintstz: $(ALINTST) $(ZLINTST) $(DZLNTST) $(TMGLIB) $(VARLIB) $(LAPACKLIB) $(XBLASLIB) $(BLASLIB) + $(FC) $(FFLAGS) $(LDFLAGS) -o $@ $^ -xlintstds: $(DSLINTST) $(TMGLIB) $(VARLIB) ../$(LAPACKLIB) $(BLASLIB) - $(LOADER) $(FFLAGS) $(LDFLAGS) -o $@ $^ +xlintstds: $(DSLINTST) $(TMGLIB) $(VARLIB) $(LAPACKLIB) $(BLASLIB) + $(FC) $(FFLAGS) $(LDFLAGS) -o $@ $^ -xlintstzc: $(ZCLINTST) $(TMGLIB) $(VARLIB) ../$(LAPACKLIB) $(BLASLIB) - $(LOADER) $(FFLAGS) $(LDFLAGS) -o $@ $^ +xlintstzc: $(ZCLINTST) $(TMGLIB) $(VARLIB) $(LAPACKLIB) $(BLASLIB) + $(FC) $(FFLAGS) $(LDFLAGS) -o $@ $^ -xlintstrfs: $(SLINTSTRFP) $(TMGLIB) $(VARLIB) ../$(LAPACKLIB) $(BLASLIB) - $(LOADER) $(FFLAGS) $(LDFLAGS) -o $@ $^ +xlintstrfs: $(SLINTSTRFP) $(TMGLIB) $(VARLIB) $(LAPACKLIB) $(BLASLIB) + $(FC) $(FFLAGS) $(LDFLAGS) -o $@ $^ -xlintstrfd: $(DLINTSTRFP) $(TMGLIB) $(VARLIB) ../$(LAPACKLIB) $(BLASLIB) - $(LOADER) $(FFLAGS) $(LDFLAGS) -o $@ $^ +xlintstrfd: $(DLINTSTRFP) $(TMGLIB) $(VARLIB) $(LAPACKLIB) $(BLASLIB) + $(FC) $(FFLAGS) $(LDFLAGS) -o $@ $^ -xlintstrfc: $(CLINTSTRFP) $(TMGLIB) $(VARLIB) ../$(LAPACKLIB) $(BLASLIB) - $(LOADER) $(FFLAGS) $(LDFLAGS) -o $@ $^ +xlintstrfc: $(CLINTSTRFP) $(TMGLIB) $(VARLIB) $(LAPACKLIB) $(BLASLIB) + $(FC) $(FFLAGS) $(LDFLAGS) -o $@ $^ -xlintstrfz: $(ZLINTSTRFP) $(TMGLIB) $(VARLIB) ../$(LAPACKLIB) $(BLASLIB) - $(LOADER) $(FFLAGS) $(LDFLAGS) -o $@ $^ +xlintstrfz: $(ZLINTSTRFP) $(TMGLIB) $(VARLIB) $(LAPACKLIB) $(BLASLIB) + $(FC) $(FFLAGS) $(LDFLAGS) -o $@ $^ $(ALINTST): $(FRC) $(SCLNTST): $(FRC) diff --git a/lapack-netlib/TESTING/LIN/alaerh.f b/lapack-netlib/TESTING/LIN/alaerh.f index 1845888a6..6c8a47f1e 100644 --- a/lapack-netlib/TESTING/LIN/alaerh.f +++ b/lapack-netlib/TESTING/LIN/alaerh.f @@ -797,6 +797,18 @@ WRITE( NOUT, FMT = 9978 ) $ SUBNAM(1:LEN_TRIM( SUBNAM )), INFO, M, N, IMAT 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 * @@ -1147,6 +1159,11 @@ * What we do next * 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 * diff --git a/lapack-netlib/TESTING/LIN/alahd.f b/lapack-netlib/TESTING/LIN/alahd.f index dd75394b3..8f966c584 100644 --- a/lapack-netlib/TESTING/LIN/alahd.f +++ b/lapack-netlib/TESTING/LIN/alahd.f @@ -584,13 +584,27 @@ * * QR decomposition with column pivoting * - WRITE( IOUNIT, FMT = 9986 )PATH + WRITE( IOUNIT, FMT = 8006 )PATH WRITE( IOUNIT, FMT = 9969 ) WRITE( IOUNIT, FMT = '( '' Test ratios:'' )' ) WRITE( IOUNIT, FMT = 9940 )1 WRITE( IOUNIT, FMT = 9939 )2 WRITE( IOUNIT, FMT = 9938 )3 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 * @@ -779,6 +793,8 @@ $ 'tall-skinny or short-wide matrices' ) 8005 FORMAT( / 1X, A3, ': Householder reconstruction from TSQR', $ ' factorization output ', /,' for tall-skinny matrices.' ) + 8006 FORMAT( / 1X, A3, ': truncated QR factorization', + $ ' with column pivoting' ) * * GE matrix types * @@ -922,6 +938,36 @@ $ / 4X, '3. Geometric distribution', 10X, $ '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 * 9968 FORMAT( ' Matrix types (2-3 have condition 1/EPS):', / 4X, @@ -1030,9 +1076,8 @@ $ ' * norm(C) * EPS )' ) 9940 FORMAT( 3X, I2, ': norm(svd(A) - svd(R)) / ', $ '( M * norm(svd(R)) * EPS )' ) - 9939 FORMAT( 3X, I2, ': norm( A*P - Q*R ) / ( M * norm(A) * EPS )' - $ ) - 9938 FORMAT( 3X, I2, ': norm( I - Q''*Q ) / ( M * EPS )' ) + 9939 FORMAT( 3X, I2, ': norm( A*P - Q*R ) / ( M * norm(A) * EPS )') + 9938 FORMAT( 3X, I2, ': norm( I - Q''*Q ) / ( M * EPS )' ) 9937 FORMAT( 3X, I2, ': norm( A - R*Q ) / ( M * norm(A) * EPS )' $ ) 9935 FORMAT( 3X, I2, ': norm( B - A * X ) / ', @@ -1105,6 +1150,15 @@ 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 )') + 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 * diff --git a/lapack-netlib/TESTING/LIN/alareq.f b/lapack-netlib/TESTING/LIN/alareq.f index db18775eb..3f057fa48 100644 --- a/lapack-netlib/TESTING/LIN/alareq.f +++ b/lapack-netlib/TESTING/LIN/alareq.f @@ -28,12 +28,12 @@ *> to evaluate the input line which requested NMATS matrix types for *> PATH. The flow of control is as follows: *> -*> If NMATS = NTYPES then +*> IF NMATS = NTYPES THEN *> DOTYPE(1:NTYPES) = .TRUE. -*> else +*> ELSE *> Read the next input line for NMATS matrix types *> Set DOTYPE(I) = .TRUE. for each valid type I -*> endif +*> END IF *> \endverbatim * * Arguments: diff --git a/lapack-netlib/TESTING/LIN/cchkaa.F b/lapack-netlib/TESTING/LIN/cchkaa.F index ec1534ed4..474454a51 100644 --- a/lapack-netlib/TESTING/LIN/cchkaa.F +++ b/lapack-netlib/TESTING/LIN/cchkaa.F @@ -69,6 +69,7 @@ *> CLQ 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 +*> ZQK 19 List types on next line if 0 < NTYPES < 19 *> CTZ 3 List types on next line if 0 < NTYPES < 3 *> CLS 6 List types on next line if 0 < NTYPES < 6 *> CEQ @@ -153,12 +154,11 @@ $ NBVAL( MAXIN ), NBVAL2( MAXIN ), $ NSVAL( MAXIN ), NVAL( MAXIN ), NXVAL( MAXIN ), $ RANKVAL( MAXIN ), PIV( NMAX ) - REAL S( 2*NMAX ) - COMPLEX E( NMAX ) * .. * .. Allocatable Arrays .. INTEGER AllocateStatus - REAL, DIMENSION(:), ALLOCATABLE :: RWORK + REAL, DIMENSION(:), ALLOCATABLE :: RWORK, S + COMPLEX, DIMENSION(:), ALLOCATABLE :: E COMPLEX, DIMENSION(:,:), ALLOCATABLE :: A, B, WORK * .. * .. External Functions .. @@ -170,14 +170,14 @@ EXTERNAL ALAREQ, CCHKEQ, CCHKGB, CCHKGE, CCHKGT, CCHKHE, $ CCHKHE_ROOK, CCHKHE_RK, CCHKHE_AA, CCHKHP, $ CCHKLQ, CCHKUNHR_COL, CCHKPB, CCHKPO, CCHKPS, - $ CCHKPP, CCHKPT, CCHKQ3, CCHKQL, CCHKQR, CCHKRQ, - $ CCHKSP, CCHKSY, CCHKSY_ROOK, CCHKSY_RK, - $ CCHKSY_AA, CCHKTB, CCHKTP, CCHKTR, CCHKTZ, - $ CDRVGB, CDRVGE, CDRVGT, CDRVHE, CDRVHE_ROOK, - $ CDRVHE_RK, CDRVHE_AA, CDRVHP, CDRVLS, CDRVPB, - $ CDRVPO, CDRVPP, CDRVPT, CDRVSP, CDRVSY, - $ CDRVSY_ROOK, CDRVSY_RK, CDRVSY_AA, ILAVER, - $ CCHKQRT, CCHKQRTP + $ CCHKPP, CCHKPT, CCHKQ3, CCHKQP3RK, CCHKQL, + $ CCHKQR, CCHKRQ, CCHKSP, CCHKSY, CCHKSY_ROOK, + $ CCHKSY_RK, CCHKSY_AA, CCHKTB, CCHKTP, CCHKTR, + $ CCHKTZ, CDRVGB, CDRVGE, CDRVGT, CDRVHE, + $ CDRVHE_ROOK, CDRVHE_RK, CDRVHE_AA, CDRVHP, + $ CDRVLS, CDRVPB, CDRVPO, CDRVPP, CDRVPT, CDRVSP, + $ CDRVSY, CDRVSY_ROOK, CDRVSY_RK, CDRVSY_AA, + $ ILAVER, CCHKQRT, CCHKQRTP * .. * .. Scalars in Common .. LOGICAL LERR, OK @@ -203,6 +203,10 @@ IF (AllocateStatus /= 0) STOP "*** Not enough memory ***" ALLOCATE ( WORK( NMAX, NMAX+MAXRHS+10 ), STAT = AllocateStatus ) 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 ) IF (AllocateStatus /= 0) STOP "*** Not enough memory ***" * .. @@ -1109,6 +1113,23 @@ 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 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 * diff --git a/lapack-netlib/TESTING/LIN/cchkqp3rk.f b/lapack-netlib/TESTING/LIN/cchkqp3rk.f new file mode 100644 index 000000000..79d6add72 --- /dev/null +++ b/lapack-netlib/TESTING/LIN/cchkqp3rk.f @@ -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 diff --git a/lapack-netlib/TESTING/LIN/clatb4.f b/lapack-netlib/TESTING/LIN/clatb4.f index eeb0f03a9..233a8631a 100644 --- a/lapack-netlib/TESTING/LIN/clatb4.f +++ b/lapack-netlib/TESTING/LIN/clatb4.f @@ -154,9 +154,6 @@ * .. Intrinsic Functions .. INTRINSIC ABS, MAX, SQRT * .. -* .. External Subroutines .. - EXTERNAL SLABAD -* .. * .. Save statement .. SAVE EPS, SMALL, LARGE, BADC1, BADC2, FIRST * .. @@ -174,11 +171,6 @@ BADC1 = SQRT( BADC2 ) SMALL = SLAMCH( 'Safe minimum' ) 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 ) LARGE = ONE / SMALL END IF @@ -233,6 +225,110 @@ ELSE ANORM = ONE 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 * @@ -517,17 +613,18 @@ * * 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 - ELSE IF( IMAT.EQ.3 .OR. IMAT.EQ.9 ) THEN + ELSE IF( MAT.EQ.3 .OR. MAT.EQ.9 ) THEN CNDNUM = BADC2 ELSE CNDNUM = TWO END IF * - IF( IMAT.EQ.4 ) THEN + IF( MAT.EQ.4 ) THEN ANORM = SMALL - ELSE IF( IMAT.EQ.5 ) THEN + ELSE IF( MAT.EQ.5 ) THEN ANORM = LARGE ELSE ANORM = ONE diff --git a/lapack-netlib/TESTING/LIN/cqpt01.f b/lapack-netlib/TESTING/LIN/cqpt01.f index 79fc2dc66..149c5bb7c 100644 --- a/lapack-netlib/TESTING/LIN/cqpt01.f +++ b/lapack-netlib/TESTING/LIN/cqpt01.f @@ -33,7 +33,8 @@ *> Householder vectors, and the rest of AF contains a partially updated *> 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 * * Arguments: @@ -172,28 +173,28 @@ * NORMA = CLANGE( 'One-norm', M, N, A, LDA, RWORK ) * - DO 30 J = 1, K - DO 10 I = 1, MIN( J, M ) + DO J = 1, K + DO I = 1, MIN( J, M ) WORK( ( J-1 )*M+I ) = AF( I, J ) - 10 CONTINUE - DO 20 I = J + 1, M + END DO + DO I = J + 1, M WORK( ( J-1 )*M+I ) = ZERO - 20 CONTINUE - 30 CONTINUE - DO 40 J = K + 1, N + END DO + END DO + DO J = K + 1, N 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, $ 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 * CALL CAXPY( M, CMPLX( -ONE ), A( 1, JPVT( J ) ), 1, $ WORK( ( J-1 )*M+1 ), 1 ) - 50 CONTINUE + END DO * CQPT01 = CLANGE( 'One-norm', M, N, WORK, M, RWORK ) / $ ( REAL( MAX( M, N ) )*SLAMCH( 'Epsilon' ) ) diff --git a/lapack-netlib/TESTING/LIN/cqrt11.f b/lapack-netlib/TESTING/LIN/cqrt11.f index 494d5e9cd..a52084973 100644 --- a/lapack-netlib/TESTING/LIN/cqrt11.f +++ b/lapack-netlib/TESTING/LIN/cqrt11.f @@ -157,9 +157,9 @@ CALL CUNM2R( 'Left', 'Conjugate transpose', M, M, K, A, LDA, TAU, $ 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 - 10 CONTINUE + END DO * CQRT11 = CLANGE( 'One-norm', M, M, WORK, M, RDUMMY ) / $ ( REAL( M )*SLAMCH( 'Epsilon' ) ) diff --git a/lapack-netlib/TESTING/LIN/cqrt12.f b/lapack-netlib/TESTING/LIN/cqrt12.f index 4c29423ae..0df2d833b 100644 --- a/lapack-netlib/TESTING/LIN/cqrt12.f +++ b/lapack-netlib/TESTING/LIN/cqrt12.f @@ -28,7 +28,7 @@ *> CQRT12 computes the singular values `svlues' of the upper trapezoid *> 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 * * Arguments: @@ -125,8 +125,8 @@ EXTERNAL CLANGE, SASUM, SLAMCH, SNRM2 * .. * .. External Subroutines .. - EXTERNAL CGEBD2, CLASCL, CLASET, SAXPY, SBDSQR, SLABAD, - $ SLASCL, XERBLA + EXTERNAL CGEBD2, CLASCL, CLASET, SAXPY, SBDSQR, SLASCL, + $ XERBLA * .. * .. Intrinsic Functions .. INTRINSIC CMPLX, MAX, MIN, REAL @@ -153,17 +153,16 @@ * Copy upper triangle of A into work * CALL CLASET( 'Full', M, N, CMPLX( ZERO ), CMPLX( ZERO ), WORK, M ) - DO 20 J = 1, N - DO 10 I = 1, MIN( J, M ) + DO J = 1, N + DO I = 1, MIN( J, M ) WORK( ( J-1 )*M+I ) = A( I, J ) - 10 CONTINUE - 20 CONTINUE + END DO + END DO * * Get machine parameters * SMLNUM = SLAMCH( 'S' ) / SLAMCH( 'P' ) BIGNUM = ONE / SMLNUM - CALL SLABAD( SMLNUM, BIGNUM ) * * Scale work if max entry outside range [SMLNUM,BIGNUM] * @@ -207,9 +206,9 @@ * ELSE * - DO 30 I = 1, MN + DO I = 1, MN RWORK( I ) = ZERO - 30 CONTINUE + END DO END IF * * Compare s and singular values of work diff --git a/lapack-netlib/TESTING/LIN/dchkaa.F b/lapack-netlib/TESTING/LIN/dchkaa.F index ef9d7808c..74077eb94 100644 --- a/lapack-netlib/TESTING/LIN/dchkaa.F +++ b/lapack-netlib/TESTING/LIN/dchkaa.F @@ -63,6 +63,7 @@ *> DLQ 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 +*> DQK 19 List types on next line if 0 < NTYPES < 19 *> DTZ 3 List types on next line if 0 < NTYPES < 3 *> DLS 6 List types on next line if 0 < NTYPES < 6 *> DEQ @@ -149,12 +150,12 @@ $ NBVAL( MAXIN ), NBVAL2( MAXIN ), $ NSVAL( MAXIN ), NVAL( MAXIN ), NXVAL( MAXIN ), $ RANKVAL( MAXIN ), PIV( NMAX ) - DOUBLE PRECISION E( NMAX ), S( 2*NMAX ) * .. * .. Allocatable Arrays .. INTEGER AllocateStatus - DOUBLE PRECISION, DIMENSION(:), ALLOCATABLE :: RWORK - DOUBLE PRECISION, DIMENSION(:,:), ALLOCATABLE :: A, B, WORK + DOUBLE PRECISION, DIMENSION(:), ALLOCATABLE :: RWORK, S + DOUBLE PRECISION, DIMENSION(:), ALLOCATABLE :: E + DOUBLE PRECISION, DIMENSION(:,:), ALLOCATABLE :: A, B, WORK * .. * .. External Functions .. LOGICAL LSAME, LSAMEN @@ -164,13 +165,13 @@ * .. External Subroutines .. EXTERNAL ALAREQ, DCHKEQ, DCHKGB, DCHKGE, DCHKGT, DCHKLQ, $ DCHKORHR_COL, DCHKPB, DCHKPO, DCHKPS, DCHKPP, - $ DCHKPT, DCHKQ3, DCHKQL, DCHKQR, DCHKRQ, DCHKSP, - $ DCHKSY, DCHKSY_ROOK, DCHKSY_RK, DCHKSY_AA, - $ DCHKTB, DCHKTP, DCHKTR, DCHKTZ, DDRVGB, DDRVGE, - $ DDRVGT, DDRVLS, DDRVPB, DDRVPO, DDRVPP, DDRVPT, - $ DDRVSP, DDRVSY, DDRVSY_ROOK, DDRVSY_RK, - $ DDRVSY_AA, ILAVER, DCHKLQTP, DCHKQRT, DCHKQRTP, - $ DCHKLQT,DCHKTSQR + $ DCHKPT, DCHKQ3, DCHKQP3RK, DCHKQL, DCHKQR, + $ DCHKRQ, DCHKSP, DCHKSY, DCHKSY_ROOK, DCHKSY_RK, + $ DCHKSY_AA, DCHKTB, DCHKTP, DCHKTR, DCHKTZ, + $ DDRVGB, DDRVGE, DDRVGT, DDRVLS, DDRVPB, DDRVPO, + $ DDRVPP, DDRVPT, DDRVSP, DDRVSY, DDRVSY_ROOK, + $ DDRVSY_RK, DDRVSY_AA, ILAVER, DCHKLQTP, DCHKQRT, + $ DCHKQRTP, DCHKLQT,DCHKTSQR * .. * .. Scalars in Common .. LOGICAL LERR, OK @@ -197,6 +198,10 @@ IF (AllocateStatus /= 0) STOP "*** Not enough memory ***" ALLOCATE ( WORK( NMAX, 3*NMAX+MAXRHS+30 ), STAT = AllocateStatus ) 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 ) IF (AllocateStatus /= 0) STOP "*** Not enough memory ***" * @@ -919,9 +924,26 @@ CALL ALAREQ( PATH, NMATS, DOTYPE, NTYPES, NIN, NOUT ) * IF( TSTCHK ) THEN - CALL DCHKQ3( DOTYPE, NM, MVAL, NN, NVAL, NNB, NBVAL, NXVAL, - $ THRESH, A( 1, 1 ), A( 1, 2 ), B( 1, 1 ), - $ B( 1, 3 ), WORK, IWORK, NOUT ) + CALL DCHKQ3( DOTYPE, NM, MVAL, NN, NVAL, NNB, NBVAL, + $ NXVAL, THRESH, A( 1, 1 ), A( 1, 2 ), + $ 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 WRITE( NOUT, FMT = 9989 )PATH END IF diff --git a/lapack-netlib/TESTING/LIN/dchkq3.f b/lapack-netlib/TESTING/LIN/dchkq3.f index 1fdf07252..494008fa8 100644 --- a/lapack-netlib/TESTING/LIN/dchkq3.f +++ b/lapack-netlib/TESTING/LIN/dchkq3.f @@ -30,7 +30,7 @@ *> *> \verbatim *> -*> DCHKQ3 tests DGEQP3. +*> DCHKQ3 tests DGEQP3. *> \endverbatim * * Arguments: diff --git a/lapack-netlib/TESTING/LIN/dchkqp3rk.f b/lapack-netlib/TESTING/LIN/dchkqp3rk.f new file mode 100644 index 000000000..434d2067e --- /dev/null +++ b/lapack-netlib/TESTING/LIN/dchkqp3rk.f @@ -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 diff --git a/lapack-netlib/TESTING/LIN/dlatb4.f b/lapack-netlib/TESTING/LIN/dlatb4.f index 8825d13e7..f3bccd45b 100644 --- a/lapack-netlib/TESTING/LIN/dlatb4.f +++ b/lapack-netlib/TESTING/LIN/dlatb4.f @@ -133,7 +133,7 @@ * * .. Parameters .. DOUBLE PRECISION SHRINK, TENTH - PARAMETER ( SHRINK = 0.25D0, TENTH = 0.1D+0 ) + PARAMETER ( SHRINK = 0.25D+0, TENTH = 0.1D+0 ) DOUBLE PRECISION ONE PARAMETER ( ONE = 1.0D+0 ) DOUBLE PRECISION TWO @@ -153,9 +153,6 @@ * .. Intrinsic Functions .. INTRINSIC ABS, MAX, SQRT * .. -* .. External Subroutines .. - EXTERNAL DLABAD -* .. * .. Save statement .. SAVE EPS, SMALL, LARGE, BADC1, BADC2, FIRST * .. @@ -173,11 +170,6 @@ BADC1 = SQRT( BADC2 ) SMALL = DLAMCH( 'Safe minimum' ) 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 ) LARGE = ONE / SMALL END IF @@ -232,6 +224,110 @@ ELSE ANORM = ONE 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 * @@ -518,17 +614,18 @@ * * 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 - ELSE IF( IMAT.EQ.3 .OR. IMAT.EQ.9 ) THEN + ELSE IF( MAT.EQ.3 .OR. MAT.EQ.9 ) THEN CNDNUM = BADC2 ELSE CNDNUM = TWO END IF * - IF( IMAT.EQ.4 ) THEN + IF( MAT.EQ.4 ) THEN ANORM = SMALL - ELSE IF( IMAT.EQ.5 ) THEN + ELSE IF( MAT.EQ.5 ) THEN ANORM = LARGE ELSE ANORM = ONE diff --git a/lapack-netlib/TESTING/LIN/dqpt01.f b/lapack-netlib/TESTING/LIN/dqpt01.f index 8efbdc774..af3f5dd36 100644 --- a/lapack-netlib/TESTING/LIN/dqpt01.f +++ b/lapack-netlib/TESTING/LIN/dqpt01.f @@ -28,12 +28,13 @@ *> *> DQPT01 tests the QR-factorization with pivoting of a matrix A. The *> 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 entries below the diagonal in the first k columns are the +*> 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 *> Householder vectors, and the rest of AF contains a partially updated *> 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 * * Arguments: @@ -172,28 +173,41 @@ * NORMA = DLANGE( 'One-norm', M, N, A, LDA, RWORK ) * - DO 30 J = 1, K - DO 10 I = 1, MIN( J, M ) + DO J = 1, K +* +* 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 ) - 10 CONTINUE - DO 20 I = J + 1, M + END DO +* +* Zero out the elements below the diagonal in the work array. +* + DO I = J + 1, M WORK( ( J-1 )*M+I ) = ZERO - 20 CONTINUE - 30 CONTINUE - DO 40 J = K + 1, N + END DO + END DO +* +* 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 ) - 40 CONTINUE + END DO * CALL DORMQR( 'Left', 'No transpose', M, N, K, AF, LDA, TAU, WORK, $ 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 ), $ 1 ) - 50 CONTINUE + END DO * DQPT01 = DLANGE( 'One-norm', M, N, WORK, M, RWORK ) / $ ( DBLE( MAX( M, N ) )*DLAMCH( 'Epsilon' ) ) diff --git a/lapack-netlib/TESTING/LIN/dqrt11.f b/lapack-netlib/TESTING/LIN/dqrt11.f index 33c7fab37..38bbeb822 100644 --- a/lapack-netlib/TESTING/LIN/dqrt11.f +++ b/lapack-netlib/TESTING/LIN/dqrt11.f @@ -157,9 +157,9 @@ CALL DORM2R( 'Left', 'Transpose', M, M, K, A, LDA, TAU, 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 - 10 CONTINUE + END DO * DQRT11 = DLANGE( 'One-norm', M, M, WORK, M, RDUMMY ) / $ ( DBLE( M )*DLAMCH( 'Epsilon' ) ) diff --git a/lapack-netlib/TESTING/LIN/dqrt12.f b/lapack-netlib/TESTING/LIN/dqrt12.f index 278e01bf0..b8a124c59 100644 --- a/lapack-netlib/TESTING/LIN/dqrt12.f +++ b/lapack-netlib/TESTING/LIN/dqrt12.f @@ -26,7 +26,7 @@ *> DQRT12 computes the singular values `svlues' of the upper trapezoid *> 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 * * Arguments: @@ -113,8 +113,7 @@ EXTERNAL DASUM, DLAMCH, DLANGE, DNRM2 * .. * .. External Subroutines .. - EXTERNAL DAXPY, DBDSQR, DGEBD2, DLABAD, DLASCL, DLASET, - $ XERBLA + EXTERNAL DAXPY, DBDSQR, DGEBD2, DLASCL, DLASET, XERBLA * .. * .. Intrinsic Functions .. INTRINSIC DBLE, MAX, MIN @@ -145,17 +144,16 @@ * Copy upper triangle of A into work * CALL DLASET( 'Full', M, N, ZERO, ZERO, WORK, M ) - DO 20 J = 1, N - DO 10 I = 1, MIN( J, M ) + DO J = 1, N + DO I = 1, MIN( J, M ) WORK( ( J-1 )*M+I ) = A( I, J ) - 10 CONTINUE - 20 CONTINUE + END DO + END DO * * Get machine parameters * SMLNUM = DLAMCH( 'S' ) / DLAMCH( 'P' ) BIGNUM = ONE / SMLNUM - CALL DLABAD( SMLNUM, BIGNUM ) * * Scale work if max entry outside range [SMLNUM,BIGNUM] * @@ -199,16 +197,18 @@ * ELSE * - DO 30 I = 1, MN + DO I = 1, MN WORK( M*N+I ) = ZERO - 30 CONTINUE + END DO END IF * * Compare s and singular values of work * CALL DAXPY( MN, -ONE, S, 1, 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 ) $ DQRT12 = DQRT12 / NRMSVL * diff --git a/lapack-netlib/TESTING/LIN/schkaa.F b/lapack-netlib/TESTING/LIN/schkaa.F index a5b826d06..2b9f2ea45 100644 --- a/lapack-netlib/TESTING/LIN/schkaa.F +++ b/lapack-netlib/TESTING/LIN/schkaa.F @@ -63,6 +63,7 @@ *> SLQ 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 +*> DQK 19 List types on next line if 0 < NTYPES < 19 *> STZ 3 List types on next line if 0 < NTYPES < 3 *> SLS 6 List types on next line if 0 < NTYPES < 6 *> SEQ @@ -147,11 +148,11 @@ $ NBVAL( MAXIN ), NBVAL2( MAXIN ), $ NSVAL( MAXIN ), NVAL( MAXIN ), NXVAL( MAXIN ), $ RANKVAL( MAXIN ), PIV( NMAX ) - REAL E( NMAX ), S( 2*NMAX ) * .. * .. Allocatable Arrays .. INTEGER AllocateStatus - REAL, DIMENSION(:), ALLOCATABLE :: RWORK + REAL, DIMENSION(:), ALLOCATABLE :: RWORK, S + REAL, DIMENSION(:), ALLOCATABLE :: E REAL, DIMENSION(:,:), ALLOCATABLE :: A, B, WORK * .. * .. External Functions .. @@ -162,13 +163,13 @@ * .. External Subroutines .. EXTERNAL ALAREQ, SCHKEQ, SCHKGB, SCHKGE, SCHKGT, SCHKLQ, $ SCHKORHR_COL, SCHKPB, SCHKPO, SCHKPS, SCHKPP, - $ SCHKPT, SCHKQ3, SCHKQL, SCHKQR, SCHKRQ, SCHKSP, - $ SCHKSY, SCHKSY_ROOK, SCHKSY_RK, SCHKSY_AA, - $ SCHKTB, SCHKTP, SCHKTR, SCHKTZ, SDRVGB, SDRVGE, - $ SDRVGT, SDRVLS, SDRVPB, SDRVPO, SDRVPP, SDRVPT, - $ SDRVSP, SDRVSY, SDRVSY_ROOK, SDRVSY_RK, - $ SDRVSY_AA, ILAVER, SCHKLQTP, SCHKQRT, SCHKQRTP, - $ SCHKLQT, SCHKTSQR + $ SCHKPT, SCHKQ3, SCHKQP3RK, SCHKQL, SCHKQR, + $ SCHKRQ, SCHKSP, SCHKSY, SCHKSY_ROOK, SCHKSY_RK, + $ SCHKSY_AA, SCHKTB, SCHKTP, SCHKTR, SCHKTZ, + $ SDRVGB, SDRVGE, SDRVGT, SDRVLS, SDRVPB, SDRVPO, + $ SDRVPP, SDRVPT, SDRVSP, SDRVSY, SDRVSY_ROOK, + $ SDRVSY_RK, SDRVSY_AA, ILAVER, SCHKLQTP, SCHKQRT, + $ SCHKQRTP, SCHKLQT, SCHKTSQR * .. * .. Scalars in Common .. LOGICAL LERR, OK @@ -188,13 +189,17 @@ * .. * .. 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 ***" - ALLOCATE (B( NMAX*MAXRHS, 4 ), STAT = AllocateStatus ) + ALLOCATE ( B( NMAX*MAXRHS, 4 ), STAT = AllocateStatus ) 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 ***" - 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 ***" * .. * .. Executable Statements .. @@ -920,6 +925,23 @@ 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 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 * diff --git a/lapack-netlib/TESTING/LIN/schkqp3rk.f b/lapack-netlib/TESTING/LIN/schkqp3rk.f new file mode 100644 index 000000000..36cf9370e --- /dev/null +++ b/lapack-netlib/TESTING/LIN/schkqp3rk.f @@ -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 diff --git a/lapack-netlib/TESTING/LIN/slatb4.f b/lapack-netlib/TESTING/LIN/slatb4.f index 94d29db40..72a310727 100644 --- a/lapack-netlib/TESTING/LIN/slatb4.f +++ b/lapack-netlib/TESTING/LIN/slatb4.f @@ -153,9 +153,6 @@ * .. Intrinsic Functions .. INTRINSIC ABS, MAX, SQRT * .. -* .. External Subroutines .. - EXTERNAL SLABAD -* .. * .. Save statement .. SAVE EPS, SMALL, LARGE, BADC1, BADC2, FIRST * .. @@ -173,11 +170,6 @@ BADC1 = SQRT( BADC2 ) SMALL = SLAMCH( 'Safe minimum' ) 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 ) LARGE = ONE / SMALL END IF @@ -232,6 +224,110 @@ ELSE ANORM = ONE 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 * @@ -518,17 +614,18 @@ * * 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 - ELSE IF( IMAT.EQ.3 .OR. IMAT.EQ.9 ) THEN + ELSE IF( MAT.EQ.3 .OR. MAT.EQ.9 ) THEN CNDNUM = BADC2 ELSE CNDNUM = TWO END IF * - IF( IMAT.EQ.4 ) THEN + IF( MAT.EQ.4 ) THEN ANORM = SMALL - ELSE IF( IMAT.EQ.5 ) THEN + ELSE IF( MAT.EQ.5 ) THEN ANORM = LARGE ELSE ANORM = ONE diff --git a/lapack-netlib/TESTING/LIN/sqpt01.f b/lapack-netlib/TESTING/LIN/sqpt01.f index de0c80e53..f53686a65 100644 --- a/lapack-netlib/TESTING/LIN/sqpt01.f +++ b/lapack-netlib/TESTING/LIN/sqpt01.f @@ -33,7 +33,8 @@ *> Householder vectors, and the rest of AF contains a partially updated *> 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 * * Arguments: @@ -172,28 +173,28 @@ * NORMA = SLANGE( 'One-norm', M, N, A, LDA, RWORK ) * - DO 30 J = 1, K - DO 10 I = 1, MIN( J, M ) + DO J = 1, K + DO I = 1, MIN( J, M ) WORK( ( J-1 )*M+I ) = AF( I, J ) - 10 CONTINUE - DO 20 I = J + 1, M + END DO + DO I = J + 1, M WORK( ( J-1 )*M+I ) = ZERO - 20 CONTINUE - 30 CONTINUE - DO 40 J = K + 1, N + END DO + END DO + DO J = K + 1, N 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, $ 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 * CALL SAXPY( M, -ONE, A( 1, JPVT( J ) ), 1, WORK( ( J-1 )*M+1 ), $ 1 ) - 50 CONTINUE + END DO * SQPT01 = SLANGE( 'One-norm', M, N, WORK, M, RWORK ) / $ ( REAL( MAX( M, N ) )*SLAMCH( 'Epsilon' ) ) diff --git a/lapack-netlib/TESTING/LIN/sqrt11.f b/lapack-netlib/TESTING/LIN/sqrt11.f index d4422dacb..a3753adcf 100644 --- a/lapack-netlib/TESTING/LIN/sqrt11.f +++ b/lapack-netlib/TESTING/LIN/sqrt11.f @@ -157,9 +157,9 @@ CALL SORM2R( 'Left', 'Transpose', M, M, K, A, LDA, TAU, 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 - 10 CONTINUE + END DO * SQRT11 = SLANGE( 'One-norm', M, M, WORK, M, RDUMMY ) / $ ( REAL( M )*SLAMCH( 'Epsilon' ) ) diff --git a/lapack-netlib/TESTING/LIN/sqrt12.f b/lapack-netlib/TESTING/LIN/sqrt12.f index 2eab0ee0d..46b359e07 100644 --- a/lapack-netlib/TESTING/LIN/sqrt12.f +++ b/lapack-netlib/TESTING/LIN/sqrt12.f @@ -26,7 +26,7 @@ *> SQRT12 computes the singular values `svlues' of the upper trapezoid *> 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 * * Arguments: @@ -113,8 +113,7 @@ EXTERNAL SASUM, SLAMCH, SLANGE, SNRM2 * .. * .. External Subroutines .. - EXTERNAL SAXPY, SBDSQR, SGEBD2, SLABAD, SLASCL, SLASET, - $ XERBLA + EXTERNAL SAXPY, SBDSQR, SGEBD2, SLASCL, SLASET, XERBLA * .. * .. Intrinsic Functions .. INTRINSIC MAX, MIN, REAL @@ -145,17 +144,16 @@ * Copy upper triangle of A into work * CALL SLASET( 'Full', M, N, ZERO, ZERO, WORK, M ) - DO 20 J = 1, N - DO 10 I = 1, MIN( J, M ) + DO J = 1, N + DO I = 1, MIN( J, M ) WORK( ( J-1 )*M+I ) = A( I, J ) - 10 CONTINUE - 20 CONTINUE + END DO + END DO * * Get machine parameters * SMLNUM = SLAMCH( 'S' ) / SLAMCH( 'P' ) BIGNUM = ONE / SMLNUM - CALL SLABAD( SMLNUM, BIGNUM ) * * Scale work if max entry outside range [SMLNUM,BIGNUM] * @@ -199,9 +197,9 @@ * ELSE * - DO 30 I = 1, MN + DO I = 1, MN WORK( M*N+I ) = ZERO - 30 CONTINUE + END DO END IF * * Compare s and singular values of work diff --git a/lapack-netlib/TESTING/LIN/zchkaa.F b/lapack-netlib/TESTING/LIN/zchkaa.F index a118515a5..57d71833f 100644 --- a/lapack-netlib/TESTING/LIN/zchkaa.F +++ b/lapack-netlib/TESTING/LIN/zchkaa.F @@ -69,6 +69,7 @@ *> ZLQ 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 +*> ZQK 19 List types on next line if 0 < NTYPES < 19 *> ZTZ 3 List types on next line if 0 < NTYPES < 3 *> ZLS 6 List types on next line if 0 < NTYPES < 6 *> ZEQ @@ -153,12 +154,11 @@ $ NBVAL( MAXIN ), NBVAL2( MAXIN ), $ NSVAL( MAXIN ), NVAL( MAXIN ), NXVAL( MAXIN ), $ RANKVAL( MAXIN ), PIV( NMAX ) - DOUBLE PRECISION S( 2*NMAX ) - COMPLEX*16 E( NMAX ) -* -* .. Allocatable Arrays .. +* .. +* .. Allocatable Arrays .. 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 * .. * .. External Functions .. @@ -170,15 +170,16 @@ EXTERNAL ALAREQ, ZCHKEQ, ZCHKGB, ZCHKGE, ZCHKGT, ZCHKHE, $ ZCHKHE_ROOK, ZCHKHE_RK, ZCHKHE_AA, ZCHKHP, $ ZCHKLQ, ZCHKUNHR_COL, ZCHKPB, ZCHKPO, ZCHKPS, - $ ZCHKPP, ZCHKPT, ZCHKQ3, ZCHKQL, ZCHKQR, ZCHKRQ, - $ ZCHKSP, ZCHKSY, ZCHKSY_ROOK, ZCHKSY_RK, - $ ZCHKSY_AA, ZCHKTB, ZCHKTP, ZCHKTR, ZCHKTZ, - $ ZDRVGB, ZDRVGE, ZDRVGT, ZDRVHE, ZDRVHE_ROOK, - $ ZDRVHE_RK, ZDRVHE_AA, ZDRVHE_AA_2STAGE, ZDRVHP, - $ ZDRVLS, ZDRVPB, ZDRVPO, ZDRVPP, ZDRVPT, - $ ZDRVSP, ZDRVSY, ZDRVSY_ROOK, ZDRVSY_RK, - $ ZDRVSY_AA, ZDRVSY_AA_2STAGE, ILAVER, ZCHKQRT, - $ ZCHKQRTP, ZCHKLQT, ZCHKLQTP, ZCHKTSQR + $ ZCHKPP, ZCHKPT, ZCHKQ3, ZCHKQP3RK, ZCHKQL, + $ ZCHKQR, ZCHKRQ, ZCHKSP, ZCHKSY, ZCHKSY_ROOK, + $ ZCHKSY_RK, ZCHKSY_AA, ZCHKTB, ZCHKTP, ZCHKTR, + $ ZCHKTZ, ZDRVGB, ZDRVGE, ZDRVGT, ZDRVHE, + $ ZDRVHE_ROOK, ZDRVHE_RK, ZDRVHE_AA, + $ ZDRVHE_AA_2STAGE, ZDRVHP, ZDRVLS, ZDRVPB, + $ ZDRVPO, ZDRVPP, ZDRVPT, ZDRVSP, ZDRVSY, + $ ZDRVSY_ROOK, ZDRVSY_RK, ZDRVSY_AA, + $ ZDRVSY_AA_2STAGE, ILAVER, ZCHKQRT, ZCHKQRTP, + $ ZCHKLQT, ZCHKLQTP, ZCHKTSQR * .. * .. Scalars in Common .. LOGICAL LERR, OK @@ -197,13 +198,18 @@ DATA THREQ / 2.0D0 / , INTSTR / '0123456789' / * * .. 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 ***" - 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 ***" - ALLOCATE (B (NMAX * MAXRHS, 4), STAT = AllocateStatus) + ALLOCATE ( E( NMAX ), STAT = AllocateStatus ) 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 ***" * .. * .. Executable Statements .. @@ -1109,6 +1115,23 @@ 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 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 * diff --git a/lapack-netlib/TESTING/LIN/zchkqp3rk.f b/lapack-netlib/TESTING/LIN/zchkqp3rk.f new file mode 100644 index 000000000..302c7b1a8 --- /dev/null +++ b/lapack-netlib/TESTING/LIN/zchkqp3rk.f @@ -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 diff --git a/lapack-netlib/TESTING/LIN/zlatb4.f b/lapack-netlib/TESTING/LIN/zlatb4.f index a6977f4e9..a2b19f83d 100644 --- a/lapack-netlib/TESTING/LIN/zlatb4.f +++ b/lapack-netlib/TESTING/LIN/zlatb4.f @@ -154,9 +154,6 @@ * .. Intrinsic Functions .. INTRINSIC ABS, MAX, SQRT * .. -* .. External Subroutines .. - EXTERNAL DLABAD -* .. * .. Save statement .. SAVE EPS, SMALL, LARGE, BADC1, BADC2, FIRST * .. @@ -174,11 +171,6 @@ BADC1 = SQRT( BADC2 ) SMALL = DLAMCH( 'Safe minimum' ) 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 ) LARGE = ONE / SMALL END IF @@ -233,6 +225,110 @@ ELSE ANORM = ONE 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 * @@ -517,17 +613,18 @@ * * 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 - ELSE IF( IMAT.EQ.3 .OR. IMAT.EQ.9 ) THEN + ELSE IF( MAT.EQ.3 .OR. MAT.EQ.9 ) THEN CNDNUM = BADC2 ELSE CNDNUM = TWO END IF * - IF( IMAT.EQ.4 ) THEN + IF( MAT.EQ.4 ) THEN ANORM = SMALL - ELSE IF( IMAT.EQ.5 ) THEN + ELSE IF( MAT.EQ.5 ) THEN ANORM = LARGE ELSE ANORM = ONE diff --git a/lapack-netlib/TESTING/LIN/zqpt01.f b/lapack-netlib/TESTING/LIN/zqpt01.f index 4e53f92c8..c69eb658f 100644 --- a/lapack-netlib/TESTING/LIN/zqpt01.f +++ b/lapack-netlib/TESTING/LIN/zqpt01.f @@ -33,7 +33,7 @@ *> Householder vectors, and the rest of AF contains a partially updated *> 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 * * Arguments: @@ -172,28 +172,28 @@ * NORMA = ZLANGE( 'One-norm', M, N, A, LDA, RWORK ) * - DO 30 J = 1, K - DO 10 I = 1, MIN( J, M ) + DO J = 1, K + DO I = 1, MIN( J, M ) WORK( ( J-1 )*M+I ) = AF( I, J ) - 10 CONTINUE - DO 20 I = J + 1, M + END DO + DO I = J + 1, M WORK( ( J-1 )*M+I ) = ZERO - 20 CONTINUE - 30 CONTINUE - DO 40 J = K + 1, N + END DO + END DO + DO J = K + 1, N 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, $ 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 * CALL ZAXPY( M, DCMPLX( -ONE ), A( 1, JPVT( J ) ), 1, $ WORK( ( J-1 )*M+1 ), 1 ) - 50 CONTINUE + END DO * ZQPT01 = ZLANGE( 'One-norm', M, N, WORK, M, RWORK ) / $ ( DBLE( MAX( M, N ) )*DLAMCH( 'Epsilon' ) ) diff --git a/lapack-netlib/TESTING/LIN/zqrt11.f b/lapack-netlib/TESTING/LIN/zqrt11.f index c3be59c36..dc4af744f 100644 --- a/lapack-netlib/TESTING/LIN/zqrt11.f +++ b/lapack-netlib/TESTING/LIN/zqrt11.f @@ -158,9 +158,9 @@ CALL ZUNM2R( 'Left', 'Conjugate transpose', M, M, K, A, LDA, TAU, $ 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 - 10 CONTINUE + END DO * ZQRT11 = ZLANGE( 'One-norm', M, M, WORK, M, RDUMMY ) / $ ( DBLE( M )*DLAMCH( 'Epsilon' ) ) diff --git a/lapack-netlib/TESTING/LIN/zqrt12.f b/lapack-netlib/TESTING/LIN/zqrt12.f index 0da6be157..91477b5ea 100644 --- a/lapack-netlib/TESTING/LIN/zqrt12.f +++ b/lapack-netlib/TESTING/LIN/zqrt12.f @@ -28,7 +28,7 @@ *> ZQRT12 computes the singular values `svlues' of the upper trapezoid *> 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 * * Arguments: @@ -125,8 +125,8 @@ EXTERNAL DASUM, DLAMCH, DNRM2, ZLANGE * .. * .. External Subroutines .. - EXTERNAL DAXPY, DBDSQR, DLABAD, DLASCL, XERBLA, ZGEBD2, - $ ZLASCL, ZLASET + EXTERNAL DAXPY, DBDSQR, DLASCL, XERBLA, ZGEBD2, ZLASCL, + $ ZLASET * .. * .. Intrinsic Functions .. INTRINSIC DBLE, DCMPLX, MAX, MIN @@ -154,17 +154,16 @@ * CALL ZLASET( 'Full', M, N, DCMPLX( ZERO ), DCMPLX( ZERO ), WORK, $ M ) - DO 20 J = 1, N - DO 10 I = 1, MIN( J, M ) + DO J = 1, N + DO I = 1, MIN( J, M ) WORK( ( J-1 )*M+I ) = A( I, J ) - 10 CONTINUE - 20 CONTINUE + END DO + END DO * * Get machine parameters * SMLNUM = DLAMCH( 'S' ) / DLAMCH( 'P' ) BIGNUM = ONE / SMLNUM - CALL DLABAD( SMLNUM, BIGNUM ) * * Scale work if max entry outside range [SMLNUM,BIGNUM] * @@ -208,9 +207,9 @@ * ELSE * - DO 30 I = 1, MN + DO I = 1, MN RWORK( I ) = ZERO - 30 CONTINUE + END DO END IF * * Compare s and singular values of work @@ -218,6 +217,7 @@ CALL DAXPY( MN, -ONE, S, 1, RWORK( 1 ), 1 ) ZQRT12 = DASUM( MN, RWORK( 1 ), 1 ) / $ ( DLAMCH( 'Epsilon' )*DBLE( MAX( M, N ) ) ) +* IF( NRMSVL.NE.ZERO ) $ ZQRT12 = ZQRT12 / NRMSVL *