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