LAPACK  3.6.1
LAPACK: Linear Algebra PACKage
subroutine ddrvbd ( integer  NSIZES,
integer, dimension( * )  MM,
integer, dimension( * )  NN,
integer  NTYPES,
logical, dimension( * )  DOTYPE,
integer, dimension( 4 )  ISEED,
double precision  THRESH,
double precision, dimension( lda, * )  A,
integer  LDA,
double precision, dimension( ldu, * )  U,
integer  LDU,
double precision, dimension( ldvt, * )  VT,
integer  LDVT,
double precision, dimension( lda, * )  ASAV,
double precision, dimension( ldu, * )  USAV,
double precision, dimension( ldvt, * )  VTSAV,
double precision, dimension( * )  S,
double precision, dimension( * )  SSAV,
double precision, dimension( * )  E,
double precision, dimension( * )  WORK,
integer  LWORK,
integer, dimension( * )  IWORK,
integer  NOUT,
integer  INFO 
)

DDRVBD

Purpose: 
 DDRVBD checks the singular value decomposition (SVD) drivers
 DGESVD, DGESDD, DGESVJ, and DGEJSV.

 Both DGESVD and DGESDD factor A = U diag(S) VT, where U and VT are
 orthogonal and diag(S) is diagonal with the entries of the array S
 on its diagonal. The entries of S are the singular values,
 nonnegative and stored in decreasing order.  U and VT can be
 optionally not computed, overwritten on A, or computed partially.

 A is M by N. Let MNMIN = min( M, N ). S has dimension MNMIN.
 U can be M by M or M by MNMIN. VT can be N by N or MNMIN by N.

 When DDRVBD is called, a number of matrix "sizes" (M's and N's)
 and a number of matrix "types" are specified.  For each size (M,N)
 and each type of matrix, and for the minimal workspace as well as
 workspace adequate to permit blocking, an  M x N  matrix "A" will be
 generated and used to test the SVD routines.  For each matrix, A will
 be factored as A = U diag(S) VT and the following 12 tests computed:

 Test for DGESVD:

 (1)    | A - U diag(S) VT | / ( |A| max(M,N) ulp )

 (2)    | I - U'U | / ( M ulp )

 (3)    | I - VT VT' | / ( N ulp )

 (4)    S contains MNMIN nonnegative values in decreasing order.
        (Return 0 if true, 1/ULP if false.)

 (5)    | U - Upartial | / ( M ulp ) where Upartial is a partially
        computed U.

 (6)    | VT - VTpartial | / ( N ulp ) where VTpartial is a partially
        computed VT.

 (7)    | S - Spartial | / ( MNMIN ulp |S| ) where Spartial is the
        vector of singular values from the partial SVD

 Test for DGESDD:

 (8)    | A - U diag(S) VT | / ( |A| max(M,N) ulp )

 (9)    | I - U'U | / ( M ulp )

 (10)   | I - VT VT' | / ( N ulp )

 (11)   S contains MNMIN nonnegative values in decreasing order.
        (Return 0 if true, 1/ULP if false.)

 (12)   | U - Upartial | / ( M ulp ) where Upartial is a partially
        computed U.

 (13)   | VT - VTpartial | / ( N ulp ) where VTpartial is a partially
        computed VT.

 (14)   | S - Spartial | / ( MNMIN ulp |S| ) where Spartial is the
        vector of singular values from the partial SVD

 Test for DGESVJ:

 (15)   | A - U diag(S) VT | / ( |A| max(M,N) ulp )

 (16)   | I - U'U | / ( M ulp )

 (17)   | I - VT VT' | / ( N ulp )

 (18)   S contains MNMIN nonnegative values in decreasing order.
        (Return 0 if true, 1/ULP if false.)

 Test for DGEJSV:

 (19)   | A - U diag(S) VT | / ( |A| max(M,N) ulp )

 (20)   | I - U'U | / ( M ulp )

 (21)   | I - VT VT' | / ( N ulp )

 (22)   S contains MNMIN nonnegative values in decreasing order.
        (Return 0 if true, 1/ULP if false.)

 Test for DGESVDX( 'V', 'V', 'A' )/DGESVDX( 'N', 'N', 'A' )

 (23)   | A - U diag(S) VT | / ( |A| max(M,N) ulp )

 (24)   | I - U'U | / ( M ulp )

 (25)   | I - VT VT' | / ( N ulp )

 (26)   S contains MNMIN nonnegative values in decreasing order.
        (Return 0 if true, 1/ULP if false.)

 (27)   | U - Upartial | / ( M ulp ) where Upartial is a partially
        computed U.

 (28)   | VT - VTpartial | / ( N ulp ) where VTpartial is a partially
        computed VT.

 (29)   | S - Spartial | / ( MNMIN ulp |S| ) where Spartial is the
        vector of singular values from the partial SVD
      
 Test for DGESVDX( 'V', 'V', 'I' )

 (30)   | U' A VT''' - diag(S) | / ( |A| max(M,N) ulp )

 (31)   | I - U'U | / ( M ulp )

 (32)   | I - VT VT' | / ( N ulp )
      
 Test for DGESVDX( 'V', 'V', 'V' )

 (33)   | U' A VT''' - diag(S) | / ( |A| max(M,N) ulp )

 (34)   | I - U'U | / ( M ulp )

 (35)   | I - VT VT' | / ( N ulp )

 The "sizes" are specified by the arrays MM(1:NSIZES) and
 NN(1:NSIZES); the value of each element pair (MM(j),NN(j))
 specifies one size.  The "types" are specified by a logical array
 DOTYPE( 1:NTYPES ); if DOTYPE(j) is .TRUE., then matrix type "j"
 will be generated.
 Currently, the list of possible types is:

 (1)  The zero matrix.
 (2)  The identity matrix.
 (3)  A matrix of the form  U D V, where U and V are orthogonal and
      D has evenly spaced entries 1, ..., ULP with random signs
      on the diagonal.
 (4)  Same as (3), but multiplied by the underflow-threshold / ULP.
 (5)  Same as (3), but multiplied by the overflow-threshold * ULP.
Parameters
[in]NSIZES
          NSIZES is INTEGER
          The number of matrix sizes (M,N) contained in the vectors
          MM and NN.
[in]MM
          MM is INTEGER array, dimension (NSIZES)
          The values of the matrix row dimension M.
[in]NN
          NN is INTEGER array, dimension (NSIZES)
          The values of the matrix column dimension N.
[in]NTYPES
          NTYPES is INTEGER
          The number of elements in DOTYPE.   If it is zero, DDRVBD
          does nothing.  It must be at least zero.  If it is MAXTYP+1
          and NSIZES is 1, then an additional type, MAXTYP+1 is
          defined, which is to use whatever matrices are in A and B.
          This is only useful if DOTYPE(1:MAXTYP) is .FALSE. and
          DOTYPE(MAXTYP+1) is .TRUE. .
[in]DOTYPE
          DOTYPE is LOGICAL array, dimension (NTYPES)
          If DOTYPE(j) is .TRUE., then for each size (m,n), a matrix
          of type j will be generated.  If NTYPES is smaller than the
          maximum number of types defined (PARAMETER MAXTYP), then
          types NTYPES+1 through MAXTYP will not be generated.  If
          NTYPES is larger than MAXTYP, DOTYPE(MAXTYP+1) through
          DOTYPE(NTYPES) will be ignored.
[in,out]ISEED
          ISEED is INTEGER array, dimension (4)
          On entry, the seed of the random number generator.  The array
          elements should be between 0 and 4095; if not they will be
          reduced mod 4096.  Also, ISEED(4) must be odd.
          On exit, ISEED is changed and can be used in the next call to
          DDRVBD to continue the same random number sequence.
[in]THRESH
          THRESH is DOUBLE PRECISION
          The threshold value for the test ratios.  A result is
          included in the output file if RESULT >= THRESH.  The test
          ratios are scaled to be O(1), so THRESH should be a small
          multiple of 1, e.g., 10 or 100.  To have every test ratio
          printed, use THRESH = 0.
[out]A
          A is DOUBLE PRECISION array, dimension (LDA,NMAX)
          where NMAX is the maximum value of N in NN.
[in]LDA
          LDA is INTEGER
          The leading dimension of the array A.  LDA >= max(1,MMAX),
          where MMAX is the maximum value of M in MM.
[out]U
          U is DOUBLE PRECISION array, dimension (LDU,MMAX)
[in]LDU
          LDU is INTEGER
          The leading dimension of the array U.  LDU >= max(1,MMAX).
[out]VT
          VT is DOUBLE PRECISION array, dimension (LDVT,NMAX)
[in]LDVT
          LDVT is INTEGER
          The leading dimension of the array VT.  LDVT >= max(1,NMAX).
[out]ASAV
          ASAV is DOUBLE PRECISION array, dimension (LDA,NMAX)
[out]USAV
          USAV is DOUBLE PRECISION array, dimension (LDU,MMAX)
[out]VTSAV
          VTSAV is DOUBLE PRECISION array, dimension (LDVT,NMAX)
[out]S
          S is DOUBLE PRECISION array, dimension
                      (max(min(MM,NN)))
[out]SSAV
          SSAV is DOUBLE PRECISION array, dimension
                      (max(min(MM,NN)))
[out]E
          E is DOUBLE PRECISION array, dimension
                      (max(min(MM,NN)))
[out]WORK
          WORK is DOUBLE PRECISION array, dimension (LWORK)
[in]LWORK
          LWORK is INTEGER
          The number of entries in WORK.  This must be at least
          max(3*MN+MX,5*MN-4)+2*MN**2 for all pairs
          pairs  (MN,MX)=( min(MM(j),NN(j), max(MM(j),NN(j)) )
[out]IWORK
          IWORK is INTEGER array, dimension at least 8*min(M,N)
[in]NOUT
          NOUT is INTEGER
          The FORTRAN unit number for printing out error messages
          (e.g., if a routine returns IINFO not equal to 0.)
[out]INFO
          INFO is INTEGER
          If 0, then everything ran OK.
           -1: NSIZES < 0
           -2: Some MM(j) < 0
           -3: Some NN(j) < 0
           -4: NTYPES < 0
           -7: THRESH < 0
          -10: LDA < 1 or LDA < MMAX, where MMAX is max( MM(j) ).
          -12: LDU < 1 or LDU < MMAX.
          -14: LDVT < 1 or LDVT < NMAX, where NMAX is max( NN(j) ).
          -21: LWORK too small.
          If  DLATMS, or DGESVD returns an error code, the
              absolute value of it is returned.
Author
Univ. of Tennessee
Univ. of California Berkeley
Univ. of Colorado Denver
NAG Ltd.
Date
June 2016

Definition at line 357 of file ddrvbd.f.

357 *
358 * -- LAPACK test routine (version 3.6.1) --
359 * -- LAPACK is a software package provided by Univ. of Tennessee, --
360 * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
361 * June 2016
362 *
363 * .. Scalar Arguments ..
364  INTEGER info, lda, ldu, ldvt, lwork, nout, nsizes,
365  $ ntypes
366  DOUBLE PRECISION thresh
367 * ..
368 * .. Array Arguments ..
369  LOGICAL dotype( * )
370  INTEGER iseed( 4 ), iwork( * ), mm( * ), nn( * )
371  DOUBLE PRECISION a( lda, * ), asav( lda, * ), e( * ), s( * ),
372  $ ssav( * ), u( ldu, * ), usav( ldu, * ),
373  $ vt( ldvt, * ), vtsav( ldvt, * ), work( * )
374 * ..
375 *
376 * =====================================================================
377 *
378 * .. Parameters ..
379  DOUBLE PRECISION zero, one, two, half
380  parameter ( zero = 0.0d0, one = 1.0d0, two = 2.0d0,
381  $ half = 0.5d0 )
382  INTEGER maxtyp
383  parameter ( maxtyp = 5 )
384 * ..
385 * .. Local Scalars ..
386  LOGICAL badmm, badnn
387  CHARACTER jobq, jobu, jobvt, range
388  CHARACTER*3 path
389  INTEGER i, iinfo, ijq, iju, ijvt, il,iu, iws, iwtmp,
390  $ itemp, j, jsize, jtype, lswork, m, minwrk,
391  $ mmax, mnmax, mnmin, mtypes, n, nfail,
392  $ nmax, ns, nsi, nsv, ntest
393  DOUBLE PRECISION anorm, dif, div, ovfl, rtunfl, ulp,
394  $ ulpinv, unfl, vl, vu
395 * ..
396 * .. Local Arrays ..
397  CHARACTER cjob( 4 ), cjobr( 3 ), cjobv( 2 )
398  INTEGER ioldsd( 4 ), iseed2( 4 )
399  DOUBLE PRECISION result( 40 )
400 * ..
401 * .. External Functions ..
402  DOUBLE PRECISION dlamch, dlarnd
403  EXTERNAL dlamch, dlarnd
404 * ..
405 * .. External Subroutines ..
406  EXTERNAL alasvm, dbdt01, dgejsv, dgesdd, dgesvd,
407  $ dgesvdx, dgesvj, dlabad, dlacpy, dlaset,
408  $ dlatms, dort01, dort03, xerbla
409 * ..
410 * .. Intrinsic Functions ..
411  INTRINSIC abs, dble, int, max, min
412 * ..
413 * .. Scalars in Common ..
414  LOGICAL lerr, ok
415  CHARACTER*32 srnamt
416  INTEGER infot, nunit
417 * ..
418 * .. Common blocks ..
419  COMMON / infoc / infot, nunit, ok, lerr
420  COMMON / srnamc / srnamt
421 * ..
422 * .. Data statements ..
423  DATA cjob / 'N', 'O', 'S', 'A' /
424  DATA cjobr / 'A', 'V', 'I' /
425  DATA cjobv / 'N', 'V' /
426 * ..
427 * .. Executable Statements ..
428 *
429 * Check for errors
430 *
431  info = 0
432  badmm = .false.
433  badnn = .false.
434  mmax = 1
435  nmax = 1
436  mnmax = 1
437  minwrk = 1
438  DO 10 j = 1, nsizes
439  mmax = max( mmax, mm( j ) )
440  IF( mm( j ).LT.0 )
441  $ badmm = .true.
442  nmax = max( nmax, nn( j ) )
443  IF( nn( j ).LT.0 )
444  $ badnn = .true.
445  mnmax = max( mnmax, min( mm( j ), nn( j ) ) )
446  minwrk = max( minwrk, max( 3*min( mm( j ),
447  $ nn( j ) )+max( mm( j ), nn( j ) ), 5*min( mm( j ),
448  $ nn( j )-4 ) )+2*min( mm( j ), nn( j ) )**2 )
449  10 CONTINUE
450 *
451 * Check for errors
452 *
453  IF( nsizes.LT.0 ) THEN
454  info = -1
455  ELSE IF( badmm ) THEN
456  info = -2
457  ELSE IF( badnn ) THEN
458  info = -3
459  ELSE IF( ntypes.LT.0 ) THEN
460  info = -4
461  ELSE IF( lda.LT.max( 1, mmax ) ) THEN
462  info = -10
463  ELSE IF( ldu.LT.max( 1, mmax ) ) THEN
464  info = -12
465  ELSE IF( ldvt.LT.max( 1, nmax ) ) THEN
466  info = -14
467  ELSE IF( minwrk.GT.lwork ) THEN
468  info = -21
469  END IF
470 *
471  IF( info.NE.0 ) THEN
472  CALL xerbla( 'DDRVBD', -info )
473  RETURN
474  END IF
475 *
476 * Initialize constants
477 *
478  path( 1: 1 ) = 'Double precision'
479  path( 2: 3 ) = 'BD'
480  nfail = 0
481  ntest = 0
482  unfl = dlamch( 'Safe minimum' )
483  ovfl = one / unfl
484  CALL dlabad( unfl, ovfl )
485  ulp = dlamch( 'Precision' )
486  rtunfl = sqrt( unfl )
487  ulpinv = one / ulp
488  infot = 0
489 *
490 * Loop over sizes, types
491 *
492  DO 240 jsize = 1, nsizes
493  m = mm( jsize )
494  n = nn( jsize )
495  mnmin = min( m, n )
496 *
497  IF( nsizes.NE.1 ) THEN
498  mtypes = min( maxtyp, ntypes )
499  ELSE
500  mtypes = min( maxtyp+1, ntypes )
501  END IF
502 *
503  DO 230 jtype = 1, mtypes
504  IF( .NOT.dotype( jtype ) )
505  $ GO TO 230
506 *
507  DO 20 j = 1, 4
508  ioldsd( j ) = iseed( j )
509  20 CONTINUE
510 *
511 * Compute "A"
512 *
513  IF( mtypes.GT.maxtyp )
514  $ GO TO 30
515 *
516  IF( jtype.EQ.1 ) THEN
517 *
518 * Zero matrix
519 *
520  CALL dlaset( 'Full', m, n, zero, zero, a, lda )
521 *
522  ELSE IF( jtype.EQ.2 ) THEN
523 *
524 * Identity matrix
525 *
526  CALL dlaset( 'Full', m, n, zero, one, a, lda )
527 *
528  ELSE
529 *
530 * (Scaled) random matrix
531 *
532  IF( jtype.EQ.3 )
533  $ anorm = one
534  IF( jtype.EQ.4 )
535  $ anorm = unfl / ulp
536  IF( jtype.EQ.5 )
537  $ anorm = ovfl*ulp
538  CALL dlatms( m, n, 'U', iseed, 'N', s, 4, dble( mnmin ),
539  $ anorm, m-1, n-1, 'N', a, lda, work, iinfo )
540  IF( iinfo.NE.0 ) THEN
541  WRITE( nout, fmt = 9996 )'Generator', iinfo, m, n,
542  $ jtype, ioldsd
543  info = abs( iinfo )
544  RETURN
545  END IF
546  END IF
547 *
548  30 CONTINUE
549  CALL dlacpy( 'F', m, n, a, lda, asav, lda )
550 *
551 * Do for minimal and adequate (for blocking) workspace
552 *
553  DO 220 iws = 1, 4
554 *
555  DO 40 j = 1, 32
556  result( j ) = -one
557  40 CONTINUE
558 *
559 * Test DGESVD: Factorize A
560 *
561  iwtmp = max( 3*min( m, n )+max( m, n ), 5*min( m, n ) )
562  lswork = iwtmp + ( iws-1 )*( lwork-iwtmp ) / 3
563  lswork = min( lswork, lwork )
564  lswork = max( lswork, 1 )
565  IF( iws.EQ.4 )
566  $ lswork = lwork
567 *
568  IF( iws.GT.1 )
569  $ CALL dlacpy( 'F', m, n, asav, lda, a, lda )
570  srnamt = 'DGESVD'
571  CALL dgesvd( 'A', 'A', m, n, a, lda, ssav, usav, ldu,
572  $ vtsav, ldvt, work, lswork, iinfo )
573  IF( iinfo.NE.0 ) THEN
574  WRITE( nout, fmt = 9995 )'GESVD', iinfo, m, n, jtype,
575  $ lswork, ioldsd
576  info = abs( iinfo )
577  RETURN
578  END IF
579 *
580 * Do tests 1--4
581 *
582  CALL dbdt01( m, n, 0, asav, lda, usav, ldu, ssav, e,
583  $ vtsav, ldvt, work, result( 1 ) )
584  IF( m.NE.0 .AND. n.NE.0 ) THEN
585  CALL dort01( 'Columns', m, m, usav, ldu, work, lwork,
586  $ result( 2 ) )
587  CALL dort01( 'Rows', n, n, vtsav, ldvt, work, lwork,
588  $ result( 3 ) )
589  END IF
590  result( 4 ) = zero
591  DO 50 i = 1, mnmin - 1
592  IF( ssav( i ).LT.ssav( i+1 ) )
593  $ result( 4 ) = ulpinv
594  IF( ssav( i ).LT.zero )
595  $ result( 4 ) = ulpinv
596  50 CONTINUE
597  IF( mnmin.GE.1 ) THEN
598  IF( ssav( mnmin ).LT.zero )
599  $ result( 4 ) = ulpinv
600  END IF
601 *
602 * Do partial SVDs, comparing to SSAV, USAV, and VTSAV
603 *
604  result( 5 ) = zero
605  result( 6 ) = zero
606  result( 7 ) = zero
607  DO 80 iju = 0, 3
608  DO 70 ijvt = 0, 3
609  IF( ( iju.EQ.3 .AND. ijvt.EQ.3 ) .OR.
610  $ ( iju.EQ.1 .AND. ijvt.EQ.1 ) )GO TO 70
611  jobu = cjob( iju+1 )
612  jobvt = cjob( ijvt+1 )
613  CALL dlacpy( 'F', m, n, asav, lda, a, lda )
614  srnamt = 'DGESVD'
615  CALL dgesvd( jobu, jobvt, m, n, a, lda, s, u, ldu,
616  $ vt, ldvt, work, lswork, iinfo )
617 *
618 * Compare U
619 *
620  dif = zero
621  IF( m.GT.0 .AND. n.GT.0 ) THEN
622  IF( iju.EQ.1 ) THEN
623  CALL dort03( 'C', m, mnmin, m, mnmin, usav,
624  $ ldu, a, lda, work, lwork, dif,
625  $ iinfo )
626  ELSE IF( iju.EQ.2 ) THEN
627  CALL dort03( 'C', m, mnmin, m, mnmin, usav,
628  $ ldu, u, ldu, work, lwork, dif,
629  $ iinfo )
630  ELSE IF( iju.EQ.3 ) THEN
631  CALL dort03( 'C', m, m, m, mnmin, usav, ldu,
632  $ u, ldu, work, lwork, dif,
633  $ iinfo )
634  END IF
635  END IF
636  result( 5 ) = max( result( 5 ), dif )
637 *
638 * Compare VT
639 *
640  dif = zero
641  IF( m.GT.0 .AND. n.GT.0 ) THEN
642  IF( ijvt.EQ.1 ) THEN
643  CALL dort03( 'R', n, mnmin, n, mnmin, vtsav,
644  $ ldvt, a, lda, work, lwork, dif,
645  $ iinfo )
646  ELSE IF( ijvt.EQ.2 ) THEN
647  CALL dort03( 'R', n, mnmin, n, mnmin, vtsav,
648  $ ldvt, vt, ldvt, work, lwork,
649  $ dif, iinfo )
650  ELSE IF( ijvt.EQ.3 ) THEN
651  CALL dort03( 'R', n, n, n, mnmin, vtsav,
652  $ ldvt, vt, ldvt, work, lwork,
653  $ dif, iinfo )
654  END IF
655  END IF
656  result( 6 ) = max( result( 6 ), dif )
657 *
658 * Compare S
659 *
660  dif = zero
661  div = max( mnmin*ulp*s( 1 ), unfl )
662  DO 60 i = 1, mnmin - 1
663  IF( ssav( i ).LT.ssav( i+1 ) )
664  $ dif = ulpinv
665  IF( ssav( i ).LT.zero )
666  $ dif = ulpinv
667  dif = max( dif, abs( ssav( i )-s( i ) ) / div )
668  60 CONTINUE
669  result( 7 ) = max( result( 7 ), dif )
670  70 CONTINUE
671  80 CONTINUE
672 *
673 * Test DGESDD: Factorize A
674 *
675  iwtmp = 5*mnmin*mnmin + 9*mnmin + max( m, n )
676  lswork = iwtmp + ( iws-1 )*( lwork-iwtmp ) / 3
677  lswork = min( lswork, lwork )
678  lswork = max( lswork, 1 )
679  IF( iws.EQ.4 )
680  $ lswork = lwork
681 *
682  CALL dlacpy( 'F', m, n, asav, lda, a, lda )
683  srnamt = 'DGESDD'
684  CALL dgesdd( 'A', m, n, a, lda, ssav, usav, ldu, vtsav,
685  $ ldvt, work, lswork, iwork, iinfo )
686  IF( iinfo.NE.0 ) THEN
687  WRITE( nout, fmt = 9995 )'GESDD', iinfo, m, n, jtype,
688  $ lswork, ioldsd
689  info = abs( iinfo )
690  RETURN
691  END IF
692 *
693 * Do tests 8--11
694 *
695  CALL dbdt01( m, n, 0, asav, lda, usav, ldu, ssav, e,
696  $ vtsav, ldvt, work, result( 8 ) )
697  IF( m.NE.0 .AND. n.NE.0 ) THEN
698  CALL dort01( 'Columns', m, m, usav, ldu, work, lwork,
699  $ result( 9 ) )
700  CALL dort01( 'Rows', n, n, vtsav, ldvt, work, lwork,
701  $ result( 10 ) )
702  END IF
703  result( 11 ) = zero
704  DO 90 i = 1, mnmin - 1
705  IF( ssav( i ).LT.ssav( i+1 ) )
706  $ result( 11 ) = ulpinv
707  IF( ssav( i ).LT.zero )
708  $ result( 11 ) = ulpinv
709  90 CONTINUE
710  IF( mnmin.GE.1 ) THEN
711  IF( ssav( mnmin ).LT.zero )
712  $ result( 11 ) = ulpinv
713  END IF
714 *
715 * Do partial SVDs, comparing to SSAV, USAV, and VTSAV
716 *
717  result( 12 ) = zero
718  result( 13 ) = zero
719  result( 14 ) = zero
720  DO 110 ijq = 0, 2
721  jobq = cjob( ijq+1 )
722  CALL dlacpy( 'F', m, n, asav, lda, a, lda )
723  srnamt = 'DGESDD'
724  CALL dgesdd( jobq, m, n, a, lda, s, u, ldu, vt, ldvt,
725  $ work, lswork, iwork, iinfo )
726 *
727 * Compare U
728 *
729  dif = zero
730  IF( m.GT.0 .AND. n.GT.0 ) THEN
731  IF( ijq.EQ.1 ) THEN
732  IF( m.GE.n ) THEN
733  CALL dort03( 'C', m, mnmin, m, mnmin, usav,
734  $ ldu, a, lda, work, lwork, dif,
735  $ info )
736  ELSE
737  CALL dort03( 'C', m, mnmin, m, mnmin, usav,
738  $ ldu, u, ldu, work, lwork, dif,
739  $ info )
740  END IF
741  ELSE IF( ijq.EQ.2 ) THEN
742  CALL dort03( 'C', m, mnmin, m, mnmin, usav, ldu,
743  $ u, ldu, work, lwork, dif, info )
744  END IF
745  END IF
746  result( 12 ) = max( result( 12 ), dif )
747 *
748 * Compare VT
749 *
750  dif = zero
751  IF( m.GT.0 .AND. n.GT.0 ) THEN
752  IF( ijq.EQ.1 ) THEN
753  IF( m.GE.n ) THEN
754  CALL dort03( 'R', n, mnmin, n, mnmin, vtsav,
755  $ ldvt, vt, ldvt, work, lwork,
756  $ dif, info )
757  ELSE
758  CALL dort03( 'R', n, mnmin, n, mnmin, vtsav,
759  $ ldvt, a, lda, work, lwork, dif,
760  $ info )
761  END IF
762  ELSE IF( ijq.EQ.2 ) THEN
763  CALL dort03( 'R', n, mnmin, n, mnmin, vtsav,
764  $ ldvt, vt, ldvt, work, lwork, dif,
765  $ info )
766  END IF
767  END IF
768  result( 13 ) = max( result( 13 ), dif )
769 *
770 * Compare S
771 *
772  dif = zero
773  div = max( mnmin*ulp*s( 1 ), unfl )
774  DO 100 i = 1, mnmin - 1
775  IF( ssav( i ).LT.ssav( i+1 ) )
776  $ dif = ulpinv
777  IF( ssav( i ).LT.zero )
778  $ dif = ulpinv
779  dif = max( dif, abs( ssav( i )-s( i ) ) / div )
780  100 CONTINUE
781  result( 14 ) = max( result( 14 ), dif )
782  110 CONTINUE
783 *
784 * Test DGESVJ: Factorize A
785 * Note: DGESVJ does not work for M < N
786 *
787  result( 15 ) = zero
788  result( 16 ) = zero
789  result( 17 ) = zero
790  result( 18 ) = zero
791 *
792  IF( m.GE.n ) THEN
793  iwtmp = 5*mnmin*mnmin + 9*mnmin + max( m, n )
794  lswork = iwtmp + ( iws-1 )*( lwork-iwtmp ) / 3
795  lswork = min( lswork, lwork )
796  lswork = max( lswork, 1 )
797  IF( iws.EQ.4 )
798  $ lswork = lwork
799 *
800  CALL dlacpy( 'F', m, n, asav, lda, usav, lda )
801  srnamt = 'DGESVJ'
802  CALL dgesvj( 'G', 'U', 'V', m, n, usav, lda, ssav,
803  & 0, a, ldvt, work, lwork, info )
804 *
805 * DGESVJ retuns V not VT, so we transpose to use the same
806 * test suite.
807 *
808  DO j=1,n
809  DO i=1,n
810  vtsav(j,i) = a(i,j)
811  END DO
812  END DO
813 *
814  IF( iinfo.NE.0 ) THEN
815  WRITE( nout, fmt = 9995 )'GESVJ', iinfo, m, n,
816  $ jtype, lswork, ioldsd
817  info = abs( iinfo )
818  RETURN
819  END IF
820 *
821 * Do tests 15--18
822 *
823  CALL dbdt01( m, n, 0, asav, lda, usav, ldu, ssav, e,
824  $ vtsav, ldvt, work, result( 15 ) )
825  IF( m.NE.0 .AND. n.NE.0 ) THEN
826  CALL dort01( 'Columns', m, m, usav, ldu, work,
827  $ lwork, result( 16 ) )
828  CALL dort01( 'Rows', n, n, vtsav, ldvt, work,
829  $ lwork, result( 17 ) )
830  END IF
831  result( 18 ) = zero
832  DO 120 i = 1, mnmin - 1
833  IF( ssav( i ).LT.ssav( i+1 ) )
834  $ result( 18 ) = ulpinv
835  IF( ssav( i ).LT.zero )
836  $ result( 18 ) = ulpinv
837  120 CONTINUE
838  IF( mnmin.GE.1 ) THEN
839  IF( ssav( mnmin ).LT.zero )
840  $ result( 18 ) = ulpinv
841  END IF
842  END IF
843 *
844 * Test DGEJSV: Factorize A
845 * Note: DGEJSV does not work for M < N
846 *
847  result( 19 ) = zero
848  result( 20 ) = zero
849  result( 21 ) = zero
850  result( 22 ) = zero
851  IF( m.GE.n ) THEN
852  iwtmp = 5*mnmin*mnmin + 9*mnmin + max( m, n )
853  lswork = iwtmp + ( iws-1 )*( lwork-iwtmp ) / 3
854  lswork = min( lswork, lwork )
855  lswork = max( lswork, 1 )
856  IF( iws.EQ.4 )
857  $ lswork = lwork
858 *
859  CALL dlacpy( 'F', m, n, asav, lda, vtsav, lda )
860  srnamt = 'DGEJSV'
861  CALL dgejsv( 'G', 'U', 'V', 'R', 'N', 'N',
862  & m, n, vtsav, lda, ssav, usav, ldu, a, ldvt,
863  & work, lwork, iwork, info )
864 *
865 * DGEJSV retuns V not VT, so we transpose to use the same
866 * test suite.
867 *
868  DO 140 j=1,n
869  DO 130 i=1,n
870  vtsav(j,i) = a(i,j)
871  130 END DO
872  140 END DO
873 *
874  IF( iinfo.NE.0 ) THEN
875  WRITE( nout, fmt = 9995 )'GESVJ', iinfo, m, n,
876  $ jtype, lswork, ioldsd
877  info = abs( iinfo )
878  RETURN
879  END IF
880 *
881 * Do tests 19--22
882 *
883  CALL dbdt01( m, n, 0, asav, lda, usav, ldu, ssav, e,
884  $ vtsav, ldvt, work, result( 19 ) )
885  IF( m.NE.0 .AND. n.NE.0 ) THEN
886  CALL dort01( 'Columns', m, m, usav, ldu, work,
887  $ lwork, result( 20 ) )
888  CALL dort01( 'Rows', n, n, vtsav, ldvt, work,
889  $ lwork, result( 21 ) )
890  END IF
891  result( 22 ) = zero
892  DO 150 i = 1, mnmin - 1
893  IF( ssav( i ).LT.ssav( i+1 ) )
894  $ result( 22 ) = ulpinv
895  IF( ssav( i ).LT.zero )
896  $ result( 22 ) = ulpinv
897  150 CONTINUE
898  IF( mnmin.GE.1 ) THEN
899  IF( ssav( mnmin ).LT.zero )
900  $ result( 22 ) = ulpinv
901  END IF
902  END IF
903 *
904 * Test DGESVDX
905 *
906  CALL dlacpy( 'F', m, n, asav, lda, a, lda )
907  CALL dgesvdx( 'V', 'V', 'A', m, n, a, lda,
908  $ vl, vu, il, iu, ns, ssav, usav, ldu,
909  $ vtsav, ldvt, work, lwork, iwork,
910  $ iinfo )
911  IF( iinfo.NE.0 ) THEN
912  WRITE( nout, fmt = 9995 )'GESVDX', iinfo, m, n,
913  $ jtype, lswork, ioldsd
914  info = abs( iinfo )
915  RETURN
916  END IF
917 *
918 * Do tests 23--29
919 *
920  result( 23 ) = zero
921  result( 24 ) = zero
922  result( 25 ) = zero
923  CALL dbdt01( m, n, 0, asav, lda, usav, ldu, ssav, e,
924  $ vtsav, ldvt, work, result( 23 ) )
925  IF( m.NE.0 .AND. n.NE.0 ) THEN
926  CALL dort01( 'Columns', m, m, usav, ldu, work, lwork,
927  $ result( 24 ) )
928  CALL dort01( 'Rows', n, n, vtsav, ldvt, work, lwork,
929  $ result( 25 ) )
930  END IF
931  result( 26 ) = zero
932  DO 160 i = 1, mnmin - 1
933  IF( ssav( i ).LT.ssav( i+1 ) )
934  $ result( 26 ) = ulpinv
935  IF( ssav( i ).LT.zero )
936  $ result( 26 ) = ulpinv
937  160 CONTINUE
938  IF( mnmin.GE.1 ) THEN
939  IF( ssav( mnmin ).LT.zero )
940  $ result( 26 ) = ulpinv
941  END IF
942 *
943 * Do partial SVDs, comparing to SSAV, USAV, and VTSAV
944 *
945  result( 27 ) = zero
946  result( 28 ) = zero
947  result( 29 ) = zero
948  DO 180 iju = 0, 1
949  DO 170 ijvt = 0, 1
950  IF( ( iju.EQ.0 .AND. ijvt.EQ.0 ) .OR.
951  $ ( iju.EQ.1 .AND. ijvt.EQ.1 ) )GO TO 170
952  jobu = cjobv( iju+1 )
953  jobvt = cjobv( ijvt+1 )
954  range = cjobr( 1 )
955  CALL dlacpy( 'F', m, n, asav, lda, a, lda )
956  CALL dgesvdx( jobu, jobvt, range, m, n, a, lda,
957  $ vl, vu, il, iu, ns, s, u, ldu,
958  $ vt, ldvt, work, lwork, iwork,
959  $ iinfo )
960 *
961 * Compare U
962 *
963  dif = zero
964  IF( m.GT.0 .AND. n.GT.0 ) THEN
965  IF( iju.EQ.1 ) THEN
966  CALL dort03( 'C', m, mnmin, m, mnmin, usav,
967  $ ldu, u, ldu, work, lwork, dif,
968  $ iinfo )
969  END IF
970  END IF
971  result( 27 ) = max( result( 27 ), dif )
972 *
973 * Compare VT
974 *
975  dif = zero
976  IF( m.GT.0 .AND. n.GT.0 ) THEN
977  IF( ijvt.EQ.1 ) THEN
978  CALL dort03( 'R', n, mnmin, n, mnmin, vtsav,
979  $ ldvt, vt, ldvt, work, lwork,
980  $ dif, iinfo )
981  END IF
982  END IF
983  result( 28 ) = max( result( 28 ), dif )
984 *
985 * Compare S
986 *
987  dif = zero
988  div = max( mnmin*ulp*s( 1 ), unfl )
989  DO 190 i = 1, mnmin - 1
990  IF( ssav( i ).LT.ssav( i+1 ) )
991  $ dif = ulpinv
992  IF( ssav( i ).LT.zero )
993  $ dif = ulpinv
994  dif = max( dif, abs( ssav( i )-s( i ) ) / div )
995  190 CONTINUE
996  result( 29 ) = max( result( 29 ), dif )
997  170 CONTINUE
998  180 CONTINUE
999 *
1000 * Do tests 30--32: DGESVDX( 'V', 'V', 'I' )
1001 *
1002  DO 200 i = 1, 4
1003  iseed2( i ) = iseed( i )
1004  200 CONTINUE
1005  IF( mnmin.LE.1 ) THEN
1006  il = 1
1007  iu = max( 1, mnmin )
1008  ELSE
1009  il = 1 + int( ( mnmin-1 )*dlarnd( 1, iseed2 ) )
1010  iu = 1 + int( ( mnmin-1 )*dlarnd( 1, iseed2 ) )
1011  IF( iu.LT.il ) THEN
1012  itemp = iu
1013  iu = il
1014  il = itemp
1015  END IF
1016  END IF
1017  CALL dlacpy( 'F', m, n, asav, lda, a, lda )
1018  CALL dgesvdx( 'V', 'V', 'I', m, n, a, lda,
1019  $ vl, vu, il, iu, nsi, s, u, ldu,
1020  $ vt, ldvt, work, lwork, iwork,
1021  $ iinfo )
1022  IF( iinfo.NE.0 ) THEN
1023  WRITE( nout, fmt = 9995 )'GESVDX', iinfo, m, n,
1024  $ jtype, lswork, ioldsd
1025  info = abs( iinfo )
1026  RETURN
1027  END IF
1028 *
1029  result( 30 ) = zero
1030  result( 31 ) = zero
1031  result( 32 ) = zero
1032  CALL dbdt05( m, n, asav, lda, s, nsi, u, ldu,
1033  $ vt, ldvt, work, result( 30 ) )
1034  CALL dort01( 'Columns', m, nsi, u, ldu, work, lwork,
1035  $ result( 31 ) )
1036  CALL dort01( 'Rows', nsi, n, vt, ldvt, work, lwork,
1037  $ result( 32 ) )
1038 *
1039 * Do tests 33--35: DGESVDX( 'V', 'V', 'V' )
1040 *
1041  IF( mnmin.GT.0 .AND. nsi.GT.1 ) THEN
1042  IF( il.NE.1 ) THEN
1043  vu = ssav( il ) +
1044  $ max( half*abs( ssav( il )-ssav( il-1 ) ),
1045  $ ulp*anorm, two*rtunfl )
1046  ELSE
1047  vu = ssav( 1 ) +
1048  $ max( half*abs( ssav( ns )-ssav( 1 ) ),
1049  $ ulp*anorm, two*rtunfl )
1050  END IF
1051  IF( iu.NE.ns ) THEN
1052  vl = ssav( iu ) - max( ulp*anorm, two*rtunfl,
1053  $ half*abs( ssav( iu+1 )-ssav( iu ) ) )
1054  ELSE
1055  vl = ssav( ns ) - max( ulp*anorm, two*rtunfl,
1056  $ half*abs( ssav( ns )-ssav( 1 ) ) )
1057  END IF
1058  vl = max( vl,zero )
1059  vu = max( vu,zero )
1060  IF( vl.GE.vu ) vu = max( vu*2, vu+vl+half )
1061  ELSE
1062  vl = zero
1063  vu = one
1064  END IF
1065  CALL dlacpy( 'F', m, n, asav, lda, a, lda )
1066  CALL dgesvdx( 'V', 'V', 'V', m, n, a, lda,
1067  $ vl, vu, il, iu, nsv, s, u, ldu,
1068  $ vt, ldvt, work, lwork, iwork,
1069  $ iinfo )
1070  IF( iinfo.NE.0 ) THEN
1071  WRITE( nout, fmt = 9995 )'GESVDX', iinfo, m, n,
1072  $ jtype, lswork, ioldsd
1073  info = abs( iinfo )
1074  RETURN
1075  END IF
1076 *
1077  result( 33 ) = zero
1078  result( 34 ) = zero
1079  result( 35 ) = zero
1080  CALL dbdt05( m, n, asav, lda, s, nsv, u, ldu,
1081  $ vt, ldvt, work, result( 33 ) )
1082  CALL dort01( 'Columns', m, nsv, u, ldu, work, lwork,
1083  $ result( 34 ) )
1084  CALL dort01( 'Rows', nsv, n, vt, ldvt, work, lwork,
1085  $ result( 35 ) )
1086 *
1087 * End of Loop -- Check for RESULT(j) > THRESH
1088 *
1089  DO 210 j = 1, 35
1090  IF( result( j ).GE.thresh ) THEN
1091  IF( nfail.EQ.0 ) THEN
1092  WRITE( nout, fmt = 9999 )
1093  WRITE( nout, fmt = 9998 )
1094  END IF
1095  WRITE( nout, fmt = 9997 )m, n, jtype, iws, ioldsd,
1096  $ j, result( j )
1097  nfail = nfail + 1
1098  END IF
1099  210 CONTINUE
1100  ntest = ntest + 35
1101  220 CONTINUE
1102  230 CONTINUE
1103  240 CONTINUE
1104 *
1105 * Summary
1106 *
1107  CALL alasvm( path, nout, nfail, ntest, 0 )
1108 *
1109  9999 FORMAT( ' SVD -- Real Singular Value Decomposition Driver ',
1110  $ / ' Matrix types (see DDRVBD for details):',
1111  $ / / ' 1 = Zero matrix', / ' 2 = Identity matrix',
1112  $ / ' 3 = Evenly spaced singular values near 1',
1113  $ / ' 4 = Evenly spaced singular values near underflow',
1114  $ / ' 5 = Evenly spaced singular values near overflow', / /
1115  $ ' Tests performed: ( A is dense, U and V are orthogonal,',
1116  $ / 19x, ' S is an array, and Upartial, VTpartial, and',
1117  $ / 19x, ' Spartial are partially computed U, VT and S),', / )
1118  9998 FORMAT( ' 1 = | A - U diag(S) VT | / ( |A| max(M,N) ulp ) ',
1119  $ / ' 2 = | I - U**T U | / ( M ulp ) ',
1120  $ / ' 3 = | I - VT VT**T | / ( N ulp ) ',
1121  $ / ' 4 = 0 if S contains min(M,N) nonnegative values in',
1122  $ ' decreasing order, else 1/ulp',
1123  $ / ' 5 = | U - Upartial | / ( M ulp )',
1124  $ / ' 6 = | VT - VTpartial | / ( N ulp )',
1125  $ / ' 7 = | S - Spartial | / ( min(M,N) ulp |S| )',
1126  $ / ' 8 = | A - U diag(S) VT | / ( |A| max(M,N) ulp ) ',
1127  $ / ' 9 = | I - U**T U | / ( M ulp ) ',
1128  $ / '10 = | I - VT VT**T | / ( N ulp ) ',
1129  $ / '11 = 0 if S contains min(M,N) nonnegative values in',
1130  $ ' decreasing order, else 1/ulp',
1131  $ / '12 = | U - Upartial | / ( M ulp )',
1132  $ / '13 = | VT - VTpartial | / ( N ulp )',
1133  $ / '14 = | S - Spartial | / ( min(M,N) ulp |S| )',
1134  $ / '15 = | A - U diag(S) VT | / ( |A| max(M,N) ulp ) ',
1135  $ / '16 = | I - U**T U | / ( M ulp ) ',
1136  $ / '17 = | I - VT VT**T | / ( N ulp ) ',
1137  $ / '18 = 0 if S contains min(M,N) nonnegative values in',
1138  $ ' decreasing order, else 1/ulp',
1139  $ / '19 = | U - Upartial | / ( M ulp )',
1140  $ / '20 = | VT - VTpartial | / ( N ulp )',
1141  $ / '21 = | S - Spartial | / ( min(M,N) ulp |S| )',
1142  $ / '22 = 0 if S contains min(M,N) nonnegative values in',
1143  $ ' decreasing order, else 1/ulp',
1144  $ / '23 = | A - U diag(S) VT | / ( |A| max(M,N) ulp ),',
1145  $ ' DGESVDX(V,V,A) ',
1146  $ / '24 = | I - U**T U | / ( M ulp ) ',
1147  $ / '25 = | I - VT VT**T | / ( N ulp ) ',
1148  $ / '26 = 0 if S contains min(M,N) nonnegative values in',
1149  $ ' decreasing order, else 1/ulp',
1150  $ / '27 = | U - Upartial | / ( M ulp )',
1151  $ / '28 = | VT - VTpartial | / ( N ulp )',
1152  $ / '29 = | S - Spartial | / ( min(M,N) ulp |S| )',
1153  $ / '30 = | U**T A VT**T - diag(S) | / ( |A| max(M,N) ulp ),',
1154  $ ' DGESVDX(V,V,I) ',
1155  $ / '31 = | I - U**T U | / ( M ulp ) ',
1156  $ / '32 = | I - VT VT**T | / ( N ulp ) ',
1157  $ / '33 = | U**T A VT**T - diag(S) | / ( |A| max(M,N) ulp ),',
1158  $ ' DGESVDX(V,V,V) ',
1159  $ / '34 = | I - U**T U | / ( M ulp ) ',
1160  $ / '35 = | I - VT VT**T | / ( N ulp ) ',
1161  $ / / )
1162  9997 FORMAT( ' M=', i5, ', N=', i5, ', type ', i1, ', IWS=', i1,
1163  $ ', seed=', 4( i4, ',' ), ' test(', i2, ')=', g11.4 )
1164  9996 FORMAT( ' DDRVBD: ', a, ' returned INFO=', i6, '.', / 9x, 'M=',
1165  $ i6, ', N=', i6, ', JTYPE=', i6, ', ISEED=(', 3( i5, ',' ),
1166  $ i5, ')' )
1167  9995 FORMAT( ' DDRVBD: ', a, ' returned INFO=', i6, '.', / 9x, 'M=',
1168  $ i6, ', N=', i6, ', JTYPE=', i6, ', LSWORK=', i6, / 9x,
1169  $ 'ISEED=(', 3( i5, ',' ), i5, ')' )
1170 *
1171  RETURN
1172 *
1173 * End of DDRVBD
1174 *
alasvm
subroutine alasvm(TYPE, NOUT, NFAIL, NRUN, NERRS)
ALASVM
Definition: alasvm.f:75
dlaset
subroutine dlaset(UPLO, M, N, ALPHA, BETA, A, LDA)
DLASET initializes the off-diagonal elements and the diagonal elements of a matrix to given values...
Definition: dlaset.f:112
dlamch
double precision function dlamch(CMACH)
DLAMCH
Definition: dlamch.f:65
dgesvj
subroutine dgesvj(JOBA, JOBU, JOBV, M, N, A, LDA, SVA, MV, V, LDV, WORK, LWORK, INFO)
DGESVJ
Definition: dgesvj.f:339
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
dbdt01
subroutine dbdt01(M, N, KD, A, LDA, Q, LDQ, D, E, PT, LDPT, WORK, RESID)
DBDT01
Definition: dbdt01.f:142
xerbla
subroutine xerbla(SRNAME, INFO)
XERBLA
Definition: xerbla.f:62
dlabad
subroutine dlabad(SMALL, LARGE)
DLABAD
Definition: dlabad.f:76
dgesdd
subroutine dgesdd(JOBZ, M, N, A, LDA, S, U, LDU, VT, LDVT, WORK, LWORK, IWORK, INFO)
DGESDD
Definition: dgesdd.f:221
dlarnd
double precision function dlarnd(IDIST, ISEED)
DLARND
Definition: dlarnd.f:75
dbdt05
subroutine dbdt05(M, N, A, LDA, S, NS, U, LDU, VT, LDVT, WORK, RESID)
Definition: dbdt05.f:126
dgesvdx
subroutine dgesvdx(JOBU, JOBVT, RANGE, M, N, A, LDA, VL, VU, IL, IU, NS, S, U, LDU, VT, LDVT, WORK, LWORK, IWORK, INFO)
DGESVDX computes the singular value decomposition (SVD) for GE matrices
Definition: dgesvdx.f:265
dlatms
subroutine dlatms(M, N, DIST, ISEED, SYM, D, MODE, COND, DMAX, KL, KU, PACK, A, LDA, WORK, INFO)
DLATMS
Definition: dlatms.f:323
dort01
subroutine dort01(ROWCOL, M, N, U, LDU, WORK, LWORK, RESID)
DORT01
Definition: dort01.f:118
dgesvd
subroutine dgesvd(JOBU, JOBVT, M, N, A, LDA, S, U, LDU, VT, LDVT, WORK, LWORK, INFO)
DGESVD computes the singular value decomposition (SVD) for GE matrices
Definition: dgesvd.f:213
dort03
subroutine dort03(RC, MU, MV, N, K, U, LDU, V, LDV, WORK, LWORK, RESULT, INFO)
DORT03
Definition: dort03.f:158
dgejsv
subroutine dgejsv(JOBA, JOBU, JOBV, JOBR, JOBT, JOBP, M, N, A, LDA, SVA, U, LDU, V, LDV, WORK, LWORK, IWORK, INFO)
DGEJSV
Definition: dgejsv.f:478

Here is the call graph for this function:

Here is the caller graph for this function: