LAPACK 3.11.0 LAPACK: Linear Algebra PACKage
Searching...
No Matches

## ◆ ddrvls()

 subroutine ddrvls ( logical, dimension( * ) DOTYPE, integer NM, integer, dimension( * ) MVAL, integer NN, integer, dimension( * ) NVAL, integer NNS, integer, dimension( * ) NSVAL, integer NNB, integer, dimension( * ) NBVAL, integer, dimension( * ) NXVAL, double precision THRESH, logical TSTERR, double precision, dimension( * ) A, double precision, dimension( * ) COPYA, double precision, dimension( * ) B, double precision, dimension( * ) COPYB, double precision, dimension( * ) C, double precision, dimension( * ) S, double precision, dimension( * ) COPYS, integer NOUT )

DDRVLS

Purpose:
``` DDRVLS tests the least squares driver routines DGELS, DGELST,
DGETSLS, DGELSS, DGELSY, and DGELSD.```
Parameters
 [in] DOTYPE ``` DOTYPE is LOGICAL array, dimension (NTYPES) The matrix types to be used for testing. Matrices of type j (for 1 <= j <= NTYPES) are used for testing if DOTYPE(j) = .TRUE.; if DOTYPE(j) = .FALSE., then type j is not used. The matrix of type j is generated as follows: j=1: A = U*D*V where U and V are random orthogonal matrices and D has random entries (> 0.1) taken from a uniform distribution (0,1). A is full rank. j=2: The same of 1, but A is scaled up. j=3: The same of 1, but A is scaled down. j=4: A = U*D*V where U and V are random orthogonal matrices and D has 3*min(M,N)/4 random entries (> 0.1) taken from a uniform distribution (0,1) and the remaining entries set to 0. A is rank-deficient. j=5: The same of 4, but A is scaled up. j=6: The same of 5, but A is scaled down.``` [in] NM ``` NM is INTEGER The number of values of M contained in the vector MVAL.``` [in] MVAL ``` MVAL is INTEGER array, dimension (NM) The values of the matrix row dimension M.``` [in] NN ``` NN is INTEGER The number of values of N contained in the vector NVAL.``` [in] NVAL ``` NVAL is INTEGER array, dimension (NN) The values of the matrix column dimension N.``` [in] NNS ``` NNS is INTEGER The number of values of NRHS contained in the vector NSVAL.``` [in] NSVAL ``` NSVAL is INTEGER array, dimension (NNS) The values of the number of right hand sides NRHS.``` [in] NNB ``` NNB is INTEGER The number of values of NB and NX contained in the vectors NBVAL and NXVAL. The blocking parameters are used in pairs (NB,NX).``` [in] NBVAL ``` NBVAL is INTEGER array, dimension (NNB) The values of the blocksize NB.``` [in] NXVAL ``` NXVAL is INTEGER array, dimension (NNB) The values of the crossover point NX.``` [in] THRESH ``` THRESH is DOUBLE PRECISION The threshold value for the test ratios. A result is included in the output file if RESULT >= THRESH. To have every test ratio printed, use THRESH = 0.``` [in] TSTERR ``` TSTERR is LOGICAL Flag that indicates whether error exits are to be tested.``` [out] A ``` A is DOUBLE PRECISION array, dimension (MMAX*NMAX) where MMAX is the maximum value of M in MVAL and NMAX is the maximum value of N in NVAL.``` [out] COPYA ` COPYA is DOUBLE PRECISION array, dimension (MMAX*NMAX)` [out] B ``` B is DOUBLE PRECISION array, dimension (MMAX*NSMAX) where MMAX is the maximum value of M in MVAL and NSMAX is the maximum value of NRHS in NSVAL.``` [out] COPYB ` COPYB is DOUBLE PRECISION array, dimension (MMAX*NSMAX)` [out] C ` C is DOUBLE PRECISION array, dimension (MMAX*NSMAX)` [out] S ``` S is DOUBLE PRECISION array, dimension (min(MMAX,NMAX))``` [out] COPYS ``` COPYS is DOUBLE PRECISION array, dimension (min(MMAX,NMAX))``` [in] NOUT ``` NOUT is INTEGER The unit number for output.```

Definition at line 189 of file ddrvls.f.

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 = 18 )
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_DGELST, LWORK_DGETSLS,
229 \$ LWORK_DGELSS, 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,
249* ..
250* .. Intrinsic Functions ..
251 INTRINSIC dble, int, 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, ?GELST, ?GETSLS, ?GELSY, ?GELSS and ?GELSD
334* routines.
335*
336 DO im = 1, nm
337 m = mval( im )
338 lda = max( 1, m )
339 DO in = 1, nn
340 n = nval( in )
341 mnmin = max(min( m, n ),1)
342 ldb = max( 1, m, n )
343 DO ins = 1, nns
344 nrhs = nsval( ins )
345 DO irank = 1, 2
346 DO iscale = 1, 3
347 itype = ( irank-1 )*3 + iscale
348 IF( dotype( itype ) ) THEN
349 IF( irank.EQ.1 ) THEN
350 DO itran = 1, 2
351 IF( itran.EQ.1 ) THEN
352 trans = 'N'
353 ELSE
354 trans = 'T'
355 END IF
356*
357* Compute workspace needed for DGELS
358 CALL dgels( trans, m, n, nrhs, a, lda,
359 \$ b, ldb, wq, -1, info )
360 lwork_dgels = int( wq( 1 ) )
361* Compute workspace needed for DGELST
362 CALL dgelst( trans, m, n, nrhs, a, lda,
363 \$ b, ldb, wq, -1, info )
364 lwork_dgelst = int( wq( 1 ) )
365* Compute workspace needed for DGETSLS
366 CALL dgetsls( trans, m, n, nrhs, a, lda,
367 \$ b, ldb, wq, -1, info )
368 lwork_dgetsls = int( wq( 1 ) )
369 ENDDO
370 END IF
371* Compute workspace needed for DGELSY
372 CALL dgelsy( m, n, nrhs, a, lda, b, ldb, iwq,
373 \$ rcond, crank, wq, -1, info )
374 lwork_dgelsy = int( wq( 1 ) )
375* Compute workspace needed for DGELSS
376 CALL dgelss( m, n, nrhs, a, lda, b, ldb, s,
377 \$ rcond, crank, wq, -1 , info )
378 lwork_dgelss = int( wq( 1 ) )
379* Compute workspace needed for DGELSD
380 CALL dgelsd( m, n, nrhs, a, lda, b, ldb, s,
381 \$ rcond, crank, wq, -1, iwq, info )
382 lwork_dgelsd = int( wq( 1 ) )
383* Compute LIWORK workspace needed for DGELSY and DGELSD
384 liwork = max( liwork, n, iwq( 1 ) )
385* Compute LWORK workspace needed for all functions
386 lwork = max( lwork, lwork_dgels, lwork_dgelst,
387 \$ lwork_dgetsls, lwork_dgelsy,
388 \$ lwork_dgelss, lwork_dgelsd )
389 END IF
390 ENDDO
391 ENDDO
392 ENDDO
393 ENDDO
394 ENDDO
395*
396 lwlsy = lwork
397*
398 ALLOCATE( work( lwork ) )
399 ALLOCATE( iwork( liwork ) )
400*
401 DO 150 im = 1, nm
402 m = mval( im )
403 lda = max( 1, m )
404*
405 DO 140 in = 1, nn
406 n = nval( in )
407 mnmin = max(min( m, n ),1)
408 ldb = max( 1, m, n )
409 mb = (mnmin+1)
410*
411 DO 130 ins = 1, nns
412 nrhs = nsval( ins )
413*
414 DO 120 irank = 1, 2
415 DO 110 iscale = 1, 3
416 itype = ( irank-1 )*3 + iscale
417 IF( .NOT.dotype( itype ) )
418 \$ GO TO 110
419* =====================================================
420* Begin test DGELS
421* =====================================================
422 IF( irank.EQ.1 ) THEN
423*
424* Generate a matrix of scaling type ISCALE
425*
426 CALL dqrt13( iscale, m, n, copya, lda, norma,
427 \$ iseed )
428*
429* Loop for testing different block sizes.
430*
431 DO inb = 1, nnb
432 nb = nbval( inb )
433 CALL xlaenv( 1, nb )
434 CALL xlaenv( 3, nxval( inb ) )
435*
436* Loop for testing non-transposed and transposed.
437*
438 DO itran = 1, 2
439 IF( itran.EQ.1 ) THEN
440 trans = 'N'
441 nrows = m
442 ncols = n
443 ELSE
444 trans = 'T'
445 nrows = n
446 ncols = m
447 END IF
448 ldwork = max( 1, ncols )
449*
450* Set up a consistent rhs
451*
452 IF( ncols.GT.0 ) THEN
453 CALL dlarnv( 2, iseed, ncols*nrhs,
454 \$ work )
455 CALL dscal( ncols*nrhs,
456 \$ one / dble( ncols ), work,
457 \$ 1 )
458 END IF
459 CALL dgemm( trans, 'No transpose', nrows,
460 \$ nrhs, ncols, one, copya, lda,
461 \$ work, ldwork, zero, b, ldb )
462 CALL dlacpy( 'Full', nrows, nrhs, b, ldb,
463 \$ copyb, ldb )
464*
465* Solve LS or overdetermined system
466*
467 IF( m.GT.0 .AND. n.GT.0 ) THEN
468 CALL dlacpy( 'Full', m, n, copya, lda,
469 \$ a, lda )
470 CALL dlacpy( 'Full', nrows, nrhs,
471 \$ copyb, ldb, b, ldb )
472 END IF
473 srnamt = 'DGELS '
474 CALL dgels( trans, m, n, nrhs, a, lda, b,
475 \$ ldb, work, lwork, info )
476 IF( info.NE.0 )
477 \$ CALL alaerh( path, 'DGELS ', info, 0,
478 \$ trans, m, n, nrhs, -1, nb,
479 \$ itype, nfail, nerrs,
480 \$ nout )
481*
482* Test 1: Check correctness of results
483* for DGELS, compute the residual:
484* RESID = norm(B - A*X) /
485* / ( max(m,n) * norm(A) * norm(X) * EPS )
486*
487 IF( nrows.GT.0 .AND. nrhs.GT.0 )
488 \$ CALL dlacpy( 'Full', nrows, nrhs,
489 \$ copyb, ldb, c, ldb )
490 CALL dqrt16( trans, m, n, nrhs, copya,
491 \$ lda, b, ldb, c, ldb, work,
492 \$ result( 1 ) )
493*
494* Test 2: Check correctness of results
495* for DGELS.
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, compute:
501* r = norm((B- A*X)**T * A) /
502* / (norm(A)*norm(B)*max(M,N,NRHS)*EPS)
503*
504 result( 2 ) = dqrt17( trans, 1, m, n,
505 \$ nrhs, copya, lda, b, ldb,
506 \$ copyb, ldb, c, work,
507 \$ lwork )
508 ELSE
509*
510* Solving overdetermined system
511*
512 result( 2 ) = dqrt14( trans, m, n,
513 \$ nrhs, copya, lda, b, ldb,
514 \$ work, lwork )
515 END IF
516*
517* Print information about the tests that
518* did not pass the threshold.
519*
520 DO k = 1, 2
521 IF( result( k ).GE.thresh ) THEN
522 IF( nfail.EQ.0 .AND. nerrs.EQ.0 )
523 \$ CALL alahd( nout, path )
524 WRITE( nout, fmt = 9999 ) trans, m,
525 \$ n, nrhs, nb, itype, k,
526 \$ result( k )
527 nfail = nfail + 1
528 END IF
529 END DO
530 nrun = nrun + 2
531 END DO
532 END DO
533 END IF
534* =====================================================
535* End test DGELS
536* =====================================================
537* =====================================================
538* Begin test DGELST
539* =====================================================
540 IF( irank.EQ.1 ) THEN
541*
542* Generate a matrix of scaling type ISCALE
543*
544 CALL dqrt13( iscale, m, n, copya, lda, norma,
545 \$ iseed )
546*
547* Loop for testing different block sizes.
548*
549 DO inb = 1, nnb
550 nb = nbval( inb )
551 CALL xlaenv( 1, nb )
552*
553* Loop for testing non-transposed and transposed.
554*
555 DO itran = 1, 2
556 IF( itran.EQ.1 ) THEN
557 trans = 'N'
558 nrows = m
559 ncols = n
560 ELSE
561 trans = 'T'
562 nrows = n
563 ncols = m
564 END IF
565 ldwork = max( 1, ncols )
566*
567* Set up a consistent rhs
568*
569 IF( ncols.GT.0 ) THEN
570 CALL dlarnv( 2, iseed, ncols*nrhs,
571 \$ work )
572 CALL dscal( ncols*nrhs,
573 \$ one / dble( ncols ), work,
574 \$ 1 )
575 END IF
576 CALL dgemm( trans, 'No transpose', nrows,
577 \$ nrhs, ncols, one, copya, lda,
578 \$ work, ldwork, zero, b, ldb )
579 CALL dlacpy( 'Full', nrows, nrhs, b, ldb,
580 \$ copyb, ldb )
581*
582* Solve LS or overdetermined system
583*
584 IF( m.GT.0 .AND. n.GT.0 ) THEN
585 CALL dlacpy( 'Full', m, n, copya, lda,
586 \$ a, lda )
587 CALL dlacpy( 'Full', nrows, nrhs,
588 \$ copyb, ldb, b, ldb )
589 END IF
590 srnamt = 'DGELST'
591 CALL dgelst( trans, m, n, nrhs, a, lda, b,
592 \$ ldb, work, lwork, info )
593 IF( info.NE.0 )
594 \$ CALL alaerh( path, 'DGELST', info, 0,
595 \$ trans, m, n, nrhs, -1, nb,
596 \$ itype, nfail, nerrs,
597 \$ nout )
598*
599* Test 3: Check correctness of results
600* for DGELST, compute the residual:
601* RESID = norm(B - A*X) /
602* / ( max(m,n) * norm(A) * norm(X) * EPS )
603*
604 IF( nrows.GT.0 .AND. nrhs.GT.0 )
605 \$ CALL dlacpy( 'Full', nrows, nrhs,
606 \$ copyb, ldb, c, ldb )
607 CALL dqrt16( trans, m, n, nrhs, copya,
608 \$ lda, b, ldb, c, ldb, work,
609 \$ result( 3 ) )
610*
611* Test 4: Check correctness of results
612* for DGELST.
613*
614 IF( ( itran.EQ.1 .AND. m.GE.n ) .OR.
615 \$ ( itran.EQ.2 .AND. m.LT.n ) ) THEN
616*
617* Solving LS system, compute:
618* r = norm((B- A*X)**T * A) /
619* / (norm(A)*norm(B)*max(M,N,NRHS)*EPS)
620*
621 result( 4 ) = dqrt17( trans, 1, m, n,
622 \$ nrhs, copya, lda, b, ldb,
623 \$ copyb, ldb, c, work,
624 \$ lwork )
625 ELSE
626*
627* Solving overdetermined system
628*
629 result( 4 ) = dqrt14( trans, m, n,
630 \$ nrhs, copya, lda, b, ldb,
631 \$ work, lwork )
632 END IF
633*
634* Print information about the tests that
635* did not pass the threshold.
636*
637 DO k = 3, 4
638 IF( result( k ).GE.thresh ) THEN
639 IF( nfail.EQ.0 .AND. nerrs.EQ.0 )
640 \$ CALL alahd( nout, path )
641 WRITE( nout, fmt = 9999 ) trans, m,
642 \$ n, nrhs, nb, itype, k,
643 \$ result( k )
644 nfail = nfail + 1
645 END IF
646 END DO
647 nrun = nrun + 2
648 END DO
649 END DO
650 END IF
651* =====================================================
652* End test DGELST
653* =====================================================
654* =====================================================
655* Begin test DGETSLS
656* =====================================================
657 IF( irank.EQ.1 ) THEN
658*
659* Generate a matrix of scaling type ISCALE
660*
661 CALL dqrt13( iscale, m, n, copya, lda, norma,
662 \$ iseed )
663*
664* Loop for testing different block sizes MB.
665*
666 DO imb = 1, nnb
667 mb = nbval( imb )
668 CALL xlaenv( 1, mb )
669*
670* Loop for testing different block sizes NB.
671*
672 DO inb = 1, nnb
673 nb = nbval( inb )
674 CALL xlaenv( 2, nb )
675*
676* Loop for testing non-transposed
677* and transposed.
678*
679 DO itran = 1, 2
680 IF( itran.EQ.1 ) THEN
681 trans = 'N'
682 nrows = m
683 ncols = n
684 ELSE
685 trans = 'T'
686 nrows = n
687 ncols = m
688 END IF
689 ldwork = max( 1, ncols )
690*
691* Set up a consistent rhs
692*
693 IF( ncols.GT.0 ) THEN
694 CALL dlarnv( 2, iseed, ncols*nrhs,
695 \$ work )
696 CALL dscal( ncols*nrhs,
697 \$ one / dble( ncols ),
698 \$ work, 1 )
699 END IF
700 CALL dgemm( trans, 'No transpose',
701 \$ nrows, nrhs, ncols, one,
702 \$ copya, lda, work, ldwork,
703 \$ zero, b, ldb )
704 CALL dlacpy( 'Full', nrows, nrhs,
705 \$ b, ldb, copyb, ldb )
706*
707* Solve LS or overdetermined system
708*
709 IF( m.GT.0 .AND. n.GT.0 ) THEN
710 CALL dlacpy( 'Full', m, n,
711 \$ copya, lda, a, lda )
712 CALL dlacpy( 'Full', nrows, nrhs,
713 \$ copyb, ldb, b, ldb )
714 END IF
715 srnamt = 'DGETSLS'
716 CALL dgetsls( trans, m, n, nrhs,
717 \$ a, lda, b, ldb, work, lwork,
718 \$ info )
719 IF( info.NE.0 )
720 \$ CALL alaerh( path, 'DGETSLS', info,
721 \$ 0, trans, m, n, nrhs,
722 \$ -1, nb, itype, nfail,
723 \$ nerrs, nout )
724*
725* Test 5: Check correctness of results
726* for DGETSLS, compute the residual:
727* RESID = norm(B - A*X) /
728* / ( max(m,n) * norm(A) * norm(X) * EPS )
729*
730 IF( nrows.GT.0 .AND. nrhs.GT.0 )
731 \$ CALL dlacpy( 'Full', nrows, nrhs,
732 \$ copyb, ldb, c, ldb )
733 CALL dqrt16( trans, m, n, nrhs,
734 \$ copya, lda, b, ldb,
735 \$ c, ldb, work,
736 \$ result( 5 ) )
737*
738* Test 6: Check correctness of results
739* for DGETSLS.
740*
741 IF( ( itran.EQ.1 .AND. m.GE.n ) .OR.
742 \$ ( itran.EQ.2 .AND. m.LT.n ) ) THEN
743*
744* Solving LS system, compute:
745* r = norm((B- A*X)**T * A) /
746* / (norm(A)*norm(B)*max(M,N,NRHS)*EPS)
747*
748 result( 6 ) = dqrt17( trans, 1, m,
749 \$ n, nrhs, copya, lda,
750 \$ b, ldb, copyb, ldb,
751 \$ c, work, lwork )
752 ELSE
753*
754* Solving overdetermined system
755*
756 result( 6 ) = dqrt14( trans, m, n,
757 \$ nrhs, copya, lda,
758 \$ b, ldb, work, lwork )
759 END IF
760*
761* Print information about the tests that
762* did not pass the threshold.
763*
764 DO k = 5, 6
765 IF( result( k ).GE.thresh ) THEN
766 IF( nfail.EQ.0 .AND. nerrs.EQ.0 )
767 \$ CALL alahd( nout, path )
768 WRITE( nout, fmt = 9997 ) trans,
769 \$ m, n, nrhs, mb, nb, itype,
770 \$ k, result( k )
771 nfail = nfail + 1
772 END IF
773 END DO
774 nrun = nrun + 2
775 END DO
776 END DO
777 END DO
778 END IF
779* =====================================================
780* End test DGETSLS
781* =====================================================
782*
783* Generate a matrix of scaling type ISCALE and rank
784* type IRANK.
785*
786 CALL dqrt15( iscale, irank, m, n, nrhs, copya, lda,
787 \$ copyb, ldb, copys, rank, norma, normb,
788 \$ iseed, work, lwork )
789*
790* workspace used: MAX(M+MIN(M,N),NRHS*MIN(M,N),2*N+M)
791*
792 ldwork = max( 1, m )
793*
794* Loop for testing different block sizes.
795*
796 DO 100 inb = 1, nnb
797 nb = nbval( inb )
798 CALL xlaenv( 1, nb )
799 CALL xlaenv( 3, nxval( inb ) )
800*
801* Test DGELSY
802*
803* DGELSY: Compute the minimum-norm solution X
804* to min( norm( A * X - B ) )
805* using the rank-revealing orthogonal
806* factorization.
807*
808* Initialize vector IWORK.
809*
810 DO 70 j = 1, n
811 iwork( j ) = 0
812 70 CONTINUE
813*
814 CALL dlacpy( 'Full', m, n, copya, lda, a, lda )
815 CALL dlacpy( 'Full', m, nrhs, copyb, ldb, b,
816 \$ ldb )
817*
818 srnamt = 'DGELSY'
819 CALL dgelsy( m, n, nrhs, a, lda, b, ldb, iwork,
820 \$ rcond, crank, work, lwlsy, info )
821 IF( info.NE.0 )
822 \$ CALL alaerh( path, 'DGELSY', info, 0, ' ', m,
823 \$ n, nrhs, -1, nb, itype, nfail,
824 \$ nerrs, nout )
825*
826* Test 7: Compute relative error in svd
827* workspace: M*N + 4*MIN(M,N) + MAX(M,N)
828*
829 result( 7 ) = dqrt12( crank, crank, a, lda,
830 \$ copys, work, lwork )
831*
832* Test 8: Compute error in solution
833* workspace: M*NRHS + M
834*
835 CALL dlacpy( 'Full', m, nrhs, copyb, ldb, work,
836 \$ ldwork )
837 CALL dqrt16( 'No transpose', m, n, nrhs, copya,
838 \$ lda, b, ldb, work, ldwork,
839 \$ work( m*nrhs+1 ), result( 8 ) )
840*
841* Test 9: Check norm of r'*A
842* workspace: NRHS*(M+N)
843*
844 result( 9 ) = zero
845 IF( m.GT.crank )
846 \$ result( 9 ) = dqrt17( 'No transpose', 1, m,
847 \$ n, nrhs, copya, lda, b, ldb,
848 \$ copyb, ldb, c, work, lwork )
849*
850* Test 10: Check if x is in the rowspace of A
851* workspace: (M+NRHS)*(N+2)
852*
853 result( 10 ) = zero
854*
855 IF( n.GT.crank )
856 \$ result( 10 ) = dqrt14( 'No transpose', m, n,
857 \$ nrhs, copya, lda, b, ldb,
858 \$ work, lwork )
859*
860* Test DGELSS
861*
862* DGELSS: Compute the minimum-norm solution X
863* to min( norm( A * X - B ) )
864* using the SVD.
865*
866 CALL dlacpy( 'Full', m, n, copya, lda, a, lda )
867 CALL dlacpy( 'Full', m, nrhs, copyb, ldb, b,
868 \$ ldb )
869 srnamt = 'DGELSS'
870 CALL dgelss( m, n, nrhs, a, lda, b, ldb, s,
871 \$ rcond, crank, work, lwork, info )
872 IF( info.NE.0 )
873 \$ CALL alaerh( path, 'DGELSS', info, 0, ' ', m,
874 \$ n, nrhs, -1, nb, itype, nfail,
875 \$ nerrs, nout )
876*
877* workspace used: 3*min(m,n) +
878* max(2*min(m,n),nrhs,max(m,n))
879*
880* Test 11: Compute relative error in svd
881*
882 IF( rank.GT.0 ) THEN
883 CALL daxpy( mnmin, -one, copys, 1, s, 1 )
884 result( 11 ) = dasum( mnmin, s, 1 ) /
885 \$ dasum( mnmin, copys, 1 ) /
886 \$ ( eps*dble( mnmin ) )
887 ELSE
888 result( 11 ) = zero
889 END IF
890*
891* Test 12: Compute error in solution
892*
893 CALL dlacpy( 'Full', m, nrhs, copyb, ldb, work,
894 \$ ldwork )
895 CALL dqrt16( 'No transpose', m, n, nrhs, copya,
896 \$ lda, b, ldb, work, ldwork,
897 \$ work( m*nrhs+1 ), result( 12 ) )
898*
899* Test 13: Check norm of r'*A
900*
901 result( 13 ) = zero
902 IF( m.GT.crank )
903 \$ result( 13 ) = dqrt17( 'No transpose', 1, m,
904 \$ n, nrhs, copya, lda, b, ldb,
905 \$ copyb, ldb, c, work, lwork )
906*
907* Test 14: Check if x is in the rowspace of A
908*
909 result( 14 ) = zero
910 IF( n.GT.crank )
911 \$ result( 14 ) = dqrt14( 'No transpose', m, n,
912 \$ nrhs, copya, lda, b, ldb,
913 \$ work, lwork )
914*
915* Test DGELSD
916*
917* DGELSD: Compute the minimum-norm solution X
918* to min( norm( A * X - B ) ) using a
919* divide and conquer SVD.
920*
921* Initialize vector IWORK.
922*
923 DO 80 j = 1, n
924 iwork( j ) = 0
925 80 CONTINUE
926*
927 CALL dlacpy( 'Full', m, n, copya, lda, a, lda )
928 CALL dlacpy( 'Full', m, nrhs, copyb, ldb, b,
929 \$ ldb )
930*
931 srnamt = 'DGELSD'
932 CALL dgelsd( m, n, nrhs, a, lda, b, ldb, s,
933 \$ rcond, crank, work, lwork, iwork,
934 \$ info )
935 IF( info.NE.0 )
936 \$ CALL alaerh( path, 'DGELSD', info, 0, ' ', m,
937 \$ n, nrhs, -1, nb, itype, nfail,
938 \$ nerrs, nout )
939*
940* Test 15: Compute relative error in svd
941*
942 IF( rank.GT.0 ) THEN
943 CALL daxpy( mnmin, -one, copys, 1, s, 1 )
944 result( 15 ) = dasum( mnmin, s, 1 ) /
945 \$ dasum( mnmin, copys, 1 ) /
946 \$ ( eps*dble( mnmin ) )
947 ELSE
948 result( 15 ) = zero
949 END IF
950*
951* Test 16: Compute error in solution
952*
953 CALL dlacpy( 'Full', m, nrhs, copyb, ldb, work,
954 \$ ldwork )
955 CALL dqrt16( 'No transpose', m, n, nrhs, copya,
956 \$ lda, b, ldb, work, ldwork,
957 \$ work( m*nrhs+1 ), result( 16 ) )
958*
959* Test 17: Check norm of r'*A
960*
961 result( 17 ) = zero
962 IF( m.GT.crank )
963 \$ result( 17 ) = dqrt17( 'No transpose', 1, m,
964 \$ n, nrhs, copya, lda, b, ldb,
965 \$ copyb, ldb, c, work, lwork )
966*
967* Test 18: Check if x is in the rowspace of A
968*
969 result( 18 ) = zero
970 IF( n.GT.crank )
971 \$ result( 18 ) = dqrt14( 'No transpose', m, n,
972 \$ nrhs, copya, lda, b, ldb,
973 \$ work, lwork )
974*
975* Print information about the tests that did not
976* pass the threshold.
977*
978 DO 90 k = 7, 18
979 IF( result( k ).GE.thresh ) THEN
980 IF( nfail.EQ.0 .AND. nerrs.EQ.0 )
981 \$ CALL alahd( nout, path )
982 WRITE( nout, fmt = 9998 )m, n, nrhs, nb,
983 \$ itype, k, result( k )
984 nfail = nfail + 1
985 END IF
986 90 CONTINUE
987 nrun = nrun + 12
988*
989 100 CONTINUE
990
991
992
993
994
995
996 110 CONTINUE
997 120 CONTINUE
998 130 CONTINUE
999 140 CONTINUE
1000 150 CONTINUE
1001*
1002* Print a summary of the results.
1003*
1004 CALL alasvm( path, nout, nfail, nrun, nerrs )
1005*
1006 9999 FORMAT( ' TRANS=''', a1, ''', M=', i5, ', N=', i5, ', NRHS=', i4,
1007 \$ ', NB=', i4, ', type', i2, ', test(', i2, ')=', g12.5 )
1008 9998 FORMAT( ' M=', i5, ', N=', i5, ', NRHS=', i4, ', NB=', i4,
1009 \$ ', type', i2, ', test(', i2, ')=', g12.5 )
1010 9997 FORMAT( ' TRANS=''', a1,' M=', i5, ', N=', i5, ', NRHS=', i4,
1011 \$ ', MB=', i4,', NB=', i4,', type', i2,
1012 \$ ', test(', i2, ')=', g12.5 )
1013*
1014 DEALLOCATE( work )
1015 DEALLOCATE( iwork )
1016 RETURN
1017*
1018* End of DDRVLS
1019*
double precision function dlamch(CMACH)
DLAMCH
Definition: dlamch.f:69
subroutine dlarnv(IDIST, ISEED, N, X)
DLARNV returns a vector of random numbers from a uniform or normal distribution.
Definition: dlarnv.f:97
subroutine dlacpy(UPLO, M, N, A, LDA, B, LDB)
DLACPY copies all or part of one two-dimensional array to another.
Definition: dlacpy.f:103
subroutine alasvm(TYPE, NOUT, NFAIL, NRUN, NERRS)
ALASVM
Definition: alasvm.f:73
subroutine xlaenv(ISPEC, NVALUE)
XLAENV
Definition: xlaenv.f:81
subroutine alahd(IOUNIT, PATH)
ALAHD
Definition: alahd.f:107
subroutine alaerh(PATH, SUBNAM, INFO, INFOE, OPTS, M, N, KL, KU, N5, IMAT, NFAIL, NERRS, NOUT)
ALAERH
Definition: alaerh.f:147
double precision function dasum(N, DX, INCX)
DASUM
Definition: dasum.f:71
subroutine dscal(N, DA, DX, INCX)
DSCAL
Definition: dscal.f:79
subroutine daxpy(N, DA, DX, INCX, DY, INCY)
DAXPY
Definition: daxpy.f:89
subroutine dgemm(TRANSA, TRANSB, M, N, K, ALPHA, A, LDA, B, LDB, BETA, C, LDC)
DGEMM
Definition: dgemm.f:187
subroutine derrls(PATH, NUNIT)
DERRLS
Definition: derrls.f:55
subroutine dqrt15(SCALE, RKSEL, M, N, NRHS, A, LDA, B, LDB, S, RANK, NORMA, NORMB, ISEED, WORK, LWORK)
DQRT15
Definition: dqrt15.f:148
subroutine dqrt13(SCALE, M, N, A, LDA, NORMA, ISEED)
DQRT13
Definition: dqrt13.f:91
subroutine dqrt16(TRANS, M, N, NRHS, A, LDA, X, LDX, B, LDB, RWORK, RESID)
DQRT16
Definition: dqrt16.f:133
double precision function dqrt12(M, N, A, LDA, S, WORK, LWORK)
DQRT12
Definition: dqrt12.f:89
double precision function dqrt14(TRANS, M, N, NRHS, A, LDA, X, LDX, WORK, LWORK)
DQRT14
Definition: dqrt14.f:116
double precision function dqrt17(TRANS, IRESID, M, N, NRHS, A, LDA, X, LDX, B, LDB, C, WORK, LWORK)
DQRT17
Definition: dqrt17.f:153
subroutine dgels(TRANS, M, N, NRHS, A, LDA, B, LDB, WORK, LWORK, INFO)
DGELS solves overdetermined or underdetermined systems for GE matrices
Definition: dgels.f:183
subroutine dgetsls(TRANS, M, N, NRHS, A, LDA, B, LDB, WORK, LWORK, INFO)
DGETSLS
Definition: dgetsls.f:162
subroutine dgelst(TRANS, M, N, NRHS, A, LDA, B, LDB, WORK, LWORK, INFO)
DGELST solves overdetermined or underdetermined systems for GE matrices using QR or LQ factorization ...
Definition: dgelst.f:194
subroutine dgelsy(M, N, NRHS, A, LDA, B, LDB, JPVT, RCOND, RANK, WORK, LWORK, INFO)
DGELSY solves overdetermined or underdetermined systems for GE matrices
Definition: dgelsy.f:204
subroutine dgelsd(M, N, NRHS, A, LDA, B, LDB, S, RCOND, RANK, WORK, LWORK, IWORK, INFO)
DGELSD computes the minimum-norm solution to a linear least squares problem for GE matrices
Definition: dgelsd.f:209
subroutine dgelss(M, N, NRHS, A, LDA, B, LDB, S, RCOND, RANK, WORK, LWORK, INFO)
DGELSS solves overdetermined or underdetermined systems for GE matrices
Definition: dgelss.f:172
Here is the call graph for this function:
Here is the caller graph for this function: