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, WORK, IWORK, 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 IWORK( * ), MVAL( * ), NBVAL( * ), NSVAL( * ), 23* $ NVAL( * ), NXVAL( * ) 24* DOUBLE PRECISION A( * ), B( * ), C( * ), COPYA( * ), COPYB( * ), 25* $ COPYS( * ), S( * ), WORK( * ) 26* .. 27* 28* 29*> \par Purpose: 30* ============= 31*> 32*> \verbatim 33*> 34*> DDRVLS tests the least squares driver routines DGELS, DGELSS, DGELSX, 35*> DGELSY 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[out] WORK 173*> \verbatim 174*> WORK is DOUBLE PRECISION array, 175*> dimension (MMAX*NMAX + 4*NMAX + MMAX). 176*> \endverbatim 177*> 178*> \param[out] IWORK 179*> \verbatim 180*> IWORK is INTEGER array, dimension (15*NMAX) 181*> \endverbatim 182*> 183*> \param[in] NOUT 184*> \verbatim 185*> NOUT is INTEGER 186*> The unit number for output. 187*> \endverbatim 188* 189* Authors: 190* ======== 191* 192*> \author Univ. of Tennessee 193*> \author Univ. of California Berkeley 194*> \author Univ. of Colorado Denver 195*> \author NAG Ltd. 196* 197*> \date November 2011 198* 199*> \ingroup double_lin 200* 201* ===================================================================== 202 SUBROUTINE DDRVLS( DOTYPE, NM, MVAL, NN, NVAL, NNS, NSVAL, NNB, 203 $ NBVAL, NXVAL, THRESH, TSTERR, A, COPYA, B, 204 $ COPYB, C, S, COPYS, WORK, IWORK, NOUT ) 205* 206* -- LAPACK test routine (version 3.4.0) -- 207* -- LAPACK is a software package provided by Univ. of Tennessee, -- 208* -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..-- 209* November 2011 210* 211* .. Scalar Arguments .. 212 LOGICAL TSTERR 213 INTEGER NM, NN, NNB, NNS, NOUT 214 DOUBLE PRECISION THRESH 215* .. 216* .. Array Arguments .. 217 LOGICAL DOTYPE( * ) 218 INTEGER IWORK( * ), MVAL( * ), NBVAL( * ), NSVAL( * ), 219 $ NVAL( * ), NXVAL( * ) 220 DOUBLE PRECISION A( * ), B( * ), C( * ), COPYA( * ), COPYB( * ), 221 $ COPYS( * ), S( * ), WORK( * ) 222* .. 223* 224* ===================================================================== 225* 226* .. Parameters .. 227 INTEGER NTESTS 228 PARAMETER ( NTESTS = 18 ) 229 INTEGER SMLSIZ 230 PARAMETER ( SMLSIZ = 25 ) 231 DOUBLE PRECISION ONE, TWO, ZERO 232 PARAMETER ( ONE = 1.0D0, TWO = 2.0D0, ZERO = 0.0D0 ) 233* .. 234* .. Local Scalars .. 235 CHARACTER TRANS 236 CHARACTER*3 PATH 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 242* .. 243* .. Local Arrays .. 244 INTEGER ISEED( 4 ), ISEEDY( 4 ) 245 DOUBLE PRECISION RESULT( NTESTS ) 246* .. 247* .. External Functions .. 248 DOUBLE PRECISION DASUM, DLAMCH, DQRT12, DQRT14, DQRT17 249 EXTERNAL DASUM, DLAMCH, DQRT12, DQRT14, DQRT17 250* .. 251* .. External Subroutines .. 252 EXTERNAL ALAERH, ALAHD, ALASVM, DAXPY, DERRLS, DGELS, 253 $ DGELSD, DGELSS, DGELSX, DGELSY, DGEMM, DLACPY, 254 $ DLARNV, DLASRT, DQRT13, DQRT15, DQRT16, DSCAL, 255 $ XLAENV 256* .. 257* .. Intrinsic Functions .. 258 INTRINSIC DBLE, INT, LOG, MAX, MIN, SQRT 259* .. 260* .. Scalars in Common .. 261 LOGICAL LERR, OK 262 CHARACTER*32 SRNAMT 263 INTEGER INFOT, IOUNIT 264* .. 265* .. Common blocks .. 266 COMMON / INFOC / INFOT, IOUNIT, OK, LERR 267 COMMON / SRNAMC / SRNAMT 268* .. 269* .. Data statements .. 270 DATA ISEEDY / 1988, 1989, 1990, 1991 / 271* .. 272* .. Executable Statements .. 273* 274* Initialize constants and the random number seed. 275* 276 PATH( 1: 1 ) = 'Double precision' 277 PATH( 2: 3 ) = 'LS' 278 NRUN = 0 279 NFAIL = 0 280 NERRS = 0 281 DO 10 I = 1, 4 282 ISEED( I ) = ISEEDY( I ) 283 10 CONTINUE 284 EPS = DLAMCH( 'Epsilon' ) 285* 286* Threshold for rank estimation 287* 288 RCOND = SQRT( EPS ) - ( SQRT( EPS )-EPS ) / 2 289* 290* Test the error exits 291* 292 CALL XLAENV( 2, 2 ) 293 CALL XLAENV( 9, SMLSIZ ) 294 IF( TSTERR ) 295 $ CALL DERRLS( PATH, NOUT ) 296* 297* Print the header if NM = 0 or NN = 0 and THRESH = 0. 298* 299 IF( ( NM.EQ.0 .OR. NN.EQ.0 ) .AND. THRESH.EQ.ZERO ) 300 $ CALL ALAHD( NOUT, PATH ) 301 INFOT = 0 302 CALL XLAENV( 2, 2 ) 303 CALL XLAENV( 9, SMLSIZ ) 304* 305 DO 150 IM = 1, NM 306 M = MVAL( IM ) 307 LDA = MAX( 1, M ) 308* 309 DO 140 IN = 1, NN 310 N = NVAL( IN ) 311 MNMIN = MIN( M, N ) 312 LDB = MAX( 1, M, N ) 313* 314 DO 130 INS = 1, NNS 315 NRHS = NSVAL( INS ) 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 ) 321* 322 DO 120 IRANK = 1, 2 323 DO 110 ISCALE = 1, 3 324 ITYPE = ( IRANK-1 )*3 + ISCALE 325 IF( .NOT.DOTYPE( ITYPE ) ) 326 $ GO TO 110 327* 328 IF( IRANK.EQ.1 ) THEN 329* 330* Test DGELS 331* 332* Generate a matrix of scaling type ISCALE 333* 334 CALL DQRT13( ISCALE, M, N, COPYA, LDA, NORMA, 335 $ ISEED ) 336 DO 40 INB = 1, NNB 337 NB = NBVAL( INB ) 338 CALL XLAENV( 1, NB ) 339 CALL XLAENV( 3, NXVAL( INB ) ) 340* 341 DO 30 ITRAN = 1, 2 342 IF( ITRAN.EQ.1 ) THEN 343 TRANS = 'N' 344 NROWS = M 345 NCOLS = N 346 ELSE 347 TRANS = 'T' 348 NROWS = N 349 NCOLS = M 350 END IF 351 LDWORK = MAX( 1, NCOLS ) 352* 353* Set up a consistent rhs 354* 355 IF( NCOLS.GT.0 ) THEN 356 CALL DLARNV( 2, ISEED, NCOLS*NRHS, 357 $ WORK ) 358 CALL DSCAL( NCOLS*NRHS, 359 $ ONE / DBLE( NCOLS ), WORK, 360 $ 1 ) 361 END IF 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, 366 $ COPYB, LDB ) 367* 368* Solve LS or overdetermined system 369* 370 IF( M.GT.0 .AND. N.GT.0 ) THEN 371 CALL DLACPY( 'Full', M, N, COPYA, LDA, 372 $ A, LDA ) 373 CALL DLACPY( 'Full', NROWS, NRHS, 374 $ COPYB, LDB, B, LDB ) 375 END IF 376 SRNAMT = 'DGELS ' 377 CALL DGELS( TRANS, M, N, NRHS, A, LDA, B, 378 $ LDB, WORK, LWORK, INFO ) 379 IF( INFO.NE.0 ) 380 $ CALL ALAERH( PATH, 'DGELS ', INFO, 0, 381 $ TRANS, M, N, NRHS, -1, NB, 382 $ ITYPE, NFAIL, NERRS, 383 $ NOUT ) 384* 385* Check correctness of results 386* 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, 393 $ RESULT( 1 ) ) 394* 395 IF( ( ITRAN.EQ.1 .AND. M.GE.N ) .OR. 396 $ ( ITRAN.EQ.2 .AND. M.LT.N ) ) THEN 397* 398* Solving LS system 399* 400 RESULT( 2 ) = DQRT17( TRANS, 1, M, N, 401 $ NRHS, COPYA, LDA, B, LDB, 402 $ COPYB, LDB, C, WORK, 403 $ LWORK ) 404 ELSE 405* 406* Solving overdetermined system 407* 408 RESULT( 2 ) = DQRT14( TRANS, M, N, 409 $ NRHS, COPYA, LDA, B, LDB, 410 $ WORK, LWORK ) 411 END IF 412* 413* Print information about the tests that 414* did not pass the threshold. 415* 416 DO 20 K = 1, 2 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, 422 $ RESULT( K ) 423 NFAIL = NFAIL + 1 424 END IF 425 20 CONTINUE 426 NRUN = NRUN + 2 427 30 CONTINUE 428 40 CONTINUE 429 END IF 430* 431* Generate a matrix of scaling type ISCALE and rank 432* type IRANK. 433* 434 CALL DQRT15( ISCALE, IRANK, M, N, NRHS, COPYA, LDA, 435 $ COPYB, LDB, COPYS, RANK, NORMA, NORMB, 436 $ ISEED, WORK, LWORK ) 437* 438* workspace used: MAX(M+MIN(M,N),NRHS*MIN(M,N),2*N+M) 439* 440* Initialize vector IWORK. 441* 442 DO 50 J = 1, N 443 IWORK( J ) = 0 444 50 CONTINUE 445 LDWORK = MAX( 1, M ) 446* 447* Test DGELSX 448* 449* DGELSX: Compute the minimum-norm solution X 450* to min( norm( A * X - B ) ) using a complete 451* orthogonal factorization. 452* 453 CALL DLACPY( 'Full', M, N, COPYA, LDA, A, LDA ) 454 CALL DLACPY( 'Full', M, NRHS, COPYB, LDB, B, LDB ) 455* 456 SRNAMT = 'DGELSX' 457 CALL DGELSX( M, N, NRHS, A, LDA, B, LDB, IWORK, 458 $ RCOND, CRANK, WORK, INFO ) 459 IF( INFO.NE.0 ) 460 $ CALL ALAERH( PATH, 'DGELSX', INFO, 0, ' ', M, N, 461 $ NRHS, -1, NB, ITYPE, NFAIL, NERRS, 462 $ NOUT ) 463* 464* workspace used: MAX( MNMIN+3*N, 2*MNMIN+NRHS ) 465* 466* Test 3: Compute relative error in svd 467* workspace: M*N + 4*MIN(M,N) + MAX(M,N) 468* 469 RESULT( 3 ) = DQRT12( CRANK, CRANK, A, LDA, COPYS, 470 $ WORK, LWORK ) 471* 472* Test 4: Compute error in solution 473* workspace: M*NRHS + M 474* 475 CALL DLACPY( 'Full', M, NRHS, COPYB, LDB, WORK, 476 $ LDWORK ) 477 CALL DQRT16( 'No transpose', M, N, NRHS, COPYA, 478 $ LDA, B, LDB, WORK, LDWORK, 479 $ WORK( M*NRHS+1 ), RESULT( 4 ) ) 480* 481* Test 5: Check norm of r'*A 482* workspace: NRHS*(M+N) 483* 484 RESULT( 5 ) = ZERO 485 IF( M.GT.CRANK ) 486 $ RESULT( 5 ) = DQRT17( 'No transpose', 1, M, N, 487 $ NRHS, COPYA, LDA, B, LDB, COPYB, 488 $ LDB, C, WORK, LWORK ) 489* 490* Test 6: Check if x is in the rowspace of A 491* workspace: (M+NRHS)*(N+2) 492* 493 RESULT( 6 ) = ZERO 494* 495 IF( N.GT.CRANK ) 496 $ RESULT( 6 ) = DQRT14( 'No transpose', M, N, 497 $ NRHS, COPYA, LDA, B, LDB, WORK, 498 $ LWORK ) 499* 500* Print information about the tests that did not 501* pass the threshold. 502* 503 DO 60 K = 3, 6 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 ) 509 NFAIL = NFAIL + 1 510 END IF 511 60 CONTINUE 512 NRUN = NRUN + 4 513* 514* Loop for testing different block sizes. 515* 516 DO 100 INB = 1, NNB 517 NB = NBVAL( INB ) 518 CALL XLAENV( 1, NB ) 519 CALL XLAENV( 3, NXVAL( INB ) ) 520* 521* Test DGELSY 522* 523* DGELSY: Compute the minimum-norm solution X 524* to min( norm( A * X - B ) ) 525* using the rank-revealing orthogonal 526* factorization. 527* 528* Initialize vector IWORK. 529* 530 DO 70 J = 1, N 531 IWORK( J ) = 0 532 70 CONTINUE 533* 534* Set LWLSY to the adequate value. 535* 536 LWLSY = MAX( 1, MNMIN+2*N+NB*( N+1 ), 537 $ 2*MNMIN+NB*NRHS ) 538* 539 CALL DLACPY( 'Full', M, N, COPYA, LDA, A, LDA ) 540 CALL DLACPY( 'Full', M, NRHS, COPYB, LDB, B, 541 $ LDB ) 542* 543 SRNAMT = 'DGELSY' 544 CALL DGELSY( M, N, NRHS, A, LDA, B, LDB, IWORK, 545 $ RCOND, CRANK, WORK, LWLSY, INFO ) 546 IF( INFO.NE.0 ) 547 $ CALL ALAERH( PATH, 'DGELSY', INFO, 0, ' ', M, 548 $ N, NRHS, -1, NB, ITYPE, NFAIL, 549 $ NERRS, NOUT ) 550* 551* Test 7: Compute relative error in svd 552* workspace: M*N + 4*MIN(M,N) + MAX(M,N) 553* 554 RESULT( 7 ) = DQRT12( CRANK, CRANK, A, LDA, 555 $ COPYS, WORK, LWORK ) 556* 557* Test 8: Compute error in solution 558* workspace: M*NRHS + M 559* 560 CALL DLACPY( 'Full', M, NRHS, COPYB, LDB, WORK, 561 $ LDWORK ) 562 CALL DQRT16( 'No transpose', M, N, NRHS, COPYA, 563 $ LDA, B, LDB, WORK, LDWORK, 564 $ WORK( M*NRHS+1 ), RESULT( 8 ) ) 565* 566* Test 9: Check norm of r'*A 567* workspace: NRHS*(M+N) 568* 569 RESULT( 9 ) = ZERO 570 IF( M.GT.CRANK ) 571 $ RESULT( 9 ) = DQRT17( 'No transpose', 1, M, 572 $ N, NRHS, COPYA, LDA, B, LDB, 573 $ COPYB, LDB, C, WORK, LWORK ) 574* 575* Test 10: Check if x is in the rowspace of A 576* workspace: (M+NRHS)*(N+2) 577* 578 RESULT( 10 ) = ZERO 579* 580 IF( N.GT.CRANK ) 581 $ RESULT( 10 ) = DQRT14( 'No transpose', M, N, 582 $ NRHS, COPYA, LDA, B, LDB, 583 $ WORK, LWORK ) 584* 585* Test DGELSS 586* 587* DGELSS: Compute the minimum-norm solution X 588* to min( norm( A * X - B ) ) 589* using the SVD. 590* 591 CALL DLACPY( 'Full', M, N, COPYA, LDA, A, LDA ) 592 CALL DLACPY( 'Full', M, NRHS, COPYB, LDB, B, 593 $ LDB ) 594 SRNAMT = 'DGELSS' 595 CALL DGELSS( M, N, NRHS, A, LDA, B, LDB, S, 596 $ RCOND, CRANK, WORK, LWORK, INFO ) 597 IF( INFO.NE.0 ) 598 $ CALL ALAERH( PATH, 'DGELSS', INFO, 0, ' ', M, 599 $ N, NRHS, -1, NB, ITYPE, NFAIL, 600 $ NERRS, NOUT ) 601* 602* workspace used: 3*min(m,n) + 603* max(2*min(m,n),nrhs,max(m,n)) 604* 605* Test 11: Compute relative error in svd 606* 607 IF( RANK.GT.0 ) THEN 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 ) ) 612 ELSE 613 RESULT( 11 ) = ZERO 614 END IF 615* 616* Test 12: Compute error in solution 617* 618 CALL DLACPY( 'Full', M, NRHS, COPYB, LDB, WORK, 619 $ LDWORK ) 620 CALL DQRT16( 'No transpose', M, N, NRHS, COPYA, 621 $ LDA, B, LDB, WORK, LDWORK, 622 $ WORK( M*NRHS+1 ), RESULT( 12 ) ) 623* 624* Test 13: Check norm of r'*A 625* 626 RESULT( 13 ) = ZERO 627 IF( M.GT.CRANK ) 628 $ RESULT( 13 ) = DQRT17( 'No transpose', 1, M, 629 $ N, NRHS, COPYA, LDA, B, LDB, 630 $ COPYB, LDB, C, WORK, LWORK ) 631* 632* Test 14: Check if x is in the rowspace of A 633* 634 RESULT( 14 ) = ZERO 635 IF( N.GT.CRANK ) 636 $ RESULT( 14 ) = DQRT14( 'No transpose', M, N, 637 $ NRHS, COPYA, LDA, B, LDB, 638 $ WORK, LWORK ) 639* 640* Test DGELSD 641* 642* DGELSD: Compute the minimum-norm solution X 643* to min( norm( A * X - B ) ) using a 644* divide and conquer SVD. 645* 646* Initialize vector IWORK. 647* 648 DO 80 J = 1, N 649 IWORK( J ) = 0 650 80 CONTINUE 651* 652 CALL DLACPY( 'Full', M, N, COPYA, LDA, A, LDA ) 653 CALL DLACPY( 'Full', M, NRHS, COPYB, LDB, B, 654 $ LDB ) 655* 656 SRNAMT = 'DGELSD' 657 CALL DGELSD( M, N, NRHS, A, LDA, B, LDB, S, 658 $ RCOND, CRANK, WORK, LWORK, IWORK, 659 $ INFO ) 660 IF( INFO.NE.0 ) 661 $ CALL ALAERH( PATH, 'DGELSD', INFO, 0, ' ', M, 662 $ N, NRHS, -1, NB, ITYPE, NFAIL, 663 $ NERRS, NOUT ) 664* 665* Test 15: Compute relative error in svd 666* 667 IF( RANK.GT.0 ) THEN 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 ) ) 672 ELSE 673 RESULT( 15 ) = ZERO 674 END IF 675* 676* Test 16: Compute error in solution 677* 678 CALL DLACPY( 'Full', M, NRHS, COPYB, LDB, WORK, 679 $ LDWORK ) 680 CALL DQRT16( 'No transpose', M, N, NRHS, COPYA, 681 $ LDA, B, LDB, WORK, LDWORK, 682 $ WORK( M*NRHS+1 ), RESULT( 16 ) ) 683* 684* Test 17: Check norm of r'*A 685* 686 RESULT( 17 ) = ZERO 687 IF( M.GT.CRANK ) 688 $ RESULT( 17 ) = DQRT17( 'No transpose', 1, M, 689 $ N, NRHS, COPYA, LDA, B, LDB, 690 $ COPYB, LDB, C, WORK, LWORK ) 691* 692* Test 18: Check if x is in the rowspace of A 693* 694 RESULT( 18 ) = ZERO 695 IF( N.GT.CRANK ) 696 $ RESULT( 18 ) = DQRT14( 'No transpose', M, N, 697 $ NRHS, COPYA, LDA, B, LDB, 698 $ WORK, LWORK ) 699* 700* Print information about the tests that did not 701* pass the threshold. 702* 703 DO 90 K = 7, NTESTS 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 ) 709 NFAIL = NFAIL + 1 710 END IF 711 90 CONTINUE 712 NRUN = NRUN + 12 713* 714 100 CONTINUE 715 110 CONTINUE 716 120 CONTINUE 717 130 CONTINUE 718 140 CONTINUE 719 150 CONTINUE 720* 721* Print a summary of the results. 722* 723 CALL ALASVM( PATH, NOUT, NFAIL, NRUN, NERRS ) 724* 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 ) 729 RETURN 730* 731* End of DDRVLS 732* 733 END 734