LAPACK 3.12.0
LAPACK: Linear Algebra PACKage
Loading...
Searching...
No Matches

◆ zchkbd()

subroutine zchkbd ( integer  nsizes,
integer, dimension( * )  mval,
integer, dimension( * )  nval,
integer  ntypes,
logical, dimension( * )  dotype,
integer  nrhs,
integer, dimension( 4 )  iseed,
double precision  thresh,
complex*16, dimension( lda, * )  a,
integer  lda,
double precision, dimension( * )  bd,
double precision, dimension( * )  be,
double precision, dimension( * )  s1,
double precision, dimension( * )  s2,
complex*16, dimension( ldx, * )  x,
integer  ldx,
complex*16, dimension( ldx, * )  y,
complex*16, dimension( ldx, * )  z,
complex*16, dimension( ldq, * )  q,
integer  ldq,
complex*16, dimension( ldpt, * )  pt,
integer  ldpt,
complex*16, dimension( ldpt, * )  u,
complex*16, dimension( ldpt, * )  vt,
complex*16, dimension( * )  work,
integer  lwork,
double precision, dimension( * )  rwork,
integer  nout,
integer  info 
)

ZCHKBD

Purpose:
 ZCHKBD checks the singular value decomposition (SVD) routines.

 ZGEBRD reduces a complex general m by n matrix A to real upper or
 lower bidiagonal form by an orthogonal transformation: Q' * A * P = B
 (or A = Q * B * P').  The matrix B is upper bidiagonal if m >= n
 and lower bidiagonal if m < n.

 ZUNGBR generates the orthogonal matrices Q and P' from ZGEBRD.
 Note that Q and P are not necessarily square.

 ZBDSQR computes the singular value decomposition of the bidiagonal
 matrix B as B = U S V'.  It is called three times to compute
    1)  B = U S1 V', where S1 is the diagonal matrix of singular
        values and the columns of the matrices U and V are the left
        and right singular vectors, respectively, of B.
    2)  Same as 1), but the singular values are stored in S2 and the
        singular vectors are not computed.
    3)  A = (UQ) S (P'V'), the SVD of the original matrix A.
 In addition, ZBDSQR has an option to apply the left orthogonal matrix
 U to a matrix X, useful in least squares applications.

 For each pair of matrix dimensions (M,N) and each selected matrix
 type, an M by N matrix A and an M by NRHS matrix X are generated.
 The problem dimensions are as follows
    A:          M x N
    Q:          M x min(M,N) (but M x M if NRHS > 0)
    P:          min(M,N) x N
    B:          min(M,N) x min(M,N)
    U, V:       min(M,N) x min(M,N)
    S1, S2      diagonal, order min(M,N)
    X:          M x NRHS

 For each generated matrix, 14 tests are performed:

 Test ZGEBRD and ZUNGBR

 (1)   | A - Q B PT | / ( |A| max(M,N) ulp ), PT = P'

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

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

 Test ZBDSQR on bidiagonal matrix B

 (4)   | B - U S1 VT | / ( |B| min(M,N) ulp ), VT = V'

 (5)   | Y - U Z | / ( |Y| max(min(M,N),k) ulp ), where Y = Q' X
                                                  and   Z = U' Y.
 (6)   | I - U' U | / ( min(M,N) ulp )

 (7)   | I - VT VT' | / ( min(M,N) ulp )

 (8)   S1 contains min(M,N) nonnegative values in decreasing order.
       (Return 0 if true, 1/ULP if false.)

 (9)   0 if the true singular values of B are within THRESH of
       those in S1.  2*THRESH if they are not.  (Tested using
       DSVDCH)

 (10)  | S1 - S2 | / ( |S1| ulp ), where S2 is computed without
                                   computing U and V.

 Test ZBDSQR on matrix A

 (11)  | A - (QU) S (VT PT) | / ( |A| max(M,N) ulp )

 (12)  | X - (QU) Z | / ( |X| max(M,k) ulp )

 (13)  | I - (QU)'(QU) | / ( M ulp )

 (14)  | I - (VT PT) (PT'VT') | / ( N ulp )

 The possible matrix types are

 (1)  The zero matrix.
 (2)  The identity matrix.

 (3)  A diagonal matrix with evenly spaced entries
      1, ..., ULP  and random signs.
      (ULP = (first number larger than 1) - 1 )
 (4)  A diagonal matrix with geometrically spaced entries
      1, ..., ULP  and random signs.
 (5)  A diagonal matrix with "clustered" entries 1, ULP, ..., ULP
      and random signs.

 (6)  Same as (3), but multiplied by SQRT( overflow threshold )
 (7)  Same as (3), but multiplied by SQRT( underflow threshold )

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

 (9)  A matrix of the form  U D V, where U and V are orthogonal and
      D has geometrically spaced entries 1, ..., ULP with random
      signs on the diagonal.

 (10) A matrix of the form  U D V, where U and V are orthogonal and
      D has "clustered" entries 1, ULP,..., ULP with random
      signs on the diagonal.

 (11) Same as (8), but multiplied by SQRT( overflow threshold )
 (12) Same as (8), but multiplied by SQRT( underflow threshold )

 (13) Rectangular matrix with random entries chosen from (-1,1).
 (14) Same as (13), but multiplied by SQRT( overflow threshold )
 (15) Same as (13), but multiplied by SQRT( underflow threshold )

 Special case:
 (16) A bidiagonal matrix with random entries chosen from a
      logarithmic distribution on [ulp^2,ulp^(-2)]  (I.e., each
      entry is  e^x, where x is chosen uniformly on
      [ 2 log(ulp), -2 log(ulp) ] .)  For *this* type:
      (a) ZGEBRD is not called to reduce it to bidiagonal form.
      (b) the bidiagonal is  min(M,N) x min(M,N); if M<N, the
          matrix will be lower bidiagonal, otherwise upper.
      (c) only tests 5--8 and 14 are performed.

 A subset of the full set of matrix types may be selected through
 the logical array DOTYPE.
Parameters
[in]NSIZES
          NSIZES is INTEGER
          The number of values of M and N contained in the vectors
          MVAL and NVAL.  The matrix sizes are used in pairs (M,N).
[in]MVAL
          MVAL is INTEGER array, dimension (NM)
          The values of the matrix row dimension M.
[in]NVAL
          NVAL is INTEGER array, dimension (NM)
          The values of the matrix column dimension N.
[in]NTYPES
          NTYPES is INTEGER
          The number of elements in DOTYPE.   If it is zero, ZCHKBD
          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]NRHS
          NRHS is INTEGER
          The number of columns in the "right-hand side" matrices X, Y,
          and Z, used in testing ZBDSQR.  If NRHS = 0, then the
          operations on the right-hand side will not be tested.
          NRHS must be at least 0.
[in,out]ISEED
          ISEED is INTEGER array, dimension (4)
          On entry ISEED specifies 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.  The values of ISEED are changed on exit, and can be
          used in the next call to ZCHKBD 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.  To have
          every test ratio printed, use THRESH = 0.  Note that the
          expected value of the test ratios is O(1), so THRESH should
          be a reasonably small multiple of 1, e.g., 10 or 100.
[out]A
          A is COMPLEX*16 array, dimension (LDA,NMAX)
          where NMAX is the maximum value of N in NVAL.
[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 MVAL.
[out]BD
          BD is DOUBLE PRECISION array, dimension
                      (max(min(MVAL(j),NVAL(j))))
[out]BE
          BE is DOUBLE PRECISION array, dimension
                      (max(min(MVAL(j),NVAL(j))))
[out]S1
          S1 is DOUBLE PRECISION array, dimension
                      (max(min(MVAL(j),NVAL(j))))
[out]S2
          S2 is DOUBLE PRECISION array, dimension
                      (max(min(MVAL(j),NVAL(j))))
[out]X
          X is COMPLEX*16 array, dimension (LDX,NRHS)
[in]LDX
          LDX is INTEGER
          The leading dimension of the arrays X, Y, and Z.
          LDX >= max(1,MMAX).
[out]Y
          Y is COMPLEX*16 array, dimension (LDX,NRHS)
[out]Z
          Z is COMPLEX*16 array, dimension (LDX,NRHS)
[out]Q
          Q is COMPLEX*16 array, dimension (LDQ,MMAX)
[in]LDQ
          LDQ is INTEGER
          The leading dimension of the array Q.  LDQ >= max(1,MMAX).
[out]PT
          PT is COMPLEX*16 array, dimension (LDPT,NMAX)
[in]LDPT
          LDPT is INTEGER
          The leading dimension of the arrays PT, U, and V.
          LDPT >= max(1, max(min(MVAL(j),NVAL(j)))).
[out]U
          U is COMPLEX*16 array, dimension
                      (LDPT,max(min(MVAL(j),NVAL(j))))
[out]VT
          VT is COMPLEX*16 array, dimension
                      (LDPT,max(min(MVAL(j),NVAL(j))))
[out]WORK
          WORK is COMPLEX*16 array, dimension (LWORK)
[in]LWORK
          LWORK is INTEGER
          The number of entries in WORK.  This must be at least
          3(M+N) and  M(M + max(M,N,k) + 1) + N*min(M,N)  for all
          pairs  (M,N)=(MM(j),NN(j))
[out]RWORK
          RWORK is DOUBLE PRECISION array, dimension
                      (5*max(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
           -6: NRHS  < 0
           -8: THRESH < 0
          -11: LDA < 1 or LDA < MMAX, where MMAX is max( MM(j) ).
          -17: LDB < 1 or LDB < MMAX.
          -21: LDQ < 1 or LDQ < MMAX.
          -23: LDP < 1 or LDP < MNMAX.
          -27: LWORK too small.
          If  ZLATMR, CLATMS, ZGEBRD, ZUNGBR, or ZBDSQR,
              returns an error code, the
              absolute value of it is returned.

-----------------------------------------------------------------------

     Some Local Variables and Parameters:
     ---- ----- --------- --- ----------

     ZERO, ONE       Real 0 and 1.
     MAXTYP          The number of types defined.
     NTEST           The number of tests performed, or which can
                     be performed so far, for the current matrix.
     MMAX            Largest value in NN.
     NMAX            Largest value in NN.
     MNMIN           min(MM(j), NN(j)) (the dimension of the bidiagonal
                     matrix.)
     MNMAX           The maximum value of MNMIN for j=1,...,NSIZES.
     NFAIL           The number of tests which have exceeded THRESH
     COND, IMODE     Values to be passed to the matrix generators.
     ANORM           Norm of A; passed to matrix generators.

     OVFL, UNFL      Overflow and underflow thresholds.
     RTOVFL, RTUNFL  Square roots of the previous 2 values.
     ULP, ULPINV     Finest relative precision and its inverse.

             The following four arrays decode JTYPE:
     KTYPE(j)        The general type (1-10) for type "j".
     KMODE(j)        The MODE value to be passed to the matrix
                     generator for type "j".
     KMAGN(j)        The order of magnitude ( O(1),
                     O(overflow^(1/2) ), O(underflow^(1/2) )
Author
Univ. of Tennessee
Univ. of California Berkeley
Univ. of Colorado Denver
NAG Ltd.

Definition at line 411 of file zchkbd.f.

415*
416* -- LAPACK test routine --
417* -- LAPACK is a software package provided by Univ. of Tennessee, --
418* -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
419*
420* .. Scalar Arguments ..
421 INTEGER INFO, LDA, LDPT, LDQ, LDX, LWORK, NOUT, NRHS,
422 $ NSIZES, NTYPES
423 DOUBLE PRECISION THRESH
424* ..
425* .. Array Arguments ..
426 LOGICAL DOTYPE( * )
427 INTEGER ISEED( 4 ), MVAL( * ), NVAL( * )
428 DOUBLE PRECISION BD( * ), BE( * ), RWORK( * ), S1( * ), S2( * )
429 COMPLEX*16 A( LDA, * ), PT( LDPT, * ), Q( LDQ, * ),
430 $ U( LDPT, * ), VT( LDPT, * ), WORK( * ),
431 $ X( LDX, * ), Y( LDX, * ), Z( LDX, * )
432* ..
433*
434* ======================================================================
435*
436* .. Parameters ..
437 DOUBLE PRECISION ZERO, ONE, TWO, HALF
438 parameter( zero = 0.0d0, one = 1.0d0, two = 2.0d0,
439 $ half = 0.5d0 )
440 COMPLEX*16 CZERO, CONE
441 parameter( czero = ( 0.0d+0, 0.0d+0 ),
442 $ cone = ( 1.0d+0, 0.0d+0 ) )
443 INTEGER MAXTYP
444 parameter( maxtyp = 16 )
445* ..
446* .. Local Scalars ..
447 LOGICAL BADMM, BADNN, BIDIAG
448 CHARACTER UPLO
449 CHARACTER*3 PATH
450 INTEGER I, IINFO, IMODE, ITYPE, J, JCOL, JSIZE, JTYPE,
451 $ LOG2UI, M, MINWRK, MMAX, MNMAX, MNMIN, MQ,
452 $ MTYPES, N, NFAIL, NMAX, NTEST
453 DOUBLE PRECISION AMNINV, ANORM, COND, OVFL, RTOVFL, RTUNFL,
454 $ TEMP1, TEMP2, ULP, ULPINV, UNFL
455* ..
456* .. Local Arrays ..
457 INTEGER IOLDSD( 4 ), IWORK( 1 ), KMAGN( MAXTYP ),
458 $ KMODE( MAXTYP ), KTYPE( MAXTYP )
459 DOUBLE PRECISION DUMMA( 1 ), RESULT( 14 )
460* ..
461* .. External Functions ..
462 DOUBLE PRECISION DLAMCH, DLARND
463 EXTERNAL dlamch, dlarnd
464* ..
465* .. External Subroutines ..
466 EXTERNAL alasum, dcopy, dlahd2, dsvdch, xerbla, zbdsqr,
469* ..
470* .. Intrinsic Functions ..
471 INTRINSIC abs, exp, int, log, max, min, sqrt
472* ..
473* .. Scalars in Common ..
474 LOGICAL LERR, OK
475 CHARACTER*32 SRNAMT
476 INTEGER INFOT, NUNIT
477* ..
478* .. Common blocks ..
479 COMMON / infoc / infot, nunit, ok, lerr
480 COMMON / srnamc / srnamt
481* ..
482* .. Data statements ..
483 DATA ktype / 1, 2, 5*4, 5*6, 3*9, 10 /
484 DATA kmagn / 2*1, 3*1, 2, 3, 3*1, 2, 3, 1, 2, 3, 0 /
485 DATA kmode / 2*0, 4, 3, 1, 4, 4, 4, 3, 1, 4, 4, 0,
486 $ 0, 0, 0 /
487* ..
488* .. Executable Statements ..
489*
490* Check for errors
491*
492 info = 0
493*
494 badmm = .false.
495 badnn = .false.
496 mmax = 1
497 nmax = 1
498 mnmax = 1
499 minwrk = 1
500 DO 10 j = 1, nsizes
501 mmax = max( mmax, mval( j ) )
502 IF( mval( j ).LT.0 )
503 $ badmm = .true.
504 nmax = max( nmax, nval( j ) )
505 IF( nval( j ).LT.0 )
506 $ badnn = .true.
507 mnmax = max( mnmax, min( mval( j ), nval( j ) ) )
508 minwrk = max( minwrk, 3*( mval( j )+nval( j ) ),
509 $ mval( j )*( mval( j )+max( mval( j ), nval( j ),
510 $ nrhs )+1 )+nval( j )*min( nval( j ), mval( j ) ) )
511 10 CONTINUE
512*
513* Check for errors
514*
515 IF( nsizes.LT.0 ) THEN
516 info = -1
517 ELSE IF( badmm ) THEN
518 info = -2
519 ELSE IF( badnn ) THEN
520 info = -3
521 ELSE IF( ntypes.LT.0 ) THEN
522 info = -4
523 ELSE IF( nrhs.LT.0 ) THEN
524 info = -6
525 ELSE IF( lda.LT.mmax ) THEN
526 info = -11
527 ELSE IF( ldx.LT.mmax ) THEN
528 info = -17
529 ELSE IF( ldq.LT.mmax ) THEN
530 info = -21
531 ELSE IF( ldpt.LT.mnmax ) THEN
532 info = -23
533 ELSE IF( minwrk.GT.lwork ) THEN
534 info = -27
535 END IF
536*
537 IF( info.NE.0 ) THEN
538 CALL xerbla( 'ZCHKBD', -info )
539 RETURN
540 END IF
541*
542* Initialize constants
543*
544 path( 1: 1 ) = 'Zomplex precision'
545 path( 2: 3 ) = 'BD'
546 nfail = 0
547 ntest = 0
548 unfl = dlamch( 'Safe minimum' )
549 ovfl = dlamch( 'Overflow' )
550 ulp = dlamch( 'Precision' )
551 ulpinv = one / ulp
552 log2ui = int( log( ulpinv ) / log( two ) )
553 rtunfl = sqrt( unfl )
554 rtovfl = sqrt( ovfl )
555 infot = 0
556*
557* Loop over sizes, types
558*
559 DO 180 jsize = 1, nsizes
560 m = mval( jsize )
561 n = nval( jsize )
562 mnmin = min( m, n )
563 amninv = one / max( m, n, 1 )
564*
565 IF( nsizes.NE.1 ) THEN
566 mtypes = min( maxtyp, ntypes )
567 ELSE
568 mtypes = min( maxtyp+1, ntypes )
569 END IF
570*
571 DO 170 jtype = 1, mtypes
572 IF( .NOT.dotype( jtype ) )
573 $ GO TO 170
574*
575 DO 20 j = 1, 4
576 ioldsd( j ) = iseed( j )
577 20 CONTINUE
578*
579 DO 30 j = 1, 14
580 result( j ) = -one
581 30 CONTINUE
582*
583 uplo = ' '
584*
585* Compute "A"
586*
587* Control parameters:
588*
589* KMAGN KMODE KTYPE
590* =1 O(1) clustered 1 zero
591* =2 large clustered 2 identity
592* =3 small exponential (none)
593* =4 arithmetic diagonal, (w/ eigenvalues)
594* =5 random symmetric, w/ eigenvalues
595* =6 nonsymmetric, w/ singular values
596* =7 random diagonal
597* =8 random symmetric
598* =9 random nonsymmetric
599* =10 random bidiagonal (log. distrib.)
600*
601 IF( mtypes.GT.maxtyp )
602 $ GO TO 100
603*
604 itype = ktype( jtype )
605 imode = kmode( jtype )
606*
607* Compute norm
608*
609 GO TO ( 40, 50, 60 )kmagn( jtype )
610*
611 40 CONTINUE
612 anorm = one
613 GO TO 70
614*
615 50 CONTINUE
616 anorm = ( rtovfl*ulp )*amninv
617 GO TO 70
618*
619 60 CONTINUE
620 anorm = rtunfl*max( m, n )*ulpinv
621 GO TO 70
622*
623 70 CONTINUE
624*
625 CALL zlaset( 'Full', lda, n, czero, czero, a, lda )
626 iinfo = 0
627 cond = ulpinv
628*
629 bidiag = .false.
630 IF( itype.EQ.1 ) THEN
631*
632* Zero matrix
633*
634 iinfo = 0
635*
636 ELSE IF( itype.EQ.2 ) THEN
637*
638* Identity
639*
640 DO 80 jcol = 1, mnmin
641 a( jcol, jcol ) = anorm
642 80 CONTINUE
643*
644 ELSE IF( itype.EQ.4 ) THEN
645*
646* Diagonal Matrix, [Eigen]values Specified
647*
648 CALL zlatms( mnmin, mnmin, 'S', iseed, 'N', rwork, imode,
649 $ cond, anorm, 0, 0, 'N', a, lda, work,
650 $ iinfo )
651*
652 ELSE IF( itype.EQ.5 ) THEN
653*
654* Symmetric, eigenvalues specified
655*
656 CALL zlatms( mnmin, mnmin, 'S', iseed, 'S', rwork, imode,
657 $ cond, anorm, m, n, 'N', a, lda, work,
658 $ iinfo )
659*
660 ELSE IF( itype.EQ.6 ) THEN
661*
662* Nonsymmetric, singular values specified
663*
664 CALL zlatms( m, n, 'S', iseed, 'N', rwork, imode, cond,
665 $ anorm, m, n, 'N', a, lda, work, iinfo )
666*
667 ELSE IF( itype.EQ.7 ) THEN
668*
669* Diagonal, random entries
670*
671 CALL zlatmr( mnmin, mnmin, 'S', iseed, 'N', work, 6, one,
672 $ cone, 'T', 'N', work( mnmin+1 ), 1, one,
673 $ work( 2*mnmin+1 ), 1, one, 'N', iwork, 0, 0,
674 $ zero, anorm, 'NO', a, lda, iwork, iinfo )
675*
676 ELSE IF( itype.EQ.8 ) THEN
677*
678* Symmetric, random entries
679*
680 CALL zlatmr( mnmin, mnmin, 'S', iseed, 'S', work, 6, one,
681 $ cone, 'T', 'N', work( mnmin+1 ), 1, one,
682 $ work( m+mnmin+1 ), 1, one, 'N', iwork, m, n,
683 $ zero, anorm, 'NO', a, lda, iwork, iinfo )
684*
685 ELSE IF( itype.EQ.9 ) THEN
686*
687* Nonsymmetric, random entries
688*
689 CALL zlatmr( m, n, 'S', iseed, 'N', work, 6, one, cone,
690 $ 'T', 'N', work( mnmin+1 ), 1, one,
691 $ work( m+mnmin+1 ), 1, one, 'N', iwork, m, n,
692 $ zero, anorm, 'NO', a, lda, iwork, iinfo )
693*
694 ELSE IF( itype.EQ.10 ) THEN
695*
696* Bidiagonal, random entries
697*
698 temp1 = -two*log( ulp )
699 DO 90 j = 1, mnmin
700 bd( j ) = exp( temp1*dlarnd( 2, iseed ) )
701 IF( j.LT.mnmin )
702 $ be( j ) = exp( temp1*dlarnd( 2, iseed ) )
703 90 CONTINUE
704*
705 iinfo = 0
706 bidiag = .true.
707 IF( m.GE.n ) THEN
708 uplo = 'U'
709 ELSE
710 uplo = 'L'
711 END IF
712 ELSE
713 iinfo = 1
714 END IF
715*
716 IF( iinfo.EQ.0 ) THEN
717*
718* Generate Right-Hand Side
719*
720 IF( bidiag ) THEN
721 CALL zlatmr( mnmin, nrhs, 'S', iseed, 'N', work, 6,
722 $ one, cone, 'T', 'N', work( mnmin+1 ), 1,
723 $ one, work( 2*mnmin+1 ), 1, one, 'N',
724 $ iwork, mnmin, nrhs, zero, one, 'NO', y,
725 $ ldx, iwork, iinfo )
726 ELSE
727 CALL zlatmr( m, nrhs, 'S', iseed, 'N', work, 6, one,
728 $ cone, 'T', 'N', work( m+1 ), 1, one,
729 $ work( 2*m+1 ), 1, one, 'N', iwork, m,
730 $ nrhs, zero, one, 'NO', x, ldx, iwork,
731 $ iinfo )
732 END IF
733 END IF
734*
735* Error Exit
736*
737 IF( iinfo.NE.0 ) THEN
738 WRITE( nout, fmt = 9998 )'Generator', iinfo, m, n,
739 $ jtype, ioldsd
740 info = abs( iinfo )
741 RETURN
742 END IF
743*
744 100 CONTINUE
745*
746* Call ZGEBRD and ZUNGBR to compute B, Q, and P, do tests.
747*
748 IF( .NOT.bidiag ) THEN
749*
750* Compute transformations to reduce A to bidiagonal form:
751* B := Q' * A * P.
752*
753 CALL zlacpy( ' ', m, n, a, lda, q, ldq )
754 CALL zgebrd( m, n, q, ldq, bd, be, work, work( mnmin+1 ),
755 $ work( 2*mnmin+1 ), lwork-2*mnmin, iinfo )
756*
757* Check error code from ZGEBRD.
758*
759 IF( iinfo.NE.0 ) THEN
760 WRITE( nout, fmt = 9998 )'ZGEBRD', iinfo, m, n,
761 $ jtype, ioldsd
762 info = abs( iinfo )
763 RETURN
764 END IF
765*
766 CALL zlacpy( ' ', m, n, q, ldq, pt, ldpt )
767 IF( m.GE.n ) THEN
768 uplo = 'U'
769 ELSE
770 uplo = 'L'
771 END IF
772*
773* Generate Q
774*
775 mq = m
776 IF( nrhs.LE.0 )
777 $ mq = mnmin
778 CALL zungbr( 'Q', m, mq, n, q, ldq, work,
779 $ work( 2*mnmin+1 ), lwork-2*mnmin, iinfo )
780*
781* Check error code from ZUNGBR.
782*
783 IF( iinfo.NE.0 ) THEN
784 WRITE( nout, fmt = 9998 )'ZUNGBR(Q)', iinfo, m, n,
785 $ jtype, ioldsd
786 info = abs( iinfo )
787 RETURN
788 END IF
789*
790* Generate P'
791*
792 CALL zungbr( 'P', mnmin, n, m, pt, ldpt, work( mnmin+1 ),
793 $ work( 2*mnmin+1 ), lwork-2*mnmin, iinfo )
794*
795* Check error code from ZUNGBR.
796*
797 IF( iinfo.NE.0 ) THEN
798 WRITE( nout, fmt = 9998 )'ZUNGBR(P)', iinfo, m, n,
799 $ jtype, ioldsd
800 info = abs( iinfo )
801 RETURN
802 END IF
803*
804* Apply Q' to an M by NRHS matrix X: Y := Q' * X.
805*
806 CALL zgemm( 'Conjugate transpose', 'No transpose', m,
807 $ nrhs, m, cone, q, ldq, x, ldx, czero, y,
808 $ ldx )
809*
810* Test 1: Check the decomposition A := Q * B * PT
811* 2: Check the orthogonality of Q
812* 3: Check the orthogonality of PT
813*
814 CALL zbdt01( m, n, 1, a, lda, q, ldq, bd, be, pt, ldpt,
815 $ work, rwork, result( 1 ) )
816 CALL zunt01( 'Columns', m, mq, q, ldq, work, lwork,
817 $ rwork, result( 2 ) )
818 CALL zunt01( 'Rows', mnmin, n, pt, ldpt, work, lwork,
819 $ rwork, result( 3 ) )
820 END IF
821*
822* Use ZBDSQR to form the SVD of the bidiagonal matrix B:
823* B := U * S1 * VT, and compute Z = U' * Y.
824*
825 CALL dcopy( mnmin, bd, 1, s1, 1 )
826 IF( mnmin.GT.0 )
827 $ CALL dcopy( mnmin-1, be, 1, rwork, 1 )
828 CALL zlacpy( ' ', m, nrhs, y, ldx, z, ldx )
829 CALL zlaset( 'Full', mnmin, mnmin, czero, cone, u, ldpt )
830 CALL zlaset( 'Full', mnmin, mnmin, czero, cone, vt, ldpt )
831*
832 CALL zbdsqr( uplo, mnmin, mnmin, mnmin, nrhs, s1, rwork, vt,
833 $ ldpt, u, ldpt, z, ldx, rwork( mnmin+1 ),
834 $ iinfo )
835*
836* Check error code from ZBDSQR.
837*
838 IF( iinfo.NE.0 ) THEN
839 WRITE( nout, fmt = 9998 )'ZBDSQR(vects)', iinfo, m, n,
840 $ jtype, ioldsd
841 info = abs( iinfo )
842 IF( iinfo.LT.0 ) THEN
843 RETURN
844 ELSE
845 result( 4 ) = ulpinv
846 GO TO 150
847 END IF
848 END IF
849*
850* Use ZBDSQR to compute only the singular values of the
851* bidiagonal matrix B; U, VT, and Z should not be modified.
852*
853 CALL dcopy( mnmin, bd, 1, s2, 1 )
854 IF( mnmin.GT.0 )
855 $ CALL dcopy( mnmin-1, be, 1, rwork, 1 )
856*
857 CALL zbdsqr( uplo, mnmin, 0, 0, 0, s2, rwork, vt, ldpt, u,
858 $ ldpt, z, ldx, rwork( mnmin+1 ), iinfo )
859*
860* Check error code from ZBDSQR.
861*
862 IF( iinfo.NE.0 ) THEN
863 WRITE( nout, fmt = 9998 )'ZBDSQR(values)', iinfo, m, n,
864 $ jtype, ioldsd
865 info = abs( iinfo )
866 IF( iinfo.LT.0 ) THEN
867 RETURN
868 ELSE
869 result( 9 ) = ulpinv
870 GO TO 150
871 END IF
872 END IF
873*
874* Test 4: Check the decomposition B := U * S1 * VT
875* 5: Check the computation Z := U' * Y
876* 6: Check the orthogonality of U
877* 7: Check the orthogonality of VT
878*
879 CALL zbdt03( uplo, mnmin, 1, bd, be, u, ldpt, s1, vt, ldpt,
880 $ work, result( 4 ) )
881 CALL zbdt02( mnmin, nrhs, y, ldx, z, ldx, u, ldpt, work,
882 $ rwork, result( 5 ) )
883 CALL zunt01( 'Columns', mnmin, mnmin, u, ldpt, work, lwork,
884 $ rwork, result( 6 ) )
885 CALL zunt01( 'Rows', mnmin, mnmin, vt, ldpt, work, lwork,
886 $ rwork, result( 7 ) )
887*
888* Test 8: Check that the singular values are sorted in
889* non-increasing order and are non-negative
890*
891 result( 8 ) = zero
892 DO 110 i = 1, mnmin - 1
893 IF( s1( i ).LT.s1( i+1 ) )
894 $ result( 8 ) = ulpinv
895 IF( s1( i ).LT.zero )
896 $ result( 8 ) = ulpinv
897 110 CONTINUE
898 IF( mnmin.GE.1 ) THEN
899 IF( s1( mnmin ).LT.zero )
900 $ result( 8 ) = ulpinv
901 END IF
902*
903* Test 9: Compare ZBDSQR with and without singular vectors
904*
905 temp2 = zero
906*
907 DO 120 j = 1, mnmin
908 temp1 = abs( s1( j )-s2( j ) ) /
909 $ max( sqrt( unfl )*max( s1( 1 ), one ),
910 $ ulp*max( abs( s1( j ) ), abs( s2( j ) ) ) )
911 temp2 = max( temp1, temp2 )
912 120 CONTINUE
913*
914 result( 9 ) = temp2
915*
916* Test 10: Sturm sequence test of singular values
917* Go up by factors of two until it succeeds
918*
919 temp1 = thresh*( half-ulp )
920*
921 DO 130 j = 0, log2ui
922 CALL dsvdch( mnmin, bd, be, s1, temp1, iinfo )
923 IF( iinfo.EQ.0 )
924 $ GO TO 140
925 temp1 = temp1*two
926 130 CONTINUE
927*
928 140 CONTINUE
929 result( 10 ) = temp1
930*
931* Use ZBDSQR to form the decomposition A := (QU) S (VT PT)
932* from the bidiagonal form A := Q B PT.
933*
934 IF( .NOT.bidiag ) THEN
935 CALL dcopy( mnmin, bd, 1, s2, 1 )
936 IF( mnmin.GT.0 )
937 $ CALL dcopy( mnmin-1, be, 1, rwork, 1 )
938*
939 CALL zbdsqr( uplo, mnmin, n, m, nrhs, s2, rwork, pt,
940 $ ldpt, q, ldq, y, ldx, rwork( mnmin+1 ),
941 $ iinfo )
942*
943* Test 11: Check the decomposition A := Q*U * S2 * VT*PT
944* 12: Check the computation Z := U' * Q' * X
945* 13: Check the orthogonality of Q*U
946* 14: Check the orthogonality of VT*PT
947*
948 CALL zbdt01( m, n, 0, a, lda, q, ldq, s2, dumma, pt,
949 $ ldpt, work, rwork, result( 11 ) )
950 CALL zbdt02( m, nrhs, x, ldx, y, ldx, q, ldq, work,
951 $ rwork, result( 12 ) )
952 CALL zunt01( 'Columns', m, mq, q, ldq, work, lwork,
953 $ rwork, result( 13 ) )
954 CALL zunt01( 'Rows', mnmin, n, pt, ldpt, work, lwork,
955 $ rwork, result( 14 ) )
956 END IF
957*
958* End of Loop -- Check for RESULT(j) > THRESH
959*
960 150 CONTINUE
961 DO 160 j = 1, 14
962 IF( result( j ).GE.thresh ) THEN
963 IF( nfail.EQ.0 )
964 $ CALL dlahd2( nout, path )
965 WRITE( nout, fmt = 9999 )m, n, jtype, ioldsd, j,
966 $ result( j )
967 nfail = nfail + 1
968 END IF
969 160 CONTINUE
970 IF( .NOT.bidiag ) THEN
971 ntest = ntest + 14
972 ELSE
973 ntest = ntest + 5
974 END IF
975*
976 170 CONTINUE
977 180 CONTINUE
978*
979* Summary
980*
981 CALL alasum( path, nout, nfail, ntest, 0 )
982*
983 RETURN
984*
985* End of ZCHKBD
986*
987 9999 FORMAT( ' M=', i5, ', N=', i5, ', type ', i2, ', seed=',
988 $ 4( i4, ',' ), ' test(', i2, ')=', g11.4 )
989 9998 FORMAT( ' ZCHKBD: ', a, ' returned INFO=', i6, '.', / 9x, 'M=',
990 $ i6, ', N=', i6, ', JTYPE=', i6, ', ISEED=(', 3( i5, ',' ),
991 $ i5, ')' )
992*
subroutine alasum(type, nout, nfail, nrun, nerrs)
ALASUM
Definition alasum.f:73
subroutine xerbla(srname, info)
Definition cblat2.f:3285
subroutine dlahd2(iounit, path)
DLAHD2
Definition dlahd2.f:65
double precision function dlarnd(idist, iseed)
DLARND
Definition dlarnd.f:73
subroutine dsvdch(n, s, e, svd, tol, info)
DSVDCH
Definition dsvdch.f:97
subroutine zbdsqr(uplo, n, ncvt, nru, ncc, d, e, vt, ldvt, u, ldu, c, ldc, rwork, info)
ZBDSQR
Definition zbdsqr.f:233
subroutine dcopy(n, dx, incx, dy, incy)
DCOPY
Definition dcopy.f:82
subroutine zgebrd(m, n, a, lda, d, e, tauq, taup, work, lwork, info)
ZGEBRD
Definition zgebrd.f:205
subroutine zgemm(transa, transb, m, n, k, alpha, a, lda, b, ldb, beta, c, ldc)
ZGEMM
Definition zgemm.f:188
subroutine zlacpy(uplo, m, n, a, lda, b, ldb)
ZLACPY copies all or part of one two-dimensional array to another.
Definition zlacpy.f:103
double precision function dlamch(cmach)
DLAMCH
Definition dlamch.f:69
subroutine zlaset(uplo, m, n, alpha, beta, a, lda)
ZLASET initializes the off-diagonal elements and the diagonal elements of a matrix to given values.
Definition zlaset.f:106
subroutine zungbr(vect, m, n, k, a, lda, tau, work, lwork, info)
ZUNGBR
Definition zungbr.f:157
subroutine zbdt01(m, n, kd, a, lda, q, ldq, d, e, pt, ldpt, work, rwork, resid)
ZBDT01
Definition zbdt01.f:147
subroutine zbdt02(m, n, b, ldb, c, ldc, u, ldu, work, rwork, resid)
ZBDT02
Definition zbdt02.f:120
subroutine zbdt03(uplo, n, kd, d, e, u, ldu, s, vt, ldvt, work, resid)
ZBDT03
Definition zbdt03.f:135
subroutine zlatmr(m, n, dist, iseed, sym, d, mode, cond, dmax, rsign, grade, dl, model, condl, dr, moder, condr, pivtng, ipivot, kl, ku, sparse, anorm, pack, a, lda, iwork, info)
ZLATMR
Definition zlatmr.f:490
subroutine zlatms(m, n, dist, iseed, sym, d, mode, cond, dmax, kl, ku, pack, a, lda, work, info)
ZLATMS
Definition zlatms.f:332
subroutine zunt01(rowcol, m, n, u, ldu, work, lwork, rwork, resid)
ZUNT01
Definition zunt01.f:126
Here is the call graph for this function:
Here is the caller graph for this function: