floraachy
/
OpenBLAS
1
0
Fork 0
Code Issues Pull Requests Actions Packages Projects Releases Wiki Activity

Add tests for ?GELST (Reference-LAPACK PR739)

This commit is contained in:
Martin Kroeker 2022-11-19 22:39:16 +01:00 committed by GitHub
parent d0afbd8d29
commit 1497336b20
No known key found for this signature in database
GPG Key ID: 4AEE18F83AFDEB23
9 changed files with 1267 additions and 390 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

22
lapack-netlib/TESTING/LIN/alahd.f
View File

@ -608,17 +608,18 @@
ELSE IF( LSAMEN( 2, P2, 'LS' ) ) THEN ELSE IF( LSAMEN( 2, P2, 'LS' ) ) THEN
* *
* LS: Least Squares driver routines for * LS: Least Squares driver routines for
* LS, LSD, LSS, LSX and LSY. * LS, LST, TSLS, LSD, LSS, LSX and LSY.
* *
WRITE( IOUNIT, FMT = 9984 )PATH WRITE( IOUNIT, FMT = 9984 )PATH
WRITE( IOUNIT, FMT = 9967 ) WRITE( IOUNIT, FMT = 9967 )
WRITE( IOUNIT, FMT = 9921 )C1, C1, C1, C1 WRITE( IOUNIT, FMT = 9921 )C1, C1, C1, C1, C1, C1
WRITE( IOUNIT, FMT = 9935 )1 WRITE( IOUNIT, FMT = 9935 )1
WRITE( IOUNIT, FMT = 9931 )2 WRITE( IOUNIT, FMT = 9931 )2
WRITE( IOUNIT, FMT = 9933 )3 WRITE( IOUNIT, FMT = 9919 )
WRITE( IOUNIT, FMT = 9935 )4 WRITE( IOUNIT, FMT = 9933 )7
WRITE( IOUNIT, FMT = 9934 )5 WRITE( IOUNIT, FMT = 9935 )8
WRITE( IOUNIT, FMT = 9932 )6 WRITE( IOUNIT, FMT = 9934 )9
WRITE( IOUNIT, FMT = 9932 )10
WRITE( IOUNIT, FMT = 9920 ) WRITE( IOUNIT, FMT = 9920 )
WRITE( IOUNIT, FMT = '( '' Messages:'' )' ) WRITE( IOUNIT, FMT = '( '' Messages:'' )' )
* *
@ -1048,10 +1049,11 @@
$ 'check if X is in the row space of A or A'' ', $ 'check if X is in the row space of A or A'' ',
$ '(overdetermined case)' ) $ '(overdetermined case)' )
9929 FORMAT( ' Test ratios (1-3: ', A1, 'TZRZF):' ) 9929 FORMAT( ' Test ratios (1-3: ', A1, 'TZRZF):' )
9920 FORMAT( 3X, ' 7-10: same as 3-6', 3X, ' 11-14: same as 3-6' ) 9919 FORMAT( 3X, ' 3-4: same as 1-2', 3X, ' 5-6: same as 1-2' )
9921 FORMAT( ' Test ratios:', / ' (1-2: ', A1, 'GELS, 3-6: ', A1, 9920 FORMAT( 3X, ' 11-14: same as 7-10', 3X, ' 15-18: same as 7-10' )
$ 'GELSY, 7-10: ', A1, 'GELSS, 11-14: ', A1, 'GELSD, 15-16: ', 9921 FORMAT( ' Test ratios:', / ' (1-2: ', A1, 'GELS, 3-4: ', A1,
$ A1, 'GETSLS)') $ 'GELST, 5-6: ', A1, 'GETSLS, 7-10: ', A1, 'GELSY, 11-14: ',
$ A1, 'GETSS, 15-18: ', A1, 'GELSD)' )
9928 FORMAT( 7X, 'where ALPHA = ( 1 + SQRT( 17 ) ) / 8' ) 9928 FORMAT( 7X, 'where ALPHA = ( 1 + SQRT( 17 ) ) / 8' )
9927 FORMAT( 3X, I2, ': ABS( Largest element in L )', / 12X, 9927 FORMAT( 3X, I2, ': ABS( Largest element in L )', / 12X,
$ ' - ( 1 / ( 1 - ALPHA ) ) + THRESH' ) $ ' - ( 1 / ( 1 - ALPHA ) ) + THRESH' )

388
lapack-netlib/TESTING/LIN/cdrvls.f
View File

@ -31,7 +31,8 @@
*> *>
*> \verbatim *> \verbatim
*> *>
*> CDRVLS tests the least squares driver routines CGELS, CGETSLS, CGELSS, CGELSY *> CDRVLS tests the least squares driver routines CGELS, CGELST,
*> CGETSLS, CGELSS, CGELSY
*> and CGELSD. *> and CGELSD.
*> \endverbatim *> \endverbatim
* *
@ -211,7 +212,7 @@
* *
* .. Parameters .. * .. Parameters ..
INTEGER NTESTS INTEGER NTESTS
PARAMETER ( NTESTS = 16 ) PARAMETER ( NTESTS = 18 )
INTEGER SMLSIZ INTEGER SMLSIZ
PARAMETER ( SMLSIZ = 25 ) PARAMETER ( SMLSIZ = 25 )
REAL ONE, ZERO REAL ONE, ZERO
@ -228,8 +229,8 @@
$ LWLSY, LWORK, M, MNMIN, N, NB, NCOLS, NERRS, $ LWLSY, LWORK, M, MNMIN, N, NB, NCOLS, NERRS,
$ NFAIL, NRHS, NROWS, NRUN, RANK, MB, $ NFAIL, NRHS, NROWS, NRUN, RANK, MB,
$ MMAX, NMAX, NSMAX, LIWORK, LRWORK, $ MMAX, NMAX, NSMAX, LIWORK, LRWORK,
$ LWORK_CGELS, LWORK_CGETSLS, LWORK_CGELSS, $ LWORK_CGELS, LWORK_CGELST, LWORK_CGETSLS,
$ LWORK_CGELSY, LWORK_CGELSD, $ LWORK_CGELSS, LWORK_CGELSY, LWORK_CGELSD,
$ LRWORK_CGELSY, LRWORK_CGELSS, LRWORK_CGELSD $ LRWORK_CGELSY, LRWORK_CGELSS, LRWORK_CGELSD
REAL EPS, NORMA, NORMB, RCOND REAL EPS, NORMA, NORMB, RCOND
* .. * ..
@ -249,7 +250,7 @@
* .. * ..
* .. External Subroutines .. * .. External Subroutines ..
EXTERNAL ALAERH, ALAHD, ALASVM, CERRLS, CGELS, CGELSD, EXTERNAL ALAERH, ALAHD, ALASVM, CERRLS, CGELS, CGELSD,
$ CGELSS, CGELSY, CGEMM, CGETSLS, CLACPY, $ CGELSS, CGELST, CGELSY, CGEMM, CGETSLS, CLACPY,
$ CLARNV, CQRT13, CQRT15, CQRT16, CSSCAL, $ CLARNV, CQRT13, CQRT15, CQRT16, CSSCAL,
$ SAXPY, XLAENV $ SAXPY, XLAENV
* .. * ..
@ -334,7 +335,8 @@
LIWORK = 1 LIWORK = 1
* *
* Iterate through all test cases and compute necessary workspace * Iterate through all test cases and compute necessary workspace
* sizes for ?GELS, ?GETSLS, ?GELSY, ?GELSS and ?GELSD routines. * sizes for ?GELS, ?GELST, ?GETSLS, ?GELSY, ?GELSS and ?GELSD
* routines.
* *
DO IM = 1, NM DO IM = 1, NM
M = MVAL( IM ) M = MVAL( IM )
@ -361,6 +363,10 @@
CALL CGELS( TRANS, M, N, NRHS, A, LDA, CALL CGELS( TRANS, M, N, NRHS, A, LDA,
$ B, LDB, WQ, -1, INFO ) $ B, LDB, WQ, -1, INFO )
LWORK_CGELS = INT( WQ( 1 ) ) LWORK_CGELS = INT( WQ( 1 ) )
* Compute workspace needed for CGELST
CALL CGELST( TRANS, M, N, NRHS, A, LDA,
$ B, LDB, WQ, -1, INFO )
LWORK_CGELST = INT ( WQ ( 1 ) )
* Compute workspace needed for CGETSLS * Compute workspace needed for CGETSLS
CALL CGETSLS( TRANS, M, N, NRHS, A, LDA, CALL CGETSLS( TRANS, M, N, NRHS, A, LDA,
$ B, LDB, WQ, -1, INFO ) $ B, LDB, WQ, -1, INFO )
@ -425,21 +431,26 @@
ITYPE = ( IRANK-1 )*3 + ISCALE ITYPE = ( IRANK-1 )*3 + ISCALE
IF( .NOT.DOTYPE( ITYPE ) ) IF( .NOT.DOTYPE( ITYPE ) )
$ GO TO 100 $ GO TO 100
* * =====================================================
* Begin test CGELS
* =====================================================
IF( IRANK.EQ.1 ) THEN IF( IRANK.EQ.1 ) THEN
* *
* Test CGELS
*
* Generate a matrix of scaling type ISCALE * Generate a matrix of scaling type ISCALE
* *
CALL CQRT13( ISCALE, M, N, COPYA, LDA, NORMA, CALL CQRT13( ISCALE, M, N, COPYA, LDA, NORMA,
$ ISEED ) $ ISEED )
DO 40 INB = 1, NNB *
* Loop for testing different block sizes.
*
DO INB = 1, NNB
NB = NBVAL( INB ) NB = NBVAL( INB )
CALL XLAENV( 1, NB ) CALL XLAENV( 1, NB )
CALL XLAENV( 3, NXVAL( INB ) ) CALL XLAENV( 3, NXVAL( INB ) )
* *
DO 30 ITRAN = 1, 2 * Loop for testing non-transposed and transposed.
*
DO ITRAN = 1, 2
IF( ITRAN.EQ.1 ) THEN IF( ITRAN.EQ.1 ) THEN
TRANS = 'N' TRANS = 'N'
NROWS = M NROWS = M
@ -484,15 +495,20 @@
$ ITYPE, NFAIL, NERRS, $ ITYPE, NFAIL, NERRS,
$ NOUT ) $ NOUT )
* *
* Check correctness of results * Test 1: Check correctness of results
* for CGELS, compute the residual:
* RESID = norm(B - A*X) /
* / ( max(m,n) * norm(A) * norm(X) * EPS )
* *
LDWORK = MAX( 1, NROWS )
IF( NROWS.GT.0 .AND. NRHS.GT.0 ) IF( NROWS.GT.0 .AND. NRHS.GT.0 )
$ CALL CLACPY( 'Full', NROWS, NRHS, $ CALL CLACPY( 'Full', NROWS, NRHS,
$ COPYB, LDB, C, LDB ) $ COPYB, LDB, C, LDB )
CALL CQRT16( TRANS, M, N, NRHS, COPYA, CALL CQRT16( TRANS, M, N, NRHS, COPYA,
$ LDA, B, LDB, C, LDB, RWORK, $ LDA, B, LDB, C, LDB, RWORK,
$ RESULT( 1 ) ) $ RESULT( 1 ) )
*
* Test 2: Check correctness of results
* for CGELS.
* *
IF( ( ITRAN.EQ.1 .AND. M.GE.N ) .OR. IF( ( ITRAN.EQ.1 .AND. M.GE.N ) .OR.
$ ( ITRAN.EQ.2 .AND. M.LT.N ) ) THEN $ ( ITRAN.EQ.2 .AND. M.LT.N ) ) THEN
@ -515,7 +531,7 @@
* Print information about the tests that * Print information about the tests that
* did not pass the threshold. * did not pass the threshold.
* *
DO 20 K = 1, 2 DO K = 1, 2
IF( RESULT( K ).GE.THRESH ) THEN IF( RESULT( K ).GE.THRESH ) THEN
IF( NFAIL.EQ.0 .AND. NERRS.EQ.0 ) IF( NFAIL.EQ.0 .AND. NERRS.EQ.0 )
$ CALL ALAHD( NOUT, PATH ) $ CALL ALAHD( NOUT, PATH )
@ -524,26 +540,34 @@
$ RESULT( K ) $ RESULT( K )
NFAIL = NFAIL + 1 NFAIL = NFAIL + 1
END IF END IF
20 CONTINUE END DO
NRUN = NRUN + 2 NRUN = NRUN + 2
30 CONTINUE END DO
40 CONTINUE END DO
* END IF
* * =====================================================
* Test CGETSLS * End test CGELS
* =====================================================
* =====================================================
* Begin test CGELST
* =====================================================
IF( IRANK.EQ.1 ) THEN
* *
* Generate a matrix of scaling type ISCALE * Generate a matrix of scaling type ISCALE
* *
CALL CQRT13( ISCALE, M, N, COPYA, LDA, NORMA, CALL CQRT13( ISCALE, M, N, COPYA, LDA, NORMA,
$ ISEED ) $ ISEED )
DO 65 INB = 1, NNB
MB = NBVAL( INB )
CALL XLAENV( 1, MB )
DO 62 IMB = 1, NNB
NB = NBVAL( IMB )
CALL XLAENV( 2, NB )
* *
DO 60 ITRAN = 1, 2 * Loop for testing different block sizes.
*
DO INB = 1, NNB
NB = NBVAL( INB )
CALL XLAENV( 1, NB )
CALL XLAENV( 3, NXVAL( INB ) )
*
* Loop for testing non-transposed and transposed.
*
DO ITRAN = 1, 2
IF( ITRAN.EQ.1 ) THEN IF( ITRAN.EQ.1 ) THEN
TRANS = 'N' TRANS = 'N'
NROWS = M NROWS = M
@ -560,9 +584,9 @@
IF( NCOLS.GT.0 ) THEN IF( NCOLS.GT.0 ) THEN
CALL CLARNV( 2, ISEED, NCOLS*NRHS, CALL CLARNV( 2, ISEED, NCOLS*NRHS,
$ WORK ) $ WORK )
CALL CSCAL( NCOLS*NRHS, CALL CSSCAL( NCOLS*NRHS,
$ CONE / REAL( NCOLS ), WORK, $ ONE / REAL( NCOLS ), WORK,
$ 1 ) $ 1 )
END IF END IF
CALL CGEMM( TRANS, 'No transpose', NROWS, CALL CGEMM( TRANS, 'No transpose', NROWS,
$ NRHS, NCOLS, CONE, COPYA, LDA, $ NRHS, NCOLS, CONE, COPYA, LDA,
@ -578,31 +602,37 @@
CALL CLACPY( 'Full', NROWS, NRHS, CALL CLACPY( 'Full', NROWS, NRHS,
$ COPYB, LDB, B, LDB ) $ COPYB, LDB, B, LDB )
END IF END IF
SRNAMT = 'CGETSLS ' SRNAMT = 'CGELST'
CALL CGETSLS( TRANS, M, N, NRHS, A, CALL CGELST( TRANS, M, N, NRHS, A, LDA, B,
$ LDA, B, LDB, WORK, LWORK, INFO ) $ LDB, WORK, LWORK, INFO )
*
IF( INFO.NE.0 ) IF( INFO.NE.0 )
$ CALL ALAERH( PATH, 'CGETSLS ', INFO, 0, $ CALL ALAERH( PATH, 'CGELST', INFO, 0,
$ TRANS, M, N, NRHS, -1, NB, $ TRANS, M, N, NRHS, -1, NB,
$ ITYPE, NFAIL, NERRS, $ ITYPE, NFAIL, NERRS,
$ NOUT ) $ NOUT )
* *
* Check correctness of results * Test 3: Check correctness of results
* for CGELST, compute the residual:
* RESID = norm(B - A*X) /
* / ( max(m,n) * norm(A) * norm(X) * EPS )
* *
LDWORK = MAX( 1, NROWS )
IF( NROWS.GT.0 .AND. NRHS.GT.0 ) IF( NROWS.GT.0 .AND. NRHS.GT.0 )
$ CALL CLACPY( 'Full', NROWS, NRHS, $ CALL CLACPY( 'Full', NROWS, NRHS,
$ COPYB, LDB, C, LDB ) $ COPYB, LDB, C, LDB )
CALL CQRT16( TRANS, M, N, NRHS, COPYA, CALL CQRT16( TRANS, M, N, NRHS, COPYA,
$ LDA, B, LDB, C, LDB, WORK2, $ LDA, B, LDB, C, LDB, RWORK,
$ RESULT( 15 ) ) $ RESULT( 3 ) )
*
* Test 4: Check correctness of results
* for CGELST.
* *
IF( ( ITRAN.EQ.1 .AND. M.GE.N ) .OR. IF( ( ITRAN.EQ.1 .AND. M.GE.N ) .OR.
$ ( ITRAN.EQ.2 .AND. M.LT.N ) ) THEN $ ( ITRAN.EQ.2 .AND. M.LT.N ) ) THEN
* *
* Solving LS system * Solving LS system
* *
RESULT( 16 ) = CQRT17( TRANS, 1, M, N, RESULT( 4 ) = CQRT17( TRANS, 1, M, N,
$ NRHS, COPYA, LDA, B, LDB, $ NRHS, COPYA, LDA, B, LDB,
$ COPYB, LDB, C, WORK, $ COPYB, LDB, C, WORK,
$ LWORK ) $ LWORK )
@ -610,7 +640,7 @@
* *
* Solving overdetermined system * Solving overdetermined system
* *
RESULT( 16 ) = CQRT14( TRANS, M, N, RESULT( 4 ) = CQRT14( TRANS, M, N,
$ NRHS, COPYA, LDA, B, LDB, $ NRHS, COPYA, LDA, B, LDB,
$ WORK, LWORK ) $ WORK, LWORK )
END IF END IF
@ -618,21 +648,151 @@
* Print information about the tests that * Print information about the tests that
* did not pass the threshold. * did not pass the threshold.
* *
DO 50 K = 15, 16 DO K = 3, 4
IF( RESULT( K ).GE.THRESH ) THEN IF( RESULT( K ).GE.THRESH ) THEN
IF( NFAIL.EQ.0 .AND. NERRS.EQ.0 ) IF( NFAIL.EQ.0 .AND. NERRS.EQ.0 )
$ CALL ALAHD( NOUT, PATH ) $ CALL ALAHD( NOUT, PATH )
WRITE( NOUT, FMT = 9997 )TRANS, M, WRITE( NOUT, FMT = 9999 )TRANS, M,
$ N, NRHS, MB, NB, ITYPE, K, $ N, NRHS, NB, ITYPE, K,
$ RESULT( K ) $ RESULT( K )
NFAIL = NFAIL + 1 NFAIL = NFAIL + 1
END IF END IF
50 CONTINUE END DO
NRUN = NRUN + 2 NRUN = NRUN + 2
60 CONTINUE END DO
62 CONTINUE END DO
65 CONTINUE
END IF END IF
* =====================================================
* End test CGELST
* =====================================================
* =====================================================
* Begin test CGELSTSLS
* =====================================================
IF( IRANK.EQ.1 ) THEN
*
* Generate a matrix of scaling type ISCALE
*
CALL CQRT13( ISCALE, M, N, COPYA, LDA, NORMA,
$ ISEED )
*
* Loop for testing different block sizes MB.
*
DO INB = 1, NNB
MB = NBVAL( INB )
CALL XLAENV( 1, MB )
*
* Loop for testing different block sizes NB.
*
DO IMB = 1, NNB
NB = NBVAL( IMB )
CALL XLAENV( 2, NB )
*
* Loop for testing non-transposed
* and transposed.
*
DO ITRAN = 1, 2
IF( ITRAN.EQ.1 ) THEN
TRANS = 'N'
NROWS = M
NCOLS = N
ELSE
TRANS = 'C'
NROWS = N
NCOLS = M
END IF
LDWORK = MAX( 1, NCOLS )
*
* Set up a consistent rhs
*
IF( NCOLS.GT.0 ) THEN
CALL CLARNV( 2, ISEED, NCOLS*NRHS,
$ WORK )
CALL CSCAL( NCOLS*NRHS,
$ CONE / REAL( NCOLS ),
$ WORK, 1 )
END IF
CALL CGEMM( TRANS, 'No transpose',
$ NROWS, NRHS, NCOLS, CONE,
$ COPYA, LDA, WORK, LDWORK,
$ CZERO, B, LDB )
CALL CLACPY( 'Full', NROWS, NRHS,
$ B, LDB, COPYB, LDB )
*
* Solve LS or overdetermined system
*
IF( M.GT.0 .AND. N.GT.0 ) THEN
CALL CLACPY( 'Full', M, N,
$ COPYA, LDA, A, LDA )
CALL CLACPY( 'Full', NROWS, NRHS,
$ COPYB, LDB, B, LDB )
END IF
SRNAMT = 'CGETSLS '
CALL CGETSLS( TRANS, M, N, NRHS, A,
$ LDA, B, LDB, WORK, LWORK,
$ INFO )
IF( INFO.NE.0 )
$ CALL ALAERH( PATH, 'CGETSLS ', INFO,
$ 0, TRANS, M, N, NRHS,
$ -1, NB, ITYPE, NFAIL,
$ NERRS, NOUT )
*
* Test 5: Check correctness of results
* for CGETSLS, compute the residual:
* RESID = norm(B - A*X) /
* / ( max(m,n) * norm(A) * norm(X) * EPS )
*
IF( NROWS.GT.0 .AND. NRHS.GT.0 )
$ CALL CLACPY( 'Full', NROWS, NRHS,
$ COPYB, LDB, C, LDB )
CALL CQRT16( TRANS, M, N, NRHS,
$ COPYA, LDA, B, LDB,
$ C, LDB, WORK2,
$ RESULT( 5 ) )
*
* Test 6: Check correctness of results
* for CGETSLS.
*
IF( ( ITRAN.EQ.1 .AND. M.GE.N ) .OR.
$ ( ITRAN.EQ.2 .AND. M.LT.N ) ) THEN
*
* Solving LS system, compute:
* r = norm((B- A*X)**T * A) /
* / (norm(A)*norm(B)*max(M,N,NRHS)*EPS)
*
RESULT( 6 ) = CQRT17( TRANS, 1, M,
$ N, NRHS, COPYA, LDA,
$ B, LDB, COPYB, LDB,
$ C, WORK, LWORK )
ELSE
*
* Solving overdetermined system
*
RESULT( 6 ) = CQRT14( TRANS, M, N,
$ NRHS, COPYA, LDA, B,
$ LDB, WORK, LWORK )
END IF
*
* Print information about the tests that
* did not pass the threshold.
*
DO K = 5, 6
IF( RESULT( K ).GE.THRESH ) THEN
IF( NFAIL.EQ.0 .AND. NERRS.EQ.0 )
$ CALL ALAHD( NOUT, PATH )
WRITE( NOUT, FMT = 9997 )TRANS,
$ M, N, NRHS, MB, NB, ITYPE, K,
$ RESULT( K )
NFAIL = NFAIL + 1
END IF
END DO
NRUN = NRUN + 2
END DO
END DO
END DO
END IF
* =====================================================
* End test CGELSTSLS
* ====================================================
* *
* Generate a matrix of scaling type ISCALE and rank * Generate a matrix of scaling type ISCALE and rank
* type IRANK. * type IRANK.
@ -680,37 +840,37 @@
* *
* workspace used: 2*MNMIN+NB*NB+NB*MAX(N,NRHS) * workspace used: 2*MNMIN+NB*NB+NB*MAX(N,NRHS)
* *
* Test 3: Compute relative error in svd * Test 7: Compute relative error in svd
* workspace: M*N + 4*MIN(M,N) + MAX(M,N) * workspace: M*N + 4*MIN(M,N) + MAX(M,N)
* *
RESULT( 3 ) = CQRT12( CRANK, CRANK, A, LDA, RESULT( 7 ) = CQRT12( CRANK, CRANK, A, LDA,
$ COPYS, WORK, LWORK, RWORK ) $ COPYS, WORK, LWORK, RWORK )
* *
* Test 4: Compute error in solution * Test 8: Compute error in solution
* workspace: M*NRHS + M * workspace: M*NRHS + M
* *
CALL CLACPY( 'Full', M, NRHS, COPYB, LDB, WORK, CALL CLACPY( 'Full', M, NRHS, COPYB, LDB, WORK,
$ LDWORK ) $ LDWORK )
CALL CQRT16( 'No transpose', M, N, NRHS, COPYA, CALL CQRT16( 'No transpose', M, N, NRHS, COPYA,
$ LDA, B, LDB, WORK, LDWORK, RWORK, $ LDA, B, LDB, WORK, LDWORK, RWORK,
$ RESULT( 4 ) ) $ RESULT( 8 ) )
* *
* Test 5: Check norm of r'*A * Test 9: Check norm of r'*A
* workspace: NRHS*(M+N) * workspace: NRHS*(M+N)
* *
RESULT( 5 ) = ZERO RESULT( 9 ) = ZERO
IF( M.GT.CRANK ) IF( M.GT.CRANK )
$ RESULT( 5 ) = CQRT17( 'No transpose', 1, M, $ RESULT( 9 ) = CQRT17( 'No transpose', 1, M,
$ N, NRHS, COPYA, LDA, B, LDB, $ N, NRHS, COPYA, LDA, B, LDB,
$ COPYB, LDB, C, WORK, LWORK ) $ COPYB, LDB, C, WORK, LWORK )
* *
* Test 6: Check if x is in the rowspace of A * Test 10: Check if x is in the rowspace of A
* workspace: (M+NRHS)*(N+2) * workspace: (M+NRHS)*(N+2)
* *
RESULT( 6 ) = ZERO RESULT( 10 ) = ZERO
* *
IF( N.GT.CRANK ) IF( N.GT.CRANK )
$ RESULT( 6 ) = CQRT14( 'No transpose', M, N, $ RESULT( 10 ) = CQRT14( 'No transpose', M, N,
$ NRHS, COPYA, LDA, B, LDB, $ NRHS, COPYA, LDA, B, LDB,
$ WORK, LWORK ) $ WORK, LWORK )
* *
@ -736,62 +896,6 @@
* workspace used: 3*min(m,n) + * workspace used: 3*min(m,n) +
* max(2*min(m,n),nrhs,max(m,n)) * max(2*min(m,n),nrhs,max(m,n))
* *
* Test 7: Compute relative error in svd
*
IF( RANK.GT.0 ) THEN
CALL SAXPY( MNMIN, -ONE, COPYS, 1, S, 1 )
RESULT( 7 ) = SASUM( MNMIN, S, 1 ) /
$ SASUM( MNMIN, COPYS, 1 ) /
$ ( EPS*REAL( MNMIN ) )
ELSE
RESULT( 7 ) = ZERO
END IF
*
* Test 8: Compute error in solution
*
CALL CLACPY( 'Full', M, NRHS, COPYB, LDB, WORK,
$ LDWORK )
CALL CQRT16( 'No transpose', M, N, NRHS, COPYA,
$ LDA, B, LDB, WORK, LDWORK, RWORK,
$ RESULT( 8 ) )
*
* Test 9: Check norm of r'*A
*
RESULT( 9 ) = ZERO
IF( M.GT.CRANK )
$ RESULT( 9 ) = CQRT17( 'No transpose', 1, M,
$ N, NRHS, COPYA, LDA, B, LDB,
$ COPYB, LDB, C, WORK, LWORK )
*
* Test 10: Check if x is in the rowspace of A
*
RESULT( 10 ) = ZERO
IF( N.GT.CRANK )
$ RESULT( 10 ) = CQRT14( 'No transpose', M, N,
$ NRHS, COPYA, LDA, B, LDB,
$ WORK, LWORK )
*
* Test CGELSD
*
* CGELSD: Compute the minimum-norm solution X
* to min( norm( A * X - B ) ) using a
* divide and conquer SVD.
*
CALL XLAENV( 9, 25 )
*
CALL CLACPY( 'Full', M, N, COPYA, LDA, A, LDA )
CALL CLACPY( 'Full', M, NRHS, COPYB, LDB, B,
$ LDB )
*
SRNAMT = 'CGELSD'
CALL CGELSD( M, N, NRHS, A, LDA, B, LDB, S,
$ RCOND, CRANK, WORK, LWORK, RWORK,
$ IWORK, INFO )
IF( INFO.NE.0 )
$ CALL ALAERH( PATH, 'CGELSD', INFO, 0, ' ', M,
$ N, NRHS, -1, NB, ITYPE, NFAIL,
$ NERRS, NOUT )
*
* Test 11: Compute relative error in svd * Test 11: Compute relative error in svd
* *
IF( RANK.GT.0 ) THEN IF( RANK.GT.0 ) THEN
@ -827,10 +931,66 @@
$ NRHS, COPYA, LDA, B, LDB, $ NRHS, COPYA, LDA, B, LDB,
$ WORK, LWORK ) $ WORK, LWORK )
* *
* Test CGELSD
*
* CGELSD: Compute the minimum-norm solution X
* to min( norm( A * X - B ) ) using a
* divide and conquer SVD.
*
CALL XLAENV( 9, 25 )
*
CALL CLACPY( 'Full', M, N, COPYA, LDA, A, LDA )
CALL CLACPY( 'Full', M, NRHS, COPYB, LDB, B,
$ LDB )
*
SRNAMT = 'CGELSD'
CALL CGELSD( M, N, NRHS, A, LDA, B, LDB, S,
$ RCOND, CRANK, WORK, LWORK, RWORK,
$ IWORK, INFO )
IF( INFO.NE.0 )
$ CALL ALAERH( PATH, 'CGELSD', INFO, 0, ' ', M,
$ N, NRHS, -1, NB, ITYPE, NFAIL,
$ NERRS, NOUT )
*
* Test 15: Compute relative error in svd
*
IF( RANK.GT.0 ) THEN
CALL SAXPY( MNMIN, -ONE, COPYS, 1, S, 1 )
RESULT( 15 ) = SASUM( MNMIN, S, 1 ) /
$ SASUM( MNMIN, COPYS, 1 ) /
$ ( EPS*REAL( MNMIN ) )
ELSE
RESULT( 15 ) = ZERO
END IF
*
* Test 16: Compute error in solution
*
CALL CLACPY( 'Full', M, NRHS, COPYB, LDB, WORK,
$ LDWORK )
CALL CQRT16( 'No transpose', M, N, NRHS, COPYA,
$ LDA, B, LDB, WORK, LDWORK, RWORK,
$ RESULT( 16 ) )
*
* Test 17: Check norm of r'*A
*
RESULT( 17 ) = ZERO
IF( M.GT.CRANK )
$ RESULT( 17 ) = CQRT17( 'No transpose', 1, M,
$ N, NRHS, COPYA, LDA, B, LDB,
$ COPYB, LDB, C, WORK, LWORK )
*
* Test 18: Check if x is in the rowspace of A
*
RESULT( 18 ) = ZERO
IF( N.GT.CRANK )
$ RESULT( 18 ) = CQRT14( 'No transpose', M, N,
$ NRHS, COPYA, LDA, B, LDB,
$ WORK, LWORK )
*
* Print information about the tests that did not * Print information about the tests that did not
* pass the threshold. * pass the threshold.
* *
DO 80 K = 3, 14 DO 80 K = 7, 18
IF( RESULT( K ).GE.THRESH ) THEN IF( RESULT( K ).GE.THRESH ) THEN
IF( NFAIL.EQ.0 .AND. NERRS.EQ.0 ) IF( NFAIL.EQ.0 .AND. NERRS.EQ.0 )
$ CALL ALAHD( NOUT, PATH ) $ CALL ALAHD( NOUT, PATH )

61
lapack-netlib/TESTING/LIN/cerrls.f
View File

@ -22,7 +22,7 @@
*> \verbatim *> \verbatim
*> *>
*> CERRLS tests the error exits for the COMPLEX least squares *> CERRLS tests the error exits for the COMPLEX least squares
*> driver routines (CGELS, CGELSS, CGELSY, CGELSD). *> driver routines (CGELS, CGELST, CGETSLS, CGELSS, CGELSY, CGELSD).
*> \endverbatim *> \endverbatim
* *
* Arguments: * Arguments:
@ -83,7 +83,8 @@
EXTERNAL LSAMEN EXTERNAL LSAMEN
* .. * ..
* .. External Subroutines .. * .. External Subroutines ..
EXTERNAL ALAESM, CGELS, CGELSD, CGELSS, CGELSY, CHKXER EXTERNAL ALAESM, CHKXER, CGELS, CGELSD, CGELSS, CGELST,
$ CGELSY, CGETSLS
* .. * ..
* .. Scalars in Common .. * .. Scalars in Common ..
LOGICAL LERR, OK LOGICAL LERR, OK
@ -130,10 +131,66 @@
INFOT = 8 INFOT = 8
CALL CGELS( 'N', 2, 0, 0, A, 2, B, 1, W, 2, INFO ) CALL CGELS( 'N', 2, 0, 0, A, 2, B, 1, W, 2, INFO )
CALL CHKXER( 'CGELS ', INFOT, NOUT, LERR, OK ) CALL CHKXER( 'CGELS ', INFOT, NOUT, LERR, OK )
INFOT = 8
CALL CGELS( 'N', 0, 2, 0, A, 1, B, 1, W, 2, INFO )
CALL CHKXER( 'CGELS', INFOT, NOUT, LERR, OK )
INFOT = 10 INFOT = 10
CALL CGELS( 'N', 1, 1, 0, A, 1, B, 1, W, 1, INFO ) CALL CGELS( 'N', 1, 1, 0, A, 1, B, 1, W, 1, INFO )
CALL CHKXER( 'CGELS ', INFOT, NOUT, LERR, OK ) CALL CHKXER( 'CGELS ', INFOT, NOUT, LERR, OK )
* *
* CGELST
*
SRNAMT = 'CGELST'
INFOT = 1
CALL CGELST( '/', 0, 0, 0, A, 1, B, 1, W, 1, INFO )
CALL CHKXER( 'CGELST', INFOT, NOUT, LERR, OK )
INFOT = 2
CALL CGELST( 'N', -1, 0, 0, A, 1, B, 1, W, 1, INFO )
CALL CHKXER( 'CGELST', INFOT, NOUT, LERR, OK )
INFOT = 3
CALL CGELST( 'N', 0, -1, 0, A, 1, B, 1, W, 1, INFO )
CALL CHKXER( 'CGELST', INFOT, NOUT, LERR, OK )
INFOT = 4
CALL CGELST( 'N', 0, 0, -1, A, 1, B, 1, W, 1, INFO )
CALL CHKXER( 'CGELST', INFOT, NOUT, LERR, OK )
INFOT = 6
CALL CGELST( 'N', 2, 0, 0, A, 1, B, 2, W, 2, INFO )
CALL CHKXER( 'CGELST', INFOT, NOUT, LERR, OK )
INFOT = 8
CALL CGELST( 'N', 2, 0, 0, A, 2, B, 1, W, 2, INFO )
CALL CHKXER( 'CGELST', INFOT, NOUT, LERR, OK )
INFOT = 8
CALL CGELST( 'N', 0, 2, 0, A, 1, B, 1, W, 2, INFO )
CALL CHKXER( 'CGELST', INFOT, NOUT, LERR, OK )
INFOT = 10
CALL CGELST( 'N', 1, 1, 0, A, 1, B, 1, W, 1, INFO )
CALL CHKXER( 'CGELST', INFOT, NOUT, LERR, OK )
*
* CGETSLS
*
SRNAMT = 'CGETSLS'
INFOT = 1
CALL CGETSLS( '/', 0, 0, 0, A, 1, B, 1, W, 1, INFO )
CALL CHKXER( 'CGETSLS', INFOT, NOUT, LERR, OK )
INFOT = 2
CALL CGETSLS( 'N', -1, 0, 0, A, 1, B, 1, W, 1, INFO )
CALL CHKXER( 'CGETSLS', INFOT, NOUT, LERR, OK )
INFOT = 3
CALL CGETSLS( 'N', 0, -1, 0, A, 1, B, 1, W, 1, INFO )
CALL CHKXER( 'CGETSLS', INFOT, NOUT, LERR, OK )
INFOT = 4
CALL CGETSLS( 'N', 0, 0, -1, A, 1, B, 1, W, 1, INFO )
CALL CHKXER( 'CGETSLS', INFOT, NOUT, LERR, OK )
INFOT = 6
CALL CGETSLS( 'N', 2, 0, 0, A, 1, B, 2, W, 2, INFO )
CALL CHKXER( 'CGETSLS', INFOT, NOUT, LERR, OK )
INFOT = 8
CALL CGETSLS( 'N', 2, 0, 0, A, 2, B, 1, W, 2, INFO )
CALL CHKXER( 'CGETSLS', INFOT, NOUT, LERR, OK )
INFOT = 8
CALL CGETSLS( 'N', 0, 2, 0, A, 1, B, 1, W, 2, INFO )
CALL CHKXER( 'CGETSLS', INFOT, NOUT, LERR, OK )
*
* CGELSS * CGELSS
* *
SRNAMT = 'CGELSS' SRNAMT = 'CGELSS'

339
lapack-netlib/TESTING/LIN/ddrvls.f
View File

@ -31,8 +31,8 @@
*> *>
*> \verbatim *> \verbatim
*> *>
*> DDRVLS tests the least squares driver routines DGELS, DGETSLS, DGELSS, DGELSY, *> DDRVLS tests the least squares driver routines DGELS, DGELST,
*> and DGELSD. *> DGETSLS, DGELSS, DGELSY, and DGELSD.
*> \endverbatim *> \endverbatim
* *
* Arguments: * Arguments:
@ -211,7 +211,7 @@
* *
* .. Parameters .. * .. Parameters ..
INTEGER NTESTS INTEGER NTESTS
PARAMETER ( NTESTS = 16 ) PARAMETER ( NTESTS = 18 )
INTEGER SMLSIZ INTEGER SMLSIZ
PARAMETER ( SMLSIZ = 25 ) PARAMETER ( SMLSIZ = 25 )
DOUBLE PRECISION ONE, TWO, ZERO DOUBLE PRECISION ONE, TWO, ZERO
@ -225,8 +225,8 @@
$ LWLSY, LWORK, M, MNMIN, N, NB, NCOLS, NERRS, $ LWLSY, LWORK, M, MNMIN, N, NB, NCOLS, NERRS,
$ NFAIL, NRHS, NROWS, NRUN, RANK, MB, $ NFAIL, NRHS, NROWS, NRUN, RANK, MB,
$ MMAX, NMAX, NSMAX, LIWORK, $ MMAX, NMAX, NSMAX, LIWORK,
$ LWORK_DGELS, LWORK_DGETSLS, LWORK_DGELSS, $ LWORK_DGELS, LWORK_DGELST, LWORK_DGETSLS,
$ LWORK_DGELSY, LWORK_DGELSD $ LWORK_DGELSS, LWORK_DGELSY, LWORK_DGELSD
DOUBLE PRECISION EPS, NORMA, NORMB, RCOND DOUBLE PRECISION EPS, NORMA, NORMB, RCOND
* .. * ..
* .. Local Arrays .. * .. Local Arrays ..
@ -243,12 +243,12 @@
* .. * ..
* .. External Subroutines .. * .. External Subroutines ..
EXTERNAL ALAERH, ALAHD, ALASVM, DAXPY, DERRLS, DGELS, EXTERNAL ALAERH, ALAHD, ALASVM, DAXPY, DERRLS, DGELS,
$ DGELSD, DGELSS, DGELSY, DGEMM, DLACPY, $ DGELSD, DGELSS, DGELST, DGELSY, DGEMM,
$ DLARNV, DLASRT, DQRT13, DQRT15, DQRT16, DSCAL, $ DGETSLS, DLACPY, DLARNV, DQRT13, DQRT15,
$ XLAENV $ DQRT16, DSCAL, XLAENV
* .. * ..
* .. Intrinsic Functions .. * .. Intrinsic Functions ..
INTRINSIC DBLE, INT, LOG, MAX, MIN, SQRT INTRINSIC DBLE, INT, MAX, MIN, SQRT
* .. * ..
* .. Scalars in Common .. * .. Scalars in Common ..
LOGICAL LERR, OK LOGICAL LERR, OK
@ -330,7 +330,8 @@
LIWORK = 1 LIWORK = 1
* *
* Iterate through all test cases and compute necessary workspace * Iterate through all test cases and compute necessary workspace
* sizes for ?GELS, ?GETSLS, ?GELSY, ?GELSS and ?GELSD routines. * sizes for ?GELS, ?GELST, ?GETSLS, ?GELSY, ?GELSS and ?GELSD
* routines.
* *
DO IM = 1, NM DO IM = 1, NM
M = MVAL( IM ) M = MVAL( IM )
@ -357,6 +358,10 @@
CALL DGELS( TRANS, M, N, NRHS, A, LDA, CALL DGELS( TRANS, M, N, NRHS, A, LDA,
$ B, LDB, WQ, -1, INFO ) $ B, LDB, WQ, -1, INFO )
LWORK_DGELS = INT ( WQ ( 1 ) ) LWORK_DGELS = INT ( WQ ( 1 ) )
* Compute workspace needed for DGELST
CALL DGELST( TRANS, M, N, NRHS, A, LDA,
$ B, LDB, WQ, -1, INFO )
LWORK_DGELST = INT ( WQ ( 1 ) )
* Compute workspace needed for DGETSLS * Compute workspace needed for DGETSLS
CALL DGETSLS( TRANS, M, N, NRHS, A, LDA, CALL DGETSLS( TRANS, M, N, NRHS, A, LDA,
$ B, LDB, WQ, -1, INFO ) $ B, LDB, WQ, -1, INFO )
@ -378,9 +383,9 @@
* Compute LIWORK workspace needed for DGELSY and DGELSD * Compute LIWORK workspace needed for DGELSY and DGELSD
LIWORK = MAX( LIWORK, N, IWQ( 1 ) ) LIWORK = MAX( LIWORK, N, IWQ( 1 ) )
* Compute LWORK workspace needed for all functions * Compute LWORK workspace needed for all functions
LWORK = MAX( LWORK, LWORK_DGELS, LWORK_DGETSLS, LWORK = MAX( LWORK, LWORK_DGELS, LWORK_DGELST,
$ LWORK_DGELSY, LWORK_DGELSS, $ LWORK_DGETSLS, LWORK_DGELSY,
$ LWORK_DGELSD ) $ LWORK_DGELSS, LWORK_DGELSD )
END IF END IF
ENDDO ENDDO
ENDDO ENDDO
@ -411,21 +416,26 @@
ITYPE = ( IRANK-1 )*3 + ISCALE ITYPE = ( IRANK-1 )*3 + ISCALE
IF( .NOT.DOTYPE( ITYPE ) ) IF( .NOT.DOTYPE( ITYPE ) )
$ GO TO 110 $ GO TO 110
* * =====================================================
* Begin test DGELS
* =====================================================
IF( IRANK.EQ.1 ) THEN IF( IRANK.EQ.1 ) THEN
* *
* Test DGELS
*
* Generate a matrix of scaling type ISCALE * Generate a matrix of scaling type ISCALE
* *
CALL DQRT13( ISCALE, M, N, COPYA, LDA, NORMA, CALL DQRT13( ISCALE, M, N, COPYA, LDA, NORMA,
$ ISEED ) $ ISEED )
DO 40 INB = 1, NNB *
* Loop for testing different block sizes.
*
DO INB = 1, NNB
NB = NBVAL( INB ) NB = NBVAL( INB )
CALL XLAENV( 1, NB ) CALL XLAENV( 1, NB )
CALL XLAENV( 3, NXVAL( INB ) ) CALL XLAENV( 3, NXVAL( INB ) )
* *
DO 30 ITRAN = 1, 2 * Loop for testing non-transposed and transposed.
*
DO ITRAN = 1, 2
IF( ITRAN.EQ.1 ) THEN IF( ITRAN.EQ.1 ) THEN
TRANS = 'N' TRANS = 'N'
NROWS = M NROWS = M
@ -469,20 +479,27 @@
$ ITYPE, NFAIL, NERRS, $ ITYPE, NFAIL, NERRS,
$ NOUT ) $ NOUT )
* *
* Check correctness of results * Test 1: Check correctness of results
* for DGELS, compute the residual:
* RESID = norm(B - A*X) /
* / ( max(m,n) * norm(A) * norm(X) * EPS )
* *
LDWORK = MAX( 1, NROWS )
IF( NROWS.GT.0 .AND. NRHS.GT.0 ) IF( NROWS.GT.0 .AND. NRHS.GT.0 )
$ CALL DLACPY( 'Full', NROWS, NRHS, $ CALL DLACPY( 'Full', NROWS, NRHS,
$ COPYB, LDB, C, LDB ) $ COPYB, LDB, C, LDB )
CALL DQRT16( TRANS, M, N, NRHS, COPYA, CALL DQRT16( TRANS, M, N, NRHS, COPYA,
$ LDA, B, LDB, C, LDB, WORK, $ LDA, B, LDB, C, LDB, WORK,
$ RESULT( 1 ) ) $ RESULT( 1 ) )
*
* Test 2: Check correctness of results
* for DGELS.
* *
IF( ( ITRAN.EQ.1 .AND. M.GE.N ) .OR. IF( ( ITRAN.EQ.1 .AND. M.GE.N ) .OR.
$ ( ITRAN.EQ.2 .AND. M.LT.N ) ) THEN $ ( ITRAN.EQ.2 .AND. M.LT.N ) ) THEN
* *
* Solving LS system * Solving LS system, compute:
* r = norm((B- A*X)**T * A) /
* / (norm(A)*norm(B)*max(M,N,NRHS)*EPS)
* *
RESULT( 2 ) = DQRT17( TRANS, 1, M, N, RESULT( 2 ) = DQRT17( TRANS, 1, M, N,
$ NRHS, COPYA, LDA, B, LDB, $ NRHS, COPYA, LDA, B, LDB,
@ -500,35 +517,42 @@
* Print information about the tests that * Print information about the tests that
* did not pass the threshold. * did not pass the threshold.
* *
DO 20 K = 1, 2 DO K = 1, 2
IF( RESULT( K ).GE.THRESH ) THEN IF( RESULT( K ).GE.THRESH ) THEN
IF( NFAIL.EQ.0 .AND. NERRS.EQ.0 ) IF( NFAIL.EQ.0 .AND. NERRS.EQ.0 )
$ CALL ALAHD( NOUT, PATH ) $ CALL ALAHD( NOUT, PATH )
WRITE( NOUT, FMT = 9999 )TRANS, M, WRITE( NOUT, FMT = 9999 ) TRANS, M,
$ N, NRHS, NB, ITYPE, K, $ N, NRHS, NB, ITYPE, K,
$ RESULT( K ) $ RESULT( K )
NFAIL = NFAIL + 1 NFAIL = NFAIL + 1
END IF END IF
20 CONTINUE END DO
NRUN = NRUN + 2 NRUN = NRUN + 2
30 CONTINUE END DO
40 CONTINUE END DO
* END IF
* * =====================================================
* Test DGETSLS * End test DGELS
* =====================================================
* =====================================================
* Begin test DGELST
* =====================================================
IF( IRANK.EQ.1 ) THEN
* *
* Generate a matrix of scaling type ISCALE * Generate a matrix of scaling type ISCALE
* *
CALL DQRT13( ISCALE, M, N, COPYA, LDA, NORMA, CALL DQRT13( ISCALE, M, N, COPYA, LDA, NORMA,
$ ISEED ) $ ISEED )
DO 65 INB = 1, NNB
MB = NBVAL( INB )
CALL XLAENV( 1, MB )
DO 62 IMB = 1, NNB
NB = NBVAL( IMB )
CALL XLAENV( 2, NB )
* *
DO 60 ITRAN = 1, 2 * Loop for testing different block sizes.
*
DO INB = 1, NNB
NB = NBVAL( INB )
CALL XLAENV( 1, NB )
*
* Loop for testing non-transposed and transposed.
*
DO ITRAN = 1, 2
IF( ITRAN.EQ.1 ) THEN IF( ITRAN.EQ.1 ) THEN
TRANS = 'N' TRANS = 'N'
NROWS = M NROWS = M
@ -563,31 +587,38 @@
CALL DLACPY( 'Full', NROWS, NRHS, CALL DLACPY( 'Full', NROWS, NRHS,
$ COPYB, LDB, B, LDB ) $ COPYB, LDB, B, LDB )
END IF END IF
SRNAMT = 'DGETSLS ' SRNAMT = 'DGELST'
CALL DGETSLS( TRANS, M, N, NRHS, A, CALL DGELST( TRANS, M, N, NRHS, A, LDA, B,
$ LDA, B, LDB, WORK, LWORK, INFO ) $ LDB, WORK, LWORK, INFO )
IF( INFO.NE.0 ) IF( INFO.NE.0 )
$ CALL ALAERH( PATH, 'DGETSLS ', INFO, 0, $ CALL ALAERH( PATH, 'DGELST', INFO, 0,
$ TRANS, M, N, NRHS, -1, NB, $ TRANS, M, N, NRHS, -1, NB,
$ ITYPE, NFAIL, NERRS, $ ITYPE, NFAIL, NERRS,
$ NOUT ) $ NOUT )
* *
* Check correctness of results * Test 3: Check correctness of results
* for DGELST, compute the residual:
* RESID = norm(B - A*X) /
* / ( max(m,n) * norm(A) * norm(X) * EPS )
* *
LDWORK = MAX( 1, NROWS )
IF( NROWS.GT.0 .AND. NRHS.GT.0 ) IF( NROWS.GT.0 .AND. NRHS.GT.0 )
$ CALL DLACPY( 'Full', NROWS, NRHS, $ CALL DLACPY( 'Full', NROWS, NRHS,
$ COPYB, LDB, C, LDB ) $ COPYB, LDB, C, LDB )
CALL DQRT16( TRANS, M, N, NRHS, COPYA, CALL DQRT16( TRANS, M, N, NRHS, COPYA,
$ LDA, B, LDB, C, LDB, WORK, $ LDA, B, LDB, C, LDB, WORK,
$ RESULT( 15 ) ) $ RESULT( 3 ) )
*
* Test 4: Check correctness of results
* for DGELST.
* *
IF( ( ITRAN.EQ.1 .AND. M.GE.N ) .OR. IF( ( ITRAN.EQ.1 .AND. M.GE.N ) .OR.
$ ( ITRAN.EQ.2 .AND. M.LT.N ) ) THEN $ ( ITRAN.EQ.2 .AND. M.LT.N ) ) THEN
* *
* Solving LS system * Solving LS system, compute:
* r = norm((B- A*X)**T * A) /
* / (norm(A)*norm(B)*max(M,N,NRHS)*EPS)
* *
RESULT( 16 ) = DQRT17( TRANS, 1, M, N, RESULT( 4 ) = DQRT17( TRANS, 1, M, N,
$ NRHS, COPYA, LDA, B, LDB, $ NRHS, COPYA, LDA, B, LDB,
$ COPYB, LDB, C, WORK, $ COPYB, LDB, C, WORK,
$ LWORK ) $ LWORK )
@ -595,7 +626,7 @@
* *
* Solving overdetermined system * Solving overdetermined system
* *
RESULT( 16 ) = DQRT14( TRANS, M, N, RESULT( 4 ) = DQRT14( TRANS, M, N,
$ NRHS, COPYA, LDA, B, LDB, $ NRHS, COPYA, LDA, B, LDB,
$ WORK, LWORK ) $ WORK, LWORK )
END IF END IF
@ -603,21 +634,151 @@
* Print information about the tests that * Print information about the tests that
* did not pass the threshold. * did not pass the threshold.
* *
DO 50 K = 15, 16 DO K = 3, 4
IF( RESULT( K ).GE.THRESH ) THEN IF( RESULT( K ).GE.THRESH ) THEN
IF( NFAIL.EQ.0 .AND. NERRS.EQ.0 ) IF( NFAIL.EQ.0 .AND. NERRS.EQ.0 )
$ CALL ALAHD( NOUT, PATH ) $ CALL ALAHD( NOUT, PATH )
WRITE( NOUT, FMT = 9997 )TRANS, M, WRITE( NOUT, FMT = 9999 ) TRANS, M,
$ N, NRHS, MB, NB, ITYPE, K, $ N, NRHS, NB, ITYPE, K,
$ RESULT( K ) $ RESULT( K )
NFAIL = NFAIL + 1 NFAIL = NFAIL + 1
END IF END IF
50 CONTINUE END DO
NRUN = NRUN + 2 NRUN = NRUN + 2
60 CONTINUE END DO
62 CONTINUE END DO
65 CONTINUE
END IF END IF
* =====================================================
* End test DGELST
* =====================================================
* =====================================================
* Begin test DGETSLS
* =====================================================
IF( IRANK.EQ.1 ) THEN
*
* Generate a matrix of scaling type ISCALE
*
CALL DQRT13( ISCALE, M, N, COPYA, LDA, NORMA,
$ ISEED )
*
* Loop for testing different block sizes MB.
*
DO IMB = 1, NNB
MB = NBVAL( IMB )
CALL XLAENV( 1, MB )
*
* Loop for testing different block sizes NB.
*
DO INB = 1, NNB
NB = NBVAL( INB )
CALL XLAENV( 2, NB )
*
* Loop for testing non-transposed
* and transposed.
*
DO ITRAN = 1, 2
IF( ITRAN.EQ.1 ) THEN
TRANS = 'N'
NROWS = M
NCOLS = N
ELSE
TRANS = 'T'
NROWS = N
NCOLS = M
END IF
LDWORK = MAX( 1, NCOLS )
*
* Set up a consistent rhs
*
IF( NCOLS.GT.0 ) THEN
CALL DLARNV( 2, ISEED, NCOLS*NRHS,
$ WORK )
CALL DSCAL( NCOLS*NRHS,
$ ONE / DBLE( NCOLS ),
$ WORK, 1 )
END IF
CALL DGEMM( TRANS, 'No transpose',
$ NROWS, NRHS, NCOLS, ONE,
$ COPYA, LDA, WORK, LDWORK,
$ ZERO, B, LDB )
CALL DLACPY( 'Full', NROWS, NRHS,
$ B, LDB, COPYB, LDB )
*
* Solve LS or overdetermined system
*
IF( M.GT.0 .AND. N.GT.0 ) THEN
CALL DLACPY( 'Full', M, N,
$ COPYA, LDA, A, LDA )
CALL DLACPY( 'Full', NROWS, NRHS,
$ COPYB, LDB, B, LDB )
END IF
SRNAMT = 'DGETSLS'
CALL DGETSLS( TRANS, M, N, NRHS,
$ A, LDA, B, LDB, WORK, LWORK,
$ INFO )
IF( INFO.NE.0 )
$ CALL ALAERH( PATH, 'DGETSLS', INFO,
$ 0, TRANS, M, N, NRHS,
$ -1, NB, ITYPE, NFAIL,
$ NERRS, NOUT )
*
* Test 5: Check correctness of results
* for DGETSLS, compute the residual:
* RESID = norm(B - A*X) /
* / ( max(m,n) * norm(A) * norm(X) * EPS )
*
IF( NROWS.GT.0 .AND. NRHS.GT.0 )
$ CALL DLACPY( 'Full', NROWS, NRHS,
$ COPYB, LDB, C, LDB )
CALL DQRT16( TRANS, M, N, NRHS,
$ COPYA, LDA, B, LDB,
$ C, LDB, WORK,
$ RESULT( 5 ) )
*
* Test 6: Check correctness of results
* for DGETSLS.
*
IF( ( ITRAN.EQ.1 .AND. M.GE.N ) .OR.
$ ( ITRAN.EQ.2 .AND. M.LT.N ) ) THEN
*
* Solving LS system, compute:
* r = norm((B- A*X)**T * A) /
* / (norm(A)*norm(B)*max(M,N,NRHS)*EPS)
*
RESULT( 6 ) = DQRT17( TRANS, 1, M,
$ N, NRHS, COPYA, LDA,
$ B, LDB, COPYB, LDB,
$ C, WORK, LWORK )
ELSE
*
* Solving overdetermined system
*
RESULT( 6 ) = DQRT14( TRANS, M, N,
$ NRHS, COPYA, LDA,
$ B, LDB, WORK, LWORK )
END IF
*
* Print information about the tests that
* did not pass the threshold.
*
DO K = 5, 6
IF( RESULT( K ).GE.THRESH ) THEN
IF( NFAIL.EQ.0 .AND. NERRS.EQ.0 )
$ CALL ALAHD( NOUT, PATH )
WRITE( NOUT, FMT = 9997 ) TRANS,
$ M, N, NRHS, MB, NB, ITYPE,
$ K, RESULT( K )
NFAIL = NFAIL + 1
END IF
END DO
NRUN = NRUN + 2
END DO
END DO
END DO
END IF
* =====================================================
* End test DGETSLS
* =====================================================
* *
* Generate a matrix of scaling type ISCALE and rank * Generate a matrix of scaling type ISCALE and rank
* type IRANK. * type IRANK.
@ -662,37 +823,37 @@
$ N, NRHS, -1, NB, ITYPE, NFAIL, $ N, NRHS, -1, NB, ITYPE, NFAIL,
$ NERRS, NOUT ) $ NERRS, NOUT )
* *
* Test 3: Compute relative error in svd * Test 7: Compute relative error in svd
* workspace: M*N + 4*MIN(M,N) + MAX(M,N) * workspace: M*N + 4*MIN(M,N) + MAX(M,N)
* *
RESULT( 3 ) = DQRT12( CRANK, CRANK, A, LDA, RESULT( 7 ) = DQRT12( CRANK, CRANK, A, LDA,
$ COPYS, WORK, LWORK ) $ COPYS, WORK, LWORK )
* *
* Test 4: Compute error in solution * Test 8: Compute error in solution
* workspace: M*NRHS + M * workspace: M*NRHS + M
* *
CALL DLACPY( 'Full', M, NRHS, COPYB, LDB, WORK, CALL DLACPY( 'Full', M, NRHS, COPYB, LDB, WORK,
$ LDWORK ) $ LDWORK )
CALL DQRT16( 'No transpose', M, N, NRHS, COPYA, CALL DQRT16( 'No transpose', M, N, NRHS, COPYA,
$ LDA, B, LDB, WORK, LDWORK, $ LDA, B, LDB, WORK, LDWORK,
$ WORK( M*NRHS+1 ), RESULT( 4 ) ) $ WORK( M*NRHS+1 ), RESULT( 8 ) )
* *
* Test 5: Check norm of r'*A * Test 9: Check norm of r'*A
* workspace: NRHS*(M+N) * workspace: NRHS*(M+N)
* *
RESULT( 5 ) = ZERO RESULT( 9 ) = ZERO
IF( M.GT.CRANK ) IF( M.GT.CRANK )
$ RESULT( 5 ) = DQRT17( 'No transpose', 1, M, $ RESULT( 9 ) = DQRT17( 'No transpose', 1, M,
$ N, NRHS, COPYA, LDA, B, LDB, $ N, NRHS, COPYA, LDA, B, LDB,
$ COPYB, LDB, C, WORK, LWORK ) $ COPYB, LDB, C, WORK, LWORK )
* *
* Test 6: Check if x is in the rowspace of A * Test 10: Check if x is in the rowspace of A
* workspace: (M+NRHS)*(N+2) * workspace: (M+NRHS)*(N+2)
* *
RESULT( 6 ) = ZERO RESULT( 10 ) = ZERO
* *
IF( N.GT.CRANK ) IF( N.GT.CRANK )
$ RESULT( 6 ) = DQRT14( 'No transpose', M, N, $ RESULT( 10 ) = DQRT14( 'No transpose', M, N,
$ NRHS, COPYA, LDA, B, LDB, $ NRHS, COPYA, LDA, B, LDB,
$ WORK, LWORK ) $ WORK, LWORK )
* *
@ -716,38 +877,38 @@
* workspace used: 3*min(m,n) + * workspace used: 3*min(m,n) +
* max(2*min(m,n),nrhs,max(m,n)) * max(2*min(m,n),nrhs,max(m,n))
* *
* Test 7: Compute relative error in svd * Test 11: Compute relative error in svd
* *
IF( RANK.GT.0 ) THEN IF( RANK.GT.0 ) THEN
CALL DAXPY( MNMIN, -ONE, COPYS, 1, S, 1 ) CALL DAXPY( MNMIN, -ONE, COPYS, 1, S, 1 )
RESULT( 7 ) = DASUM( MNMIN, S, 1 ) / RESULT( 11 ) = DASUM( MNMIN, S, 1 ) /
$ DASUM( MNMIN, COPYS, 1 ) / $ DASUM( MNMIN, COPYS, 1 ) /
$ ( EPS*DBLE( MNMIN ) ) $ ( EPS*DBLE( MNMIN ) )
ELSE ELSE
RESULT( 7 ) = ZERO RESULT( 11 ) = ZERO
END IF END IF
* *
* Test 8: Compute error in solution * Test 12: Compute error in solution
* *
CALL DLACPY( 'Full', M, NRHS, COPYB, LDB, WORK, CALL DLACPY( 'Full', M, NRHS, COPYB, LDB, WORK,
$ LDWORK ) $ LDWORK )
CALL DQRT16( 'No transpose', M, N, NRHS, COPYA, CALL DQRT16( 'No transpose', M, N, NRHS, COPYA,
$ LDA, B, LDB, WORK, LDWORK, $ LDA, B, LDB, WORK, LDWORK,
$ WORK( M*NRHS+1 ), RESULT( 8 ) ) $ WORK( M*NRHS+1 ), RESULT( 12 ) )
* *
* Test 9: Check norm of r'*A * Test 13: Check norm of r'*A
* *
RESULT( 9 ) = ZERO RESULT( 13 ) = ZERO
IF( M.GT.CRANK ) IF( M.GT.CRANK )
$ RESULT( 9 ) = DQRT17( 'No transpose', 1, M, $ RESULT( 13 ) = DQRT17( 'No transpose', 1, M,
$ N, NRHS, COPYA, LDA, B, LDB, $ N, NRHS, COPYA, LDA, B, LDB,
$ COPYB, LDB, C, WORK, LWORK ) $ COPYB, LDB, C, WORK, LWORK )
* *
* Test 10: Check if x is in the rowspace of A * Test 14: Check if x is in the rowspace of A
* *
RESULT( 10 ) = ZERO RESULT( 14 ) = ZERO
IF( N.GT.CRANK ) IF( N.GT.CRANK )
$ RESULT( 10 ) = DQRT14( 'No transpose', M, N, $ RESULT( 14 ) = DQRT14( 'No transpose', M, N,
$ NRHS, COPYA, LDA, B, LDB, $ NRHS, COPYA, LDA, B, LDB,
$ WORK, LWORK ) $ WORK, LWORK )
* *
@ -776,45 +937,45 @@
$ N, NRHS, -1, NB, ITYPE, NFAIL, $ N, NRHS, -1, NB, ITYPE, NFAIL,
$ NERRS, NOUT ) $ NERRS, NOUT )
* *
* Test 11: Compute relative error in svd * Test 15: Compute relative error in svd
* *
IF( RANK.GT.0 ) THEN IF( RANK.GT.0 ) THEN
CALL DAXPY( MNMIN, -ONE, COPYS, 1, S, 1 ) CALL DAXPY( MNMIN, -ONE, COPYS, 1, S, 1 )
RESULT( 11 ) = DASUM( MNMIN, S, 1 ) / RESULT( 15 ) = DASUM( MNMIN, S, 1 ) /
$ DASUM( MNMIN, COPYS, 1 ) / $ DASUM( MNMIN, COPYS, 1 ) /
$ ( EPS*DBLE( MNMIN ) ) $ ( EPS*DBLE( MNMIN ) )
ELSE ELSE
RESULT( 11 ) = ZERO RESULT( 15 ) = ZERO
END IF END IF
* *
* Test 12: Compute error in solution * Test 16: Compute error in solution
* *
CALL DLACPY( 'Full', M, NRHS, COPYB, LDB, WORK, CALL DLACPY( 'Full', M, NRHS, COPYB, LDB, WORK,
$ LDWORK ) $ LDWORK )
CALL DQRT16( 'No transpose', M, N, NRHS, COPYA, CALL DQRT16( 'No transpose', M, N, NRHS, COPYA,
$ LDA, B, LDB, WORK, LDWORK, $ LDA, B, LDB, WORK, LDWORK,
$ WORK( M*NRHS+1 ), RESULT( 12 ) ) $ WORK( M*NRHS+1 ), RESULT( 16 ) )
* *
* Test 13: Check norm of r'*A * Test 17: Check norm of r'*A
* *
RESULT( 13 ) = ZERO RESULT( 17 ) = ZERO
IF( M.GT.CRANK ) IF( M.GT.CRANK )
$ RESULT( 13 ) = DQRT17( 'No transpose', 1, M, $ RESULT( 17 ) = DQRT17( 'No transpose', 1, M,
$ N, NRHS, COPYA, LDA, B, LDB, $ N, NRHS, COPYA, LDA, B, LDB,
$ COPYB, LDB, C, WORK, LWORK ) $ COPYB, LDB, C, WORK, LWORK )
* *
* Test 14: Check if x is in the rowspace of A * Test 18: Check if x is in the rowspace of A
* *
RESULT( 14 ) = ZERO RESULT( 18 ) = ZERO
IF( N.GT.CRANK ) IF( N.GT.CRANK )
$ RESULT( 14 ) = DQRT14( 'No transpose', M, N, $ RESULT( 18 ) = DQRT14( 'No transpose', M, N,
$ NRHS, COPYA, LDA, B, LDB, $ NRHS, COPYA, LDA, B, LDB,
$ WORK, LWORK ) $ WORK, LWORK )
* *
* Print information about the tests that did not * Print information about the tests that did not
* pass the threshold. * pass the threshold.
* *
DO 90 K = 3, 14 DO 90 K = 7, 18
IF( RESULT( K ).GE.THRESH ) THEN IF( RESULT( K ).GE.THRESH ) THEN
IF( NFAIL.EQ.0 .AND. NERRS.EQ.0 ) IF( NFAIL.EQ.0 .AND. NERRS.EQ.0 )
$ CALL ALAHD( NOUT, PATH ) $ CALL ALAHD( NOUT, PATH )
@ -826,6 +987,12 @@
NRUN = NRUN + 12 NRUN = NRUN + 12
* *
100 CONTINUE 100 CONTINUE
110 CONTINUE 110 CONTINUE
120 CONTINUE 120 CONTINUE
130 CONTINUE 130 CONTINUE

61
lapack-netlib/TESTING/LIN/derrls.f
View File

@ -22,7 +22,7 @@
*> \verbatim *> \verbatim
*> *>
*> DERRLS tests the error exits for the DOUBLE PRECISION least squares *> DERRLS tests the error exits for the DOUBLE PRECISION least squares
*> driver routines (DGELS, SGELSS, SGELSY, SGELSD). *> driver routines (DGELS, DGELST, DGETSLS, SGELSS, SGELSY, SGELSD).
*> \endverbatim *> \endverbatim
* *
* Arguments: * Arguments:
@ -83,7 +83,8 @@
EXTERNAL LSAMEN EXTERNAL LSAMEN
* .. * ..
* .. External Subroutines .. * .. External Subroutines ..
EXTERNAL ALAESM, CHKXER, DGELS, DGELSD, DGELSS, DGELSY EXTERNAL ALAESM, CHKXER, DGELS, DGELSD, DGELSS, DGELST,
$ DGELSY, DGETSLS
* .. * ..
* .. Scalars in Common .. * .. Scalars in Common ..
LOGICAL LERR, OK LOGICAL LERR, OK
@ -130,10 +131,66 @@
INFOT = 8 INFOT = 8
CALL DGELS( 'N', 2, 0, 0, A, 2, B, 1, W, 2, INFO ) CALL DGELS( 'N', 2, 0, 0, A, 2, B, 1, W, 2, INFO )
CALL CHKXER( 'DGELS ', INFOT, NOUT, LERR, OK ) CALL CHKXER( 'DGELS ', INFOT, NOUT, LERR, OK )
INFOT = 8
CALL DGELS( 'N', 0, 2, 0, A, 1, B, 1, W, 2, INFO )
CALL CHKXER( 'DGELS', INFOT, NOUT, LERR, OK )
INFOT = 10 INFOT = 10
CALL DGELS( 'N', 1, 1, 0, A, 1, B, 1, W, 1, INFO ) CALL DGELS( 'N', 1, 1, 0, A, 1, B, 1, W, 1, INFO )
CALL CHKXER( 'DGELS ', INFOT, NOUT, LERR, OK ) CALL CHKXER( 'DGELS ', INFOT, NOUT, LERR, OK )
* *
* DGELST
*
SRNAMT = 'DGELST'
INFOT = 1
CALL DGELST( '/', 0, 0, 0, A, 1, B, 1, W, 1, INFO )
CALL CHKXER( 'DGELST', INFOT, NOUT, LERR, OK )
INFOT = 2
CALL DGELST( 'N', -1, 0, 0, A, 1, B, 1, W, 1, INFO )
CALL CHKXER( 'DGELST', INFOT, NOUT, LERR, OK )
INFOT = 3
CALL DGELST( 'N', 0, -1, 0, A, 1, B, 1, W, 1, INFO )
CALL CHKXER( 'DGELST', INFOT, NOUT, LERR, OK )
INFOT = 4
CALL DGELST( 'N', 0, 0, -1, A, 1, B, 1, W, 1, INFO )
CALL CHKXER( 'DGELST', INFOT, NOUT, LERR, OK )
INFOT = 6
CALL DGELST( 'N', 2, 0, 0, A, 1, B, 2, W, 2, INFO )
CALL CHKXER( 'DGELST', INFOT, NOUT, LERR, OK )
INFOT = 8
CALL DGELST( 'N', 2, 0, 0, A, 2, B, 1, W, 2, INFO )
CALL CHKXER( 'DGELST', INFOT, NOUT, LERR, OK )
INFOT = 8
CALL DGELST( 'N', 0, 2, 0, A, 1, B, 1, W, 2, INFO )
CALL CHKXER( 'DGELST', INFOT, NOUT, LERR, OK )
INFOT = 10
CALL DGELST( 'N', 1, 1, 0, A, 1, B, 1, W, 1, INFO )
CALL CHKXER( 'DGELST', INFOT, NOUT, LERR, OK )
*
* DGETSLS
*
SRNAMT = 'DGETSLS'
INFOT = 1
CALL DGETSLS( '/', 0, 0, 0, A, 1, B, 1, W, 1, INFO )
CALL CHKXER( 'DGETSLS', INFOT, NOUT, LERR, OK )
INFOT = 2
CALL DGETSLS( 'N', -1, 0, 0, A, 1, B, 1, W, 1, INFO )
CALL CHKXER( 'DGETSLS', INFOT, NOUT, LERR, OK )
INFOT = 3
CALL DGETSLS( 'N', 0, -1, 0, A, 1, B, 1, W, 1, INFO )
CALL CHKXER( 'DGETSLS', INFOT, NOUT, LERR, OK )
INFOT = 4
CALL DGETSLS( 'N', 0, 0, -1, A, 1, B, 1, W, 1, INFO )
CALL CHKXER( 'DGETSLS', INFOT, NOUT, LERR, OK )
INFOT = 6
CALL DGETSLS( 'N', 2, 0, 0, A, 1, B, 2, W, 2, INFO )
CALL CHKXER( 'DGETSLS', INFOT, NOUT, LERR, OK )
INFOT = 8
CALL DGETSLS( 'N', 2, 0, 0, A, 2, B, 1, W, 2, INFO )
CALL CHKXER( 'DGETSLS', INFOT, NOUT, LERR, OK )
INFOT = 8
CALL DGETSLS( 'N', 0, 2, 0, A, 1, B, 1, W, 2, INFO )
CALL CHKXER( 'DGETSLS', INFOT, NOUT, LERR, OK )
*
* DGELSS * DGELSS
* *
SRNAMT = 'DGELSS' SRNAMT = 'DGELSS'

331
lapack-netlib/TESTING/LIN/sdrvls.f
View File

@ -31,8 +31,8 @@
*> *>
*> \verbatim *> \verbatim
*> *>
*> SDRVLS tests the least squares driver routines SGELS, SGETSLS, SGELSS, SGELSY, *> SDRVLS tests the least squares driver routines SGELS, SGELST,
*> and SGELSD. *> SGETSLS, SGELSS, SGELSY and SGELSD.
*> \endverbatim *> \endverbatim
* *
* Arguments: * Arguments:
@ -211,7 +211,7 @@
* *
* .. Parameters .. * .. Parameters ..
INTEGER NTESTS INTEGER NTESTS
PARAMETER ( NTESTS = 16 ) PARAMETER ( NTESTS = 18 )
INTEGER SMLSIZ INTEGER SMLSIZ
PARAMETER ( SMLSIZ = 25 ) PARAMETER ( SMLSIZ = 25 )
REAL ONE, TWO, ZERO REAL ONE, TWO, ZERO
@ -225,8 +225,8 @@
$ LWLSY, LWORK, M, MNMIN, N, NB, NCOLS, NERRS, $ LWLSY, LWORK, M, MNMIN, N, NB, NCOLS, NERRS,
$ NFAIL, NRHS, NROWS, NRUN, RANK, MB, $ NFAIL, NRHS, NROWS, NRUN, RANK, MB,
$ MMAX, NMAX, NSMAX, LIWORK, $ MMAX, NMAX, NSMAX, LIWORK,
$ LWORK_SGELS, LWORK_SGETSLS, LWORK_SGELSS, $ LWORK_SGELS, LWORK_SGELST, LWORK_SGETSLS,
$ LWORK_SGELSY, LWORK_SGELSD $ LWORK_SGELSS, LWORK_SGELSY, LWORK_SGELSD
REAL EPS, NORMA, NORMB, RCOND REAL EPS, NORMA, NORMB, RCOND
* .. * ..
* .. Local Arrays .. * .. Local Arrays ..
@ -243,12 +243,12 @@
* .. * ..
* .. External Subroutines .. * .. External Subroutines ..
EXTERNAL ALAERH, ALAHD, ALASVM, SAXPY, SERRLS, SGELS, EXTERNAL ALAERH, ALAHD, ALASVM, SAXPY, SERRLS, SGELS,
$ SGELSD, SGELSS, SGELSY, SGEMM, SLACPY, $ SGELSD, SGELSS, SGELST, SGELSY, SGEMM,
$ SLARNV, SQRT13, SQRT15, SQRT16, SSCAL, $ SGETSLS, SLACPY, SLARNV, SQRT13, SQRT15,
$ XLAENV, SGETSLS $ SQRT16, SSCAL, XLAENV
* .. * ..
* .. Intrinsic Functions .. * .. Intrinsic Functions ..
INTRINSIC INT, LOG, MAX, MIN, REAL, SQRT INTRINSIC INT, MAX, MIN, REAL, SQRT
* .. * ..
* .. Scalars in Common .. * .. Scalars in Common ..
LOGICAL LERR, OK LOGICAL LERR, OK
@ -330,7 +330,8 @@
LIWORK = 1 LIWORK = 1
* *
* Iterate through all test cases and compute necessary workspace * Iterate through all test cases and compute necessary workspace
* sizes for ?GELS, ?GETSLS, ?GELSY, ?GELSS and ?GELSD routines. * sizes for ?GELS, ?GELST, ?GETSLS, ?GELSY, ?GELSS and ?GELSD
* routines.
* *
DO IM = 1, NM DO IM = 1, NM
M = MVAL( IM ) M = MVAL( IM )
@ -357,6 +358,10 @@
CALL SGELS( TRANS, M, N, NRHS, A, LDA, CALL SGELS( TRANS, M, N, NRHS, A, LDA,
$ B, LDB, WQ( 1 ), -1, INFO ) $ B, LDB, WQ( 1 ), -1, INFO )
LWORK_SGELS = INT ( WQ( 1 ) ) LWORK_SGELS = INT ( WQ( 1 ) )
* Compute workspace needed for SGELST
CALL SGELST( TRANS, M, N, NRHS, A, LDA,
$ B, LDB, WQ, -1, INFO )
LWORK_SGELST = INT ( WQ ( 1 ) )
* Compute workspace needed for SGETSLS * Compute workspace needed for SGETSLS
CALL SGETSLS( TRANS, M, N, NRHS, A, LDA, CALL SGETSLS( TRANS, M, N, NRHS, A, LDA,
$ B, LDB, WQ( 1 ), -1, INFO ) $ B, LDB, WQ( 1 ), -1, INFO )
@ -378,9 +383,9 @@
* Compute LIWORK workspace needed for SGELSY and SGELSD * Compute LIWORK workspace needed for SGELSY and SGELSD
LIWORK = MAX( LIWORK, N, IWQ( 1 ) ) LIWORK = MAX( LIWORK, N, IWQ( 1 ) )
* Compute LWORK workspace needed for all functions * Compute LWORK workspace needed for all functions
LWORK = MAX( LWORK, LWORK_SGELS, LWORK_SGETSLS, LWORK = MAX( LWORK, LWORK_SGELS, LWORK_SGELST,
$ LWORK_SGELSY, LWORK_SGELSS, $ LWORK_SGETSLS, LWORK_SGELSY,
$ LWORK_SGELSD ) $ LWORK_SGELSS, LWORK_SGELSD )
END IF END IF
ENDDO ENDDO
ENDDO ENDDO
@ -411,21 +416,26 @@
ITYPE = ( IRANK-1 )*3 + ISCALE ITYPE = ( IRANK-1 )*3 + ISCALE
IF( .NOT.DOTYPE( ITYPE ) ) IF( .NOT.DOTYPE( ITYPE ) )
$ GO TO 110 $ GO TO 110
* * =====================================================
* Begin test SGELS
* =====================================================
IF( IRANK.EQ.1 ) THEN IF( IRANK.EQ.1 ) THEN
* *
* Test SGELS
*
* Generate a matrix of scaling type ISCALE * Generate a matrix of scaling type ISCALE
* *
CALL SQRT13( ISCALE, M, N, COPYA, LDA, NORMA, CALL SQRT13( ISCALE, M, N, COPYA, LDA, NORMA,
$ ISEED ) $ ISEED )
DO 40 INB = 1, NNB *
* Loop for testing different block sizes.
*
DO INB = 1, NNB
NB = NBVAL( INB ) NB = NBVAL( INB )
CALL XLAENV( 1, NB ) CALL XLAENV( 1, NB )
CALL XLAENV( 3, NXVAL( INB ) ) CALL XLAENV( 3, NXVAL( INB ) )
* *
DO 30 ITRAN = 1, 2 * Loop for testing non-transposed and transposed.
*
DO ITRAN = 1, 2
IF( ITRAN.EQ.1 ) THEN IF( ITRAN.EQ.1 ) THEN
TRANS = 'N' TRANS = 'N'
NROWS = M NROWS = M
@ -469,20 +479,27 @@
$ ITYPE, NFAIL, NERRS, $ ITYPE, NFAIL, NERRS,
$ NOUT ) $ NOUT )
* *
* Check correctness of results * Test 1: Check correctness of results
* for SGELS, compute the residual:
* RESID = norm(B - A*X) /
* / ( max(m,n) * norm(A) * norm(X) * EPS )
* *
LDWORK = MAX( 1, NROWS )
IF( NROWS.GT.0 .AND. NRHS.GT.0 ) IF( NROWS.GT.0 .AND. NRHS.GT.0 )
$ CALL SLACPY( 'Full', NROWS, NRHS, $ CALL SLACPY( 'Full', NROWS, NRHS,
$ COPYB, LDB, C, LDB ) $ COPYB, LDB, C, LDB )
CALL SQRT16( TRANS, M, N, NRHS, COPYA, CALL SQRT16( TRANS, M, N, NRHS, COPYA,
$ LDA, B, LDB, C, LDB, WORK, $ LDA, B, LDB, C, LDB, WORK,
$ RESULT( 1 ) ) $ RESULT( 1 ) )
*
* Test 2: Check correctness of results
* for SGELS.
* *
IF( ( ITRAN.EQ.1 .AND. M.GE.N ) .OR. IF( ( ITRAN.EQ.1 .AND. M.GE.N ) .OR.
$ ( ITRAN.EQ.2 .AND. M.LT.N ) ) THEN $ ( ITRAN.EQ.2 .AND. M.LT.N ) ) THEN
* *
* Solving LS system * Solving LS system, compute:
* r = norm((B- A*X)**T * A) /
* / (norm(A)*norm(B)*max(M,N,NRHS)*EPS)
* *
RESULT( 2 ) = SQRT17( TRANS, 1, M, N, RESULT( 2 ) = SQRT17( TRANS, 1, M, N,
$ NRHS, COPYA, LDA, B, LDB, $ NRHS, COPYA, LDA, B, LDB,
@ -500,7 +517,7 @@
* Print information about the tests that * Print information about the tests that
* did not pass the threshold. * did not pass the threshold.
* *
DO 20 K = 1, 2 DO K = 1, 2
IF( RESULT( K ).GE.THRESH ) THEN IF( RESULT( K ).GE.THRESH ) THEN
IF( NFAIL.EQ.0 .AND. NERRS.EQ.0 ) IF( NFAIL.EQ.0 .AND. NERRS.EQ.0 )
$ CALL ALAHD( NOUT, PATH ) $ CALL ALAHD( NOUT, PATH )
@ -509,26 +526,33 @@
$ RESULT( K ) $ RESULT( K )
NFAIL = NFAIL + 1 NFAIL = NFAIL + 1
END IF END IF
20 CONTINUE END DO
NRUN = NRUN + 2 NRUN = NRUN + 2
30 CONTINUE END DO
40 CONTINUE END DO
* END IF
* * =====================================================
* Test SGETSLS * End test SGELS
* =====================================================
* =====================================================
* Begin test SGELST
* =====================================================
IF( IRANK.EQ.1 ) THEN
* *
* Generate a matrix of scaling type ISCALE * Generate a matrix of scaling type ISCALE
* *
CALL SQRT13( ISCALE, M, N, COPYA, LDA, NORMA, CALL SQRT13( ISCALE, M, N, COPYA, LDA, NORMA,
$ ISEED ) $ ISEED )
DO 65 INB = 1, NNB
MB = NBVAL( INB )
CALL XLAENV( 1, MB )
DO 62 IMB = 1, NNB
NB = NBVAL( IMB )
CALL XLAENV( 2, NB )
* *
DO 60 ITRAN = 1, 2 * Loop for testing different block sizes.
*
DO INB = 1, NNB
NB = NBVAL( INB )
CALL XLAENV( 1, NB )
*
* Loop for testing non-transposed and transposed.
*
DO ITRAN = 1, 2
IF( ITRAN.EQ.1 ) THEN IF( ITRAN.EQ.1 ) THEN
TRANS = 'N' TRANS = 'N'
NROWS = M NROWS = M
@ -563,31 +587,38 @@
CALL SLACPY( 'Full', NROWS, NRHS, CALL SLACPY( 'Full', NROWS, NRHS,
$ COPYB, LDB, B, LDB ) $ COPYB, LDB, B, LDB )
END IF END IF
SRNAMT = 'SGETSLS ' SRNAMT = 'SGELST'
CALL SGETSLS( TRANS, M, N, NRHS, A, CALL SGELST( TRANS, M, N, NRHS, A, LDA, B,
$ LDA, B, LDB, WORK, LWORK, INFO ) $ LDB, WORK, LWORK, INFO )
IF( INFO.NE.0 ) IF( INFO.NE.0 )
$ CALL ALAERH( PATH, 'SGETSLS ', INFO, 0, $ CALL ALAERH( PATH, 'SGELST', INFO, 0,
$ TRANS, M, N, NRHS, -1, NB, $ TRANS, M, N, NRHS, -1, NB,
$ ITYPE, NFAIL, NERRS, $ ITYPE, NFAIL, NERRS,
$ NOUT ) $ NOUT )
* *
* Check correctness of results * Test 3: Check correctness of results
* for SGELST, compute the residual:
* RESID = norm(B - A*X) /
* / ( max(m,n) * norm(A) * norm(X) * EPS )
* *
LDWORK = MAX( 1, NROWS )
IF( NROWS.GT.0 .AND. NRHS.GT.0 ) IF( NROWS.GT.0 .AND. NRHS.GT.0 )
$ CALL SLACPY( 'Full', NROWS, NRHS, $ CALL SLACPY( 'Full', NROWS, NRHS,
$ COPYB, LDB, C, LDB ) $ COPYB, LDB, C, LDB )
CALL SQRT16( TRANS, M, N, NRHS, COPYA, CALL SQRT16( TRANS, M, N, NRHS, COPYA,
$ LDA, B, LDB, C, LDB, WORK, $ LDA, B, LDB, C, LDB, WORK,
$ RESULT( 15 ) ) $ RESULT( 3 ) )
*
* Test 4: Check correctness of results
* for SGELST.
* *
IF( ( ITRAN.EQ.1 .AND. M.GE.N ) .OR. IF( ( ITRAN.EQ.1 .AND. M.GE.N ) .OR.
$ ( ITRAN.EQ.2 .AND. M.LT.N ) ) THEN $ ( ITRAN.EQ.2 .AND. M.LT.N ) ) THEN
* *
* Solving LS system * Solving LS system, compute:
* r = norm((B- A*X)**T * A) /
* / (norm(A)*norm(B)*max(M,N,NRHS)*EPS)
* *
RESULT( 16 ) = SQRT17( TRANS, 1, M, N, RESULT( 4 ) = SQRT17( TRANS, 1, M, N,
$ NRHS, COPYA, LDA, B, LDB, $ NRHS, COPYA, LDA, B, LDB,
$ COPYB, LDB, C, WORK, $ COPYB, LDB, C, WORK,
$ LWORK ) $ LWORK )
@ -595,7 +626,7 @@
* *
* Solving overdetermined system * Solving overdetermined system
* *
RESULT( 16 ) = SQRT14( TRANS, M, N, RESULT( 4 ) = SQRT14( TRANS, M, N,
$ NRHS, COPYA, LDA, B, LDB, $ NRHS, COPYA, LDA, B, LDB,
$ WORK, LWORK ) $ WORK, LWORK )
END IF END IF
@ -603,21 +634,151 @@
* Print information about the tests that * Print information about the tests that
* did not pass the threshold. * did not pass the threshold.
* *
DO 50 K = 15, 16 DO K = 3, 4
IF( RESULT( K ).GE.THRESH ) THEN IF( RESULT( K ).GE.THRESH ) THEN
IF( NFAIL.EQ.0 .AND. NERRS.EQ.0 ) IF( NFAIL.EQ.0 .AND. NERRS.EQ.0 )
$ CALL ALAHD( NOUT, PATH ) $ CALL ALAHD( NOUT, PATH )
WRITE( NOUT, FMT = 9997 )TRANS, M, WRITE( NOUT, FMT = 9999 ) TRANS, M,
$ N, NRHS, MB, NB, ITYPE, K, $ N, NRHS, NB, ITYPE, K,
$ RESULT( K ) $ RESULT( K )
NFAIL = NFAIL + 1 NFAIL = NFAIL + 1
END IF END IF
50 CONTINUE END DO
NRUN = NRUN + 2 NRUN = NRUN + 2
60 CONTINUE END DO
62 CONTINUE END DO
65 CONTINUE
END IF END IF
* =====================================================
* End test SGELST
* =====================================================
* =====================================================
* Begin test SGETSLS
* =====================================================
IF( IRANK.EQ.1 ) THEN
*
* Generate a matrix of scaling type ISCALE
*
CALL SQRT13( ISCALE, M, N, COPYA, LDA, NORMA,
$ ISEED )
*
* Loop for testing different block sizes MB.
*
DO IMB = 1, NNB
MB = NBVAL( IMB )
CALL XLAENV( 1, MB )
*
* Loop for testing different block sizes NB.
*
DO INB = 1, NNB
NB = NBVAL( INB )
CALL XLAENV( 2, NB )
*
* Loop for testing non-transposed
* and transposed.
*
DO ITRAN = 1, 2
IF( ITRAN.EQ.1 ) THEN
TRANS = 'N'
NROWS = M
NCOLS = N
ELSE
TRANS = 'T'
NROWS = N
NCOLS = M
END IF
LDWORK = MAX( 1, NCOLS )
*
* Set up a consistent rhs
*
IF( NCOLS.GT.0 ) THEN
CALL SLARNV( 2, ISEED, NCOLS*NRHS,
$ WORK )
CALL SSCAL( NCOLS*NRHS,
$ ONE / REAL( NCOLS ),
$ WORK, 1 )
END IF
CALL SGEMM( TRANS, 'No transpose',
$ NROWS, NRHS, NCOLS, ONE,
$ COPYA, LDA, WORK, LDWORK,
$ ZERO, B, LDB )
CALL SLACPY( 'Full', NROWS, NRHS,
$ B, LDB, COPYB, LDB )
*
* Solve LS or overdetermined system
*
IF( M.GT.0 .AND. N.GT.0 ) THEN
CALL SLACPY( 'Full', M, N,
$ COPYA, LDA, A, LDA )
CALL SLACPY( 'Full', NROWS, NRHS,
$ COPYB, LDB, B, LDB )
END IF
SRNAMT = 'SGETSLS'
CALL SGETSLS( TRANS, M, N, NRHS,
$ A, LDA, B, LDB, WORK, LWORK,
$ INFO )
IF( INFO.NE.0 )
$ CALL ALAERH( PATH, 'SGETSLS', INFO,
$ 0, TRANS, M, N, NRHS,
$ -1, NB, ITYPE, NFAIL,
$ NERRS, NOUT )
*
* Test 5: Check correctness of results
* for SGETSLS, compute the residual:
* RESID = norm(B - A*X) /
* / ( max(m,n) * norm(A) * norm(X) * EPS )
*
IF( NROWS.GT.0 .AND. NRHS.GT.0 )
$ CALL SLACPY( 'Full', NROWS, NRHS,
$ COPYB, LDB, C, LDB )
CALL SQRT16( TRANS, M, N, NRHS,
$ COPYA, LDA, B, LDB,
$ C, LDB, WORK,
$ RESULT( 5 ) )
*
* Test 6: Check correctness of results
* for SGETSLS.
*
IF( ( ITRAN.EQ.1 .AND. M.GE.N ) .OR.
$ ( ITRAN.EQ.2 .AND. M.LT.N ) ) THEN
*
* Solving LS system, compute:
* r = norm((B- A*X)**T * A) /
* / (norm(A)*norm(B)*max(M,N,NRHS)*EPS)
*
RESULT( 6 ) = SQRT17( TRANS, 1, M,
$ N, NRHS, COPYA, LDA,
$ B, LDB, COPYB, LDB,
$ C, WORK, LWORK )
ELSE
*
* Solving overdetermined system
*
RESULT( 6 ) = SQRT14( TRANS, M, N,
$ NRHS, COPYA, LDA,
$ B, LDB, WORK, LWORK )
END IF
*
* Print information about the tests that
* did not pass the threshold.
*
DO K = 5, 6
IF( RESULT( K ).GE.THRESH ) THEN
IF( NFAIL.EQ.0 .AND. NERRS.EQ.0 )
$ CALL ALAHD( NOUT, PATH )
WRITE( NOUT, FMT = 9997 ) TRANS,
$ M, N, NRHS, MB, NB, ITYPE,
$ K, RESULT( K )
NFAIL = NFAIL + 1
END IF
END DO
NRUN = NRUN + 2
END DO
END DO
END DO
END IF
* =====================================================
* End test SGETSLS
* =====================================================
* *
* Generate a matrix of scaling type ISCALE and rank * Generate a matrix of scaling type ISCALE and rank
* type IRANK. * type IRANK.
@ -662,37 +823,37 @@
$ N, NRHS, -1, NB, ITYPE, NFAIL, $ N, NRHS, -1, NB, ITYPE, NFAIL,
$ NERRS, NOUT ) $ NERRS, NOUT )
* *
* Test 3: Compute relative error in svd * Test 7: Compute relative error in svd
* workspace: M*N + 4*MIN(M,N) + MAX(M,N) * workspace: M*N + 4*MIN(M,N) + MAX(M,N)
* *
RESULT( 3 ) = SQRT12( CRANK, CRANK, A, LDA, RESULT( 7 ) = SQRT12( CRANK, CRANK, A, LDA,
$ COPYS, WORK, LWORK ) $ COPYS, WORK, LWORK )
* *
* Test 4: Compute error in solution * Test 8: Compute error in solution
* workspace: M*NRHS + M * workspace: M*NRHS + M
* *
CALL SLACPY( 'Full', M, NRHS, COPYB, LDB, WORK, CALL SLACPY( 'Full', M, NRHS, COPYB, LDB, WORK,
$ LDWORK ) $ LDWORK )
CALL SQRT16( 'No transpose', M, N, NRHS, COPYA, CALL SQRT16( 'No transpose', M, N, NRHS, COPYA,
$ LDA, B, LDB, WORK, LDWORK, $ LDA, B, LDB, WORK, LDWORK,
$ WORK( M*NRHS+1 ), RESULT( 4 ) ) $ WORK( M*NRHS+1 ), RESULT( 8 ) )
* *
* Test 5: Check norm of r'*A * Test 9: Check norm of r'*A
* workspace: NRHS*(M+N) * workspace: NRHS*(M+N)
* *
RESULT( 5 ) = ZERO RESULT( 9 ) = ZERO
IF( M.GT.CRANK ) IF( M.GT.CRANK )
$ RESULT( 5 ) = SQRT17( 'No transpose', 1, M, $ RESULT( 9 ) = SQRT17( 'No transpose', 1, M,
$ N, NRHS, COPYA, LDA, B, LDB, $ N, NRHS, COPYA, LDA, B, LDB,
$ COPYB, LDB, C, WORK, LWORK ) $ COPYB, LDB, C, WORK, LWORK )
* *
* Test 6: Check if x is in the rowspace of A * Test 10: Check if x is in the rowspace of A
* workspace: (M+NRHS)*(N+2) * workspace: (M+NRHS)*(N+2)
* *
RESULT( 6 ) = ZERO RESULT( 10 ) = ZERO
* *
IF( N.GT.CRANK ) IF( N.GT.CRANK )
$ RESULT( 6 ) = SQRT14( 'No transpose', M, N, $ RESULT( 10 ) = SQRT14( 'No transpose', M, N,
$ NRHS, COPYA, LDA, B, LDB, $ NRHS, COPYA, LDA, B, LDB,
$ WORK, LWORK ) $ WORK, LWORK )
* *
@ -716,38 +877,38 @@
* workspace used: 3*min(m,n) + * workspace used: 3*min(m,n) +
* max(2*min(m,n),nrhs,max(m,n)) * max(2*min(m,n),nrhs,max(m,n))
* *
* Test 7: Compute relative error in svd * Test 11: Compute relative error in svd
* *
IF( RANK.GT.0 ) THEN IF( RANK.GT.0 ) THEN
CALL SAXPY( MNMIN, -ONE, COPYS, 1, S, 1 ) CALL SAXPY( MNMIN, -ONE, COPYS, 1, S, 1 )
RESULT( 7 ) = SASUM( MNMIN, S, 1 ) / RESULT( 11 ) = SASUM( MNMIN, S, 1 ) /
$ SASUM( MNMIN, COPYS, 1 ) / $ SASUM( MNMIN, COPYS, 1 ) /
$ ( EPS*REAL( MNMIN ) ) $ ( EPS*REAL( MNMIN ) )
ELSE ELSE
RESULT( 7 ) = ZERO RESULT( 11 ) = ZERO
END IF END IF
* *
* Test 8: Compute error in solution * Test 12: Compute error in solution
* *
CALL SLACPY( 'Full', M, NRHS, COPYB, LDB, WORK, CALL SLACPY( 'Full', M, NRHS, COPYB, LDB, WORK,
$ LDWORK ) $ LDWORK )
CALL SQRT16( 'No transpose', M, N, NRHS, COPYA, CALL SQRT16( 'No transpose', M, N, NRHS, COPYA,
$ LDA, B, LDB, WORK, LDWORK, $ LDA, B, LDB, WORK, LDWORK,
$ WORK( M*NRHS+1 ), RESULT( 8 ) ) $ WORK( M*NRHS+1 ), RESULT( 12 ) )
* *
* Test 9: Check norm of r'*A * Test 13: Check norm of r'*A
* *
RESULT( 9 ) = ZERO RESULT( 13 ) = ZERO
IF( M.GT.CRANK ) IF( M.GT.CRANK )
$ RESULT( 9 ) = SQRT17( 'No transpose', 1, M, $ RESULT( 13 ) = SQRT17( 'No transpose', 1, M,
$ N, NRHS, COPYA, LDA, B, LDB, $ N, NRHS, COPYA, LDA, B, LDB,
$ COPYB, LDB, C, WORK, LWORK ) $ COPYB, LDB, C, WORK, LWORK )
* *
* Test 10: Check if x is in the rowspace of A * Test 14: Check if x is in the rowspace of A
* *
RESULT( 10 ) = ZERO RESULT( 14 ) = ZERO
IF( N.GT.CRANK ) IF( N.GT.CRANK )
$ RESULT( 10 ) = SQRT14( 'No transpose', M, N, $ RESULT( 14 ) = SQRT14( 'No transpose', M, N,
$ NRHS, COPYA, LDA, B, LDB, $ NRHS, COPYA, LDA, B, LDB,
$ WORK, LWORK ) $ WORK, LWORK )
* *
@ -776,45 +937,45 @@
$ N, NRHS, -1, NB, ITYPE, NFAIL, $ N, NRHS, -1, NB, ITYPE, NFAIL,
$ NERRS, NOUT ) $ NERRS, NOUT )
* *
* Test 11: Compute relative error in svd * Test 15: Compute relative error in svd
* *
IF( RANK.GT.0 ) THEN IF( RANK.GT.0 ) THEN
CALL SAXPY( MNMIN, -ONE, COPYS, 1, S, 1 ) CALL SAXPY( MNMIN, -ONE, COPYS, 1, S, 1 )
RESULT( 11 ) = SASUM( MNMIN, S, 1 ) / RESULT( 15 ) = SASUM( MNMIN, S, 1 ) /
$ SASUM( MNMIN, COPYS, 1 ) / $ SASUM( MNMIN, COPYS, 1 ) /
$ ( EPS*REAL( MNMIN ) ) $ ( EPS*REAL( MNMIN ) )
ELSE ELSE
RESULT( 11 ) = ZERO RESULT( 15 ) = ZERO
END IF END IF
* *
* Test 12: Compute error in solution * Test 16: Compute error in solution
* *
CALL SLACPY( 'Full', M, NRHS, COPYB, LDB, WORK, CALL SLACPY( 'Full', M, NRHS, COPYB, LDB, WORK,
$ LDWORK ) $ LDWORK )
CALL SQRT16( 'No transpose', M, N, NRHS, COPYA, CALL SQRT16( 'No transpose', M, N, NRHS, COPYA,
$ LDA, B, LDB, WORK, LDWORK, $ LDA, B, LDB, WORK, LDWORK,
$ WORK( M*NRHS+1 ), RESULT( 12 ) ) $ WORK( M*NRHS+1 ), RESULT( 16 ) )
* *
* Test 13: Check norm of r'*A * Test 17: Check norm of r'*A
* *
RESULT( 13 ) = ZERO RESULT( 17 ) = ZERO
IF( M.GT.CRANK ) IF( M.GT.CRANK )
$ RESULT( 13 ) = SQRT17( 'No transpose', 1, M, $ RESULT( 17 ) = SQRT17( 'No transpose', 1, M,
$ N, NRHS, COPYA, LDA, B, LDB, $ N, NRHS, COPYA, LDA, B, LDB,
$ COPYB, LDB, C, WORK, LWORK ) $ COPYB, LDB, C, WORK, LWORK )
* *
* Test 14: Check if x is in the rowspace of A * Test 18: Check if x is in the rowspace of A
* *
RESULT( 14 ) = ZERO RESULT( 18 ) = ZERO
IF( N.GT.CRANK ) IF( N.GT.CRANK )
$ RESULT( 14 ) = SQRT14( 'No transpose', M, N, $ RESULT( 18 ) = SQRT14( 'No transpose', M, N,
$ NRHS, COPYA, LDA, B, LDB, $ NRHS, COPYA, LDA, B, LDB,
$ WORK, LWORK ) $ WORK, LWORK )
* *
* Print information about the tests that did not * Print information about the tests that did not
* pass the threshold. * pass the threshold.
* *
DO 90 K = 3, 14 DO 90 K = 7, 18
IF( RESULT( K ).GE.THRESH ) THEN IF( RESULT( K ).GE.THRESH ) THEN
IF( NFAIL.EQ.0 .AND. NERRS.EQ.0 ) IF( NFAIL.EQ.0 .AND. NERRS.EQ.0 )
$ CALL ALAHD( NOUT, PATH ) $ CALL ALAHD( NOUT, PATH )

61
lapack-netlib/TESTING/LIN/serrls.f
View File

@ -22,7 +22,7 @@
*> \verbatim *> \verbatim
*> *>
*> SERRLS tests the error exits for the REAL least squares *> SERRLS tests the error exits for the REAL least squares
*> driver routines (SGELS, SGELSS, SGELSY, SGELSD). *> driver routines (SGELS, SGELST, SGETSLS, SGELSS, SGELSY, SGELSD).
*> \endverbatim *> \endverbatim
* *
* Arguments: * Arguments:
@ -83,7 +83,8 @@
EXTERNAL LSAMEN EXTERNAL LSAMEN
* .. * ..
* .. External Subroutines .. * .. External Subroutines ..
EXTERNAL ALAESM, CHKXER, SGELS, SGELSD, SGELSS, SGELSY EXTERNAL ALAESM, CHKXER, SGELS, SGELSD, SGELSS, SGELST,
$ SGELSY, SGETSLS
* .. * ..
* .. Scalars in Common .. * .. Scalars in Common ..
LOGICAL LERR, OK LOGICAL LERR, OK
@ -130,10 +131,66 @@
INFOT = 8 INFOT = 8
CALL SGELS( 'N', 2, 0, 0, A, 2, B, 1, W, 2, INFO ) CALL SGELS( 'N', 2, 0, 0, A, 2, B, 1, W, 2, INFO )
CALL CHKXER( 'SGELS ', INFOT, NOUT, LERR, OK ) CALL CHKXER( 'SGELS ', INFOT, NOUT, LERR, OK )
INFOT = 8
CALL SGELS( 'N', 0, 2, 0, A, 1, B, 1, W, 2, INFO )
CALL CHKXER( 'DGELS', INFOT, NOUT, LERR, OK )
INFOT = 10 INFOT = 10
CALL SGELS( 'N', 1, 1, 0, A, 1, B, 1, W, 1, INFO ) CALL SGELS( 'N', 1, 1, 0, A, 1, B, 1, W, 1, INFO )
CALL CHKXER( 'SGELS ', INFOT, NOUT, LERR, OK ) CALL CHKXER( 'SGELS ', INFOT, NOUT, LERR, OK )
* *
* SGELST
*
SRNAMT = 'SGELST'
INFOT = 1
CALL SGELST( '/', 0, 0, 0, A, 1, B, 1, W, 1, INFO )
CALL CHKXER( 'SGELST', INFOT, NOUT, LERR, OK )
INFOT = 2
CALL SGELST( 'N', -1, 0, 0, A, 1, B, 1, W, 1, INFO )
CALL CHKXER( 'SGELST', INFOT, NOUT, LERR, OK )
INFOT = 3
CALL SGELST( 'N', 0, -1, 0, A, 1, B, 1, W, 1, INFO )
CALL CHKXER( 'SGELST', INFOT, NOUT, LERR, OK )
INFOT = 4
CALL SGELST( 'N', 0, 0, -1, A, 1, B, 1, W, 1, INFO )
CALL CHKXER( 'SGELST', INFOT, NOUT, LERR, OK )
INFOT = 6
CALL SGELST( 'N', 2, 0, 0, A, 1, B, 2, W, 2, INFO )
CALL CHKXER( 'SGELST', INFOT, NOUT, LERR, OK )
INFOT = 8
CALL SGELST( 'N', 2, 0, 0, A, 2, B, 1, W, 2, INFO )
CALL CHKXER( 'SGELST', INFOT, NOUT, LERR, OK )
INFOT = 8
CALL SGELST( 'N', 0, 2, 0, A, 1, B, 1, W, 2, INFO )
CALL CHKXER( 'SGELST', INFOT, NOUT, LERR, OK )
INFOT = 10
CALL SGELST( 'N', 1, 1, 0, A, 1, B, 1, W, 1, INFO )
CALL CHKXER( 'SGELST', INFOT, NOUT, LERR, OK )
*
* SGETSLS
*
SRNAMT = 'SGETSLS'
INFOT = 1
CALL SGETSLS( '/', 0, 0, 0, A, 1, B, 1, W, 1, INFO )
CALL CHKXER( 'SGETSLS', INFOT, NOUT, LERR, OK )
INFOT = 2
CALL SGETSLS( 'N', -1, 0, 0, A, 1, B, 1, W, 1, INFO )
CALL CHKXER( 'SGETSLS', INFOT, NOUT, LERR, OK )
INFOT = 3
CALL SGETSLS( 'N', 0, -1, 0, A, 1, B, 1, W, 1, INFO )
CALL CHKXER( 'SGETSLS', INFOT, NOUT, LERR, OK )
INFOT = 4
CALL SGETSLS( 'N', 0, 0, -1, A, 1, B, 1, W, 1, INFO )
CALL CHKXER( 'SGETSLS', INFOT, NOUT, LERR, OK )
INFOT = 6
CALL SGETSLS( 'N', 2, 0, 0, A, 1, B, 2, W, 2, INFO )
CALL CHKXER( 'SGETSLS', INFOT, NOUT, LERR, OK )
INFOT = 8
CALL SGETSLS( 'N', 2, 0, 0, A, 2, B, 1, W, 2, INFO )
CALL CHKXER( 'SGETSLS', INFOT, NOUT, LERR, OK )
INFOT = 8
CALL SGETSLS( 'N', 0, 2, 0, A, 1, B, 1, W, 2, INFO )
CALL CHKXER( 'SGETSLS', INFOT, NOUT, LERR, OK )
*
* SGELSS * SGELSS
* *
SRNAMT = 'SGELSS' SRNAMT = 'SGELSS'

333
lapack-netlib/TESTING/LIN/zdrvls.f
View File

@ -31,8 +31,8 @@
*> *>
*> \verbatim *> \verbatim
*> *>
*> ZDRVLS tests the least squares driver routines ZGELS, ZGETSLS, ZGELSS, ZGELSY *> ZDRVLS tests the least squares driver routines ZGELS, ZGELST,
*> and ZGELSD. *> ZGETSLS, ZGELSS, ZGELSY and ZGELSD.
*> \endverbatim *> \endverbatim
* *
* Arguments: * Arguments:
@ -211,7 +211,7 @@
* *
* .. Parameters .. * .. Parameters ..
INTEGER NTESTS INTEGER NTESTS
PARAMETER ( NTESTS = 16 ) PARAMETER ( NTESTS = 18 )
INTEGER SMLSIZ INTEGER SMLSIZ
PARAMETER ( SMLSIZ = 25 ) PARAMETER ( SMLSIZ = 25 )
DOUBLE PRECISION ONE, ZERO DOUBLE PRECISION ONE, ZERO
@ -228,8 +228,8 @@
$ LWLSY, LWORK, M, MNMIN, N, NB, NCOLS, NERRS, $ LWLSY, LWORK, M, MNMIN, N, NB, NCOLS, NERRS,
$ NFAIL, NRHS, NROWS, NRUN, RANK, MB, $ NFAIL, NRHS, NROWS, NRUN, RANK, MB,
$ MMAX, NMAX, NSMAX, LIWORK, LRWORK, $ MMAX, NMAX, NSMAX, LIWORK, LRWORK,
$ LWORK_ZGELS, LWORK_ZGETSLS, LWORK_ZGELSS, $ LWORK_ZGELS, LWORK_ZGELST, LWORK_ZGETSLS,
$ LWORK_ZGELSY, LWORK_ZGELSD, $ LWORK_ZGELSS, LWORK_ZGELSY, LWORK_ZGELSD,
$ LRWORK_ZGELSY, LRWORK_ZGELSS, LRWORK_ZGELSD $ LRWORK_ZGELSY, LRWORK_ZGELSS, LRWORK_ZGELSD
DOUBLE PRECISION EPS, NORMA, NORMB, RCOND DOUBLE PRECISION EPS, NORMA, NORMB, RCOND
* .. * ..
@ -248,10 +248,10 @@
EXTERNAL DASUM, DLAMCH, ZQRT12, ZQRT14, ZQRT17 EXTERNAL DASUM, DLAMCH, ZQRT12, ZQRT14, ZQRT17
* .. * ..
* .. External Subroutines .. * .. External Subroutines ..
EXTERNAL ALAERH, ALAHD, ALASVM, DAXPY, DLASRT, XLAENV, EXTERNAL ALAERH, ALAHD, ALASVM, DAXPY, ZERRLS, ZGELS,
$ ZDSCAL, ZERRLS, ZGELS, ZGELSD, ZGELSS, $ ZGELSD, ZGELSS, ZGELST, ZGELSY, ZGEMM,
$ ZGELSY, ZGEMM, ZLACPY, ZLARNV, ZQRT13, ZQRT15, $ ZGETSLS, ZLACPY, ZLARNV, ZQRT13, ZQRT15,
$ ZQRT16, ZGETSLS $ ZQRT16, ZDSCAL, XLAENV
* .. * ..
* .. Intrinsic Functions .. * .. Intrinsic Functions ..
INTRINSIC DBLE, MAX, MIN, INT, SQRT INTRINSIC DBLE, MAX, MIN, INT, SQRT
@ -334,7 +334,8 @@
LIWORK = 1 LIWORK = 1
* *
* Iterate through all test cases and compute necessary workspace * Iterate through all test cases and compute necessary workspace
* sizes for ?GELS, ?GETSLS, ?GELSY, ?GELSS and ?GELSD routines. * sizes for ?GELS, ?GELST, ?GETSLS, ?GELSY, ?GELSS and ?GELSD
* routines.
* *
DO IM = 1, NM DO IM = 1, NM
M = MVAL( IM ) M = MVAL( IM )
@ -361,6 +362,10 @@
CALL ZGELS( TRANS, M, N, NRHS, A, LDA, CALL ZGELS( TRANS, M, N, NRHS, A, LDA,
$ B, LDB, WQ, -1, INFO ) $ B, LDB, WQ, -1, INFO )
LWORK_ZGELS = INT ( WQ( 1 ) ) LWORK_ZGELS = INT ( WQ( 1 ) )
* Compute workspace needed for ZGELST
CALL ZGELST( TRANS, M, N, NRHS, A, LDA,
$ B, LDB, WQ, -1, INFO )
LWORK_ZGELST = INT ( WQ ( 1 ) )
* Compute workspace needed for ZGETSLS * Compute workspace needed for ZGETSLS
CALL ZGETSLS( TRANS, M, N, NRHS, A, LDA, CALL ZGETSLS( TRANS, M, N, NRHS, A, LDA,
$ B, LDB, WQ, -1, INFO ) $ B, LDB, WQ, -1, INFO )
@ -390,9 +395,9 @@
LRWORK = MAX( LRWORK, LRWORK_ZGELSY, LRWORK = MAX( LRWORK, LRWORK_ZGELSY,
$ LRWORK_ZGELSS, LRWORK_ZGELSD ) $ LRWORK_ZGELSS, LRWORK_ZGELSD )
* Compute LWORK workspace needed for all functions * Compute LWORK workspace needed for all functions
LWORK = MAX( LWORK, LWORK_ZGELS, LWORK_ZGETSLS, LWORK = MAX( LWORK, LWORK_ZGELS, LWORK_ZGELST,
$ LWORK_ZGELSY, LWORK_ZGELSS, $ LWORK_ZGETSLS, LWORK_ZGELSY,
$ LWORK_ZGELSD ) $ LWORK_ZGELSS, LWORK_ZGELSD )
END IF END IF
ENDDO ENDDO
ENDDO ENDDO
@ -425,21 +430,26 @@
ITYPE = ( IRANK-1 )*3 + ISCALE ITYPE = ( IRANK-1 )*3 + ISCALE
IF( .NOT.DOTYPE( ITYPE ) ) IF( .NOT.DOTYPE( ITYPE ) )
$ GO TO 100 $ GO TO 100
* * =====================================================
* Begin test ZGELS
* =====================================================
IF( IRANK.EQ.1 ) THEN IF( IRANK.EQ.1 ) THEN
* *
* Test ZGELS
*
* Generate a matrix of scaling type ISCALE * Generate a matrix of scaling type ISCALE
* *
CALL ZQRT13( ISCALE, M, N, COPYA, LDA, NORMA, CALL ZQRT13( ISCALE, M, N, COPYA, LDA, NORMA,
$ ISEED ) $ ISEED )
DO 40 INB = 1, NNB *
* Loop for testing different block sizes.
*
DO INB = 1, NNB
NB = NBVAL( INB ) NB = NBVAL( INB )
CALL XLAENV( 1, NB ) CALL XLAENV( 1, NB )
CALL XLAENV( 3, NXVAL( INB ) ) CALL XLAENV( 3, NXVAL( INB ) )
* *
DO 30 ITRAN = 1, 2 * Loop for testing non-transposed and transposed.
*
DO ITRAN = 1, 2
IF( ITRAN.EQ.1 ) THEN IF( ITRAN.EQ.1 ) THEN
TRANS = 'N' TRANS = 'N'
NROWS = M NROWS = M
@ -484,15 +494,20 @@
$ ITYPE, NFAIL, NERRS, $ ITYPE, NFAIL, NERRS,
$ NOUT ) $ NOUT )
* *
* Check correctness of results * Test 1: Check correctness of results
* for ZGELS, compute the residual:
* RESID = norm(B - A*X) /
* / ( max(m,n) * norm(A) * norm(X) * EPS )
* *
LDWORK = MAX( 1, NROWS )
IF( NROWS.GT.0 .AND. NRHS.GT.0 ) IF( NROWS.GT.0 .AND. NRHS.GT.0 )
$ CALL ZLACPY( 'Full', NROWS, NRHS, $ CALL ZLACPY( 'Full', NROWS, NRHS,
$ COPYB, LDB, C, LDB ) $ COPYB, LDB, C, LDB )
CALL ZQRT16( TRANS, M, N, NRHS, COPYA, CALL ZQRT16( TRANS, M, N, NRHS, COPYA,
$ LDA, B, LDB, C, LDB, RWORK, $ LDA, B, LDB, C, LDB, RWORK,
$ RESULT( 1 ) ) $ RESULT( 1 ) )
*
* Test 2: Check correctness of results
* for ZGELS.
* *
IF( ( ITRAN.EQ.1 .AND. M.GE.N ) .OR. IF( ( ITRAN.EQ.1 .AND. M.GE.N ) .OR.
$ ( ITRAN.EQ.2 .AND. M.LT.N ) ) THEN $ ( ITRAN.EQ.2 .AND. M.LT.N ) ) THEN
@ -515,7 +530,7 @@
* Print information about the tests that * Print information about the tests that
* did not pass the threshold. * did not pass the threshold.
* *
DO 20 K = 1, 2 DO K = 1, 2
IF( RESULT( K ).GE.THRESH ) THEN IF( RESULT( K ).GE.THRESH ) THEN
IF( NFAIL.EQ.0 .AND. NERRS.EQ.0 ) IF( NFAIL.EQ.0 .AND. NERRS.EQ.0 )
$ CALL ALAHD( NOUT, PATH ) $ CALL ALAHD( NOUT, PATH )
@ -524,26 +539,34 @@
$ RESULT( K ) $ RESULT( K )
NFAIL = NFAIL + 1 NFAIL = NFAIL + 1
END IF END IF
20 CONTINUE END DO
NRUN = NRUN + 2 NRUN = NRUN + 2
30 CONTINUE END DO
40 CONTINUE END DO
* END IF
* * =====================================================
* Test ZGETSLS * End test ZGELS
* =====================================================
* =====================================================
* Begin test ZGELST
* =====================================================
IF( IRANK.EQ.1 ) THEN
* *
* Generate a matrix of scaling type ISCALE * Generate a matrix of scaling type ISCALE
* *
CALL ZQRT13( ISCALE, M, N, COPYA, LDA, NORMA, CALL ZQRT13( ISCALE, M, N, COPYA, LDA, NORMA,
$ ISEED ) $ ISEED )
DO 65 INB = 1, NNB
MB = NBVAL( INB )
CALL XLAENV( 1, MB )
DO 62 IMB = 1, NNB
NB = NBVAL( IMB )
CALL XLAENV( 2, NB )
* *
DO 60 ITRAN = 1, 2 * Loop for testing different block sizes.
*
DO INB = 1, NNB
NB = NBVAL( INB )
CALL XLAENV( 1, NB )
CALL XLAENV( 3, NXVAL( INB ) )
*
* Loop for testing non-transposed and transposed.
*
DO ITRAN = 1, 2
IF( ITRAN.EQ.1 ) THEN IF( ITRAN.EQ.1 ) THEN
TRANS = 'N' TRANS = 'N'
NROWS = M NROWS = M
@ -560,9 +583,9 @@
IF( NCOLS.GT.0 ) THEN IF( NCOLS.GT.0 ) THEN
CALL ZLARNV( 2, ISEED, NCOLS*NRHS, CALL ZLARNV( 2, ISEED, NCOLS*NRHS,
$ WORK ) $ WORK )
CALL ZSCAL( NCOLS*NRHS, CALL ZDSCAL( NCOLS*NRHS,
$ CONE / DBLE( NCOLS ), WORK, $ ONE / DBLE( NCOLS ), WORK,
$ 1 ) $ 1 )
END IF END IF
CALL ZGEMM( TRANS, 'No transpose', NROWS, CALL ZGEMM( TRANS, 'No transpose', NROWS,
$ NRHS, NCOLS, CONE, COPYA, LDA, $ NRHS, NCOLS, CONE, COPYA, LDA,
@ -578,31 +601,37 @@
CALL ZLACPY( 'Full', NROWS, NRHS, CALL ZLACPY( 'Full', NROWS, NRHS,
$ COPYB, LDB, B, LDB ) $ COPYB, LDB, B, LDB )
END IF END IF
SRNAMT = 'ZGETSLS ' SRNAMT = 'ZGELST'
CALL ZGETSLS( TRANS, M, N, NRHS, A, CALL ZGELST( TRANS, M, N, NRHS, A, LDA, B,
$ LDA, B, LDB, WORK, LWORK, INFO ) $ LDB, WORK, LWORK, INFO )
*
IF( INFO.NE.0 ) IF( INFO.NE.0 )
$ CALL ALAERH( PATH, 'ZGETSLS ', INFO, 0, $ CALL ALAERH( PATH, 'ZGELST', INFO, 0,
$ TRANS, M, N, NRHS, -1, NB, $ TRANS, M, N, NRHS, -1, NB,
$ ITYPE, NFAIL, NERRS, $ ITYPE, NFAIL, NERRS,
$ NOUT ) $ NOUT )
* *
* Check correctness of results * Test 3: Check correctness of results
* for ZGELST, compute the residual:
* RESID = norm(B - A*X) /
* / ( max(m,n) * norm(A) * norm(X) * EPS )
* *
LDWORK = MAX( 1, NROWS )
IF( NROWS.GT.0 .AND. NRHS.GT.0 ) IF( NROWS.GT.0 .AND. NRHS.GT.0 )
$ CALL ZLACPY( 'Full', NROWS, NRHS, $ CALL ZLACPY( 'Full', NROWS, NRHS,
$ COPYB, LDB, C, LDB ) $ COPYB, LDB, C, LDB )
CALL ZQRT16( TRANS, M, N, NRHS, COPYA, CALL ZQRT16( TRANS, M, N, NRHS, COPYA,
$ LDA, B, LDB, C, LDB, WORK2, $ LDA, B, LDB, C, LDB, RWORK,
$ RESULT( 15 ) ) $ RESULT( 3 ) )
*
* Test 4: Check correctness of results
* for ZGELST.
* *
IF( ( ITRAN.EQ.1 .AND. M.GE.N ) .OR. IF( ( ITRAN.EQ.1 .AND. M.GE.N ) .OR.
$ ( ITRAN.EQ.2 .AND. M.LT.N ) ) THEN $ ( ITRAN.EQ.2 .AND. M.LT.N ) ) THEN
* *
* Solving LS system * Solving LS system
* *
RESULT( 16 ) = ZQRT17( TRANS, 1, M, N, RESULT( 4 ) = ZQRT17( TRANS, 1, M, N,
$ NRHS, COPYA, LDA, B, LDB, $ NRHS, COPYA, LDA, B, LDB,
$ COPYB, LDB, C, WORK, $ COPYB, LDB, C, WORK,
$ LWORK ) $ LWORK )
@ -610,7 +639,7 @@
* *
* Solving overdetermined system * Solving overdetermined system
* *
RESULT( 16 ) = ZQRT14( TRANS, M, N, RESULT( 4 ) = ZQRT14( TRANS, M, N,
$ NRHS, COPYA, LDA, B, LDB, $ NRHS, COPYA, LDA, B, LDB,
$ WORK, LWORK ) $ WORK, LWORK )
END IF END IF
@ -618,21 +647,151 @@
* Print information about the tests that * Print information about the tests that
* did not pass the threshold. * did not pass the threshold.
* *
DO 50 K = 15, 16 DO K = 3, 4
IF( RESULT( K ).GE.THRESH ) THEN IF( RESULT( K ).GE.THRESH ) THEN
IF( NFAIL.EQ.0 .AND. NERRS.EQ.0 ) IF( NFAIL.EQ.0 .AND. NERRS.EQ.0 )
$ CALL ALAHD( NOUT, PATH ) $ CALL ALAHD( NOUT, PATH )
WRITE( NOUT, FMT = 9997 )TRANS, M, WRITE( NOUT, FMT = 9999 )TRANS, M,
$ N, NRHS, MB, NB, ITYPE, K, $ N, NRHS, NB, ITYPE, K,
$ RESULT( K ) $ RESULT( K )
NFAIL = NFAIL + 1 NFAIL = NFAIL + 1
END IF END IF
50 CONTINUE END DO
NRUN = NRUN + 2 NRUN = NRUN + 2
60 CONTINUE END DO
62 CONTINUE END DO
65 CONTINUE
END IF END IF
* =====================================================
* End test ZGELST
* =====================================================
* =====================================================
* Begin test ZGELSTSLS
* =====================================================
IF( IRANK.EQ.1 ) THEN
*
* Generate a matrix of scaling type ISCALE
*
CALL ZQRT13( ISCALE, M, N, COPYA, LDA, NORMA,
$ ISEED )
*
* Loop for testing different block sizes MB.
*
DO INB = 1, NNB
MB = NBVAL( INB )
CALL XLAENV( 1, MB )
*
* Loop for testing different block sizes NB.
*
DO IMB = 1, NNB
NB = NBVAL( IMB )
CALL XLAENV( 2, NB )
*
* Loop for testing non-transposed
* and transposed.
*
DO ITRAN = 1, 2
IF( ITRAN.EQ.1 ) THEN
TRANS = 'N'
NROWS = M
NCOLS = N
ELSE
TRANS = 'C'
NROWS = N
NCOLS = M
END IF
LDWORK = MAX( 1, NCOLS )
*
* Set up a consistent rhs
*
IF( NCOLS.GT.0 ) THEN
CALL ZLARNV( 2, ISEED, NCOLS*NRHS,
$ WORK )
CALL ZSCAL( NCOLS*NRHS,
$ CONE / DBLE( NCOLS ),
$ WORK, 1 )
END IF
CALL ZGEMM( TRANS, 'No transpose',
$ NROWS, NRHS, NCOLS, CONE,
$ COPYA, LDA, WORK, LDWORK,
$ CZERO, B, LDB )
CALL ZLACPY( 'Full', NROWS, NRHS,
$ B, LDB, COPYB, LDB )
*
* Solve LS or overdetermined system
*
IF( M.GT.0 .AND. N.GT.0 ) THEN
CALL ZLACPY( 'Full', M, N,
$ COPYA, LDA, A, LDA )
CALL ZLACPY( 'Full', NROWS, NRHS,
$ COPYB, LDB, B, LDB )
END IF
SRNAMT = 'ZGETSLS '
CALL ZGETSLS( TRANS, M, N, NRHS, A,
$ LDA, B, LDB, WORK, LWORK,
$ INFO )
IF( INFO.NE.0 )
$ CALL ALAERH( PATH, 'ZGETSLS ', INFO,
$ 0, TRANS, M, N, NRHS,
$ -1, NB, ITYPE, NFAIL,
$ NERRS, NOUT )
*
* Test 5: Check correctness of results
* for ZGETSLS, compute the residual:
* RESID = norm(B - A*X) /
* / ( max(m,n) * norm(A) * norm(X) * EPS )
*
IF( NROWS.GT.0 .AND. NRHS.GT.0 )
$ CALL ZLACPY( 'Full', NROWS, NRHS,
$ COPYB, LDB, C, LDB )
CALL ZQRT16( TRANS, M, N, NRHS,
$ COPYA, LDA, B, LDB,
$ C, LDB, WORK2,
$ RESULT( 5 ) )
*
* Test 6: Check correctness of results
* for ZGETSLS.
*
IF( ( ITRAN.EQ.1 .AND. M.GE.N ) .OR.
$ ( ITRAN.EQ.2 .AND. M.LT.N ) ) THEN
*
* Solving LS system, compute:
* r = norm((B- A*X)**T * A) /
* / (norm(A)*norm(B)*max(M,N,NRHS)*EPS)
*
RESULT( 6 ) = ZQRT17( TRANS, 1, M,
$ N, NRHS, COPYA, LDA,
$ B, LDB, COPYB, LDB,
$ C, WORK, LWORK )
ELSE
*
* Solving overdetermined system
*
RESULT( 6 ) = ZQRT14( TRANS, M, N,
$ NRHS, COPYA, LDA, B,
$ LDB, WORK, LWORK )
END IF
*
* Print information about the tests that
* did not pass the threshold.
*
DO K = 5, 6
IF( RESULT( K ).GE.THRESH ) THEN
IF( NFAIL.EQ.0 .AND. NERRS.EQ.0 )
$ CALL ALAHD( NOUT, PATH )
WRITE( NOUT, FMT = 9997 )TRANS,
$ M, N, NRHS, MB, NB, ITYPE, K,
$ RESULT( K )
NFAIL = NFAIL + 1
END IF
END DO
NRUN = NRUN + 2
END DO
END DO
END DO
END IF
* =====================================================
* End test ZGELSTSLS
* =====================================================
* *
* Generate a matrix of scaling type ISCALE and rank * Generate a matrix of scaling type ISCALE and rank
* type IRANK. * type IRANK.
@ -680,37 +839,37 @@
* *
* workspace used: 2*MNMIN+NB*NB+NB*MAX(N,NRHS) * workspace used: 2*MNMIN+NB*NB+NB*MAX(N,NRHS)
* *
* Test 3: Compute relative error in svd * Test 7: Compute relative error in svd
* workspace: M*N + 4*MIN(M,N) + MAX(M,N) * workspace: M*N + 4*MIN(M,N) + MAX(M,N)
* *
RESULT( 3 ) = ZQRT12( CRANK, CRANK, A, LDA, RESULT( 7 ) = ZQRT12( CRANK, CRANK, A, LDA,
$ COPYS, WORK, LWORK, RWORK ) $ COPYS, WORK, LWORK, RWORK )
* *
* Test 4: Compute error in solution * Test 8: Compute error in solution
* workspace: M*NRHS + M * workspace: M*NRHS + M
* *
CALL ZLACPY( 'Full', M, NRHS, COPYB, LDB, WORK, CALL ZLACPY( 'Full', M, NRHS, COPYB, LDB, WORK,
$ LDWORK ) $ LDWORK )
CALL ZQRT16( 'No transpose', M, N, NRHS, COPYA, CALL ZQRT16( 'No transpose', M, N, NRHS, COPYA,
$ LDA, B, LDB, WORK, LDWORK, RWORK, $ LDA, B, LDB, WORK, LDWORK, RWORK,
$ RESULT( 4 ) ) $ RESULT( 8 ) )
* *
* Test 5: Check norm of r'*A * Test 9: Check norm of r'*A
* workspace: NRHS*(M+N) * workspace: NRHS*(M+N)
* *
RESULT( 5 ) = ZERO RESULT( 9 ) = ZERO
IF( M.GT.CRANK ) IF( M.GT.CRANK )
$ RESULT( 5 ) = ZQRT17( 'No transpose', 1, M, $ RESULT( 9 ) = ZQRT17( 'No transpose', 1, M,
$ N, NRHS, COPYA, LDA, B, LDB, $ N, NRHS, COPYA, LDA, B, LDB,
$ COPYB, LDB, C, WORK, LWORK ) $ COPYB, LDB, C, WORK, LWORK )
* *
* Test 6: Check if x is in the rowspace of A * Test 10: Check if x is in the rowspace of A
* workspace: (M+NRHS)*(N+2) * workspace: (M+NRHS)*(N+2)
* *
RESULT( 6 ) = ZERO RESULT( 10 ) = ZERO
* *
IF( N.GT.CRANK ) IF( N.GT.CRANK )
$ RESULT( 6 ) = ZQRT14( 'No transpose', M, N, $ RESULT( 10 ) = ZQRT14( 'No transpose', M, N,
$ NRHS, COPYA, LDA, B, LDB, $ NRHS, COPYA, LDA, B, LDB,
$ WORK, LWORK ) $ WORK, LWORK )
* *
@ -736,38 +895,38 @@
* workspace used: 3*min(m,n) + * workspace used: 3*min(m,n) +
* max(2*min(m,n),nrhs,max(m,n)) * max(2*min(m,n),nrhs,max(m,n))
* *
* Test 7: Compute relative error in svd * Test 11: Compute relative error in svd
* *
IF( RANK.GT.0 ) THEN IF( RANK.GT.0 ) THEN
CALL DAXPY( MNMIN, -ONE, COPYS, 1, S, 1 ) CALL DAXPY( MNMIN, -ONE, COPYS, 1, S, 1 )
RESULT( 7 ) = DASUM( MNMIN, S, 1 ) / RESULT( 11 ) = DASUM( MNMIN, S, 1 ) /
$ DASUM( MNMIN, COPYS, 1 ) / $ DASUM( MNMIN, COPYS, 1 ) /
$ ( EPS*DBLE( MNMIN ) ) $ ( EPS*DBLE( MNMIN ) )
ELSE ELSE
RESULT( 7 ) = ZERO RESULT( 11 ) = ZERO
END IF END IF
* *
* Test 8: Compute error in solution * Test 12: Compute error in solution
* *
CALL ZLACPY( 'Full', M, NRHS, COPYB, LDB, WORK, CALL ZLACPY( 'Full', M, NRHS, COPYB, LDB, WORK,
$ LDWORK ) $ LDWORK )
CALL ZQRT16( 'No transpose', M, N, NRHS, COPYA, CALL ZQRT16( 'No transpose', M, N, NRHS, COPYA,
$ LDA, B, LDB, WORK, LDWORK, RWORK, $ LDA, B, LDB, WORK, LDWORK, RWORK,
$ RESULT( 8 ) ) $ RESULT( 12 ) )
* *
* Test 9: Check norm of r'*A * Test 13: Check norm of r'*A
* *
RESULT( 9 ) = ZERO RESULT( 13 ) = ZERO
IF( M.GT.CRANK ) IF( M.GT.CRANK )
$ RESULT( 9 ) = ZQRT17( 'No transpose', 1, M, $ RESULT( 13 ) = ZQRT17( 'No transpose', 1, M,
$ N, NRHS, COPYA, LDA, B, LDB, $ N, NRHS, COPYA, LDA, B, LDB,
$ COPYB, LDB, C, WORK, LWORK ) $ COPYB, LDB, C, WORK, LWORK )
* *
* Test 10: Check if x is in the rowspace of A * Test 14: Check if x is in the rowspace of A
* *
RESULT( 10 ) = ZERO RESULT( 14 ) = ZERO
IF( N.GT.CRANK ) IF( N.GT.CRANK )
$ RESULT( 10 ) = ZQRT14( 'No transpose', M, N, $ RESULT( 14 ) = ZQRT14( 'No transpose', M, N,
$ NRHS, COPYA, LDA, B, LDB, $ NRHS, COPYA, LDA, B, LDB,
$ WORK, LWORK ) $ WORK, LWORK )
* *
@ -792,45 +951,45 @@
$ N, NRHS, -1, NB, ITYPE, NFAIL, $ N, NRHS, -1, NB, ITYPE, NFAIL,
$ NERRS, NOUT ) $ NERRS, NOUT )
* *
* Test 11: Compute relative error in svd * Test 15: Compute relative error in svd
* *
IF( RANK.GT.0 ) THEN IF( RANK.GT.0 ) THEN
CALL DAXPY( MNMIN, -ONE, COPYS, 1, S, 1 ) CALL DAXPY( MNMIN, -ONE, COPYS, 1, S, 1 )
RESULT( 11 ) = DASUM( MNMIN, S, 1 ) / RESULT( 15 ) = DASUM( MNMIN, S, 1 ) /
$ DASUM( MNMIN, COPYS, 1 ) / $ DASUM( MNMIN, COPYS, 1 ) /
$ ( EPS*DBLE( MNMIN ) ) $ ( EPS*DBLE( MNMIN ) )
ELSE ELSE
RESULT( 11 ) = ZERO RESULT( 15 ) = ZERO
END IF END IF
* *
* Test 12: Compute error in solution * Test 16: Compute error in solution
* *
CALL ZLACPY( 'Full', M, NRHS, COPYB, LDB, WORK, CALL ZLACPY( 'Full', M, NRHS, COPYB, LDB, WORK,
$ LDWORK ) $ LDWORK )
CALL ZQRT16( 'No transpose', M, N, NRHS, COPYA, CALL ZQRT16( 'No transpose', M, N, NRHS, COPYA,
$ LDA, B, LDB, WORK, LDWORK, RWORK, $ LDA, B, LDB, WORK, LDWORK, RWORK,
$ RESULT( 12 ) ) $ RESULT( 16 ) )
* *
* Test 13: Check norm of r'*A * Test 17: Check norm of r'*A
* *
RESULT( 13 ) = ZERO RESULT( 17 ) = ZERO
IF( M.GT.CRANK ) IF( M.GT.CRANK )
$ RESULT( 13 ) = ZQRT17( 'No transpose', 1, M, $ RESULT( 17 ) = ZQRT17( 'No transpose', 1, M,
$ N, NRHS, COPYA, LDA, B, LDB, $ N, NRHS, COPYA, LDA, B, LDB,
$ COPYB, LDB, C, WORK, LWORK ) $ COPYB, LDB, C, WORK, LWORK )
* *
* Test 14: Check if x is in the rowspace of A * Test 18: Check if x is in the rowspace of A
* *
RESULT( 14 ) = ZERO RESULT( 18 ) = ZERO
IF( N.GT.CRANK ) IF( N.GT.CRANK )
$ RESULT( 14 ) = ZQRT14( 'No transpose', M, N, $ RESULT( 18 ) = ZQRT14( 'No transpose', M, N,
$ NRHS, COPYA, LDA, B, LDB, $ NRHS, COPYA, LDA, B, LDB,
$ WORK, LWORK ) $ WORK, LWORK )
* *
* Print information about the tests that did not * Print information about the tests that did not
* pass the threshold. * pass the threshold.
* *
DO 80 K = 3, 14 DO 80 K = 7, 18
IF( RESULT( K ).GE.THRESH ) THEN IF( RESULT( K ).GE.THRESH ) THEN
IF( NFAIL.EQ.0 .AND. NERRS.EQ.0 ) IF( NFAIL.EQ.0 .AND. NERRS.EQ.0 )
$ CALL ALAHD( NOUT, PATH ) $ CALL ALAHD( NOUT, PATH )

61
lapack-netlib/TESTING/LIN/zerrls.f
View File

@ -22,7 +22,7 @@
*> \verbatim *> \verbatim
*> *>
*> ZERRLS tests the error exits for the COMPLEX*16 least squares *> ZERRLS tests the error exits for the COMPLEX*16 least squares
*> driver routines (ZGELS, CGELSS, CGELSY, CGELSD). *> driver routines (ZGELS, ZGELST, ZGETSLS, CGELSS, CGELSY, CGELSD).
*> \endverbatim *> \endverbatim
* *
* Arguments: * Arguments:
@ -83,7 +83,8 @@
EXTERNAL LSAMEN EXTERNAL LSAMEN
* .. * ..
* .. External Subroutines .. * .. External Subroutines ..
EXTERNAL ALAESM, CHKXER, ZGELS, ZGELSD, ZGELSS, ZGELSY EXTERNAL ALAESM, CHKXER, ZGELS, ZGELSD, ZGELSS, ZGELST,
$ ZGELSY, ZGETSLS
* .. * ..
* .. Scalars in Common .. * .. Scalars in Common ..
LOGICAL LERR, OK LOGICAL LERR, OK
@ -130,10 +131,66 @@
INFOT = 8 INFOT = 8
CALL ZGELS( 'N', 2, 0, 0, A, 2, B, 1, W, 2, INFO ) CALL ZGELS( 'N', 2, 0, 0, A, 2, B, 1, W, 2, INFO )
CALL CHKXER( 'ZGELS ', INFOT, NOUT, LERR, OK ) CALL CHKXER( 'ZGELS ', INFOT, NOUT, LERR, OK )
INFOT = 8
CALL ZGELS( 'N', 0, 2, 0, A, 1, B, 1, W, 2, INFO )
CALL CHKXER( 'ZGELS', INFOT, NOUT, LERR, OK )
INFOT = 10 INFOT = 10
CALL ZGELS( 'N', 1, 1, 0, A, 1, B, 1, W, 1, INFO ) CALL ZGELS( 'N', 1, 1, 0, A, 1, B, 1, W, 1, INFO )
CALL CHKXER( 'ZGELS ', INFOT, NOUT, LERR, OK ) CALL CHKXER( 'ZGELS ', INFOT, NOUT, LERR, OK )
* *
* ZGELST
*
SRNAMT = 'ZGELST'
INFOT = 1
CALL ZGELST( '/', 0, 0, 0, A, 1, B, 1, W, 1, INFO )
CALL CHKXER( 'ZGELST', INFOT, NOUT, LERR, OK )
INFOT = 2
CALL ZGELST( 'N', -1, 0, 0, A, 1, B, 1, W, 1, INFO )
CALL CHKXER( 'ZGELST', INFOT, NOUT, LERR, OK )
INFOT = 3
CALL ZGELST( 'N', 0, -1, 0, A, 1, B, 1, W, 1, INFO )
CALL CHKXER( 'ZGELST', INFOT, NOUT, LERR, OK )
INFOT = 4
CALL ZGELST( 'N', 0, 0, -1, A, 1, B, 1, W, 1, INFO )
CALL CHKXER( 'ZGELST', INFOT, NOUT, LERR, OK )
INFOT = 6
CALL ZGELST( 'N', 2, 0, 0, A, 1, B, 2, W, 2, INFO )
CALL CHKXER( 'ZGELST', INFOT, NOUT, LERR, OK )
INFOT = 8
CALL ZGELST( 'N', 2, 0, 0, A, 2, B, 1, W, 2, INFO )
CALL CHKXER( 'ZGELST', INFOT, NOUT, LERR, OK )
INFOT = 8
CALL ZGELST( 'N', 0, 2, 0, A, 1, B, 1, W, 2, INFO )
CALL CHKXER( 'ZGELST', INFOT, NOUT, LERR, OK )
INFOT = 10
CALL ZGELST( 'N', 1, 1, 0, A, 1, B, 1, W, 1, INFO )
CALL CHKXER( 'ZGELST', INFOT, NOUT, LERR, OK )
*
* ZGETSLS
*
SRNAMT = 'ZGETSLS'
INFOT = 1
CALL ZGETSLS( '/', 0, 0, 0, A, 1, B, 1, W, 1, INFO )
CALL CHKXER( 'ZGETSLS', INFOT, NOUT, LERR, OK )
INFOT = 2
CALL ZGETSLS( 'N', -1, 0, 0, A, 1, B, 1, W, 1, INFO )
CALL CHKXER( 'ZGETSLS', INFOT, NOUT, LERR, OK )
INFOT = 3
CALL ZGETSLS( 'N', 0, -1, 0, A, 1, B, 1, W, 1, INFO )
CALL CHKXER( 'ZGETSLS', INFOT, NOUT, LERR, OK )
INFOT = 4
CALL ZGETSLS( 'N', 0, 0, -1, A, 1, B, 1, W, 1, INFO )
CALL CHKXER( 'ZGETSLS', INFOT, NOUT, LERR, OK )
INFOT = 6
CALL ZGETSLS( 'N', 2, 0, 0, A, 1, B, 2, W, 2, INFO )
CALL CHKXER( 'ZGETSLS', INFOT, NOUT, LERR, OK )
INFOT = 8
CALL ZGETSLS( 'N', 2, 0, 0, A, 2, B, 1, W, 2, INFO )
CALL CHKXER( 'ZGETSLS', INFOT, NOUT, LERR, OK )
INFOT = 8
CALL ZGETSLS( 'N', 0, 2, 0, A, 1, B, 1, W, 2, INFO )
CALL CHKXER( 'ZGETSLS', INFOT, NOUT, LERR, OK )
*
* ZGELSS * ZGELSS
* *
SRNAMT = 'ZGELSS' SRNAMT = 'ZGELSS'