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