1*> \brief \b DDRVLS 2* 3* =========== DOCUMENTATION =========== 4* 5* Online html documentation available at 6* http://www.netlib.org/lapack/explore-html/ 7* 8* Definition: 9* =========== 10* 11* SUBROUTINE DDRVLS( DOTYPE, NM, MVAL, NN, NVAL, NNS, NSVAL, NNB, 12* NBVAL, NXVAL, THRESH, TSTERR, A, COPYA, B, 13* COPYB, C, S, COPYS, NOUT ) 14* 15* .. Scalar Arguments .. 16* LOGICAL TSTERR 17* INTEGER NM, NN, NNB, NNS, NOUT 18* DOUBLE PRECISION THRESH 19* .. 20* .. Array Arguments .. 21* LOGICAL DOTYPE( * ) 22* INTEGER MVAL( * ), NBVAL( * ), NSVAL( * ), 23* $ NVAL( * ), NXVAL( * ) 24* DOUBLE PRECISION A( * ), B( * ), C( * ), COPYA( * ), COPYB( * ), 25* $ COPYS( * ), S( * ) 26* .. 27* 28* 29*> \par Purpose: 30* ============= 31*> 32*> \verbatim 33*> 34*> DDRVLS tests the least squares driver routines DGELS, DGETSLS, DGELSS, DGELSY, 35*> and DGELSD. 36*> \endverbatim 37* 38* Arguments: 39* ========== 40* 41*> \param[in] DOTYPE 42*> \verbatim 43*> DOTYPE is LOGICAL array, dimension (NTYPES) 44*> The matrix types to be used for testing. Matrices of type j 45*> (for 1 <= j <= NTYPES) are used for testing if DOTYPE(j) = 46*> .TRUE.; if DOTYPE(j) = .FALSE., then type j is not used. 47*> The matrix of type j is generated as follows: 48*> j=1: A = U*D*V where U and V are random orthogonal matrices 49*> and D has random entries (> 0.1) taken from a uniform 50*> distribution (0,1). A is full rank. 51*> j=2: The same of 1, but A is scaled up. 52*> j=3: The same of 1, but A is scaled down. 53*> j=4: A = U*D*V where U and V are random orthogonal matrices 54*> and D has 3*min(M,N)/4 random entries (> 0.1) taken 55*> from a uniform distribution (0,1) and the remaining 56*> entries set to 0. A is rank-deficient. 57*> j=5: The same of 4, but A is scaled up. 58*> j=6: The same of 5, but A is scaled down. 59*> \endverbatim 60*> 61*> \param[in] NM 62*> \verbatim 63*> NM is INTEGER 64*> The number of values of M contained in the vector MVAL. 65*> \endverbatim 66*> 67*> \param[in] MVAL 68*> \verbatim 69*> MVAL is INTEGER array, dimension (NM) 70*> The values of the matrix row dimension M. 71*> \endverbatim 72*> 73*> \param[in] NN 74*> \verbatim 75*> NN is INTEGER 76*> The number of values of N contained in the vector NVAL. 77*> \endverbatim 78*> 79*> \param[in] NVAL 80*> \verbatim 81*> NVAL is INTEGER array, dimension (NN) 82*> The values of the matrix column dimension N. 83*> \endverbatim 84*> 85*> \param[in] NNS 86*> \verbatim 87*> NNS is INTEGER 88*> The number of values of NRHS contained in the vector NSVAL. 89*> \endverbatim 90*> 91*> \param[in] NSVAL 92*> \verbatim 93*> NSVAL is INTEGER array, dimension (NNS) 94*> The values of the number of right hand sides NRHS. 95*> \endverbatim 96*> 97*> \param[in] NNB 98*> \verbatim 99*> NNB is INTEGER 100*> The number of values of NB and NX contained in the 101*> vectors NBVAL and NXVAL. The blocking parameters are used 102*> in pairs (NB,NX). 103*> \endverbatim 104*> 105*> \param[in] NBVAL 106*> \verbatim 107*> NBVAL is INTEGER array, dimension (NNB) 108*> The values of the blocksize NB. 109*> \endverbatim 110*> 111*> \param[in] NXVAL 112*> \verbatim 113*> NXVAL is INTEGER array, dimension (NNB) 114*> The values of the crossover point NX. 115*> \endverbatim 116*> 117*> \param[in] THRESH 118*> \verbatim 119*> THRESH is DOUBLE PRECISION 120*> The threshold value for the test ratios. A result is 121*> included in the output file if RESULT >= THRESH. To have 122*> every test ratio printed, use THRESH = 0. 123*> \endverbatim 124*> 125*> \param[in] TSTERR 126*> \verbatim 127*> TSTERR is LOGICAL 128*> Flag that indicates whether error exits are to be tested. 129*> \endverbatim 130*> 131*> \param[out] A 132*> \verbatim 133*> A is DOUBLE PRECISION array, dimension (MMAX*NMAX) 134*> where MMAX is the maximum value of M in MVAL and NMAX is the 135*> maximum value of N in NVAL. 136*> \endverbatim 137*> 138*> \param[out] COPYA 139*> \verbatim 140*> COPYA is DOUBLE PRECISION array, dimension (MMAX*NMAX) 141*> \endverbatim 142*> 143*> \param[out] B 144*> \verbatim 145*> B is DOUBLE PRECISION array, dimension (MMAX*NSMAX) 146*> where MMAX is the maximum value of M in MVAL and NSMAX is the 147*> maximum value of NRHS in NSVAL. 148*> \endverbatim 149*> 150*> \param[out] COPYB 151*> \verbatim 152*> COPYB is DOUBLE PRECISION array, dimension (MMAX*NSMAX) 153*> \endverbatim 154*> 155*> \param[out] C 156*> \verbatim 157*> C is DOUBLE PRECISION array, dimension (MMAX*NSMAX) 158*> \endverbatim 159*> 160*> \param[out] S 161*> \verbatim 162*> S is DOUBLE PRECISION array, dimension 163*> (min(MMAX,NMAX)) 164*> \endverbatim 165*> 166*> \param[out] COPYS 167*> \verbatim 168*> COPYS is DOUBLE PRECISION array, dimension 169*> (min(MMAX,NMAX)) 170*> \endverbatim 171*> 172*> \param[in] NOUT 173*> \verbatim 174*> NOUT is INTEGER 175*> The unit number for output. 176*> \endverbatim 177* 178* Authors: 179* ======== 180* 181*> \author Univ. of Tennessee 182*> \author Univ. of California Berkeley 183*> \author Univ. of Colorado Denver 184*> \author NAG Ltd. 185* 186*> \date June 2017 187* 188*> \ingroup double_lin 189* 190* ===================================================================== 191 SUBROUTINE DDRVLS( DOTYPE, NM, MVAL, NN, NVAL, NNS, NSVAL, NNB, 192 $ NBVAL, NXVAL, THRESH, TSTERR, A, COPYA, B, 193 $ COPYB, C, S, COPYS, NOUT ) 194* 195* -- LAPACK test routine (version 3.7.1) -- 196* -- LAPACK is a software package provided by Univ. of Tennessee, -- 197* -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..-- 198* June 2017 199* 200* .. Scalar Arguments .. 201 LOGICAL TSTERR 202 INTEGER NM, NN, NNB, NNS, NOUT 203 DOUBLE PRECISION THRESH 204* .. 205* .. Array Arguments .. 206 LOGICAL DOTYPE( * ) 207 INTEGER MVAL( * ), NBVAL( * ), NSVAL( * ), 208 $ NVAL( * ), NXVAL( * ) 209 DOUBLE PRECISION A( * ), B( * ), C( * ), COPYA( * ), COPYB( * ), 210 $ COPYS( * ), S( * ) 211* .. 212* 213* ===================================================================== 214* 215* .. Parameters .. 216 INTEGER NTESTS 217 PARAMETER ( NTESTS = 16 ) 218 INTEGER SMLSIZ 219 PARAMETER ( SMLSIZ = 25 ) 220 DOUBLE PRECISION ONE, TWO, ZERO 221 PARAMETER ( ONE = 1.0D0, TWO = 2.0D0, ZERO = 0.0D0 ) 222* .. 223* .. Local Scalars .. 224 CHARACTER TRANS 225 CHARACTER*3 PATH 226 INTEGER CRANK, I, IM, IMB, IN, INB, INFO, INS, IRANK, 227 $ ISCALE, ITRAN, ITYPE, J, K, LDA, LDB, LDWORK, 228 $ LWLSY, LWORK, M, MNMIN, N, NB, NCOLS, NERRS, 229 $ NFAIL, NRHS, NROWS, NRUN, RANK, MB, 230 $ MMAX, NMAX, NSMAX, LIWORK, 231 $ LWORK_DGELS, LWORK_DGETSLS, LWORK_DGELSS, 232 $ LWORK_DGELSY, LWORK_DGELSD 233 DOUBLE PRECISION EPS, NORMA, NORMB, RCOND 234* .. 235* .. Local Arrays .. 236 INTEGER ISEED( 4 ), ISEEDY( 4 ), IWQ( 1 ) 237 DOUBLE PRECISION RESULT( NTESTS ), WQ( 1 ) 238* .. 239* .. Allocatable Arrays .. 240 DOUBLE PRECISION, ALLOCATABLE :: WORK (:) 241 INTEGER, ALLOCATABLE :: IWORK (:) 242* .. 243* .. External Functions .. 244 DOUBLE PRECISION DASUM, DLAMCH, DQRT12, DQRT14, DQRT17 245 EXTERNAL DASUM, DLAMCH, DQRT12, DQRT14, DQRT17 246* .. 247* .. External Subroutines .. 248 EXTERNAL ALAERH, ALAHD, ALASVM, DAXPY, DERRLS, DGELS, 249 $ DGELSD, DGELSS, DGELSY, DGEMM, DLACPY, 250 $ DLARNV, DLASRT, DQRT13, DQRT15, DQRT16, DSCAL, 251 $ XLAENV 252* .. 253* .. Intrinsic Functions .. 254 INTRINSIC DBLE, INT, LOG, MAX, MIN, SQRT 255* .. 256* .. Scalars in Common .. 257 LOGICAL LERR, OK 258 CHARACTER*32 SRNAMT 259 INTEGER INFOT, IOUNIT 260* .. 261* .. Common blocks .. 262 COMMON / INFOC / INFOT, IOUNIT, OK, LERR 263 COMMON / SRNAMC / SRNAMT 264* .. 265* .. Data statements .. 266 DATA ISEEDY / 1988, 1989, 1990, 1991 / 267* .. 268* .. Executable Statements .. 269* 270* Initialize constants and the random number seed. 271* 272 PATH( 1: 1 ) = 'Double precision' 273 PATH( 2: 3 ) = 'LS' 274 NRUN = 0 275 NFAIL = 0 276 NERRS = 0 277 DO 10 I = 1, 4 278 ISEED( I ) = ISEEDY( I ) 279 10 CONTINUE 280 EPS = DLAMCH( 'Epsilon' ) 281* 282* Threshold for rank estimation 283* 284 RCOND = SQRT( EPS ) - ( SQRT( EPS )-EPS ) / 2 285* 286* Test the error exits 287* 288 CALL XLAENV( 2, 2 ) 289 CALL XLAENV( 9, SMLSIZ ) 290 IF( TSTERR ) 291 $ CALL DERRLS( PATH, NOUT ) 292* 293* Print the header if NM = 0 or NN = 0 and THRESH = 0. 294* 295 IF( ( NM.EQ.0 .OR. NN.EQ.0 ) .AND. THRESH.EQ.ZERO ) 296 $ CALL ALAHD( NOUT, PATH ) 297 INFOT = 0 298 CALL XLAENV( 2, 2 ) 299 CALL XLAENV( 9, SMLSIZ ) 300* 301* Compute maximal workspace needed for all routines 302* 303 NMAX = 0 304 MMAX = 0 305 NSMAX = 0 306 DO I = 1, NM 307 IF ( MVAL( I ).GT.MMAX ) THEN 308 MMAX = MVAL( I ) 309 END IF 310 ENDDO 311 DO I = 1, NN 312 IF ( NVAL( I ).GT.NMAX ) THEN 313 NMAX = NVAL( I ) 314 END IF 315 ENDDO 316 DO I = 1, NNS 317 IF ( NSVAL( I ).GT.NSMAX ) THEN 318 NSMAX = NSVAL( I ) 319 END IF 320 ENDDO 321 M = MMAX 322 N = NMAX 323 NRHS = NSMAX 324 MNMIN = MAX( MIN( M, N ), 1 ) 325* 326* Compute workspace needed for routines 327* DQRT14, DQRT17 (two side cases), DQRT15 and DQRT12 328* 329 LWORK = MAX( 1, ( M+N )*NRHS, 330 $ ( N+NRHS )*( M+2 ), ( M+NRHS )*( N+2 ), 331 $ MAX( M+MNMIN, NRHS*MNMIN,2*N+M ), 332 $ MAX( M*N+4*MNMIN+MAX(M,N), M*N+2*MNMIN+4*N ) ) 333 LIWORK = 1 334* 335* Iterate through all test cases and compute necessary workspace 336* sizes for ?GELS, ?GETSLS, ?GELSY, ?GELSS and ?GELSD routines. 337* 338 DO IM = 1, NM 339 M = MVAL( IM ) 340 LDA = MAX( 1, M ) 341 DO IN = 1, NN 342 N = NVAL( IN ) 343 MNMIN = MAX(MIN( M, N ),1) 344 LDB = MAX( 1, M, N ) 345 DO INS = 1, NNS 346 NRHS = NSVAL( INS ) 347 DO IRANK = 1, 2 348 DO ISCALE = 1, 3 349 ITYPE = ( IRANK-1 )*3 + ISCALE 350 IF( DOTYPE( ITYPE ) ) THEN 351 IF( IRANK.EQ.1 ) THEN 352 DO ITRAN = 1, 2 353 IF( ITRAN.EQ.1 ) THEN 354 TRANS = 'N' 355 ELSE 356 TRANS = 'T' 357 END IF 358* 359* Compute workspace needed for DGELS 360 CALL DGELS( TRANS, M, N, NRHS, A, LDA, 361 $ B, LDB, WQ, -1, INFO ) 362 LWORK_DGELS = INT ( WQ ( 1 ) ) 363* Compute workspace needed for DGETSLS 364 CALL DGETSLS( TRANS, M, N, NRHS, A, LDA, 365 $ B, LDB, WQ, -1, INFO ) 366 LWORK_DGETSLS = INT( WQ ( 1 ) ) 367 ENDDO 368 END IF 369* Compute workspace needed for DGELSY 370 CALL DGELSY( M, N, NRHS, A, LDA, B, LDB, IWQ, 371 $ RCOND, CRANK, WQ, -1, INFO ) 372 LWORK_DGELSY = INT( WQ ( 1 ) ) 373* Compute workspace needed for DGELSS 374 CALL DGELSS( M, N, NRHS, A, LDA, B, LDB, S, 375 $ RCOND, CRANK, WQ, -1 , INFO ) 376 LWORK_DGELSS = INT( WQ ( 1 ) ) 377* Compute workspace needed for DGELSD 378 CALL DGELSD( M, N, NRHS, A, LDA, B, LDB, S, 379 $ RCOND, CRANK, WQ, -1, IWQ, INFO ) 380 LWORK_DGELSD = INT( WQ ( 1 ) ) 381* Compute LIWORK workspace needed for DGELSY and DGELSD 382 LIWORK = MAX( LIWORK, N, IWQ( 1 ) ) 383* Compute LWORK workspace needed for all functions 384 LWORK = MAX( LWORK, LWORK_DGELS, LWORK_DGETSLS, 385 $ LWORK_DGELSY, LWORK_DGELSS, 386 $ LWORK_DGELSD ) 387 END IF 388 ENDDO 389 ENDDO 390 ENDDO 391 ENDDO 392 ENDDO 393* 394 LWLSY = LWORK 395* 396 ALLOCATE( WORK( LWORK ) ) 397 ALLOCATE( IWORK( LIWORK ) ) 398* 399 DO 150 IM = 1, NM 400 M = MVAL( IM ) 401 LDA = MAX( 1, M ) 402* 403 DO 140 IN = 1, NN 404 N = NVAL( IN ) 405 MNMIN = MAX(MIN( M, N ),1) 406 LDB = MAX( 1, M, N ) 407 MB = (MNMIN+1) 408* 409 DO 130 INS = 1, NNS 410 NRHS = NSVAL( INS ) 411* 412 DO 120 IRANK = 1, 2 413 DO 110 ISCALE = 1, 3 414 ITYPE = ( IRANK-1 )*3 + ISCALE 415 IF( .NOT.DOTYPE( ITYPE ) ) 416 $ GO TO 110 417* 418 IF( IRANK.EQ.1 ) THEN 419* 420* Test DGELS 421* 422* Generate a matrix of scaling type ISCALE 423* 424 CALL DQRT13( ISCALE, M, N, COPYA, LDA, NORMA, 425 $ ISEED ) 426 DO 40 INB = 1, NNB 427 NB = NBVAL( INB ) 428 CALL XLAENV( 1, NB ) 429 CALL XLAENV( 3, NXVAL( INB ) ) 430* 431 DO 30 ITRAN = 1, 2 432 IF( ITRAN.EQ.1 ) THEN 433 TRANS = 'N' 434 NROWS = M 435 NCOLS = N 436 ELSE 437 TRANS = 'T' 438 NROWS = N 439 NCOLS = M 440 END IF 441 LDWORK = MAX( 1, NCOLS ) 442* 443* Set up a consistent rhs 444* 445 IF( NCOLS.GT.0 ) THEN 446 CALL DLARNV( 2, ISEED, NCOLS*NRHS, 447 $ WORK ) 448 CALL DSCAL( NCOLS*NRHS, 449 $ ONE / DBLE( NCOLS ), WORK, 450 $ 1 ) 451 END IF 452 CALL DGEMM( TRANS, 'No transpose', NROWS, 453 $ NRHS, NCOLS, ONE, COPYA, LDA, 454 $ WORK, LDWORK, ZERO, B, LDB ) 455 CALL DLACPY( 'Full', NROWS, NRHS, B, LDB, 456 $ COPYB, LDB ) 457* 458* Solve LS or overdetermined system 459* 460 IF( M.GT.0 .AND. N.GT.0 ) THEN 461 CALL DLACPY( 'Full', M, N, COPYA, LDA, 462 $ A, LDA ) 463 CALL DLACPY( 'Full', NROWS, NRHS, 464 $ COPYB, LDB, B, LDB ) 465 END IF 466 SRNAMT = 'DGELS ' 467 CALL DGELS( TRANS, M, N, NRHS, A, LDA, B, 468 $ LDB, WORK, LWORK, INFO ) 469 IF( INFO.NE.0 ) 470 $ CALL ALAERH( PATH, 'DGELS ', INFO, 0, 471 $ TRANS, M, N, NRHS, -1, NB, 472 $ ITYPE, NFAIL, NERRS, 473 $ NOUT ) 474* 475* Check correctness of results 476* 477 LDWORK = MAX( 1, NROWS ) 478 IF( NROWS.GT.0 .AND. NRHS.GT.0 ) 479 $ CALL DLACPY( 'Full', NROWS, NRHS, 480 $ COPYB, LDB, C, LDB ) 481 CALL DQRT16( TRANS, M, N, NRHS, COPYA, 482 $ LDA, B, LDB, C, LDB, WORK, 483 $ RESULT( 1 ) ) 484* 485 IF( ( ITRAN.EQ.1 .AND. M.GE.N ) .OR. 486 $ ( ITRAN.EQ.2 .AND. M.LT.N ) ) THEN 487* 488* Solving LS system 489* 490 RESULT( 2 ) = DQRT17( TRANS, 1, M, N, 491 $ NRHS, COPYA, LDA, B, LDB, 492 $ COPYB, LDB, C, WORK, 493 $ LWORK ) 494 ELSE 495* 496* Solving overdetermined system 497* 498 RESULT( 2 ) = DQRT14( TRANS, M, N, 499 $ NRHS, COPYA, LDA, B, LDB, 500 $ WORK, LWORK ) 501 END IF 502* 503* Print information about the tests that 504* did not pass the threshold. 505* 506 DO 20 K = 1, 2 507 IF( RESULT( K ).GE.THRESH ) THEN 508 IF( NFAIL.EQ.0 .AND. NERRS.EQ.0 ) 509 $ CALL ALAHD( NOUT, PATH ) 510 WRITE( NOUT, FMT = 9999 )TRANS, M, 511 $ N, NRHS, NB, ITYPE, K, 512 $ RESULT( K ) 513 NFAIL = NFAIL + 1 514 END IF 515 20 CONTINUE 516 NRUN = NRUN + 2 517 30 CONTINUE 518 40 CONTINUE 519* 520* 521* Test DGETSLS 522* 523* Generate a matrix of scaling type ISCALE 524* 525 CALL DQRT13( ISCALE, M, N, COPYA, LDA, NORMA, 526 $ ISEED ) 527 DO 65 INB = 1, NNB 528 MB = NBVAL( INB ) 529 CALL XLAENV( 1, MB ) 530 DO 62 IMB = 1, NNB 531 NB = NBVAL( IMB ) 532 CALL XLAENV( 2, NB ) 533* 534 DO 60 ITRAN = 1, 2 535 IF( ITRAN.EQ.1 ) THEN 536 TRANS = 'N' 537 NROWS = M 538 NCOLS = N 539 ELSE 540 TRANS = 'T' 541 NROWS = N 542 NCOLS = M 543 END IF 544 LDWORK = MAX( 1, NCOLS ) 545* 546* Set up a consistent rhs 547* 548 IF( NCOLS.GT.0 ) THEN 549 CALL DLARNV( 2, ISEED, NCOLS*NRHS, 550 $ WORK ) 551 CALL DSCAL( NCOLS*NRHS, 552 $ ONE / DBLE( NCOLS ), WORK, 553 $ 1 ) 554 END IF 555 CALL DGEMM( TRANS, 'No transpose', NROWS, 556 $ NRHS, NCOLS, ONE, COPYA, LDA, 557 $ WORK, LDWORK, ZERO, B, LDB ) 558 CALL DLACPY( 'Full', NROWS, NRHS, B, LDB, 559 $ COPYB, LDB ) 560* 561* Solve LS or overdetermined system 562* 563 IF( M.GT.0 .AND. N.GT.0 ) THEN 564 CALL DLACPY( 'Full', M, N, COPYA, LDA, 565 $ A, LDA ) 566 CALL DLACPY( 'Full', NROWS, NRHS, 567 $ COPYB, LDB, B, LDB ) 568 END IF 569 SRNAMT = 'DGETSLS ' 570 CALL DGETSLS( TRANS, M, N, NRHS, A, 571 $ LDA, B, LDB, WORK, LWORK, INFO ) 572 IF( INFO.NE.0 ) 573 $ CALL ALAERH( PATH, 'DGETSLS ', INFO, 0, 574 $ TRANS, M, N, NRHS, -1, NB, 575 $ ITYPE, NFAIL, NERRS, 576 $ NOUT ) 577* 578* Check correctness of results 579* 580 LDWORK = MAX( 1, NROWS ) 581 IF( NROWS.GT.0 .AND. NRHS.GT.0 ) 582 $ CALL DLACPY( 'Full', NROWS, NRHS, 583 $ COPYB, LDB, C, LDB ) 584 CALL DQRT16( TRANS, M, N, NRHS, COPYA, 585 $ LDA, B, LDB, C, LDB, WORK, 586 $ RESULT( 15 ) ) 587* 588 IF( ( ITRAN.EQ.1 .AND. M.GE.N ) .OR. 589 $ ( ITRAN.EQ.2 .AND. M.LT.N ) ) THEN 590* 591* Solving LS system 592* 593 RESULT( 16 ) = DQRT17( TRANS, 1, M, N, 594 $ NRHS, COPYA, LDA, B, LDB, 595 $ COPYB, LDB, C, WORK, 596 $ LWORK ) 597 ELSE 598* 599* Solving overdetermined system 600* 601 RESULT( 16 ) = DQRT14( TRANS, M, N, 602 $ NRHS, COPYA, LDA, B, LDB, 603 $ WORK, LWORK ) 604 END IF 605* 606* Print information about the tests that 607* did not pass the threshold. 608* 609 DO 50 K = 15, 16 610 IF( RESULT( K ).GE.THRESH ) THEN 611 IF( NFAIL.EQ.0 .AND. NERRS.EQ.0 ) 612 $ CALL ALAHD( NOUT, PATH ) 613 WRITE( NOUT, FMT = 9997 )TRANS, M, 614 $ N, NRHS, MB, NB, ITYPE, K, 615 $ RESULT( K ) 616 NFAIL = NFAIL + 1 617 END IF 618 50 CONTINUE 619 NRUN = NRUN + 2 620 60 CONTINUE 621 62 CONTINUE 622 65 CONTINUE 623 END IF 624* 625* Generate a matrix of scaling type ISCALE and rank 626* type IRANK. 627* 628 CALL DQRT15( ISCALE, IRANK, M, N, NRHS, COPYA, LDA, 629 $ COPYB, LDB, COPYS, RANK, NORMA, NORMB, 630 $ ISEED, WORK, LWORK ) 631* 632* workspace used: MAX(M+MIN(M,N),NRHS*MIN(M,N),2*N+M) 633* 634 LDWORK = MAX( 1, M ) 635* 636* Loop for testing different block sizes. 637* 638 DO 100 INB = 1, NNB 639 NB = NBVAL( INB ) 640 CALL XLAENV( 1, NB ) 641 CALL XLAENV( 3, NXVAL( INB ) ) 642* 643* Test DGELSY 644* 645* DGELSY: Compute the minimum-norm solution X 646* to min( norm( A * X - B ) ) 647* using the rank-revealing orthogonal 648* factorization. 649* 650* Initialize vector IWORK. 651* 652 DO 70 J = 1, N 653 IWORK( J ) = 0 654 70 CONTINUE 655* 656 CALL DLACPY( 'Full', M, N, COPYA, LDA, A, LDA ) 657 CALL DLACPY( 'Full', M, NRHS, COPYB, LDB, B, 658 $ LDB ) 659* 660 SRNAMT = 'DGELSY' 661 CALL DGELSY( M, N, NRHS, A, LDA, B, LDB, IWORK, 662 $ RCOND, CRANK, WORK, LWLSY, INFO ) 663 IF( INFO.NE.0 ) 664 $ CALL ALAERH( PATH, 'DGELSY', INFO, 0, ' ', M, 665 $ N, NRHS, -1, NB, ITYPE, NFAIL, 666 $ NERRS, NOUT ) 667* 668* Test 3: Compute relative error in svd 669* workspace: M*N + 4*MIN(M,N) + MAX(M,N) 670* 671 RESULT( 3 ) = DQRT12( CRANK, CRANK, A, LDA, 672 $ COPYS, WORK, LWORK ) 673* 674* Test 4: Compute error in solution 675* workspace: M*NRHS + M 676* 677 CALL DLACPY( 'Full', M, NRHS, COPYB, LDB, WORK, 678 $ LDWORK ) 679 CALL DQRT16( 'No transpose', M, N, NRHS, COPYA, 680 $ LDA, B, LDB, WORK, LDWORK, 681 $ WORK( M*NRHS+1 ), RESULT( 4 ) ) 682* 683* Test 5: Check norm of r'*A 684* workspace: NRHS*(M+N) 685* 686 RESULT( 5 ) = ZERO 687 IF( M.GT.CRANK ) 688 $ RESULT( 5 ) = DQRT17( 'No transpose', 1, M, 689 $ N, NRHS, COPYA, LDA, B, LDB, 690 $ COPYB, LDB, C, WORK, LWORK ) 691* 692* Test 6: Check if x is in the rowspace of A 693* workspace: (M+NRHS)*(N+2) 694* 695 RESULT( 6 ) = ZERO 696* 697 IF( N.GT.CRANK ) 698 $ RESULT( 6 ) = DQRT14( 'No transpose', M, N, 699 $ NRHS, COPYA, LDA, B, LDB, 700 $ WORK, LWORK ) 701* 702* Test DGELSS 703* 704* DGELSS: Compute the minimum-norm solution X 705* to min( norm( A * X - B ) ) 706* using the SVD. 707* 708 CALL DLACPY( 'Full', M, N, COPYA, LDA, A, LDA ) 709 CALL DLACPY( 'Full', M, NRHS, COPYB, LDB, B, 710 $ LDB ) 711 SRNAMT = 'DGELSS' 712 CALL DGELSS( M, N, NRHS, A, LDA, B, LDB, S, 713 $ RCOND, CRANK, WORK, LWORK, INFO ) 714 IF( INFO.NE.0 ) 715 $ CALL ALAERH( PATH, 'DGELSS', INFO, 0, ' ', M, 716 $ N, NRHS, -1, NB, ITYPE, NFAIL, 717 $ NERRS, NOUT ) 718* 719* workspace used: 3*min(m,n) + 720* max(2*min(m,n),nrhs,max(m,n)) 721* 722* Test 7: Compute relative error in svd 723* 724 IF( RANK.GT.0 ) THEN 725 CALL DAXPY( MNMIN, -ONE, COPYS, 1, S, 1 ) 726 RESULT( 7 ) = DASUM( MNMIN, S, 1 ) / 727 $ DASUM( MNMIN, COPYS, 1 ) / 728 $ ( EPS*DBLE( MNMIN ) ) 729 ELSE 730 RESULT( 7 ) = ZERO 731 END IF 732* 733* Test 8: Compute error in solution 734* 735 CALL DLACPY( 'Full', M, NRHS, COPYB, LDB, WORK, 736 $ LDWORK ) 737 CALL DQRT16( 'No transpose', M, N, NRHS, COPYA, 738 $ LDA, B, LDB, WORK, LDWORK, 739 $ WORK( M*NRHS+1 ), RESULT( 8 ) ) 740* 741* Test 9: Check norm of r'*A 742* 743 RESULT( 9 ) = ZERO 744 IF( M.GT.CRANK ) 745 $ RESULT( 9 ) = DQRT17( 'No transpose', 1, M, 746 $ N, NRHS, COPYA, LDA, B, LDB, 747 $ COPYB, LDB, C, WORK, LWORK ) 748* 749* Test 10: Check if x is in the rowspace of A 750* 751 RESULT( 10 ) = ZERO 752 IF( N.GT.CRANK ) 753 $ RESULT( 10 ) = DQRT14( 'No transpose', M, N, 754 $ NRHS, COPYA, LDA, B, LDB, 755 $ WORK, LWORK ) 756* 757* Test DGELSD 758* 759* DGELSD: Compute the minimum-norm solution X 760* to min( norm( A * X - B ) ) using a 761* divide and conquer SVD. 762* 763* Initialize vector IWORK. 764* 765 DO 80 J = 1, N 766 IWORK( J ) = 0 767 80 CONTINUE 768* 769 CALL DLACPY( 'Full', M, N, COPYA, LDA, A, LDA ) 770 CALL DLACPY( 'Full', M, NRHS, COPYB, LDB, B, 771 $ LDB ) 772* 773 SRNAMT = 'DGELSD' 774 CALL DGELSD( M, N, NRHS, A, LDA, B, LDB, S, 775 $ RCOND, CRANK, WORK, LWORK, IWORK, 776 $ INFO ) 777 IF( INFO.NE.0 ) 778 $ CALL ALAERH( PATH, 'DGELSD', INFO, 0, ' ', M, 779 $ N, NRHS, -1, NB, ITYPE, NFAIL, 780 $ NERRS, NOUT ) 781* 782* Test 11: Compute relative error in svd 783* 784 IF( RANK.GT.0 ) THEN 785 CALL DAXPY( MNMIN, -ONE, COPYS, 1, S, 1 ) 786 RESULT( 11 ) = DASUM( MNMIN, S, 1 ) / 787 $ DASUM( MNMIN, COPYS, 1 ) / 788 $ ( EPS*DBLE( MNMIN ) ) 789 ELSE 790 RESULT( 11 ) = ZERO 791 END IF 792* 793* Test 12: Compute error in solution 794* 795 CALL DLACPY( 'Full', M, NRHS, COPYB, LDB, WORK, 796 $ LDWORK ) 797 CALL DQRT16( 'No transpose', M, N, NRHS, COPYA, 798 $ LDA, B, LDB, WORK, LDWORK, 799 $ WORK( M*NRHS+1 ), RESULT( 12 ) ) 800* 801* Test 13: Check norm of r'*A 802* 803 RESULT( 13 ) = ZERO 804 IF( M.GT.CRANK ) 805 $ RESULT( 13 ) = DQRT17( 'No transpose', 1, M, 806 $ N, NRHS, COPYA, LDA, B, LDB, 807 $ COPYB, LDB, C, WORK, LWORK ) 808* 809* Test 14: Check if x is in the rowspace of A 810* 811 RESULT( 14 ) = ZERO 812 IF( N.GT.CRANK ) 813 $ RESULT( 14 ) = DQRT14( 'No transpose', M, N, 814 $ NRHS, COPYA, LDA, B, LDB, 815 $ WORK, LWORK ) 816* 817* Print information about the tests that did not 818* pass the threshold. 819* 820 DO 90 K = 3, 14 821 IF( RESULT( K ).GE.THRESH ) THEN 822 IF( NFAIL.EQ.0 .AND. NERRS.EQ.0 ) 823 $ CALL ALAHD( NOUT, PATH ) 824 WRITE( NOUT, FMT = 9998 )M, N, NRHS, NB, 825 $ ITYPE, K, RESULT( K ) 826 NFAIL = NFAIL + 1 827 END IF 828 90 CONTINUE 829 NRUN = NRUN + 12 830* 831 100 CONTINUE 832 110 CONTINUE 833 120 CONTINUE 834 130 CONTINUE 835 140 CONTINUE 836 150 CONTINUE 837* 838* Print a summary of the results. 839* 840 CALL ALASVM( PATH, NOUT, NFAIL, NRUN, NERRS ) 841* 842 9999 FORMAT( ' TRANS=''', A1, ''', M=', I5, ', N=', I5, ', NRHS=', I4, 843 $ ', NB=', I4, ', type', I2, ', test(', I2, ')=', G12.5 ) 844 9998 FORMAT( ' M=', I5, ', N=', I5, ', NRHS=', I4, ', NB=', I4, 845 $ ', type', I2, ', test(', I2, ')=', G12.5 ) 846 9997 FORMAT( ' TRANS=''', A1,' M=', I5, ', N=', I5, ', NRHS=', I4, 847 $ ', MB=', I4,', NB=', I4,', type', I2, 848 $ ', test(', I2, ')=', G12.5 ) 849* 850 DEALLOCATE( WORK ) 851 DEALLOCATE( IWORK ) 852 RETURN 853* 854* End of DDRVLS 855* 856 END 857