LAPACK  3.5.0
LAPACK: Linear Algebra PACKage
ddrvbd.f File Reference

Go to the source code of this file.

Functions/Subroutines

subroutine ddrvbd (NSIZES, MM, NN, NTYPES, DOTYPE, ISEED, THRESH, A, LDA, U, LDU, VT, LDVT, ASAV, USAV, VTSAV, S, SSAV, E, WORK, LWORK, IWORK, NOUT, INFO)
 DDRVBD More...
 

Function/Subroutine Documentation

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 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.)

 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
November 2011

Definition at line 321 of file ddrvbd.f.

321 *
322 * -- LAPACK test routine (version 3.4.0) --
323 * -- LAPACK is a software package provided by Univ. of Tennessee, --
324 * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
325 * November 2011
326 *
327 * .. Scalar Arguments ..
328  INTEGER info, lda, ldu, ldvt, lwork, nout, nsizes,
329  $ ntypes
330  DOUBLE PRECISION thresh
331 * ..
332 * .. Array Arguments ..
333  LOGICAL dotype( * )
334  INTEGER iseed( 4 ), iwork( * ), mm( * ), nn( * )
335  DOUBLE PRECISION a( lda, * ), asav( lda, * ), e( * ), s( * ),
336  $ ssav( * ), u( ldu, * ), usav( ldu, * ),
337  $ vt( ldvt, * ), vtsav( ldvt, * ), work( * )
338 * ..
339 *
340 * =====================================================================
341 *
342 * .. Parameters ..
343  DOUBLE PRECISION zero, one
344  parameter( zero = 0.0d0, one = 1.0d0 )
345  INTEGER maxtyp
346  parameter( maxtyp = 5 )
347 * ..
348 * .. Local Scalars ..
349  LOGICAL badmm, badnn
350  CHARACTER jobq, jobu, jobvt
351  CHARACTER*3 path
352  INTEGER i, iinfo, ijq, iju, ijvt, iws, iwtmp, j, jsize,
353  $ jtype, lswork, m, minwrk, mmax, mnmax, mnmin,
354  $ mtypes, n, nfail, nmax, ntest
355  DOUBLE PRECISION anorm, dif, div, ovfl, ulp, ulpinv, unfl
356 * ..
357 * .. Local Arrays ..
358  CHARACTER cjob( 4 )
359  INTEGER ioldsd( 4 )
360  DOUBLE PRECISION result( 22 )
361 * ..
362 * .. External Functions ..
363  DOUBLE PRECISION dlamch
364  EXTERNAL dlamch
365 * ..
366 * .. External Subroutines ..
367  EXTERNAL alasvm, dbdt01, dgesdd, dgesvd, dlabad, dlacpy,
369  $ dgejsv
370 * ..
371 * .. Intrinsic Functions ..
372  INTRINSIC abs, dble, max, min
373 * ..
374 * .. Scalars in Common ..
375  LOGICAL lerr, ok
376  CHARACTER*32 srnamt
377  INTEGER infot, nunit
378 * ..
379 * .. Common blocks ..
380  COMMON / infoc / infot, nunit, ok, lerr
381  COMMON / srnamc / srnamt
382 * ..
383 * .. Data statements ..
384  DATA cjob / 'N', 'O', 'S', 'A' /
385 * ..
386 * .. Executable Statements ..
387 *
388 * Check for errors
389 *
390  info = 0
391  badmm = .false.
392  badnn = .false.
393  mmax = 1
394  nmax = 1
395  mnmax = 1
396  minwrk = 1
397  DO 10 j = 1, nsizes
398  mmax = max( mmax, mm( j ) )
399  IF( mm( j ).LT.0 )
400  $ badmm = .true.
401  nmax = max( nmax, nn( j ) )
402  IF( nn( j ).LT.0 )
403  $ badnn = .true.
404  mnmax = max( mnmax, min( mm( j ), nn( j ) ) )
405  minwrk = max( minwrk, max( 3*min( mm( j ),
406  $ nn( j ) )+max( mm( j ), nn( j ) ), 5*min( mm( j ),
407  $ nn( j )-4 ) )+2*min( mm( j ), nn( j ) )**2 )
408  10 CONTINUE
409 *
410 * Check for errors
411 *
412  IF( nsizes.LT.0 ) THEN
413  info = -1
414  ELSE IF( badmm ) THEN
415  info = -2
416  ELSE IF( badnn ) THEN
417  info = -3
418  ELSE IF( ntypes.LT.0 ) THEN
419  info = -4
420  ELSE IF( lda.LT.max( 1, mmax ) ) THEN
421  info = -10
422  ELSE IF( ldu.LT.max( 1, mmax ) ) THEN
423  info = -12
424  ELSE IF( ldvt.LT.max( 1, nmax ) ) THEN
425  info = -14
426  ELSE IF( minwrk.GT.lwork ) THEN
427  info = -21
428  END IF
429 *
430  IF( info.NE.0 ) THEN
431  CALL xerbla( 'DDRVBD', -info )
432  RETURN
433  END IF
434 *
435 * Initialize constants
436 *
437  path( 1: 1 ) = 'Double precision'
438  path( 2: 3 ) = 'BD'
439  nfail = 0
440  ntest = 0
441  unfl = dlamch( 'Safe minimum' )
442  ovfl = one / unfl
443  CALL dlabad( unfl, ovfl )
444  ulp = dlamch( 'Precision' )
445  ulpinv = one / ulp
446  infot = 0
447 *
448 * Loop over sizes, types
449 *
450  DO 150 jsize = 1, nsizes
451  m = mm( jsize )
452  n = nn( jsize )
453  mnmin = min( m, n )
454 *
455  IF( nsizes.NE.1 ) THEN
456  mtypes = min( maxtyp, ntypes )
457  ELSE
458  mtypes = min( maxtyp+1, ntypes )
459  END IF
460 *
461  DO 140 jtype = 1, mtypes
462  IF( .NOT.dotype( jtype ) )
463  $ GO TO 140
464 *
465  DO 20 j = 1, 4
466  ioldsd( j ) = iseed( j )
467  20 CONTINUE
468 *
469 * Compute "A"
470 *
471  IF( mtypes.GT.maxtyp )
472  $ GO TO 30
473 *
474  IF( jtype.EQ.1 ) THEN
475 *
476 * Zero matrix
477 *
478  CALL dlaset( 'Full', m, n, zero, zero, a, lda )
479 *
480  ELSE IF( jtype.EQ.2 ) THEN
481 *
482 * Identity matrix
483 *
484  CALL dlaset( 'Full', m, n, zero, one, a, lda )
485 *
486  ELSE
487 *
488 * (Scaled) random matrix
489 *
490  IF( jtype.EQ.3 )
491  $ anorm = one
492  IF( jtype.EQ.4 )
493  $ anorm = unfl / ulp
494  IF( jtype.EQ.5 )
495  $ anorm = ovfl*ulp
496  CALL dlatms( m, n, 'U', iseed, 'N', s, 4, dble( mnmin ),
497  $ anorm, m-1, n-1, 'N', a, lda, work, iinfo )
498  IF( iinfo.NE.0 ) THEN
499  WRITE( nout, fmt = 9996 )'Generator', iinfo, m, n,
500  $ jtype, ioldsd
501  info = abs( iinfo )
502  RETURN
503  END IF
504  END IF
505 *
506  30 CONTINUE
507  CALL dlacpy( 'F', m, n, a, lda, asav, lda )
508 *
509 * Do for minimal and adequate (for blocking) workspace
510 *
511  DO 130 iws = 1, 4
512 *
513  DO 40 j = 1, 14
514  result( j ) = -one
515  40 CONTINUE
516 *
517 * Test DGESVD: Factorize A
518 *
519  iwtmp = max( 3*min( m, n )+max( m, n ), 5*min( m, n ) )
520  lswork = iwtmp + ( iws-1 )*( lwork-iwtmp ) / 3
521  lswork = min( lswork, lwork )
522  lswork = max( lswork, 1 )
523  IF( iws.EQ.4 )
524  $ lswork = lwork
525 *
526  IF( iws.GT.1 )
527  $ CALL dlacpy( 'F', m, n, asav, lda, a, lda )
528  srnamt = 'DGESVD'
529  CALL dgesvd( 'A', 'A', m, n, a, lda, ssav, usav, ldu,
530  $ vtsav, ldvt, work, lswork, iinfo )
531  IF( iinfo.NE.0 ) THEN
532  WRITE( nout, fmt = 9995 )'GESVD', iinfo, m, n, jtype,
533  $ lswork, ioldsd
534  info = abs( iinfo )
535  RETURN
536  END IF
537 *
538 * Do tests 1--4
539 *
540  CALL dbdt01( m, n, 0, asav, lda, usav, ldu, ssav, e,
541  $ vtsav, ldvt, work, result( 1 ) )
542  IF( m.NE.0 .AND. n.NE.0 ) THEN
543  CALL dort01( 'Columns', m, m, usav, ldu, work, lwork,
544  $ result( 2 ) )
545  CALL dort01( 'Rows', n, n, vtsav, ldvt, work, lwork,
546  $ result( 3 ) )
547  END IF
548  result( 4 ) = zero
549  DO 50 i = 1, mnmin - 1
550  IF( ssav( i ).LT.ssav( i+1 ) )
551  $ result( 4 ) = ulpinv
552  IF( ssav( i ).LT.zero )
553  $ result( 4 ) = ulpinv
554  50 CONTINUE
555  IF( mnmin.GE.1 ) THEN
556  IF( ssav( mnmin ).LT.zero )
557  $ result( 4 ) = ulpinv
558  END IF
559 *
560 * Do partial SVDs, comparing to SSAV, USAV, and VTSAV
561 *
562  result( 5 ) = zero
563  result( 6 ) = zero
564  result( 7 ) = zero
565  DO 80 iju = 0, 3
566  DO 70 ijvt = 0, 3
567  IF( ( iju.EQ.3 .AND. ijvt.EQ.3 ) .OR.
568  $ ( iju.EQ.1 .AND. ijvt.EQ.1 ) )GO TO 70
569  jobu = cjob( iju+1 )
570  jobvt = cjob( ijvt+1 )
571  CALL dlacpy( 'F', m, n, asav, lda, a, lda )
572  srnamt = 'DGESVD'
573  CALL dgesvd( jobu, jobvt, m, n, a, lda, s, u, ldu,
574  $ vt, ldvt, work, lswork, iinfo )
575 *
576 * Compare U
577 *
578  dif = zero
579  IF( m.GT.0 .AND. n.GT.0 ) THEN
580  IF( iju.EQ.1 ) THEN
581  CALL dort03( 'C', m, mnmin, m, mnmin, usav,
582  $ ldu, a, lda, work, lwork, dif,
583  $ iinfo )
584  ELSE IF( iju.EQ.2 ) THEN
585  CALL dort03( 'C', m, mnmin, m, mnmin, usav,
586  $ ldu, u, ldu, work, lwork, dif,
587  $ iinfo )
588  ELSE IF( iju.EQ.3 ) THEN
589  CALL dort03( 'C', m, m, m, mnmin, usav, ldu,
590  $ u, ldu, work, lwork, dif,
591  $ iinfo )
592  END IF
593  END IF
594  result( 5 ) = max( result( 5 ), dif )
595 *
596 * Compare VT
597 *
598  dif = zero
599  IF( m.GT.0 .AND. n.GT.0 ) THEN
600  IF( ijvt.EQ.1 ) THEN
601  CALL dort03( 'R', n, mnmin, n, mnmin, vtsav,
602  $ ldvt, a, lda, work, lwork, dif,
603  $ iinfo )
604  ELSE IF( ijvt.EQ.2 ) THEN
605  CALL dort03( 'R', n, mnmin, n, mnmin, vtsav,
606  $ ldvt, vt, ldvt, work, lwork,
607  $ dif, iinfo )
608  ELSE IF( ijvt.EQ.3 ) THEN
609  CALL dort03( 'R', n, n, n, mnmin, vtsav,
610  $ ldvt, vt, ldvt, work, lwork,
611  $ dif, iinfo )
612  END IF
613  END IF
614  result( 6 ) = max( result( 6 ), dif )
615 *
616 * Compare S
617 *
618  dif = zero
619  div = max( dble( mnmin )*ulp*s( 1 ), unfl )
620  DO 60 i = 1, mnmin - 1
621  IF( ssav( i ).LT.ssav( i+1 ) )
622  $ dif = ulpinv
623  IF( ssav( i ).LT.zero )
624  $ dif = ulpinv
625  dif = max( dif, abs( ssav( i )-s( i ) ) / div )
626  60 CONTINUE
627  result( 7 ) = max( result( 7 ), dif )
628  70 CONTINUE
629  80 CONTINUE
630 *
631 * Test DGESDD: Factorize A
632 *
633  iwtmp = 5*mnmin*mnmin + 9*mnmin + max( m, n )
634  lswork = iwtmp + ( iws-1 )*( lwork-iwtmp ) / 3
635  lswork = min( lswork, lwork )
636  lswork = max( lswork, 1 )
637  IF( iws.EQ.4 )
638  $ lswork = lwork
639 *
640  CALL dlacpy( 'F', m, n, asav, lda, a, lda )
641  srnamt = 'DGESDD'
642  CALL dgesdd( 'A', m, n, a, lda, ssav, usav, ldu, vtsav,
643  $ ldvt, work, lswork, iwork, iinfo )
644  IF( iinfo.NE.0 ) THEN
645  WRITE( nout, fmt = 9995 )'GESDD', iinfo, m, n, jtype,
646  $ lswork, ioldsd
647  info = abs( iinfo )
648  RETURN
649  END IF
650 *
651 * Do tests 8--11
652 *
653  CALL dbdt01( m, n, 0, asav, lda, usav, ldu, ssav, e,
654  $ vtsav, ldvt, work, result( 8 ) )
655  IF( m.NE.0 .AND. n.NE.0 ) THEN
656  CALL dort01( 'Columns', m, m, usav, ldu, work, lwork,
657  $ result( 9 ) )
658  CALL dort01( 'Rows', n, n, vtsav, ldvt, work, lwork,
659  $ result( 10 ) )
660  END IF
661  result( 11 ) = zero
662  DO 90 i = 1, mnmin - 1
663  IF( ssav( i ).LT.ssav( i+1 ) )
664  $ result( 11 ) = ulpinv
665  IF( ssav( i ).LT.zero )
666  $ result( 11 ) = ulpinv
667  90 CONTINUE
668  IF( mnmin.GE.1 ) THEN
669  IF( ssav( mnmin ).LT.zero )
670  $ result( 11 ) = ulpinv
671  END IF
672 *
673 * Do partial SVDs, comparing to SSAV, USAV, and VTSAV
674 *
675  result( 12 ) = zero
676  result( 13 ) = zero
677  result( 14 ) = zero
678  DO 110 ijq = 0, 2
679  jobq = cjob( ijq+1 )
680  CALL dlacpy( 'F', m, n, asav, lda, a, lda )
681  srnamt = 'DGESDD'
682  CALL dgesdd( jobq, m, n, a, lda, s, u, ldu, vt, ldvt,
683  $ work, lswork, iwork, iinfo )
684 *
685 * Compare U
686 *
687  dif = zero
688  IF( m.GT.0 .AND. n.GT.0 ) THEN
689  IF( ijq.EQ.1 ) THEN
690  IF( m.GE.n ) THEN
691  CALL dort03( 'C', m, mnmin, m, mnmin, usav,
692  $ ldu, a, lda, work, lwork, dif,
693  $ info )
694  ELSE
695  CALL dort03( 'C', m, mnmin, m, mnmin, usav,
696  $ ldu, u, ldu, work, lwork, dif,
697  $ info )
698  END IF
699  ELSE IF( ijq.EQ.2 ) THEN
700  CALL dort03( 'C', m, mnmin, m, mnmin, usav, ldu,
701  $ u, ldu, work, lwork, dif, info )
702  END IF
703  END IF
704  result( 12 ) = max( result( 12 ), dif )
705 *
706 * Compare VT
707 *
708  dif = zero
709  IF( m.GT.0 .AND. n.GT.0 ) THEN
710  IF( ijq.EQ.1 ) THEN
711  IF( m.GE.n ) THEN
712  CALL dort03( 'R', n, mnmin, n, mnmin, vtsav,
713  $ ldvt, vt, ldvt, work, lwork,
714  $ dif, info )
715  ELSE
716  CALL dort03( 'R', n, mnmin, n, mnmin, vtsav,
717  $ ldvt, a, lda, work, lwork, dif,
718  $ info )
719  END IF
720  ELSE IF( ijq.EQ.2 ) THEN
721  CALL dort03( 'R', n, mnmin, n, mnmin, vtsav,
722  $ ldvt, vt, ldvt, work, lwork, dif,
723  $ info )
724  END IF
725  END IF
726  result( 13 ) = max( result( 13 ), dif )
727 *
728 * Compare S
729 *
730  dif = zero
731  div = max( dble( mnmin )*ulp*s( 1 ), unfl )
732  DO 100 i = 1, mnmin - 1
733  IF( ssav( i ).LT.ssav( i+1 ) )
734  $ dif = ulpinv
735  IF( ssav( i ).LT.zero )
736  $ dif = ulpinv
737  dif = max( dif, abs( ssav( i )-s( i ) ) / div )
738  100 CONTINUE
739  result( 14 ) = max( result( 14 ), dif )
740  110 CONTINUE
741 *
742 * Test DGESVJ: Factorize A
743 * Note: DGESVJ does not work for M < N
744 *
745  result( 15 ) = zero
746  result( 16 ) = zero
747  result( 17 ) = zero
748  result( 18 ) = zero
749 *
750  IF( m.GE.n ) THEN
751  iwtmp = 5*mnmin*mnmin + 9*mnmin + max( m, n )
752  lswork = iwtmp + ( iws-1 )*( lwork-iwtmp ) / 3
753  lswork = min( lswork, lwork )
754  lswork = max( lswork, 1 )
755  IF( iws.EQ.4 )
756  $ lswork = lwork
757 *
758  CALL dlacpy( 'F', m, n, asav, lda, usav, lda )
759  srnamt = 'DGESVJ'
760  CALL dgesvj( 'G', 'U', 'V', m, n, usav, lda, ssav,
761  & 0, a, ldvt, work, lwork, info )
762 *
763 * DGESVJ retuns V not VT, so we transpose to use the same
764 * test suite.
765 *
766  DO j=1,n
767  DO i=1,n
768  vtsav(j,i) = a(i,j)
769  END DO
770  END DO
771 *
772  IF( iinfo.NE.0 ) THEN
773  WRITE( nout, fmt = 9995 )'GESVJ', iinfo, m, n,
774  $ jtype, lswork, ioldsd
775  info = abs( iinfo )
776  RETURN
777  END IF
778 *
779 * Do tests 15--18
780 *
781  CALL dbdt01( m, n, 0, asav, lda, usav, ldu, ssav, e,
782  $ vtsav, ldvt, work, result( 15 ) )
783  IF( m.NE.0 .AND. n.NE.0 ) THEN
784  CALL dort01( 'Columns', m, m, usav, ldu, work,
785  $ lwork, result( 16 ) )
786  CALL dort01( 'Rows', n, n, vtsav, ldvt, work,
787  $ lwork, result( 17 ) )
788  END IF
789  result( 18 ) = zero
790  DO 200 i = 1, mnmin - 1
791  IF( ssav( i ).LT.ssav( i+1 ) )
792  $ result( 18 ) = ulpinv
793  IF( ssav( i ).LT.zero )
794  $ result( 18 ) = ulpinv
795  200 CONTINUE
796  IF( mnmin.GE.1 ) THEN
797  IF( ssav( mnmin ).LT.zero )
798  $ result( 18 ) = ulpinv
799  END IF
800  END IF
801 *
802 * Test DGEJSV: Factorize A
803 * Note: DGEJSV does not work for M < N
804 *
805  result( 19 ) = zero
806  result( 20 ) = zero
807  result( 21 ) = zero
808  result( 22 ) = zero
809  IF( m.GE.n ) THEN
810  iwtmp = 5*mnmin*mnmin + 9*mnmin + max( m, n )
811  lswork = iwtmp + ( iws-1 )*( lwork-iwtmp ) / 3
812  lswork = min( lswork, lwork )
813  lswork = max( lswork, 1 )
814  IF( iws.EQ.4 )
815  $ lswork = lwork
816 *
817  CALL dlacpy( 'F', m, n, asav, lda, vtsav, lda )
818  srnamt = 'DGEJSV'
819  CALL dgejsv( 'G', 'U', 'V', 'R', 'N', 'N',
820  & m, n, vtsav, lda, ssav, usav, ldu, a, ldvt,
821  & work, lwork, iwork, info )
822 *
823 * DGEJSV retuns V not VT, so we transpose to use the same
824 * test suite.
825 *
826  DO j=1,n
827  DO i=1,n
828  vtsav(j,i) = a(i,j)
829  END DO
830  END DO
831 *
832  IF( iinfo.NE.0 ) THEN
833  WRITE( nout, fmt = 9995 )'GESVJ', iinfo, m, n,
834  $ jtype, lswork, ioldsd
835  info = abs( iinfo )
836  RETURN
837  END IF
838 *
839 * Do tests 19--22
840 *
841  CALL dbdt01( m, n, 0, asav, lda, usav, ldu, ssav, e,
842  $ vtsav, ldvt, work, result( 19 ) )
843  IF( m.NE.0 .AND. n.NE.0 ) THEN
844  CALL dort01( 'Columns', m, m, usav, ldu, work,
845  $ lwork, result( 20 ) )
846  CALL dort01( 'Rows', n, n, vtsav, ldvt, work,
847  $ lwork, result( 21 ) )
848  END IF
849  result( 22 ) = zero
850  DO 300 i = 1, mnmin - 1
851  IF( ssav( i ).LT.ssav( i+1 ) )
852  $ result( 22 ) = ulpinv
853  IF( ssav( i ).LT.zero )
854  $ result( 22 ) = ulpinv
855  300 CONTINUE
856  IF( mnmin.GE.1 ) THEN
857  IF( ssav( mnmin ).LT.zero )
858  $ result( 22 ) = ulpinv
859  END IF
860  END IF
861 *
862 * End of Loop -- Check for RESULT(j) > THRESH
863 *
864  DO 120 j = 1, 22
865  IF( result( j ).GE.thresh ) THEN
866  IF( nfail.EQ.0 ) THEN
867  WRITE( nout, fmt = 9999 )
868  WRITE( nout, fmt = 9998 )
869  END IF
870  WRITE( nout, fmt = 9997 )m, n, jtype, iws, ioldsd,
871  $ j, result( j )
872  nfail = nfail + 1
873  END IF
874  120 CONTINUE
875  ntest = ntest + 22
876 *
877  130 CONTINUE
878  140 CONTINUE
879  150 CONTINUE
880 *
881 * Summary
882 *
883  CALL alasvm( path, nout, nfail, ntest, 0 )
884 *
885  9999 FORMAT( ' SVD -- Real Singular Value Decomposition Driver ',
886  $ / ' Matrix types (see DDRVBD for details):',
887  $ / / ' 1 = Zero matrix', / ' 2 = Identity matrix',
888  $ / ' 3 = Evenly spaced singular values near 1',
889  $ / ' 4 = Evenly spaced singular values near underflow',
890  $ / ' 5 = Evenly spaced singular values near overflow', / /
891  $ ' Tests performed: ( A is dense, U and V are orthogonal,',
892  $ / 19x, ' S is an array, and Upartial, VTpartial, and',
893  $ / 19x, ' Spartial are partially computed U, VT and S),', / )
894  9998 FORMAT( ' 1 = | A - U diag(S) VT | / ( |A| max(M,N) ulp ) ',
895  $ / ' 2 = | I - U**T U | / ( M ulp ) ',
896  $ / ' 3 = | I - VT VT**T | / ( N ulp ) ',
897  $ / ' 4 = 0 if S contains min(M,N) nonnegative values in',
898  $ ' decreasing order, else 1/ulp',
899  $ / ' 5 = | U - Upartial | / ( M ulp )',
900  $ / ' 6 = | VT - VTpartial | / ( N ulp )',
901  $ / ' 7 = | S - Spartial | / ( min(M,N) ulp |S| )',
902  $ / ' 8 = | A - U diag(S) VT | / ( |A| max(M,N) ulp ) ',
903  $ / ' 9 = | I - U**T U | / ( M ulp ) ',
904  $ / '10 = | I - VT VT**T | / ( N ulp ) ',
905  $ / '11 = 0 if S contains min(M,N) nonnegative values in',
906  $ ' decreasing order, else 1/ulp',
907  $ / '12 = | U - Upartial | / ( M ulp )',
908  $ / '13 = | VT - VTpartial | / ( N ulp )',
909  $ / '14 = | S - Spartial | / ( min(M,N) ulp |S| )',
910  $ / '15 = | A - U diag(S) VT | / ( |A| max(M,N) ulp ) ',
911  $ / '16 = | I - U**T U | / ( M ulp ) ',
912  $ / '17 = | I - VT VT**T | / ( N ulp ) ',
913  $ / '18 = 0 if S contains min(M,N) nonnegative values in',
914  $ ' decreasing order, else 1/ulp',
915  $ / '19 = | U - Upartial | / ( M ulp )',
916  $ / '20 = | VT - VTpartial | / ( N ulp )',
917  $ / '21 = | S - Spartial | / ( min(M,N) ulp |S| )', / / )
918  9997 FORMAT( ' M=', i5, ', N=', i5, ', type ', i1, ', IWS=', i1,
919  $ ', seed=', 4( i4, ',' ), ' test(', i2, ')=', g11.4 )
920  9996 FORMAT( ' DDRVBD: ', a, ' returned INFO=', i6, '.', / 9x, 'M=',
921  $ i6, ', N=', i6, ', JTYPE=', i6, ', ISEED=(', 3( i5, ',' ),
922  $ i5, ')' )
923  9995 FORMAT( ' DDRVBD: ', a, ' returned INFO=', i6, '.', / 9x, 'M=',
924  $ i6, ', N=', i6, ', JTYPE=', i6, ', LSWORK=', i6, / 9x,
925  $ 'ISEED=(', 3( i5, ',' ), i5, ')' )
926 *
927  RETURN
928 *
929 * End of DDRVBD
930 *
subroutine dort03(RC, MU, MV, N, K, U, LDU, V, LDV, WORK, LWORK, RESULT, INFO)
DORT03
Definition: dort03.f:158
subroutine dlatms(M, N, DIST, ISEED, SYM, D, MODE, COND, DMAX, KL, KU, PACK, A, LDA, WORK, INFO)
DLATMS
Definition: dlatms.f:323
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:476
subroutine dgesvj(JOBA, JOBU, JOBV, M, N, A, LDA, SVA, MV, V, LDV, WORK, LWORK, INFO)
DGESVJ
Definition: dgesvj.f:337
subroutine alasvm(TYPE, NOUT, NFAIL, NRUN, NERRS)
ALASVM
Definition: alasvm.f:75
subroutine dbdt01(M, N, KD, A, LDA, Q, LDQ, D, E, PT, LDPT, WORK, RESID)
DBDT01
Definition: dbdt01.f:142
subroutine dgesdd(JOBZ, M, N, A, LDA, S, U, LDU, VT, LDVT, WORK, LWORK, IWORK, INFO)
DGESDD
Definition: dgesdd.f:218
subroutine xerbla(SRNAME, INFO)
XERBLA
Definition: xerbla.f:62
subroutine dlabad(SMALL, LARGE)
DLABAD
Definition: dlabad.f:76
double precision function dlamch(CMACH)
DLAMCH
Definition: dlamch.f:65
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
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
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
subroutine dort01(ROWCOL, M, N, U, LDU, WORK, LWORK, RESID)
DORT01
Definition: dort01.f:118

Here is the call graph for this function:

Here is the caller graph for this function: