LAPACK  3.7.0
LAPACK: Linear Algebra PACKage
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
December 2016

Definition at line 194 of file ddrvls.f.

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

Here is the call graph for this function:

Here is the caller graph for this function: