213 INTEGER nm, nn, nnb, nns, nout
214 DOUBLE PRECISION thresh
218 INTEGER iwork( * ), mval( * ), nbval( * ), nsval( * ),
219 $ nval( * ), nxval( * )
220 DOUBLE PRECISION a( * ), b( * ), c( * ), copya( * ), copyb( * ),
221 $ copys( * ), s( * ), work( * )
228 parameter( ntests = 18 )
230 parameter( smlsiz = 25 )
231 DOUBLE PRECISION one, two, zero
232 parameter( one = 1.0d0, two = 2.0d0, zero = 0.0d0 )
237 INTEGER crank, i, im, in, inb, info, ins, irank,
238 $ iscale, itran, itype, j, k, lda, ldb, ldwork,
239 $ lwlsy, lwork, m, mnmin, n, nb, ncols, nerrs,
240 $ nfail, nlvl, nrhs, nrows, nrun, rank
241 DOUBLE PRECISION eps, norma, normb, rcond
244 INTEGER iseed( 4 ), iseedy( 4 )
245 DOUBLE PRECISION result( ntests )
258 INTRINSIC dble, int, log, max, min, sqrt
263 INTEGER infot, iounit
266 COMMON / infoc / infot, iounit, ok, lerr
267 COMMON / srnamc / srnamt
270 DATA iseedy / 1988, 1989, 1990, 1991 /
276 path( 1: 1 ) =
'Double precision'
282 iseed( i ) = iseedy( i )
288 rcond = sqrt( eps ) - ( sqrt( eps )-eps ) / 2
295 $
CALL derrls( path, nout )
299 IF( ( nm.EQ.0 .OR. nn.EQ.0 ) .AND. thresh.EQ.zero )
300 $
CALL alahd( nout, path )
316 nlvl = max( int( log( max( one, dble( mnmin ) ) /
317 $ dble( smlsiz+1 ) ) / log( two ) ) + 1, 0 )
318 lwork = max( 1, ( m+nrhs )*( n+2 ), ( n+nrhs )*( m+2 ),
319 $ m*n+4*mnmin+max( m, n ), 12*mnmin+2*mnmin*smlsiz+
320 $ 8*mnmin*nlvl+mnmin*nrhs+(smlsiz+1)**2 )
324 itype = ( irank-1 )*3 + iscale
325 IF( .NOT.dotype( itype ) )
328 IF( irank.EQ.1 )
THEN
334 CALL dqrt13( iscale, m, n, copya, lda, norma,
339 CALL xlaenv( 3, nxval( inb ) )
342 IF( itran.EQ.1 )
THEN
351 ldwork = max( 1, ncols )
355 IF( ncols.GT.0 )
THEN
356 CALL dlarnv( 2, iseed, ncols*nrhs,
358 CALL dscal( ncols*nrhs,
359 $ one / dble( ncols ), work,
362 CALL dgemm( trans,
'No transpose', nrows,
363 $ nrhs, ncols, one, copya, lda,
364 $ work, ldwork, zero, b, ldb )
365 CALL dlacpy(
'Full', nrows, nrhs, b, ldb,
370 IF( m.GT.0 .AND. n.GT.0 )
THEN
371 CALL dlacpy(
'Full', m, n, copya, lda,
373 CALL dlacpy(
'Full', nrows, nrhs,
374 $ copyb, ldb, b, ldb )
377 CALL dgels( trans, m, n, nrhs, a, lda, b,
378 $ ldb, work, lwork, info )
380 $
CALL alaerh( path,
'DGELS ', info, 0,
381 $ trans, m, n, nrhs, -1, nb,
382 $ itype, nfail, nerrs,
387 ldwork = max( 1, nrows )
388 IF( nrows.GT.0 .AND. nrhs.GT.0 )
389 $
CALL dlacpy(
'Full', nrows, nrhs,
390 $ copyb, ldb, c, ldb )
391 CALL dqrt16( trans, m, n, nrhs, copya,
392 $ lda, b, ldb, c, ldb, work,
395 IF( ( itran.EQ.1 .AND. m.GE.n ) .OR.
396 $ ( itran.EQ.2 .AND. m.LT.n ) )
THEN
400 result( 2 ) =
dqrt17( trans, 1, m, n,
401 $ nrhs, copya, lda, b, ldb,
402 $ copyb, ldb, c, work,
408 result( 2 ) =
dqrt14( trans, m, n,
409 $ nrhs, copya, lda, b, ldb,
417 IF( result( k ).GE.thresh )
THEN
418 IF( nfail.EQ.0 .AND. nerrs.EQ.0 )
419 $
CALL alahd( nout, path )
420 WRITE( nout, fmt = 9999 )trans, m,
421 $ n, nrhs, nb, itype, k,
434 CALL dqrt15( iscale, irank, m, n, nrhs, copya, lda,
435 $ copyb, ldb, copys, rank, norma, normb,
436 $ iseed, work, lwork )
453 CALL dlacpy(
'Full', m, n, copya, lda, a, lda )
454 CALL dlacpy(
'Full', m, nrhs, copyb, ldb, b, ldb )
457 CALL dgelsx( m, n, nrhs, a, lda, b, ldb, iwork,
458 $ rcond, crank, work, info )
460 $
CALL alaerh( path,
'DGELSX', info, 0,
' ', m, n,
461 $ nrhs, -1, nb, itype, nfail, nerrs,
469 result( 3 ) =
dqrt12( crank, crank, a, lda, copys,
475 CALL dlacpy(
'Full', m, nrhs, copyb, ldb, work,
477 CALL dqrt16(
'No transpose', m, n, nrhs, copya,
478 $ lda, b, ldb, work, ldwork,
479 $ work( m*nrhs+1 ), result( 4 ) )
486 $ result( 5 ) =
dqrt17(
'No transpose', 1, m, n,
487 $ nrhs, copya, lda, b, ldb, copyb,
488 $ ldb, c, work, lwork )
496 $ result( 6 ) =
dqrt14(
'No transpose', m, n,
497 $ nrhs, copya, lda, b, ldb, work,
504 IF( result( k ).GE.thresh )
THEN
505 IF( nfail.EQ.0 .AND. nerrs.EQ.0 )
506 $
CALL alahd( nout, path )
507 WRITE( nout, fmt = 9998 )m, n, nrhs, nb,
508 $ itype, k, result( k )
519 CALL xlaenv( 3, nxval( inb ) )
536 lwlsy = max( 1, mnmin+2*n+nb*( n+1 ),
539 CALL dlacpy(
'Full', m, n, copya, lda, a, lda )
540 CALL dlacpy(
'Full', m, nrhs, copyb, ldb, b,
544 CALL dgelsy( m, n, nrhs, a, lda, b, ldb, iwork,
545 $ rcond, crank, work, lwlsy, info )
547 $
CALL alaerh( path,
'DGELSY', info, 0,
' ', m,
548 $ n, nrhs, -1, nb, itype, nfail,
554 result( 7 ) =
dqrt12( crank, crank, a, lda,
555 $ copys, work, lwork )
560 CALL dlacpy(
'Full', m, nrhs, copyb, ldb, work,
562 CALL dqrt16(
'No transpose', m, n, nrhs, copya,
563 $ lda, b, ldb, work, ldwork,
564 $ work( m*nrhs+1 ), result( 8 ) )
571 $ result( 9 ) =
dqrt17(
'No transpose', 1, m,
572 $ n, nrhs, copya, lda, b, ldb,
573 $ copyb, ldb, c, work, lwork )
581 $ result( 10 ) =
dqrt14(
'No transpose', m, n,
582 $ nrhs, copya, lda, b, ldb,
591 CALL dlacpy(
'Full', m, n, copya, lda, a, lda )
592 CALL dlacpy(
'Full', m, nrhs, copyb, ldb, b,
595 CALL dgelss( m, n, nrhs, a, lda, b, ldb, s,
596 $ rcond, crank, work, lwork, info )
598 $
CALL alaerh( path,
'DGELSS', info, 0,
' ', m,
599 $ n, nrhs, -1, nb, itype, nfail,
608 CALL daxpy( mnmin, -one, copys, 1, s, 1 )
609 result( 11 ) =
dasum( mnmin, s, 1 ) /
610 $
dasum( mnmin, copys, 1 ) /
611 $ ( eps*dble( mnmin ) )
618 CALL dlacpy(
'Full', m, nrhs, copyb, ldb, work,
620 CALL dqrt16(
'No transpose', m, n, nrhs, copya,
621 $ lda, b, ldb, work, ldwork,
622 $ work( m*nrhs+1 ), result( 12 ) )
628 $ result( 13 ) =
dqrt17(
'No transpose', 1, m,
629 $ n, nrhs, copya, lda, b, ldb,
630 $ copyb, ldb, c, work, lwork )
636 $ result( 14 ) =
dqrt14(
'No transpose', m, n,
637 $ nrhs, copya, lda, b, ldb,
652 CALL dlacpy(
'Full', m, n, copya, lda, a, lda )
653 CALL dlacpy(
'Full', m, nrhs, copyb, ldb, b,
657 CALL dgelsd( m, n, nrhs, a, lda, b, ldb, s,
658 $ rcond, crank, work, lwork, iwork,
661 $
CALL alaerh( path,
'DGELSD', info, 0,
' ', m,
662 $ n, nrhs, -1, nb, itype, nfail,
668 CALL daxpy( mnmin, -one, copys, 1, s, 1 )
669 result( 15 ) =
dasum( mnmin, s, 1 ) /
670 $
dasum( mnmin, copys, 1 ) /
671 $ ( eps*dble( mnmin ) )
678 CALL dlacpy(
'Full', m, nrhs, copyb, ldb, work,
680 CALL dqrt16(
'No transpose', m, n, nrhs, copya,
681 $ lda, b, ldb, work, ldwork,
682 $ work( m*nrhs+1 ), result( 16 ) )
688 $ result( 17 ) =
dqrt17(
'No transpose', 1, m,
689 $ n, nrhs, copya, lda, b, ldb,
690 $ copyb, ldb, c, work, lwork )
696 $ result( 18 ) =
dqrt14(
'No transpose', m, n,
697 $ nrhs, copya, lda, b, ldb,
704 IF( result( k ).GE.thresh )
THEN
705 IF( nfail.EQ.0 .AND. nerrs.EQ.0 )
706 $
CALL alahd( nout, path )
707 WRITE( nout, fmt = 9998 )m, n, nrhs, nb,
708 $ itype, k, result( k )
723 CALL alasvm( path, nout, nfail, nrun, nerrs )
725 9999
FORMAT(
' TRANS=''', a1,
''', M=', i5,
', N=', i5,
', NRHS=', i4,
726 $
', NB=', i4,
', type', i2,
', test(', i2,
')=', g12.5 )
727 9998
FORMAT(
' M=', i5,
', N=', i5,
', NRHS=', i4,
', NB=', i4,
728 $
', type', i2,
', test(', i2,
')=', g12.5 )
subroutine dlasrt(ID, N, D, INFO)
DLASRT sorts numbers in increasing or decreasing order.
subroutine alahd(IOUNIT, PATH)
ALAHD
double precision function dqrt17(TRANS, IRESID, M, N, NRHS, A, LDA, X, LDX, B, LDB, C, WORK, LWORK)
DQRT17
subroutine alasvm(TYPE, NOUT, NFAIL, NRUN, NERRS)
ALASVM
double precision function dqrt12(M, N, A, LDA, S, WORK, LWORK)
DQRT12
subroutine dqrt16(TRANS, M, N, NRHS, A, LDA, X, LDX, B, LDB, RWORK, RESID)
DQRT16
subroutine dgemm(TRANSA, TRANSB, M, N, K, ALPHA, A, LDA, B, LDB, BETA, C, LDC)
DGEMM
subroutine dgelsx(M, N, NRHS, A, LDA, B, LDB, JPVT, RCOND, RANK, WORK, INFO)
DGELSX solves overdetermined or underdetermined systems for GE matrices
subroutine dgels(TRANS, M, N, NRHS, A, LDA, B, LDB, WORK, LWORK, INFO)
DGELS solves overdetermined or underdetermined systems for GE matrices
subroutine alaerh(PATH, SUBNAM, INFO, INFOE, OPTS, M, N, KL, KU, N5, IMAT, NFAIL, NERRS, NOUT)
ALAERH
subroutine dlarnv(IDIST, ISEED, N, X)
DLARNV returns a vector of random numbers from a uniform or normal distribution.
subroutine dgelsy(M, N, NRHS, A, LDA, B, LDB, JPVT, RCOND, RANK, WORK, LWORK, INFO)
DGELSY solves overdetermined or underdetermined systems for GE matrices
subroutine dqrt13(SCALE, M, N, A, LDA, NORMA, ISEED)
DQRT13
subroutine dscal(N, DA, DX, INCX)
DSCAL
subroutine dgelss(M, N, NRHS, A, LDA, B, LDB, S, RCOND, RANK, WORK, LWORK, INFO)
DGELSS solves overdetermined or underdetermined systems for GE matrices
subroutine derrls(PATH, NUNIT)
DERRLS
subroutine xlaenv(ISPEC, NVALUE)
XLAENV
double precision function dlamch(CMACH)
DLAMCH
subroutine daxpy(N, DA, DX, INCX, DY, INCY)
DAXPY
subroutine dlacpy(UPLO, M, N, A, LDA, B, LDB)
DLACPY copies all or part of one two-dimensional array to another.
double precision function dqrt14(TRANS, M, N, NRHS, A, LDA, X, LDX, WORK, LWORK)
DQRT14
double precision function dasum(N, DX, INCX)
DASUM
subroutine dgelsd(M, N, NRHS, A, LDA, B, LDB, S, RCOND, RANK, WORK, LWORK, IWORK, INFO)
DGELSD computes the minimum-norm solution to a linear least squares problem for GE matrices ...
subroutine dqrt15(SCALE, RKSEL, M, N, NRHS, A, LDA, B, LDB, S, RANK, NORMA, NORMB, ISEED, WORK, LWORK)
DQRT15