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*> \ingroup double_lin 187* 188* ===================================================================== 189 SUBROUTINE DDRVLS( DOTYPE, NM, MVAL, NN, NVAL, NNS, NSVAL, NNB, 190 $ NBVAL, NXVAL, THRESH, TSTERR, A, COPYA, B, 191 $ COPYB, C, S, COPYS, NOUT ) 192* 193* -- LAPACK test routine -- 194* -- LAPACK is a software package provided by Univ. of Tennessee, -- 195* -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..-- 196* 197* .. Scalar Arguments .. 198 LOGICAL TSTERR 199 INTEGER NM, NN, NNB, NNS, NOUT 200 DOUBLE PRECISION THRESH 201* .. 202* .. Array Arguments .. 203 LOGICAL DOTYPE( * ) 204 INTEGER MVAL( * ), NBVAL( * ), NSVAL( * ), 205 $ NVAL( * ), NXVAL( * ) 206 DOUBLE PRECISION A( * ), B( * ), C( * ), COPYA( * ), COPYB( * ), 207 $ COPYS( * ), S( * ) 208* .. 209* 210* ===================================================================== 211* 212* .. Parameters .. 213 INTEGER NTESTS 214 PARAMETER ( NTESTS = 16 ) 215 INTEGER SMLSIZ 216 PARAMETER ( SMLSIZ = 25 ) 217 DOUBLE PRECISION ONE, TWO, ZERO 218 PARAMETER ( ONE = 1.0D0, TWO = 2.0D0, ZERO = 0.0D0 ) 219* .. 220* .. Local Scalars .. 221 CHARACTER TRANS 222 CHARACTER*3 PATH 223 INTEGER CRANK, I, IM, IMB, IN, INB, INFO, INS, IRANK, 224 $ ISCALE, ITRAN, ITYPE, J, K, LDA, LDB, LDWORK, 225 $ LWLSY, LWORK, M, MNMIN, N, NB, NCOLS, NERRS, 226 $ NFAIL, NRHS, NROWS, NRUN, RANK, MB, 227 $ MMAX, NMAX, NSMAX, LIWORK, 228 $ LWORK_DGELS, LWORK_DGETSLS, LWORK_DGELSS, 229 $ LWORK_DGELSY, LWORK_DGELSD 230 DOUBLE PRECISION EPS, NORMA, NORMB, RCOND 231* .. 232* .. Local Arrays .. 233 INTEGER ISEED( 4 ), ISEEDY( 4 ), IWQ( 1 ) 234 DOUBLE PRECISION RESULT( NTESTS ), WQ( 1 ) 235* .. 236* .. Allocatable Arrays .. 237 DOUBLE PRECISION, ALLOCATABLE :: WORK (:) 238 INTEGER, ALLOCATABLE :: IWORK (:) 239* .. 240* .. External Functions .. 241 DOUBLE PRECISION DASUM, DLAMCH, DQRT12, DQRT14, DQRT17 242 EXTERNAL DASUM, DLAMCH, DQRT12, DQRT14, DQRT17 243* .. 244* .. External Subroutines .. 245 EXTERNAL ALAERH, ALAHD, ALASVM, DAXPY, DERRLS, DGELS, 246 $ DGELSD, DGELSS, DGELSY, DGEMM, DLACPY, 247 $ DLARNV, DLASRT, DQRT13, DQRT15, DQRT16, DSCAL, 248 $ XLAENV 249* .. 250* .. Intrinsic Functions .. 251 INTRINSIC DBLE, INT, LOG, MAX, MIN, SQRT 252* .. 253* .. Scalars in Common .. 254 LOGICAL LERR, OK 255 CHARACTER*32 SRNAMT 256 INTEGER INFOT, IOUNIT 257* .. 258* .. Common blocks .. 259 COMMON / INFOC / INFOT, IOUNIT, OK, LERR 260 COMMON / SRNAMC / SRNAMT 261* .. 262* .. Data statements .. 263 DATA ISEEDY / 1988, 1989, 1990, 1991 / 264* .. 265* .. Executable Statements .. 266* 267* Initialize constants and the random number seed. 268* 269 PATH( 1: 1 ) = 'Double precision' 270 PATH( 2: 3 ) = 'LS' 271 NRUN = 0 272 NFAIL = 0 273 NERRS = 0 274 DO 10 I = 1, 4 275 ISEED( I ) = ISEEDY( I ) 276 10 CONTINUE 277 EPS = DLAMCH( 'Epsilon' ) 278* 279* Threshold for rank estimation 280* 281 RCOND = SQRT( EPS ) - ( SQRT( EPS )-EPS ) / 2 282* 283* Test the error exits 284* 285 CALL XLAENV( 2, 2 ) 286 CALL XLAENV( 9, SMLSIZ ) 287 IF( TSTERR ) 288 $ CALL DERRLS( PATH, NOUT ) 289* 290* Print the header if NM = 0 or NN = 0 and THRESH = 0. 291* 292 IF( ( NM.EQ.0 .OR. NN.EQ.0 ) .AND. THRESH.EQ.ZERO ) 293 $ CALL ALAHD( NOUT, PATH ) 294 INFOT = 0 295 CALL XLAENV( 2, 2 ) 296 CALL XLAENV( 9, SMLSIZ ) 297* 298* Compute maximal workspace needed for all routines 299* 300 NMAX = 0 301 MMAX = 0 302 NSMAX = 0 303 DO I = 1, NM 304 IF ( MVAL( I ).GT.MMAX ) THEN 305 MMAX = MVAL( I ) 306 END IF 307 ENDDO 308 DO I = 1, NN 309 IF ( NVAL( I ).GT.NMAX ) THEN 310 NMAX = NVAL( I ) 311 END IF 312 ENDDO 313 DO I = 1, NNS 314 IF ( NSVAL( I ).GT.NSMAX ) THEN 315 NSMAX = NSVAL( I ) 316 END IF 317 ENDDO 318 M = MMAX 319 N = NMAX 320 NRHS = NSMAX 321 MNMIN = MAX( MIN( M, N ), 1 ) 322* 323* Compute workspace needed for routines 324* DQRT14, DQRT17 (two side cases), DQRT15 and DQRT12 325* 326 LWORK = MAX( 1, ( M+N )*NRHS, 327 $ ( N+NRHS )*( M+2 ), ( M+NRHS )*( N+2 ), 328 $ MAX( M+MNMIN, NRHS*MNMIN,2*N+M ), 329 $ MAX( M*N+4*MNMIN+MAX(M,N), M*N+2*MNMIN+4*N ) ) 330 LIWORK = 1 331* 332* Iterate through all test cases and compute necessary workspace 333* sizes for ?GELS, ?GETSLS, ?GELSY, ?GELSS and ?GELSD routines. 334* 335 DO IM = 1, NM 336 M = MVAL( IM ) 337 LDA = MAX( 1, M ) 338 DO IN = 1, NN 339 N = NVAL( IN ) 340 MNMIN = MAX(MIN( M, N ),1) 341 LDB = MAX( 1, M, N ) 342 DO INS = 1, NNS 343 NRHS = NSVAL( INS ) 344 DO IRANK = 1, 2 345 DO ISCALE = 1, 3 346 ITYPE = ( IRANK-1 )*3 + ISCALE 347 IF( DOTYPE( ITYPE ) ) THEN 348 IF( IRANK.EQ.1 ) THEN 349 DO ITRAN = 1, 2 350 IF( ITRAN.EQ.1 ) THEN 351 TRANS = 'N' 352 ELSE 353 TRANS = 'T' 354 END IF 355* 356* Compute workspace needed for DGELS 357 CALL DGELS( TRANS, M, N, NRHS, A, LDA, 358 $ B, LDB, WQ, -1, INFO ) 359 LWORK_DGELS = INT ( WQ ( 1 ) ) 360* Compute workspace needed for DGETSLS 361 CALL DGETSLS( TRANS, M, N, NRHS, A, LDA, 362 $ B, LDB, WQ, -1, INFO ) 363 LWORK_DGETSLS = INT( WQ ( 1 ) ) 364 ENDDO 365 END IF 366* Compute workspace needed for DGELSY 367 CALL DGELSY( M, N, NRHS, A, LDA, B, LDB, IWQ, 368 $ RCOND, CRANK, WQ, -1, INFO ) 369 LWORK_DGELSY = INT( WQ ( 1 ) ) 370* Compute workspace needed for DGELSS 371 CALL DGELSS( M, N, NRHS, A, LDA, B, LDB, S, 372 $ RCOND, CRANK, WQ, -1 , INFO ) 373 LWORK_DGELSS = INT( WQ ( 1 ) ) 374* Compute workspace needed for DGELSD 375 CALL DGELSD( M, N, NRHS, A, LDA, B, LDB, S, 376 $ RCOND, CRANK, WQ, -1, IWQ, INFO ) 377 LWORK_DGELSD = INT( WQ ( 1 ) ) 378* Compute LIWORK workspace needed for DGELSY and DGELSD 379 LIWORK = MAX( LIWORK, N, IWQ( 1 ) ) 380* Compute LWORK workspace needed for all functions 381 LWORK = MAX( LWORK, LWORK_DGELS, LWORK_DGETSLS, 382 $ LWORK_DGELSY, LWORK_DGELSS, 383 $ LWORK_DGELSD ) 384 END IF 385 ENDDO 386 ENDDO 387 ENDDO 388 ENDDO 389 ENDDO 390* 391 LWLSY = LWORK 392* 393 ALLOCATE( WORK( LWORK ) ) 394 ALLOCATE( IWORK( LIWORK ) ) 395* 396 DO 150 IM = 1, NM 397 M = MVAL( IM ) 398 LDA = MAX( 1, M ) 399* 400 DO 140 IN = 1, NN 401 N = NVAL( IN ) 402 MNMIN = MAX(MIN( M, N ),1) 403 LDB = MAX( 1, M, N ) 404 MB = (MNMIN+1) 405* 406 DO 130 INS = 1, NNS 407 NRHS = NSVAL( INS ) 408* 409 DO 120 IRANK = 1, 2 410 DO 110 ISCALE = 1, 3 411 ITYPE = ( IRANK-1 )*3 + ISCALE 412 IF( .NOT.DOTYPE( ITYPE ) ) 413 $ GO TO 110 414* 415 IF( IRANK.EQ.1 ) THEN 416* 417* Test DGELS 418* 419* Generate a matrix of scaling type ISCALE 420* 421 CALL DQRT13( ISCALE, M, N, COPYA, LDA, NORMA, 422 $ ISEED ) 423 DO 40 INB = 1, NNB 424 NB = NBVAL( INB ) 425 CALL XLAENV( 1, NB ) 426 CALL XLAENV( 3, NXVAL( INB ) ) 427* 428 DO 30 ITRAN = 1, 2 429 IF( ITRAN.EQ.1 ) THEN 430 TRANS = 'N' 431 NROWS = M 432 NCOLS = N 433 ELSE 434 TRANS = 'T' 435 NROWS = N 436 NCOLS = M 437 END IF 438 LDWORK = MAX( 1, NCOLS ) 439* 440* Set up a consistent rhs 441* 442 IF( NCOLS.GT.0 ) THEN 443 CALL DLARNV( 2, ISEED, NCOLS*NRHS, 444 $ WORK ) 445 CALL DSCAL( NCOLS*NRHS, 446 $ ONE / DBLE( NCOLS ), WORK, 447 $ 1 ) 448 END IF 449 CALL DGEMM( TRANS, 'No transpose', NROWS, 450 $ NRHS, NCOLS, ONE, COPYA, LDA, 451 $ WORK, LDWORK, ZERO, B, LDB ) 452 CALL DLACPY( 'Full', NROWS, NRHS, B, LDB, 453 $ COPYB, LDB ) 454* 455* Solve LS or overdetermined system 456* 457 IF( M.GT.0 .AND. N.GT.0 ) THEN 458 CALL DLACPY( 'Full', M, N, COPYA, LDA, 459 $ A, LDA ) 460 CALL DLACPY( 'Full', NROWS, NRHS, 461 $ COPYB, LDB, B, LDB ) 462 END IF 463 SRNAMT = 'DGELS ' 464 CALL DGELS( TRANS, M, N, NRHS, A, LDA, B, 465 $ LDB, WORK, LWORK, INFO ) 466 IF( INFO.NE.0 ) 467 $ CALL ALAERH( PATH, 'DGELS ', INFO, 0, 468 $ TRANS, M, N, NRHS, -1, NB, 469 $ ITYPE, NFAIL, NERRS, 470 $ NOUT ) 471* 472* Check correctness of results 473* 474 LDWORK = MAX( 1, NROWS ) 475 IF( NROWS.GT.0 .AND. NRHS.GT.0 ) 476 $ CALL DLACPY( 'Full', NROWS, NRHS, 477 $ COPYB, LDB, C, LDB ) 478 CALL DQRT16( TRANS, M, N, NRHS, COPYA, 479 $ LDA, B, LDB, C, LDB, WORK, 480 $ RESULT( 1 ) ) 481* 482 IF( ( ITRAN.EQ.1 .AND. M.GE.N ) .OR. 483 $ ( ITRAN.EQ.2 .AND. M.LT.N ) ) THEN 484* 485* Solving LS system 486* 487 RESULT( 2 ) = DQRT17( TRANS, 1, M, N, 488 $ NRHS, COPYA, LDA, B, LDB, 489 $ COPYB, LDB, C, WORK, 490 $ LWORK ) 491 ELSE 492* 493* Solving overdetermined system 494* 495 RESULT( 2 ) = DQRT14( TRANS, M, N, 496 $ NRHS, COPYA, LDA, B, LDB, 497 $ WORK, LWORK ) 498 END IF 499* 500* Print information about the tests that 501* did not pass the threshold. 502* 503 DO 20 K = 1, 2 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 = 9999 )TRANS, M, 508 $ N, NRHS, NB, ITYPE, K, 509 $ RESULT( K ) 510 NFAIL = NFAIL + 1 511 END IF 512 20 CONTINUE 513 NRUN = NRUN + 2 514 30 CONTINUE 515 40 CONTINUE 516* 517* 518* Test DGETSLS 519* 520* Generate a matrix of scaling type ISCALE 521* 522 CALL DQRT13( ISCALE, M, N, COPYA, LDA, NORMA, 523 $ ISEED ) 524 DO 65 INB = 1, NNB 525 MB = NBVAL( INB ) 526 CALL XLAENV( 1, MB ) 527 DO 62 IMB = 1, NNB 528 NB = NBVAL( IMB ) 529 CALL XLAENV( 2, NB ) 530* 531 DO 60 ITRAN = 1, 2 532 IF( ITRAN.EQ.1 ) THEN 533 TRANS = 'N' 534 NROWS = M 535 NCOLS = N 536 ELSE 537 TRANS = 'T' 538 NROWS = N 539 NCOLS = M 540 END IF 541 LDWORK = MAX( 1, NCOLS ) 542* 543* Set up a consistent rhs 544* 545 IF( NCOLS.GT.0 ) THEN 546 CALL DLARNV( 2, ISEED, NCOLS*NRHS, 547 $ WORK ) 548 CALL DSCAL( NCOLS*NRHS, 549 $ ONE / DBLE( NCOLS ), WORK, 550 $ 1 ) 551 END IF 552 CALL DGEMM( TRANS, 'No transpose', NROWS, 553 $ NRHS, NCOLS, ONE, COPYA, LDA, 554 $ WORK, LDWORK, ZERO, B, LDB ) 555 CALL DLACPY( 'Full', NROWS, NRHS, B, LDB, 556 $ COPYB, LDB ) 557* 558* Solve LS or overdetermined system 559* 560 IF( M.GT.0 .AND. N.GT.0 ) THEN 561 CALL DLACPY( 'Full', M, N, COPYA, LDA, 562 $ A, LDA ) 563 CALL DLACPY( 'Full', NROWS, NRHS, 564 $ COPYB, LDB, B, LDB ) 565 END IF 566 SRNAMT = 'DGETSLS ' 567 CALL DGETSLS( TRANS, M, N, NRHS, A, 568 $ LDA, B, LDB, WORK, LWORK, INFO ) 569 IF( INFO.NE.0 ) 570 $ CALL ALAERH( PATH, 'DGETSLS ', INFO, 0, 571 $ TRANS, M, N, NRHS, -1, NB, 572 $ ITYPE, NFAIL, NERRS, 573 $ NOUT ) 574* 575* Check correctness of results 576* 577 LDWORK = MAX( 1, NROWS ) 578 IF( NROWS.GT.0 .AND. NRHS.GT.0 ) 579 $ CALL DLACPY( 'Full', NROWS, NRHS, 580 $ COPYB, LDB, C, LDB ) 581 CALL DQRT16( TRANS, M, N, NRHS, COPYA, 582 $ LDA, B, LDB, C, LDB, WORK, 583 $ RESULT( 15 ) ) 584* 585 IF( ( ITRAN.EQ.1 .AND. M.GE.N ) .OR. 586 $ ( ITRAN.EQ.2 .AND. M.LT.N ) ) THEN 587* 588* Solving LS system 589* 590 RESULT( 16 ) = DQRT17( TRANS, 1, M, N, 591 $ NRHS, COPYA, LDA, B, LDB, 592 $ COPYB, LDB, C, WORK, 593 $ LWORK ) 594 ELSE 595* 596* Solving overdetermined system 597* 598 RESULT( 16 ) = DQRT14( TRANS, M, N, 599 $ NRHS, COPYA, LDA, B, LDB, 600 $ WORK, LWORK ) 601 END IF 602* 603* Print information about the tests that 604* did not pass the threshold. 605* 606 DO 50 K = 15, 16 607 IF( RESULT( K ).GE.THRESH ) THEN 608 IF( NFAIL.EQ.0 .AND. NERRS.EQ.0 ) 609 $ CALL ALAHD( NOUT, PATH ) 610 WRITE( NOUT, FMT = 9997 )TRANS, M, 611 $ N, NRHS, MB, NB, ITYPE, K, 612 $ RESULT( K ) 613 NFAIL = NFAIL + 1 614 END IF 615 50 CONTINUE 616 NRUN = NRUN + 2 617 60 CONTINUE 618 62 CONTINUE 619 65 CONTINUE 620 END IF 621* 622* Generate a matrix of scaling type ISCALE and rank 623* type IRANK. 624* 625 CALL DQRT15( ISCALE, IRANK, M, N, NRHS, COPYA, LDA, 626 $ COPYB, LDB, COPYS, RANK, NORMA, NORMB, 627 $ ISEED, WORK, LWORK ) 628* 629* workspace used: MAX(M+MIN(M,N),NRHS*MIN(M,N),2*N+M) 630* 631 LDWORK = MAX( 1, M ) 632* 633* Loop for testing different block sizes. 634* 635 DO 100 INB = 1, NNB 636 NB = NBVAL( INB ) 637 CALL XLAENV( 1, NB ) 638 CALL XLAENV( 3, NXVAL( INB ) ) 639* 640* Test DGELSY 641* 642* DGELSY: Compute the minimum-norm solution X 643* to min( norm( A * X - B ) ) 644* using the rank-revealing orthogonal 645* factorization. 646* 647* Initialize vector IWORK. 648* 649 DO 70 J = 1, N 650 IWORK( J ) = 0 651 70 CONTINUE 652* 653 CALL DLACPY( 'Full', M, N, COPYA, LDA, A, LDA ) 654 CALL DLACPY( 'Full', M, NRHS, COPYB, LDB, B, 655 $ LDB ) 656* 657 SRNAMT = 'DGELSY' 658 CALL DGELSY( M, N, NRHS, A, LDA, B, LDB, IWORK, 659 $ RCOND, CRANK, WORK, LWLSY, INFO ) 660 IF( INFO.NE.0 ) 661 $ CALL ALAERH( PATH, 'DGELSY', INFO, 0, ' ', M, 662 $ N, NRHS, -1, NB, ITYPE, NFAIL, 663 $ NERRS, NOUT ) 664* 665* Test 3: Compute relative error in svd 666* workspace: M*N + 4*MIN(M,N) + MAX(M,N) 667* 668 RESULT( 3 ) = DQRT12( CRANK, CRANK, A, LDA, 669 $ COPYS, WORK, LWORK ) 670* 671* Test 4: Compute error in solution 672* workspace: M*NRHS + M 673* 674 CALL DLACPY( 'Full', M, NRHS, COPYB, LDB, WORK, 675 $ LDWORK ) 676 CALL DQRT16( 'No transpose', M, N, NRHS, COPYA, 677 $ LDA, B, LDB, WORK, LDWORK, 678 $ WORK( M*NRHS+1 ), RESULT( 4 ) ) 679* 680* Test 5: Check norm of r'*A 681* workspace: NRHS*(M+N) 682* 683 RESULT( 5 ) = ZERO 684 IF( M.GT.CRANK ) 685 $ RESULT( 5 ) = DQRT17( 'No transpose', 1, M, 686 $ N, NRHS, COPYA, LDA, B, LDB, 687 $ COPYB, LDB, C, WORK, LWORK ) 688* 689* Test 6: Check if x is in the rowspace of A 690* workspace: (M+NRHS)*(N+2) 691* 692 RESULT( 6 ) = ZERO 693* 694 IF( N.GT.CRANK ) 695 $ RESULT( 6 ) = DQRT14( 'No transpose', M, N, 696 $ NRHS, COPYA, LDA, B, LDB, 697 $ WORK, LWORK ) 698* 699* Test DGELSS 700* 701* DGELSS: Compute the minimum-norm solution X 702* to min( norm( A * X - B ) ) 703* using the SVD. 704* 705 CALL DLACPY( 'Full', M, N, COPYA, LDA, A, LDA ) 706 CALL DLACPY( 'Full', M, NRHS, COPYB, LDB, B, 707 $ LDB ) 708 SRNAMT = 'DGELSS' 709 CALL DGELSS( M, N, NRHS, A, LDA, B, LDB, S, 710 $ RCOND, CRANK, WORK, LWORK, INFO ) 711 IF( INFO.NE.0 ) 712 $ CALL ALAERH( PATH, 'DGELSS', INFO, 0, ' ', M, 713 $ N, NRHS, -1, NB, ITYPE, NFAIL, 714 $ NERRS, NOUT ) 715* 716* workspace used: 3*min(m,n) + 717* max(2*min(m,n),nrhs,max(m,n)) 718* 719* Test 7: Compute relative error in svd 720* 721 IF( RANK.GT.0 ) THEN 722 CALL DAXPY( MNMIN, -ONE, COPYS, 1, S, 1 ) 723 RESULT( 7 ) = DASUM( MNMIN, S, 1 ) / 724 $ DASUM( MNMIN, COPYS, 1 ) / 725 $ ( EPS*DBLE( MNMIN ) ) 726 ELSE 727 RESULT( 7 ) = ZERO 728 END IF 729* 730* Test 8: Compute error in solution 731* 732 CALL DLACPY( 'Full', M, NRHS, COPYB, LDB, WORK, 733 $ LDWORK ) 734 CALL DQRT16( 'No transpose', M, N, NRHS, COPYA, 735 $ LDA, B, LDB, WORK, LDWORK, 736 $ WORK( M*NRHS+1 ), RESULT( 8 ) ) 737* 738* Test 9: Check norm of r'*A 739* 740 RESULT( 9 ) = ZERO 741 IF( M.GT.CRANK ) 742 $ RESULT( 9 ) = DQRT17( 'No transpose', 1, M, 743 $ N, NRHS, COPYA, LDA, B, LDB, 744 $ COPYB, LDB, C, WORK, LWORK ) 745* 746* Test 10: Check if x is in the rowspace of A 747* 748 RESULT( 10 ) = ZERO 749 IF( N.GT.CRANK ) 750 $ RESULT( 10 ) = DQRT14( 'No transpose', M, N, 751 $ NRHS, COPYA, LDA, B, LDB, 752 $ WORK, LWORK ) 753* 754* Test DGELSD 755* 756* DGELSD: Compute the minimum-norm solution X 757* to min( norm( A * X - B ) ) using a 758* divide and conquer SVD. 759* 760* Initialize vector IWORK. 761* 762 DO 80 J = 1, N 763 IWORK( J ) = 0 764 80 CONTINUE 765* 766 CALL DLACPY( 'Full', M, N, COPYA, LDA, A, LDA ) 767 CALL DLACPY( 'Full', M, NRHS, COPYB, LDB, B, 768 $ LDB ) 769* 770 SRNAMT = 'DGELSD' 771 CALL DGELSD( M, N, NRHS, A, LDA, B, LDB, S, 772 $ RCOND, CRANK, WORK, LWORK, IWORK, 773 $ INFO ) 774 IF( INFO.NE.0 ) 775 $ CALL ALAERH( PATH, 'DGELSD', INFO, 0, ' ', M, 776 $ N, NRHS, -1, NB, ITYPE, NFAIL, 777 $ NERRS, NOUT ) 778* 779* Test 11: Compute relative error in svd 780* 781 IF( RANK.GT.0 ) THEN 782 CALL DAXPY( MNMIN, -ONE, COPYS, 1, S, 1 ) 783 RESULT( 11 ) = DASUM( MNMIN, S, 1 ) / 784 $ DASUM( MNMIN, COPYS, 1 ) / 785 $ ( EPS*DBLE( MNMIN ) ) 786 ELSE 787 RESULT( 11 ) = ZERO 788 END IF 789* 790* Test 12: Compute error in solution 791* 792 CALL DLACPY( 'Full', M, NRHS, COPYB, LDB, WORK, 793 $ LDWORK ) 794 CALL DQRT16( 'No transpose', M, N, NRHS, COPYA, 795 $ LDA, B, LDB, WORK, LDWORK, 796 $ WORK( M*NRHS+1 ), RESULT( 12 ) ) 797* 798* Test 13: Check norm of r'*A 799* 800 RESULT( 13 ) = ZERO 801 IF( M.GT.CRANK ) 802 $ RESULT( 13 ) = DQRT17( 'No transpose', 1, M, 803 $ N, NRHS, COPYA, LDA, B, LDB, 804 $ COPYB, LDB, C, WORK, LWORK ) 805* 806* Test 14: Check if x is in the rowspace of A 807* 808 RESULT( 14 ) = ZERO 809 IF( N.GT.CRANK ) 810 $ RESULT( 14 ) = DQRT14( 'No transpose', M, N, 811 $ NRHS, COPYA, LDA, B, LDB, 812 $ WORK, LWORK ) 813* 814* Print information about the tests that did not 815* pass the threshold. 816* 817 DO 90 K = 3, 14 818 IF( RESULT( K ).GE.THRESH ) THEN 819 IF( NFAIL.EQ.0 .AND. NERRS.EQ.0 ) 820 $ CALL ALAHD( NOUT, PATH ) 821 WRITE( NOUT, FMT = 9998 )M, N, NRHS, NB, 822 $ ITYPE, K, RESULT( K ) 823 NFAIL = NFAIL + 1 824 END IF 825 90 CONTINUE 826 NRUN = NRUN + 12 827* 828 100 CONTINUE 829 110 CONTINUE 830 120 CONTINUE 831 130 CONTINUE 832 140 CONTINUE 833 150 CONTINUE 834* 835* Print a summary of the results. 836* 837 CALL ALASVM( PATH, NOUT, NFAIL, NRUN, NERRS ) 838* 839 9999 FORMAT( ' TRANS=''', A1, ''', M=', I5, ', N=', I5, ', NRHS=', I4, 840 $ ', NB=', I4, ', type', I2, ', test(', I2, ')=', G12.5 ) 841 9998 FORMAT( ' M=', I5, ', N=', I5, ', NRHS=', I4, ', NB=', I4, 842 $ ', type', I2, ', test(', I2, ')=', G12.5 ) 843 9997 FORMAT( ' TRANS=''', A1,' M=', I5, ', N=', I5, ', NRHS=', I4, 844 $ ', MB=', I4,', NB=', I4,', type', I2, 845 $ ', test(', I2, ')=', G12.5 ) 846* 847 DEALLOCATE( WORK ) 848 DEALLOCATE( IWORK ) 849 RETURN 850* 851* End of DDRVLS 852* 853 END 854