LAPACK  3.10.1
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.

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 = 16 )
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_DGETSLS, LWORK_DGELSS,
229  $ 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,
248  $ xlaenv
249 * ..
250 * .. Intrinsic Functions ..
251  INTRINSIC dble, int, log, 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, ?GETSLS, ?GELSY, ?GELSS and ?GELSD routines.
334 *
335  DO im = 1, nm
336  m = mval( im )
337  lda = max( 1, m )
338  DO in = 1, nn
339  n = nval( in )
340  mnmin = max(min( m, n ),1)
341  ldb = max( 1, m, n )
342  DO ins = 1, nns
343  nrhs = nsval( ins )
344  DO irank = 1, 2
345  DO iscale = 1, 3
346  itype = ( irank-1 )*3 + iscale
347  IF( dotype( itype ) ) THEN
348  IF( irank.EQ.1 ) THEN
349  DO itran = 1, 2
350  IF( itran.EQ.1 ) THEN
351  trans = 'N'
352  ELSE
353  trans = 'T'
354  END IF
355 *
356 * Compute workspace needed for DGELS
357  CALL dgels( trans, m, n, nrhs, a, lda,
358  $ b, ldb, wq, -1, info )
359  lwork_dgels = int( wq( 1 ) )
360 * Compute workspace needed for DGETSLS
361  CALL dgetsls( trans, m, n, nrhs, a, lda,
362  $ b, ldb, wq, -1, info )
363  lwork_dgetsls = int( wq( 1 ) )
364  ENDDO
365  END IF
366 * Compute workspace needed for DGELSY
367  CALL dgelsy( m, n, nrhs, a, lda, b, ldb, iwq,
368  $ rcond, crank, wq, -1, info )
369  lwork_dgelsy = int( wq( 1 ) )
370 * Compute workspace needed for DGELSS
371  CALL dgelss( m, n, nrhs, a, lda, b, ldb, s,
372  $ rcond, crank, wq, -1 , info )
373  lwork_dgelss = int( wq( 1 ) )
374 * Compute workspace needed for DGELSD
375  CALL dgelsd( m, n, nrhs, a, lda, b, ldb, s,
376  $ rcond, crank, wq, -1, iwq, info )
377  lwork_dgelsd = int( wq( 1 ) )
378 * Compute LIWORK workspace needed for DGELSY and DGELSD
379  liwork = max( liwork, n, iwq( 1 ) )
380 * Compute LWORK workspace needed for all functions
381  lwork = max( lwork, lwork_dgels, lwork_dgetsls,
382  $ lwork_dgelsy, lwork_dgelss,
383  $ lwork_dgelsd )
384  END IF
385  ENDDO
386  ENDDO
387  ENDDO
388  ENDDO
389  ENDDO
390 *
391  lwlsy = lwork
392 *
393  ALLOCATE( work( lwork ) )
394  ALLOCATE( iwork( liwork ) )
395 *
396  DO 150 im = 1, nm
397  m = mval( im )
398  lda = max( 1, m )
399 *
400  DO 140 in = 1, nn
401  n = nval( in )
402  mnmin = max(min( m, n ),1)
403  ldb = max( 1, m, n )
404  mb = (mnmin+1)
405 *
406  DO 130 ins = 1, nns
407  nrhs = nsval( ins )
408 *
409  DO 120 irank = 1, 2
410  DO 110 iscale = 1, 3
411  itype = ( irank-1 )*3 + iscale
412  IF( .NOT.dotype( itype ) )
413  $ GO TO 110
414 *
415  IF( irank.EQ.1 ) THEN
416 *
417 * Test DGELS
418 *
419 * Generate a matrix of scaling type ISCALE
420 *
421  CALL dqrt13( iscale, m, n, copya, lda, norma,
422  $ iseed )
423  DO 40 inb = 1, nnb
424  nb = nbval( inb )
425  CALL xlaenv( 1, nb )
426  CALL xlaenv( 3, nxval( inb ) )
427 *
428  DO 30 itran = 1, 2
429  IF( itran.EQ.1 ) THEN
430  trans = 'N'
431  nrows = m
432  ncols = n
433  ELSE
434  trans = 'T'
435  nrows = n
436  ncols = m
437  END IF
438  ldwork = max( 1, ncols )
439 *
440 * Set up a consistent rhs
441 *
442  IF( ncols.GT.0 ) THEN
443  CALL dlarnv( 2, iseed, ncols*nrhs,
444  $ work )
445  CALL dscal( ncols*nrhs,
446  $ one / dble( ncols ), work,
447  $ 1 )
448  END IF
449  CALL dgemm( trans, 'No transpose', nrows,
450  $ nrhs, ncols, one, copya, lda,
451  $ work, ldwork, zero, b, ldb )
452  CALL dlacpy( 'Full', nrows, nrhs, b, ldb,
453  $ copyb, ldb )
454 *
455 * Solve LS or overdetermined system
456 *
457  IF( m.GT.0 .AND. n.GT.0 ) THEN
458  CALL dlacpy( 'Full', m, n, copya, lda,
459  $ a, lda )
460  CALL dlacpy( 'Full', nrows, nrhs,
461  $ copyb, ldb, b, ldb )
462  END IF
463  srnamt = 'DGELS '
464  CALL dgels( trans, m, n, nrhs, a, lda, b,
465  $ ldb, work, lwork, info )
466  IF( info.NE.0 )
467  $ CALL alaerh( path, 'DGELS ', info, 0,
468  $ trans, m, n, nrhs, -1, nb,
469  $ itype, nfail, nerrs,
470  $ nout )
471 *
472 * Check correctness of results
473 *
474  ldwork = max( 1, nrows )
475  IF( nrows.GT.0 .AND. nrhs.GT.0 )
476  $ CALL dlacpy( 'Full', nrows, nrhs,
477  $ copyb, ldb, c, ldb )
478  CALL dqrt16( trans, m, n, nrhs, copya,
479  $ lda, b, ldb, c, ldb, work,
480  $ result( 1 ) )
481 *
482  IF( ( itran.EQ.1 .AND. m.GE.n ) .OR.
483  $ ( itran.EQ.2 .AND. m.LT.n ) ) THEN
484 *
485 * Solving LS system
486 *
487  result( 2 ) = dqrt17( trans, 1, m, n,
488  $ nrhs, copya, lda, b, ldb,
489  $ copyb, ldb, c, work,
490  $ lwork )
491  ELSE
492 *
493 * Solving overdetermined system
494 *
495  result( 2 ) = dqrt14( trans, m, n,
496  $ nrhs, copya, lda, b, ldb,
497  $ work, lwork )
498  END IF
499 *
500 * Print information about the tests that
501 * did not pass the threshold.
502 *
503  DO 20 k = 1, 2
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 = 9999 )trans, m,
508  $ n, nrhs, nb, itype, k,
509  $ result( k )
510  nfail = nfail + 1
511  END IF
512  20 CONTINUE
513  nrun = nrun + 2
514  30 CONTINUE
515  40 CONTINUE
516 *
517 *
518 * Test DGETSLS
519 *
520 * Generate a matrix of scaling type ISCALE
521 *
522  CALL dqrt13( iscale, m, n, copya, lda, norma,
523  $ iseed )
524  DO 65 inb = 1, nnb
525  mb = nbval( inb )
526  CALL xlaenv( 1, mb )
527  DO 62 imb = 1, nnb
528  nb = nbval( imb )
529  CALL xlaenv( 2, nb )
530 *
531  DO 60 itran = 1, 2
532  IF( itran.EQ.1 ) THEN
533  trans = 'N'
534  nrows = m
535  ncols = n
536  ELSE
537  trans = 'T'
538  nrows = n
539  ncols = m
540  END IF
541  ldwork = max( 1, ncols )
542 *
543 * Set up a consistent rhs
544 *
545  IF( ncols.GT.0 ) THEN
546  CALL dlarnv( 2, iseed, ncols*nrhs,
547  $ work )
548  CALL dscal( ncols*nrhs,
549  $ one / dble( ncols ), work,
550  $ 1 )
551  END IF
552  CALL dgemm( trans, 'No transpose', nrows,
553  $ nrhs, ncols, one, copya, lda,
554  $ work, ldwork, zero, b, ldb )
555  CALL dlacpy( 'Full', nrows, nrhs, b, ldb,
556  $ copyb, ldb )
557 *
558 * Solve LS or overdetermined system
559 *
560  IF( m.GT.0 .AND. n.GT.0 ) THEN
561  CALL dlacpy( 'Full', m, n, copya, lda,
562  $ a, lda )
563  CALL dlacpy( 'Full', nrows, nrhs,
564  $ copyb, ldb, b, ldb )
565  END IF
566  srnamt = 'DGETSLS '
567  CALL dgetsls( trans, m, n, nrhs, a,
568  $ lda, b, ldb, work, lwork, info )
569  IF( info.NE.0 )
570  $ CALL alaerh( path, 'DGETSLS ', info, 0,
571  $ trans, m, n, nrhs, -1, nb,
572  $ itype, nfail, nerrs,
573  $ nout )
574 *
575 * Check correctness of results
576 *
577  ldwork = max( 1, nrows )
578  IF( nrows.GT.0 .AND. nrhs.GT.0 )
579  $ CALL dlacpy( 'Full', nrows, nrhs,
580  $ copyb, ldb, c, ldb )
581  CALL dqrt16( trans, m, n, nrhs, copya,
582  $ lda, b, ldb, c, ldb, work,
583  $ result( 15 ) )
584 *
585  IF( ( itran.EQ.1 .AND. m.GE.n ) .OR.
586  $ ( itran.EQ.2 .AND. m.LT.n ) ) THEN
587 *
588 * Solving LS system
589 *
590  result( 16 ) = dqrt17( trans, 1, m, n,
591  $ nrhs, copya, lda, b, ldb,
592  $ copyb, ldb, c, work,
593  $ lwork )
594  ELSE
595 *
596 * Solving overdetermined system
597 *
598  result( 16 ) = dqrt14( trans, m, n,
599  $ nrhs, copya, lda, b, ldb,
600  $ work, lwork )
601  END IF
602 *
603 * Print information about the tests that
604 * did not pass the threshold.
605 *
606  DO 50 k = 15, 16
607  IF( result( k ).GE.thresh ) THEN
608  IF( nfail.EQ.0 .AND. nerrs.EQ.0 )
609  $ CALL alahd( nout, path )
610  WRITE( nout, fmt = 9997 )trans, m,
611  $ n, nrhs, mb, nb, itype, k,
612  $ result( k )
613  nfail = nfail + 1
614  END IF
615  50 CONTINUE
616  nrun = nrun + 2
617  60 CONTINUE
618  62 CONTINUE
619  65 CONTINUE
620  END IF
621 *
622 * Generate a matrix of scaling type ISCALE and rank
623 * type IRANK.
624 *
625  CALL dqrt15( iscale, irank, m, n, nrhs, copya, lda,
626  $ copyb, ldb, copys, rank, norma, normb,
627  $ iseed, work, lwork )
628 *
629 * workspace used: MAX(M+MIN(M,N),NRHS*MIN(M,N),2*N+M)
630 *
631  ldwork = max( 1, m )
632 *
633 * Loop for testing different block sizes.
634 *
635  DO 100 inb = 1, nnb
636  nb = nbval( inb )
637  CALL xlaenv( 1, nb )
638  CALL xlaenv( 3, nxval( inb ) )
639 *
640 * Test DGELSY
641 *
642 * DGELSY: Compute the minimum-norm solution X
643 * to min( norm( A * X - B ) )
644 * using the rank-revealing orthogonal
645 * factorization.
646 *
647 * Initialize vector IWORK.
648 *
649  DO 70 j = 1, n
650  iwork( j ) = 0
651  70 CONTINUE
652 *
653  CALL dlacpy( 'Full', m, n, copya, lda, a, lda )
654  CALL dlacpy( 'Full', m, nrhs, copyb, ldb, b,
655  $ ldb )
656 *
657  srnamt = 'DGELSY'
658  CALL dgelsy( m, n, nrhs, a, lda, b, ldb, iwork,
659  $ rcond, crank, work, lwlsy, info )
660  IF( info.NE.0 )
661  $ CALL alaerh( path, 'DGELSY', info, 0, ' ', m,
662  $ n, nrhs, -1, nb, itype, nfail,
663  $ nerrs, nout )
664 *
665 * Test 3: Compute relative error in svd
666 * workspace: M*N + 4*MIN(M,N) + MAX(M,N)
667 *
668  result( 3 ) = dqrt12( crank, crank, a, lda,
669  $ copys, work, lwork )
670 *
671 * Test 4: Compute error in solution
672 * workspace: M*NRHS + M
673 *
674  CALL dlacpy( 'Full', m, nrhs, copyb, ldb, work,
675  $ ldwork )
676  CALL dqrt16( 'No transpose', m, n, nrhs, copya,
677  $ lda, b, ldb, work, ldwork,
678  $ work( m*nrhs+1 ), result( 4 ) )
679 *
680 * Test 5: Check norm of r'*A
681 * workspace: NRHS*(M+N)
682 *
683  result( 5 ) = zero
684  IF( m.GT.crank )
685  $ result( 5 ) = dqrt17( 'No transpose', 1, m,
686  $ n, nrhs, copya, lda, b, ldb,
687  $ copyb, ldb, c, work, lwork )
688 *
689 * Test 6: Check if x is in the rowspace of A
690 * workspace: (M+NRHS)*(N+2)
691 *
692  result( 6 ) = zero
693 *
694  IF( n.GT.crank )
695  $ result( 6 ) = dqrt14( 'No transpose', m, n,
696  $ nrhs, copya, lda, b, ldb,
697  $ work, lwork )
698 *
699 * Test DGELSS
700 *
701 * DGELSS: Compute the minimum-norm solution X
702 * to min( norm( A * X - B ) )
703 * using the SVD.
704 *
705  CALL dlacpy( 'Full', m, n, copya, lda, a, lda )
706  CALL dlacpy( 'Full', m, nrhs, copyb, ldb, b,
707  $ ldb )
708  srnamt = 'DGELSS'
709  CALL dgelss( m, n, nrhs, a, lda, b, ldb, s,
710  $ rcond, crank, work, lwork, info )
711  IF( info.NE.0 )
712  $ CALL alaerh( path, 'DGELSS', info, 0, ' ', m,
713  $ n, nrhs, -1, nb, itype, nfail,
714  $ nerrs, nout )
715 *
716 * workspace used: 3*min(m,n) +
717 * max(2*min(m,n),nrhs,max(m,n))
718 *
719 * Test 7: Compute relative error in svd
720 *
721  IF( rank.GT.0 ) THEN
722  CALL daxpy( mnmin, -one, copys, 1, s, 1 )
723  result( 7 ) = dasum( mnmin, s, 1 ) /
724  $ dasum( mnmin, copys, 1 ) /
725  $ ( eps*dble( mnmin ) )
726  ELSE
727  result( 7 ) = zero
728  END IF
729 *
730 * Test 8: Compute error in solution
731 *
732  CALL dlacpy( 'Full', m, nrhs, copyb, ldb, work,
733  $ ldwork )
734  CALL dqrt16( 'No transpose', m, n, nrhs, copya,
735  $ lda, b, ldb, work, ldwork,
736  $ work( m*nrhs+1 ), result( 8 ) )
737 *
738 * Test 9: Check norm of r'*A
739 *
740  result( 9 ) = zero
741  IF( m.GT.crank )
742  $ result( 9 ) = dqrt17( 'No transpose', 1, m,
743  $ n, nrhs, copya, lda, b, ldb,
744  $ copyb, ldb, c, work, lwork )
745 *
746 * Test 10: Check if x is in the rowspace of A
747 *
748  result( 10 ) = zero
749  IF( n.GT.crank )
750  $ result( 10 ) = dqrt14( 'No transpose', m, n,
751  $ nrhs, copya, lda, b, ldb,
752  $ work, lwork )
753 *
754 * Test DGELSD
755 *
756 * DGELSD: Compute the minimum-norm solution X
757 * to min( norm( A * X - B ) ) using a
758 * divide and conquer SVD.
759 *
760 * Initialize vector IWORK.
761 *
762  DO 80 j = 1, n
763  iwork( j ) = 0
764  80 CONTINUE
765 *
766  CALL dlacpy( 'Full', m, n, copya, lda, a, lda )
767  CALL dlacpy( 'Full', m, nrhs, copyb, ldb, b,
768  $ ldb )
769 *
770  srnamt = 'DGELSD'
771  CALL dgelsd( m, n, nrhs, a, lda, b, ldb, s,
772  $ rcond, crank, work, lwork, iwork,
773  $ info )
774  IF( info.NE.0 )
775  $ CALL alaerh( path, 'DGELSD', info, 0, ' ', m,
776  $ n, nrhs, -1, nb, itype, nfail,
777  $ nerrs, nout )
778 *
779 * Test 11: Compute relative error in svd
780 *
781  IF( rank.GT.0 ) THEN
782  CALL daxpy( mnmin, -one, copys, 1, s, 1 )
783  result( 11 ) = dasum( mnmin, s, 1 ) /
784  $ dasum( mnmin, copys, 1 ) /
785  $ ( eps*dble( mnmin ) )
786  ELSE
787  result( 11 ) = zero
788  END IF
789 *
790 * Test 12: Compute error in solution
791 *
792  CALL dlacpy( 'Full', m, nrhs, copyb, ldb, work,
793  $ ldwork )
794  CALL dqrt16( 'No transpose', m, n, nrhs, copya,
795  $ lda, b, ldb, work, ldwork,
796  $ work( m*nrhs+1 ), result( 12 ) )
797 *
798 * Test 13: Check norm of r'*A
799 *
800  result( 13 ) = zero
801  IF( m.GT.crank )
802  $ result( 13 ) = dqrt17( 'No transpose', 1, m,
803  $ n, nrhs, copya, lda, b, ldb,
804  $ copyb, ldb, c, work, lwork )
805 *
806 * Test 14: Check if x is in the rowspace of A
807 *
808  result( 14 ) = zero
809  IF( n.GT.crank )
810  $ result( 14 ) = dqrt14( 'No transpose', m, n,
811  $ nrhs, copya, lda, b, ldb,
812  $ work, lwork )
813 *
814 * Print information about the tests that did not
815 * pass the threshold.
816 *
817  DO 90 k = 3, 14
818  IF( result( k ).GE.thresh ) THEN
819  IF( nfail.EQ.0 .AND. nerrs.EQ.0 )
820  $ CALL alahd( nout, path )
821  WRITE( nout, fmt = 9998 )m, n, nrhs, nb,
822  $ itype, k, result( k )
823  nfail = nfail + 1
824  END IF
825  90 CONTINUE
826  nrun = nrun + 12
827 *
828  100 CONTINUE
829  110 CONTINUE
830  120 CONTINUE
831  130 CONTINUE
832  140 CONTINUE
833  150 CONTINUE
834 *
835 * Print a summary of the results.
836 *
837  CALL alasvm( path, nout, nfail, nrun, nerrs )
838 *
839  9999 FORMAT( ' TRANS=''', a1, ''', M=', i5, ', N=', i5, ', NRHS=', i4,
840  $ ', NB=', i4, ', type', i2, ', test(', i2, ')=', g12.5 )
841  9998 FORMAT( ' M=', i5, ', N=', i5, ', NRHS=', i4, ', NB=', i4,
842  $ ', type', i2, ', test(', i2, ')=', g12.5 )
843  9997 FORMAT( ' TRANS=''', a1,' M=', i5, ', N=', i5, ', NRHS=', i4,
844  $ ', MB=', i4,', NB=', i4,', type', i2,
845  $ ', test(', i2, ')=', g12.5 )
846 *
847  DEALLOCATE( work )
848  DEALLOCATE( iwork )
849  RETURN
850 *
851 * End of DDRVLS
852 *
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
subroutine dlasrt(ID, N, D, INFO)
DLASRT sorts numbers in increasing or decreasing order.
Definition: dlasrt.f:88
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 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: