LAPACK  3.9.0
LAPACK: Linear Algebra PACKage

◆ 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, 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.
Author
Univ. of Tennessee
Univ. of California Berkeley
Univ. of Colorado Denver
NAG Ltd.
Date
June 2017

Definition at line 194 of file ddrvls.f.

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  DOUBLE PRECISION THRESH
204 * ..
205 * .. Array Arguments ..
206  LOGICAL DOTYPE( * )
207  INTEGER MVAL( * ), NBVAL( * ), NSVAL( * ),
208  $ NVAL( * ), NXVAL( * )
209  DOUBLE PRECISION A( * ), B( * ), C( * ), COPYA( * ), COPYB( * ),
210  $ COPYS( * ), S( * )
211 * ..
212 *
213 * =====================================================================
214 *
215 * .. Parameters ..
216  INTEGER NTESTS
217  parameter( ntests = 16 )
218  INTEGER SMLSIZ
219  parameter( smlsiz = 25 )
220  DOUBLE PRECISION ONE, TWO, ZERO
221  parameter( one = 1.0d0, two = 2.0d0, zero = 0.0d0 )
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,
231  $ LWORK_DGELS, LWORK_DGETSLS, LWORK_DGELSS,
232  $ LWORK_DGELSY, LWORK_DGELSD
233  DOUBLE PRECISION EPS, NORMA, NORMB, RCOND
234 * ..
235 * .. Local Arrays ..
236  INTEGER ISEED( 4 ), ISEEDY( 4 ), IWQ( 1 )
237  DOUBLE PRECISION RESULT( NTESTS ), WQ( 1 )
238 * ..
239 * .. Allocatable Arrays ..
240  DOUBLE PRECISION, ALLOCATABLE :: WORK (:)
241  INTEGER, ALLOCATABLE :: IWORK (:)
242 * ..
243 * .. External Functions ..
244  DOUBLE PRECISION DASUM, DLAMCH, DQRT12, DQRT14, DQRT17
245  EXTERNAL dasum, dlamch, dqrt12, dqrt14, dqrt17
246 * ..
247 * .. External Subroutines ..
248  EXTERNAL alaerh, alahd, alasvm, daxpy, derrls, dgels,
249  $ dgelsd, dgelss, dgelsy, dgemm, dlacpy,
250  $ dlarnv, dlasrt, dqrt13, dqrt15, dqrt16, dscal,
251  $ xlaenv
252 * ..
253 * .. Intrinsic Functions ..
254  INTRINSIC dble, int, log, max, min, sqrt
255 * ..
256 * .. Scalars in Common ..
257  LOGICAL LERR, OK
258  CHARACTER*32 SRNAMT
259  INTEGER INFOT, IOUNIT
260 * ..
261 * .. Common blocks ..
262  COMMON / infoc / infot, iounit, ok, lerr
263  COMMON / srnamc / srnamt
264 * ..
265 * .. Data statements ..
266  DATA iseedy / 1988, 1989, 1990, 1991 /
267 * ..
268 * .. Executable Statements ..
269 *
270 * Initialize constants and the random number seed.
271 *
272  path( 1: 1 ) = 'Double precision'
273  path( 2: 3 ) = 'LS'
274  nrun = 0
275  nfail = 0
276  nerrs = 0
277  DO 10 i = 1, 4
278  iseed( i ) = iseedy( i )
279  10 CONTINUE
280  eps = dlamch( 'Epsilon' )
281 *
282 * Threshold for rank estimation
283 *
284  rcond = sqrt( eps ) - ( sqrt( eps )-eps ) / 2
285 *
286 * Test the error exits
287 *
288  CALL xlaenv( 2, 2 )
289  CALL xlaenv( 9, smlsiz )
290  IF( tsterr )
291  $ CALL derrls( path, nout )
292 *
293 * Print the header if NM = 0 or NN = 0 and THRESH = 0.
294 *
295  IF( ( nm.EQ.0 .OR. nn.EQ.0 ) .AND. thresh.EQ.zero )
296  $ CALL alahd( nout, path )
297  infot = 0
298  CALL xlaenv( 2, 2 )
299  CALL xlaenv( 9, smlsiz )
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 * DQRT14, DQRT17 (two side cases), DQRT15 and DQRT12
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  liwork = 1
334 *
335 * Iterate through all test cases and compute necessary workspace
336 * sizes for ?GELS, ?GETSLS, ?GELSY, ?GELSS and ?GELSD routines.
337 *
338  DO im = 1, nm
339  m = mval( im )
340  lda = max( 1, m )
341  DO in = 1, nn
342  n = nval( in )
343  mnmin = max(min( m, n ),1)
344  ldb = max( 1, m, n )
345  DO ins = 1, nns
346  nrhs = nsval( ins )
347  DO irank = 1, 2
348  DO iscale = 1, 3
349  itype = ( irank-1 )*3 + iscale
350  IF( dotype( itype ) ) THEN
351  IF( irank.EQ.1 ) THEN
352  DO itran = 1, 2
353  IF( itran.EQ.1 ) THEN
354  trans = 'N'
355  ELSE
356  trans = 'T'
357  END IF
358 *
359 * Compute workspace needed for DGELS
360  CALL dgels( trans, m, n, nrhs, a, lda,
361  $ b, ldb, wq, -1, info )
362  lwork_dgels = int( wq( 1 ) )
363 * Compute workspace needed for DGETSLS
364  CALL dgetsls( trans, m, n, nrhs, a, lda,
365  $ b, ldb, wq, -1, info )
366  lwork_dgetsls = int( wq( 1 ) )
367  ENDDO
368  END IF
369 * Compute workspace needed for DGELSY
370  CALL dgelsy( m, n, nrhs, a, lda, b, ldb, iwq,
371  $ rcond, crank, wq, -1, info )
372  lwork_dgelsy = int( wq( 1 ) )
373 * Compute workspace needed for DGELSS
374  CALL dgelss( m, n, nrhs, a, lda, b, ldb, s,
375  $ rcond, crank, wq, -1 , info )
376  lwork_dgelss = int( wq( 1 ) )
377 * Compute workspace needed for DGELSD
378  CALL dgelsd( m, n, nrhs, a, lda, b, ldb, s,
379  $ rcond, crank, wq, -1, iwq, info )
380  lwork_dgelsd = int( wq( 1 ) )
381 * Compute LIWORK workspace needed for DGELSY and DGELSD
382  liwork = max( liwork, n, iwq( 1 ) )
383 * Compute LWORK workspace needed for all functions
384  lwork = max( lwork, lwork_dgels, lwork_dgetsls,
385  $ lwork_dgelsy, lwork_dgelss,
386  $ lwork_dgelsd )
387  END IF
388  ENDDO
389  ENDDO
390  ENDDO
391  ENDDO
392  ENDDO
393 *
394  lwlsy = lwork
395 *
396  ALLOCATE( work( lwork ) )
397  ALLOCATE( iwork( liwork ) )
398 *
399  DO 150 im = 1, nm
400  m = mval( im )
401  lda = max( 1, m )
402 *
403  DO 140 in = 1, nn
404  n = nval( in )
405  mnmin = max(min( m, n ),1)
406  ldb = max( 1, m, n )
407  mb = (mnmin+1)
408 *
409  DO 130 ins = 1, nns
410  nrhs = nsval( ins )
411 *
412  DO 120 irank = 1, 2
413  DO 110 iscale = 1, 3
414  itype = ( irank-1 )*3 + iscale
415  IF( .NOT.dotype( itype ) )
416  $ GO TO 110
417 *
418  IF( irank.EQ.1 ) THEN
419 *
420 * Test DGELS
421 *
422 * Generate a matrix of scaling type ISCALE
423 *
424  CALL dqrt13( iscale, m, n, copya, lda, norma,
425  $ iseed )
426  DO 40 inb = 1, nnb
427  nb = nbval( inb )
428  CALL xlaenv( 1, nb )
429  CALL xlaenv( 3, nxval( inb ) )
430 *
431  DO 30 itran = 1, 2
432  IF( itran.EQ.1 ) THEN
433  trans = 'N'
434  nrows = m
435  ncols = n
436  ELSE
437  trans = 'T'
438  nrows = n
439  ncols = m
440  END IF
441  ldwork = max( 1, ncols )
442 *
443 * Set up a consistent rhs
444 *
445  IF( ncols.GT.0 ) THEN
446  CALL dlarnv( 2, iseed, ncols*nrhs,
447  $ work )
448  CALL dscal( ncols*nrhs,
449  $ one / dble( ncols ), work,
450  $ 1 )
451  END IF
452  CALL dgemm( trans, 'No transpose', nrows,
453  $ nrhs, ncols, one, copya, lda,
454  $ work, ldwork, zero, b, ldb )
455  CALL dlacpy( 'Full', nrows, nrhs, b, ldb,
456  $ copyb, ldb )
457 *
458 * Solve LS or overdetermined system
459 *
460  IF( m.GT.0 .AND. n.GT.0 ) THEN
461  CALL dlacpy( 'Full', m, n, copya, lda,
462  $ a, lda )
463  CALL dlacpy( 'Full', nrows, nrhs,
464  $ copyb, ldb, b, ldb )
465  END IF
466  srnamt = 'DGELS '
467  CALL dgels( trans, m, n, nrhs, a, lda, b,
468  $ ldb, work, lwork, info )
469  IF( info.NE.0 )
470  $ CALL alaerh( path, 'DGELS ', info, 0,
471  $ trans, m, n, nrhs, -1, nb,
472  $ itype, nfail, nerrs,
473  $ nout )
474 *
475 * Check correctness of results
476 *
477  ldwork = max( 1, nrows )
478  IF( nrows.GT.0 .AND. nrhs.GT.0 )
479  $ CALL dlacpy( 'Full', nrows, nrhs,
480  $ copyb, ldb, c, ldb )
481  CALL dqrt16( trans, m, n, nrhs, copya,
482  $ lda, b, ldb, c, ldb, work,
483  $ result( 1 ) )
484 *
485  IF( ( itran.EQ.1 .AND. m.GE.n ) .OR.
486  $ ( itran.EQ.2 .AND. m.LT.n ) ) THEN
487 *
488 * Solving LS system
489 *
490  result( 2 ) = dqrt17( trans, 1, m, n,
491  $ nrhs, copya, lda, b, ldb,
492  $ copyb, ldb, c, work,
493  $ lwork )
494  ELSE
495 *
496 * Solving overdetermined system
497 *
498  result( 2 ) = dqrt14( trans, m, n,
499  $ nrhs, copya, lda, b, ldb,
500  $ work, lwork )
501  END IF
502 *
503 * Print information about the tests that
504 * did not pass the threshold.
505 *
506  DO 20 k = 1, 2
507  IF( result( k ).GE.thresh ) THEN
508  IF( nfail.EQ.0 .AND. nerrs.EQ.0 )
509  $ CALL alahd( nout, path )
510  WRITE( nout, fmt = 9999 )trans, m,
511  $ n, nrhs, nb, itype, k,
512  $ result( k )
513  nfail = nfail + 1
514  END IF
515  20 CONTINUE
516  nrun = nrun + 2
517  30 CONTINUE
518  40 CONTINUE
519 *
520 *
521 * Test DGETSLS
522 *
523 * Generate a matrix of scaling type ISCALE
524 *
525  CALL dqrt13( iscale, m, n, copya, lda, norma,
526  $ iseed )
527  DO 65 inb = 1, nnb
528  mb = nbval( inb )
529  CALL xlaenv( 1, mb )
530  DO 62 imb = 1, nnb
531  nb = nbval( imb )
532  CALL xlaenv( 2, nb )
533 *
534  DO 60 itran = 1, 2
535  IF( itran.EQ.1 ) THEN
536  trans = 'N'
537  nrows = m
538  ncols = n
539  ELSE
540  trans = 'T'
541  nrows = n
542  ncols = m
543  END IF
544  ldwork = max( 1, ncols )
545 *
546 * Set up a consistent rhs
547 *
548  IF( ncols.GT.0 ) THEN
549  CALL dlarnv( 2, iseed, ncols*nrhs,
550  $ work )
551  CALL dscal( ncols*nrhs,
552  $ one / dble( ncols ), work,
553  $ 1 )
554  END IF
555  CALL dgemm( trans, 'No transpose', nrows,
556  $ nrhs, ncols, one, copya, lda,
557  $ work, ldwork, zero, b, ldb )
558  CALL dlacpy( 'Full', nrows, nrhs, b, ldb,
559  $ copyb, ldb )
560 *
561 * Solve LS or overdetermined system
562 *
563  IF( m.GT.0 .AND. n.GT.0 ) THEN
564  CALL dlacpy( 'Full', m, n, copya, lda,
565  $ a, lda )
566  CALL dlacpy( 'Full', nrows, nrhs,
567  $ copyb, ldb, b, ldb )
568  END IF
569  srnamt = 'DGETSLS '
570  CALL dgetsls( trans, m, n, nrhs, a,
571  $ lda, b, ldb, work, lwork, info )
572  IF( info.NE.0 )
573  $ CALL alaerh( path, 'DGETSLS ', info, 0,
574  $ trans, m, n, nrhs, -1, nb,
575  $ itype, nfail, nerrs,
576  $ nout )
577 *
578 * Check correctness of results
579 *
580  ldwork = max( 1, nrows )
581  IF( nrows.GT.0 .AND. nrhs.GT.0 )
582  $ CALL dlacpy( 'Full', nrows, nrhs,
583  $ copyb, ldb, c, ldb )
584  CALL dqrt16( trans, m, n, nrhs, copya,
585  $ lda, b, ldb, c, ldb, work,
586  $ result( 15 ) )
587 *
588  IF( ( itran.EQ.1 .AND. m.GE.n ) .OR.
589  $ ( itran.EQ.2 .AND. m.LT.n ) ) THEN
590 *
591 * Solving LS system
592 *
593  result( 16 ) = dqrt17( trans, 1, m, n,
594  $ nrhs, copya, lda, b, ldb,
595  $ copyb, ldb, c, work,
596  $ lwork )
597  ELSE
598 *
599 * Solving overdetermined system
600 *
601  result( 16 ) = dqrt14( trans, m, n,
602  $ nrhs, copya, lda, b, ldb,
603  $ work, lwork )
604  END IF
605 *
606 * Print information about the tests that
607 * did not pass the threshold.
608 *
609  DO 50 k = 15, 16
610  IF( result( k ).GE.thresh ) THEN
611  IF( nfail.EQ.0 .AND. nerrs.EQ.0 )
612  $ CALL alahd( nout, path )
613  WRITE( nout, fmt = 9997 )trans, m,
614  $ n, nrhs, mb, nb, itype, k,
615  $ result( k )
616  nfail = nfail + 1
617  END IF
618  50 CONTINUE
619  nrun = nrun + 2
620  60 CONTINUE
621  62 CONTINUE
622  65 CONTINUE
623  END IF
624 *
625 * Generate a matrix of scaling type ISCALE and rank
626 * type IRANK.
627 *
628  CALL dqrt15( iscale, irank, m, n, nrhs, copya, lda,
629  $ copyb, ldb, copys, rank, norma, normb,
630  $ iseed, work, lwork )
631 *
632 * workspace used: MAX(M+MIN(M,N),NRHS*MIN(M,N),2*N+M)
633 *
634  ldwork = max( 1, m )
635 *
636 * Loop for testing different block sizes.
637 *
638  DO 100 inb = 1, nnb
639  nb = nbval( inb )
640  CALL xlaenv( 1, nb )
641  CALL xlaenv( 3, nxval( inb ) )
642 *
643 * Test DGELSY
644 *
645 * DGELSY: Compute the minimum-norm solution X
646 * to min( norm( A * X - B ) )
647 * using the rank-revealing orthogonal
648 * factorization.
649 *
650 * Initialize vector IWORK.
651 *
652  DO 70 j = 1, n
653  iwork( j ) = 0
654  70 CONTINUE
655 *
656  CALL dlacpy( 'Full', m, n, copya, lda, a, lda )
657  CALL dlacpy( 'Full', m, nrhs, copyb, ldb, b,
658  $ ldb )
659 *
660  srnamt = 'DGELSY'
661  CALL dgelsy( m, n, nrhs, a, lda, b, ldb, iwork,
662  $ rcond, crank, work, lwlsy, info )
663  IF( info.NE.0 )
664  $ CALL alaerh( path, 'DGELSY', info, 0, ' ', m,
665  $ n, nrhs, -1, nb, itype, nfail,
666  $ nerrs, nout )
667 *
668 * Test 3: Compute relative error in svd
669 * workspace: M*N + 4*MIN(M,N) + MAX(M,N)
670 *
671  result( 3 ) = dqrt12( crank, crank, a, lda,
672  $ copys, work, lwork )
673 *
674 * Test 4: Compute error in solution
675 * workspace: M*NRHS + M
676 *
677  CALL dlacpy( 'Full', m, nrhs, copyb, ldb, work,
678  $ ldwork )
679  CALL dqrt16( 'No transpose', m, n, nrhs, copya,
680  $ lda, b, ldb, work, ldwork,
681  $ work( m*nrhs+1 ), result( 4 ) )
682 *
683 * Test 5: Check norm of r'*A
684 * workspace: NRHS*(M+N)
685 *
686  result( 5 ) = zero
687  IF( m.GT.crank )
688  $ result( 5 ) = dqrt17( 'No transpose', 1, m,
689  $ n, nrhs, copya, lda, b, ldb,
690  $ copyb, ldb, c, work, lwork )
691 *
692 * Test 6: Check if x is in the rowspace of A
693 * workspace: (M+NRHS)*(N+2)
694 *
695  result( 6 ) = zero
696 *
697  IF( n.GT.crank )
698  $ result( 6 ) = dqrt14( 'No transpose', m, n,
699  $ nrhs, copya, lda, b, ldb,
700  $ work, lwork )
701 *
702 * Test DGELSS
703 *
704 * DGELSS: Compute the minimum-norm solution X
705 * to min( norm( A * X - B ) )
706 * using the SVD.
707 *
708  CALL dlacpy( 'Full', m, n, copya, lda, a, lda )
709  CALL dlacpy( 'Full', m, nrhs, copyb, ldb, b,
710  $ ldb )
711  srnamt = 'DGELSS'
712  CALL dgelss( m, n, nrhs, a, lda, b, ldb, s,
713  $ rcond, crank, work, lwork, info )
714  IF( info.NE.0 )
715  $ CALL alaerh( path, 'DGELSS', info, 0, ' ', m,
716  $ n, nrhs, -1, nb, itype, nfail,
717  $ nerrs, nout )
718 *
719 * workspace used: 3*min(m,n) +
720 * max(2*min(m,n),nrhs,max(m,n))
721 *
722 * Test 7: Compute relative error in svd
723 *
724  IF( rank.GT.0 ) THEN
725  CALL daxpy( mnmin, -one, copys, 1, s, 1 )
726  result( 7 ) = dasum( mnmin, s, 1 ) /
727  $ dasum( mnmin, copys, 1 ) /
728  $ ( eps*dble( mnmin ) )
729  ELSE
730  result( 7 ) = zero
731  END IF
732 *
733 * Test 8: Compute error in solution
734 *
735  CALL dlacpy( 'Full', m, nrhs, copyb, ldb, work,
736  $ ldwork )
737  CALL dqrt16( 'No transpose', m, n, nrhs, copya,
738  $ lda, b, ldb, work, ldwork,
739  $ work( m*nrhs+1 ), result( 8 ) )
740 *
741 * Test 9: Check norm of r'*A
742 *
743  result( 9 ) = zero
744  IF( m.GT.crank )
745  $ result( 9 ) = dqrt17( 'No transpose', 1, m,
746  $ n, nrhs, copya, lda, b, ldb,
747  $ copyb, ldb, c, work, lwork )
748 *
749 * Test 10: Check if x is in the rowspace of A
750 *
751  result( 10 ) = zero
752  IF( n.GT.crank )
753  $ result( 10 ) = dqrt14( 'No transpose', m, n,
754  $ nrhs, copya, lda, b, ldb,
755  $ work, lwork )
756 *
757 * Test DGELSD
758 *
759 * DGELSD: Compute the minimum-norm solution X
760 * to min( norm( A * X - B ) ) using a
761 * divide and conquer SVD.
762 *
763 * Initialize vector IWORK.
764 *
765  DO 80 j = 1, n
766  iwork( j ) = 0
767  80 CONTINUE
768 *
769  CALL dlacpy( 'Full', m, n, copya, lda, a, lda )
770  CALL dlacpy( 'Full', m, nrhs, copyb, ldb, b,
771  $ ldb )
772 *
773  srnamt = 'DGELSD'
774  CALL dgelsd( m, n, nrhs, a, lda, b, ldb, s,
775  $ rcond, crank, work, lwork, iwork,
776  $ info )
777  IF( info.NE.0 )
778  $ CALL alaerh( path, 'DGELSD', info, 0, ' ', m,
779  $ n, nrhs, -1, nb, itype, nfail,
780  $ nerrs, nout )
781 *
782 * Test 11: Compute relative error in svd
783 *
784  IF( rank.GT.0 ) THEN
785  CALL daxpy( mnmin, -one, copys, 1, s, 1 )
786  result( 11 ) = dasum( mnmin, s, 1 ) /
787  $ dasum( mnmin, copys, 1 ) /
788  $ ( eps*dble( mnmin ) )
789  ELSE
790  result( 11 ) = zero
791  END IF
792 *
793 * Test 12: Compute error in solution
794 *
795  CALL dlacpy( 'Full', m, nrhs, copyb, ldb, work,
796  $ ldwork )
797  CALL dqrt16( 'No transpose', m, n, nrhs, copya,
798  $ lda, b, ldb, work, ldwork,
799  $ work( m*nrhs+1 ), result( 12 ) )
800 *
801 * Test 13: Check norm of r'*A
802 *
803  result( 13 ) = zero
804  IF( m.GT.crank )
805  $ result( 13 ) = dqrt17( 'No transpose', 1, m,
806  $ n, nrhs, copya, lda, b, ldb,
807  $ copyb, ldb, c, work, lwork )
808 *
809 * Test 14: Check if x is in the rowspace of A
810 *
811  result( 14 ) = zero
812  IF( n.GT.crank )
813  $ result( 14 ) = dqrt14( 'No transpose', m, n,
814  $ nrhs, copya, lda, b, ldb,
815  $ work, lwork )
816 *
817 * Print information about the tests that did not
818 * pass the threshold.
819 *
820  DO 90 k = 3, 14
821  IF( result( k ).GE.thresh ) THEN
822  IF( nfail.EQ.0 .AND. nerrs.EQ.0 )
823  $ CALL alahd( nout, path )
824  WRITE( nout, fmt = 9998 )m, n, nrhs, nb,
825  $ itype, k, result( k )
826  nfail = nfail + 1
827  END IF
828  90 CONTINUE
829  nrun = nrun + 12
830 *
831  100 CONTINUE
832  110 CONTINUE
833  120 CONTINUE
834  130 CONTINUE
835  140 CONTINUE
836  150 CONTINUE
837 *
838 * Print a summary of the results.
839 *
840  CALL alasvm( path, nout, nfail, nrun, nerrs )
841 *
842  9999 FORMAT( ' TRANS=''', a1, ''', M=', i5, ', N=', i5, ', NRHS=', i4,
843  $ ', NB=', i4, ', type', i2, ', test(', i2, ')=', g12.5 )
844  9998 FORMAT( ' M=', i5, ', N=', i5, ', NRHS=', i4, ', NB=', i4,
845  $ ', type', i2, ', test(', i2, ')=', g12.5 )
846  9997 FORMAT( ' TRANS=''', a1,' M=', i5, ', N=', i5, ', NRHS=', i4,
847  $ ', MB=', i4,', NB=', i4,', type', i2,
848  $ ', test(', i2, ')=', g12.5 )
849 *
850  DEALLOCATE( work )
851  DEALLOCATE( iwork )
852  RETURN
853 *
854 * End of DDRVLS
855 *
Here is the call graph for this function:
Here is the caller graph for this function:
alasvm
subroutine alasvm(TYPE, NOUT, NFAIL, NRUN, NERRS)
ALASVM
Definition: alasvm.f:75
dqrt17
double precision function dqrt17(TRANS, IRESID, M, N, NRHS, A, LDA, X, LDX, B, LDB, C, WORK, LWORK)
DQRT17
Definition: dqrt17.f:152
dgelsd
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:211
alahd
subroutine alahd(IOUNIT, PATH)
ALAHD
Definition: alahd.f:109
dqrt12
double precision function dqrt12(M, N, A, LDA, S, WORK, LWORK)
DQRT12
Definition: dqrt12.f:91
dgels
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:185
dgelss
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:174
dlarnv
subroutine dlarnv(IDIST, ISEED, N, X)
DLARNV returns a vector of random numbers from a uniform or normal distribution.
Definition: dlarnv.f:99
derrls
subroutine derrls(PATH, NUNIT)
DERRLS
Definition: derrls.f:57
dgetsls
subroutine dgetsls(TRANS, M, N, NRHS, A, LDA, B, LDB, WORK, LWORK, INFO)
DGETSLS
Definition: dgetsls.f:164
dgemm
subroutine dgemm(TRANSA, TRANSB, M, N, K, ALPHA, A, LDA, B, LDB, BETA, C, LDC)
DGEMM
Definition: dgemm.f:189
dqrt13
subroutine dqrt13(SCALE, M, N, A, LDA, NORMA, ISEED)
DQRT13
Definition: dqrt13.f:93
dqrt15
subroutine dqrt15(SCALE, RKSEL, M, N, NRHS, A, LDA, B, LDB, S, RANK, NORMA, NORMB, ISEED, WORK, LWORK)
DQRT15
Definition: dqrt15.f:150
alaerh
subroutine alaerh(PATH, SUBNAM, INFO, INFOE, OPTS, M, N, KL, KU, N5, IMAT, NFAIL, NERRS, NOUT)
ALAERH
Definition: alaerh.f:149
dqrt14
double precision function dqrt14(TRANS, M, N, NRHS, A, LDA, X, LDX, WORK, LWORK)
DQRT14
Definition: dqrt14.f:118
dasum
double precision function dasum(N, DX, INCX)
DASUM
Definition: dasum.f:73
dqrt16
subroutine dqrt16(TRANS, M, N, NRHS, A, LDA, X, LDX, B, LDB, RWORK, RESID)
DQRT16
Definition: dqrt16.f:135
dscal
subroutine dscal(N, DA, DX, INCX)
DSCAL
Definition: dscal.f:81
dlasrt
subroutine dlasrt(ID, N, D, INFO)
DLASRT sorts numbers in increasing or decreasing order.
Definition: dlasrt.f:90
xlaenv
subroutine xlaenv(ISPEC, NVALUE)
XLAENV
Definition: xlaenv.f:83
dlamch
double precision function dlamch(CMACH)
DLAMCH
Definition: dlamch.f:70
dlacpy
subroutine dlacpy(UPLO, M, N, A, LDA, B, LDB)
DLACPY copies all or part of one two-dimensional array to another.
Definition: dlacpy.f:105
daxpy
subroutine daxpy(N, DA, DX, INCX, DY, INCY)
DAXPY
Definition: daxpy.f:91
dgelsy
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:206