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

SDRVBD

Purpose:
 SDRVBD checks the singular value decomposition (SVD) drivers
 SGESVD, SGESDD, SGESVJ, and SGEJSV.

 Both SGESVD and SGESDD 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 SDRVBD 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 SGESVD:

 (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 SGESDD:

 (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 SGESVJ:

 (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 SGEJSV:

 (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 SGESVDX( 'V', 'V', 'A' )/SGESVDX( '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 SGESVDX( '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 SGESVDX( '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, SDRVBD
          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
          SDRVBD to continue the same random number sequence.
[in]THRESH
          THRESH is REAL
          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 REAL 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 REAL array, dimension (LDU,MMAX)
[in]LDU
          LDU is INTEGER
          The leading dimension of the array U.  LDU >= max(1,MMAX).
[out]VT
          VT is REAL array, dimension (LDVT,NMAX)
[in]LDVT
          LDVT is INTEGER
          The leading dimension of the array VT.  LDVT >= max(1,NMAX).
[out]ASAV
          ASAV is REAL array, dimension (LDA,NMAX)
[out]USAV
          USAV is REAL array, dimension (LDU,MMAX)
[out]VTSAV
          VTSAV is REAL array, dimension (LDVT,NMAX)
[out]S
          S is REAL array, dimension
                      (max(min(MM,NN)))
[out]SSAV
          SSAV is REAL array, dimension
                      (max(min(MM,NN)))
[out]E
          E is REAL array, dimension
                      (max(min(MM,NN)))
[out]WORK
          WORK is REAL 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  SLATMS, or SGESVD 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 sdrvbd.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  REAL thresh
367 * ..
368 * .. Array Arguments ..
369  LOGICAL dotype( * )
370  INTEGER iseed( 4 ), iwork( * ), mm( * ), nn( * )
371  REAL 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  REAL zero, one, two, half
380  parameter ( zero = 0.0e0, one = 1.0e0, two = 2.0e0,
381  $ half = 0.5e0 )
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  REAL 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  REAL result( 40 )
400 * ..
401 * .. External Functions ..
402  REAL slamch, slarnd
403  EXTERNAL slamch, slarnd
404 * ..
405 * .. External Subroutines ..
406  EXTERNAL alasvm, sbdt01, sgejsv, sgesdd, sgesvd,
409 * ..
410 * .. Intrinsic Functions ..
411  INTRINSIC abs, REAL, 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( 'SDRVBD', -info )
473  RETURN
474  END IF
475 *
476 * Initialize constants
477 *
478  path( 1: 1 ) = 'Single precision'
479  path( 2: 3 ) = 'BD'
480  nfail = 0
481  ntest = 0
482  unfl = slamch( 'Safe minimum' )
483  ovfl = one / unfl
484  CALL slabad( unfl, ovfl )
485  ulp = slamch( '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 slaset( 'Full', m, n, zero, zero, a, lda )
521 *
522  ELSE IF( jtype.EQ.2 ) THEN
523 *
524 * Identity matrix
525 *
526  CALL slaset( '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 slatms( m, n, 'U', iseed, 'N', s, 4, REAL( 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 slacpy( '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 SGESVD: 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 slacpy( 'F', m, n, asav, lda, a, lda )
570  srnamt = 'SGESVD'
571  CALL sgesvd( '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 sbdt01( 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 sort01( 'Columns', m, m, usav, ldu, work, lwork,
586  $ result( 2 ) )
587  CALL sort01( '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 slacpy( 'F', m, n, asav, lda, a, lda )
614  srnamt = 'SGESVD'
615  CALL sgesvd( 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 sort03( 'C', m, mnmin, m, mnmin, usav,
624  $ ldu, a, lda, work, lwork, dif,
625  $ iinfo )
626  ELSE IF( iju.EQ.2 ) THEN
627  CALL sort03( 'C', m, mnmin, m, mnmin, usav,
628  $ ldu, u, ldu, work, lwork, dif,
629  $ iinfo )
630  ELSE IF( iju.EQ.3 ) THEN
631  CALL sort03( '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 sort03( 'R', n, mnmin, n, mnmin, vtsav,
644  $ ldvt, a, lda, work, lwork, dif,
645  $ iinfo )
646  ELSE IF( ijvt.EQ.2 ) THEN
647  CALL sort03( 'R', n, mnmin, n, mnmin, vtsav,
648  $ ldvt, vt, ldvt, work, lwork,
649  $ dif, iinfo )
650  ELSE IF( ijvt.EQ.3 ) THEN
651  CALL sort03( '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 SGESDD: 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 slacpy( 'F', m, n, asav, lda, a, lda )
683  srnamt = 'SGESDD'
684  CALL sgesdd( '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 sbdt01( 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 sort01( 'Columns', m, m, usav, ldu, work, lwork,
699  $ result( 9 ) )
700  CALL sort01( '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 slacpy( 'F', m, n, asav, lda, a, lda )
723  srnamt = 'SGESDD'
724  CALL sgesdd( 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 sort03( 'C', m, mnmin, m, mnmin, usav,
734  $ ldu, a, lda, work, lwork, dif,
735  $ info )
736  ELSE
737  CALL sort03( '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 sort03( '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 sort03( 'R', n, mnmin, n, mnmin, vtsav,
755  $ ldvt, vt, ldvt, work, lwork,
756  $ dif, info )
757  ELSE
758  CALL sort03( '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 sort03( '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 SGESVJ: Factorize A
785 * Note: SGESVJ 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 slacpy( 'F', m, n, asav, lda, usav, lda )
801  srnamt = 'SGESVJ'
802  CALL sgesvj( 'G', 'U', 'V', m, n, usav, lda, ssav,
803  & 0, a, ldvt, work, lwork, info )
804 *
805 * SGESVJ 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 sbdt01( 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 sort01( 'Columns', m, m, usav, ldu, work,
827  $ lwork, result( 16 ) )
828  CALL sort01( '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 SGEJSV: Factorize A
845 * Note: SGEJSV 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 slacpy( 'F', m, n, asav, lda, vtsav, lda )
860  srnamt = 'SGEJSV'
861  CALL sgejsv( 'G', 'U', 'V', 'R', 'N', 'N',
862  & m, n, vtsav, lda, ssav, usav, ldu, a, ldvt,
863  & work, lwork, iwork, info )
864 *
865 * SGEJSV 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 sbdt01( 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 sort01( 'Columns', m, m, usav, ldu, work,
887  $ lwork, result( 20 ) )
888  CALL sort01( '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 SGESVDX
905 *
906  CALL slacpy( 'F', m, n, asav, lda, a, lda )
907  CALL sgesvdx( '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 sbdt01( 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 sort01( 'Columns', m, m, usav, ldu, work, lwork,
927  $ result( 24 ) )
928  CALL sort01( '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 slacpy( 'F', m, n, asav, lda, a, lda )
956  CALL sgesvdx( 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 sort03( '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 sort03( '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: SGESVDX( '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 )*slarnd( 1, iseed2 ) )
1010  iu = 1 + int( ( mnmin-1 )*slarnd( 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 slacpy( 'F', m, n, asav, lda, a, lda )
1018  CALL sgesvdx( '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 sbdt05( m, n, asav, lda, s, nsi, u, ldu,
1033  $ vt, ldvt, work, result( 30 ) )
1034  CALL sort01( 'Columns', m, nsi, u, ldu, work, lwork,
1035  $ result( 31 ) )
1036  CALL sort01( 'Rows', nsi, n, vt, ldvt, work, lwork,
1037  $ result( 32 ) )
1038 *
1039 * Do tests 33--35: SGESVDX( '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 slacpy( 'F', m, n, asav, lda, a, lda )
1066  CALL sgesvdx( '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 sbdt05( m, n, asav, lda, s, nsv, u, ldu,
1081  $ vt, ldvt, work, result( 33 ) )
1082  CALL sort01( 'Columns', m, nsv, u, ldu, work, lwork,
1083  $ result( 34 ) )
1084  CALL sort01( '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 SDRVBD 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  $ ' SGESVDX(V,V,A) ',
1145  $ / '23 = | A - U diag(S) VT | / ( |A| max(M,N) ulp ),'
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  $ ' SGESVDX(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  $ ' SGESVDX(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( ' SDRVBD: ', a, ' returned INFO=', i6, '.', / 9x, 'M=',
1165  $ i6, ', N=', i6, ', JTYPE=', i6, ', ISEED=(', 3( i5, ',' ),
1166  $ i5, ')' )
1167  9995 FORMAT( ' SDRVBD: ', 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 SDRVBD
1174 *
subroutine alasvm(TYPE, NOUT, NFAIL, NRUN, NERRS)
ALASVM
Definition: alasvm.f:75
subroutine sgesvj(JOBA, JOBU, JOBV, M, N, A, LDA, SVA, MV, V, LDV, WORK, LWORK, INFO)
SGESVJ
Definition: sgesvj.f:325
subroutine slabad(SMALL, LARGE)
SLABAD
Definition: slabad.f:76
subroutine xerbla(SRNAME, INFO)
XERBLA
Definition: xerbla.f:62
subroutine slacpy(UPLO, M, N, A, LDA, B, LDB)
SLACPY copies all or part of one two-dimensional array to another.
Definition: slacpy.f:105
subroutine sort01(ROWCOL, M, N, U, LDU, WORK, LWORK, RESID)
SORT01
Definition: sort01.f:118
subroutine slatms(M, N, DIST, ISEED, SYM, D, MODE, COND, DMAX, KL, KU, PACK, A, LDA, WORK, INFO)
SLATMS
Definition: slatms.f:323
subroutine sbdt01(M, N, KD, A, LDA, Q, LDQ, D, E, PT, LDPT, WORK, RESID)
SBDT01
Definition: sbdt01.f:142
subroutine sort03(RC, MU, MV, N, K, U, LDU, V, LDV, WORK, LWORK, RESULT, INFO)
SORT03
Definition: sort03.f:158
subroutine slaset(UPLO, M, N, ALPHA, BETA, A, LDA)
SLASET initializes the off-diagonal elements and the diagonal elements of a matrix to given values...
Definition: slaset.f:112
subroutine sgejsv(JOBA, JOBU, JOBV, JOBR, JOBT, JOBP, M, N, A, LDA, SVA, U, LDU, V, LDV, WORK, LWORK, IWORK, INFO)
SGEJSV
Definition: sgejsv.f:478
subroutine sgesvdx(JOBU, JOBVT, RANGE, M, N, A, LDA, VL, VU, IL, IU, NS, S, U, LDU, VT, LDVT, WORK, LWORK, IWORK, INFO)
SGESVDX computes the singular value decomposition (SVD) for GE matrices
Definition: sgesvdx.f:265
real function slarnd(IDIST, ISEED)
SLARND
Definition: slarnd.f:75
subroutine sgesvd(JOBU, JOBVT, M, N, A, LDA, S, U, LDU, VT, LDVT, WORK, LWORK, INFO)
SGESVD computes the singular value decomposition (SVD) for GE matrices
Definition: sgesvd.f:213
subroutine sgesdd(JOBZ, M, N, A, LDA, S, U, LDU, VT, LDVT, WORK, LWORK, IWORK, INFO)
SGESDD
Definition: sgesdd.f:220
real function slamch(CMACH)
SLAMCH
Definition: slamch.f:69
subroutine sbdt05(M, N, A, LDA, S, NS, U, LDU, VT, LDVT, WORK, RESID)
Definition: sbdt05.f:126

Here is the call graph for this function:

Here is the caller graph for this function: