LAPACK  3.6.1
LAPACK: Linear Algebra PACKage
subroutine zlatmt ( integer  M,
integer  N,
character  DIST,
integer, dimension( 4 )  ISEED,
character  SYM,
double precision, dimension( * )  D,
integer  MODE,
double precision  COND,
double precision  DMAX,
integer  RANK,
integer  KL,
integer  KU,
character  PACK,
complex*16, dimension( lda, * )  A,
integer  LDA,
complex*16, dimension( * )  WORK,
integer  INFO 
)

ZLATMT

Purpose:
    ZLATMT generates random matrices with specified singular values
    (or hermitian with specified eigenvalues)
    for testing LAPACK programs.

    ZLATMT operates by applying the following sequence of
    operations:

      Set the diagonal to D, where D may be input or
         computed according to MODE, COND, DMAX, and SYM
         as described below.

      Generate a matrix with the appropriate band structure, by one
         of two methods:

      Method A:
          Generate a dense M x N matrix by multiplying D on the left
              and the right by random unitary matrices, then:

          Reduce the bandwidth according to KL and KU, using
              Householder transformations.

      Method B:
          Convert the bandwidth-0 (i.e., diagonal) matrix to a
              bandwidth-1 matrix using Givens rotations, "chasing"
              out-of-band elements back, much as in QR; then convert
              the bandwidth-1 to a bandwidth-2 matrix, etc.  Note
              that for reasonably small bandwidths (relative to M and
              N) this requires less storage, as a dense matrix is not
              generated.  Also, for hermitian or symmetric matrices,
              only one triangle is generated.

      Method A is chosen if the bandwidth is a large fraction of the
          order of the matrix, and LDA is at least M (so a dense
          matrix can be stored.)  Method B is chosen if the bandwidth
          is small (< 1/2 N for hermitian or symmetric, < .3 N+M for
          non-symmetric), or LDA is less than M and not less than the
          bandwidth.

      Pack the matrix if desired. Options specified by PACK are:
         no packing
         zero out upper half (if hermitian)
         zero out lower half (if hermitian)
         store the upper half columnwise (if hermitian or upper
               triangular)
         store the lower half columnwise (if hermitian or lower
               triangular)
         store the lower triangle in banded format (if hermitian or
               lower triangular)
         store the upper triangle in banded format (if hermitian or
               upper triangular)
         store the entire matrix in banded format
      If Method B is chosen, and band format is specified, then the
         matrix will be generated in the band format, so no repacking
         will be necessary.
Parameters
[in]M
          M is INTEGER
           The number of rows of A. Not modified.
[in]N
          N is INTEGER
           The number of columns of A. N must equal M if the matrix
           is symmetric or hermitian (i.e., if SYM is not 'N')
           Not modified.
[in]DIST
          DIST is CHARACTER*1
           On entry, DIST specifies the type of distribution to be used
           to generate the random eigen-/singular values.
           'U' => UNIFORM( 0, 1 )  ( 'U' for uniform )
           'S' => UNIFORM( -1, 1 ) ( 'S' for symmetric )
           'N' => NORMAL( 0, 1 )   ( 'N' for normal )
           Not modified.
[in,out]ISEED
          ISEED is INTEGER array, dimension ( 4 )
           On entry ISEED specifies the seed of the random number
           generator. They should lie between 0 and 4095 inclusive,
           and ISEED(4) should be odd. The random number generator
           uses a linear congruential sequence limited to small
           integers, and so should produce machine independent
           random numbers. The values of ISEED are changed on
           exit, and can be used in the next call to ZLATMT
           to continue the same random number sequence.
           Changed on exit.
[in]SYM
          SYM is CHARACTER*1
           If SYM='H', the generated matrix is hermitian, with
             eigenvalues specified by D, COND, MODE, and DMAX; they
             may be positive, negative, or zero.
           If SYM='P', the generated matrix is hermitian, with
             eigenvalues (= singular values) specified by D, COND,
             MODE, and DMAX; they will not be negative.
           If SYM='N', the generated matrix is nonsymmetric, with
             singular values specified by D, COND, MODE, and DMAX;
             they will not be negative.
           If SYM='S', the generated matrix is (complex) symmetric,
             with singular values specified by D, COND, MODE, and
             DMAX; they will not be negative.
           Not modified.
[in,out]D
          D is DOUBLE PRECISION array, dimension ( MIN( M, N ) )
           This array is used to specify the singular values or
           eigenvalues of A (see SYM, above.)  If MODE=0, then D is
           assumed to contain the singular/eigenvalues, otherwise
           they will be computed according to MODE, COND, and DMAX,
           and placed in D.
           Modified if MODE is nonzero.
[in]MODE
          MODE is INTEGER
           On entry this describes how the singular/eigenvalues are to
           be specified:
           MODE = 0 means use D as input
           MODE = 1 sets D(1)=1 and D(2:RANK)=1.0/COND
           MODE = 2 sets D(1:RANK-1)=1 and D(RANK)=1.0/COND
           MODE = 3 sets D(I)=COND**(-(I-1)/(RANK-1))
           MODE = 4 sets D(i)=1 - (i-1)/(N-1)*(1 - 1/COND)
           MODE = 5 sets D to random numbers in the range
                    ( 1/COND , 1 ) such that their logarithms
                    are uniformly distributed.
           MODE = 6 set D to random numbers from same distribution
                    as the rest of the matrix.
           MODE < 0 has the same meaning as ABS(MODE), except that
              the order of the elements of D is reversed.
           Thus if MODE is positive, D has entries ranging from
              1 to 1/COND, if negative, from 1/COND to 1,
           If SYM='H', and MODE is neither 0, 6, nor -6, then
              the elements of D will also be multiplied by a random
              sign (i.e., +1 or -1.)
           Not modified.
[in]COND
          COND is DOUBLE PRECISION
           On entry, this is used as described under MODE above.
           If used, it must be >= 1. Not modified.
[in]DMAX
          DMAX is DOUBLE PRECISION
           If MODE is neither -6, 0 nor 6, the contents of D, as
           computed according to MODE and COND, will be scaled by
           DMAX / max(abs(D(i))); thus, the maximum absolute eigen- or
           singular value (which is to say the norm) will be abs(DMAX).
           Note that DMAX need not be positive: if DMAX is negative
           (or zero), D will be scaled by a negative number (or zero).
           Not modified.
[in]RANK
          RANK is INTEGER
           The rank of matrix to be generated for modes 1,2,3 only.
           D( RANK+1:N ) = 0.
           Not modified.
[in]KL
          KL is INTEGER
           This specifies the lower bandwidth of the  matrix. For
           example, KL=0 implies upper triangular, KL=1 implies upper
           Hessenberg, and KL being at least M-1 means that the matrix
           has full lower bandwidth.  KL must equal KU if the matrix
           is symmetric or hermitian.
           Not modified.
[in]KU
          KU is INTEGER
           This specifies the upper bandwidth of the  matrix. For
           example, KU=0 implies lower triangular, KU=1 implies lower
           Hessenberg, and KU being at least N-1 means that the matrix
           has full upper bandwidth.  KL must equal KU if the matrix
           is symmetric or hermitian.
           Not modified.
[in]PACK
          PACK is CHARACTER*1
           This specifies packing of matrix as follows:
           'N' => no packing
           'U' => zero out all subdiagonal entries (if symmetric
                  or hermitian)
           'L' => zero out all superdiagonal entries (if symmetric
                  or hermitian)
           'C' => store the upper triangle columnwise (only if the
                  matrix is symmetric, hermitian, or upper triangular)
           'R' => store the lower triangle columnwise (only if the
                  matrix is symmetric, hermitian, or lower triangular)
           'B' => store the lower triangle in band storage scheme
                  (only if the matrix is symmetric, hermitian, or
                  lower triangular)
           'Q' => store the upper triangle in band storage scheme
                  (only if the matrix is symmetric, hermitian, or
                  upper triangular)
           'Z' => store the entire matrix in band storage scheme
                      (pivoting can be provided for by using this
                      option to store A in the trailing rows of
                      the allocated storage)

           Using these options, the various LAPACK packed and banded
           storage schemes can be obtained:
           GB                    - use 'Z'
           PB, SB, HB, or TB     - use 'B' or 'Q'
           PP, SP, HB, or TP     - use 'C' or 'R'

           If two calls to ZLATMT differ only in the PACK parameter,
           they will generate mathematically equivalent matrices.
           Not modified.
[in,out]A
          A is COMPLEX*16 array, dimension ( LDA, N )
           On exit A is the desired test matrix.  A is first generated
           in full (unpacked) form, and then packed, if so specified
           by PACK.  Thus, the first M elements of the first N
           columns will always be modified.  If PACK specifies a
           packed or banded storage scheme, all LDA elements of the
           first N columns will be modified; the elements of the
           array which do not correspond to elements of the generated
           matrix are set to zero.
           Modified.
[in]LDA
          LDA is INTEGER
           LDA specifies the first dimension of A as declared in the
           calling program.  If PACK='N', 'U', 'L', 'C', or 'R', then
           LDA must be at least M.  If PACK='B' or 'Q', then LDA must
           be at least MIN( KL, M-1) (which is equal to MIN(KU,N-1)).
           If PACK='Z', LDA must be large enough to hold the packed
           array: MIN( KU, N-1) + MIN( KL, M-1) + 1.
           Not modified.
[out]WORK
          WORK is COMPLEX*16 array, dimension ( 3*MAX( N, M ) )
           Workspace.
           Modified.
[out]INFO
          INFO is INTEGER
           Error code.  On exit, INFO will be set to one of the
           following values:
             0 => normal return
            -1 => M negative or unequal to N and SYM='S', 'H', or 'P'
            -2 => N negative
            -3 => DIST illegal string
            -5 => SYM illegal string
            -7 => MODE not in range -6 to 6
            -8 => COND less than 1.0, and MODE neither -6, 0 nor 6
           -10 => KL negative
           -11 => KU negative, or SYM is not 'N' and KU is not equal to
                  KL
           -12 => PACK illegal string, or PACK='U' or 'L', and SYM='N';
                  or PACK='C' or 'Q' and SYM='N' and KL is not zero;
                  or PACK='R' or 'B' and SYM='N' and KU is not zero;
                  or PACK='U', 'L', 'C', 'R', 'B', or 'Q', and M is not
                  N.
           -14 => LDA is less than M, or PACK='Z' and LDA is less than
                  MIN(KU,N-1) + MIN(KL,M-1) + 1.
            1  => Error return from DLATM7
            2  => Cannot scale to DMAX (max. sing. value is 0)
            3  => Error return from ZLAGGE, ZLAGHE or ZLAGSY
Author
Univ. of Tennessee
Univ. of California Berkeley
Univ. of Colorado Denver
NAG Ltd.
Date
November 2011

Definition at line 342 of file zlatmt.f.

342 *
343 * -- LAPACK computational routine (version 3.4.0) --
344 * -- LAPACK is a software package provided by Univ. of Tennessee, --
345 * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
346 * November 2011
347 *
348 * .. Scalar Arguments ..
349  DOUBLE PRECISION cond, dmax
350  INTEGER info, kl, ku, lda, m, mode, n, rank
351  CHARACTER dist, pack, sym
352 * ..
353 * .. Array Arguments ..
354  COMPLEX*16 a( lda, * ), work( * )
355  DOUBLE PRECISION d( * )
356  INTEGER iseed( 4 )
357 * ..
358 *
359 * =====================================================================
360 *
361 * .. Parameters ..
362  DOUBLE PRECISION zero
363  parameter ( zero = 0.0d+0 )
364  DOUBLE PRECISION one
365  parameter ( one = 1.0d+0 )
366  COMPLEX*16 czero
367  parameter ( czero = ( 0.0d+0, 0.0d+0 ) )
368  DOUBLE PRECISION twopi
369  parameter ( twopi = 6.2831853071795864769252867663d+0 )
370 * ..
371 * .. Local Scalars ..
372  COMPLEX*16 c, ct, dummy, extra, s, st, ztemp
373  DOUBLE PRECISION alpha, angle, realc, temp
374  INTEGER i, ic, icol, idist, iendch, iinfo, il, ilda,
375  $ ioffg, ioffst, ipack, ipackg, ir, ir1, ir2,
376  $ irow, irsign, iskew, isym, isympk, j, jc, jch,
377  $ jkl, jku, jr, k, llb, minlda, mnmin, mr, nc,
378  $ uub
379  LOGICAL csym, givens, ilextr, iltemp, topdwn
380 * ..
381 * .. External Functions ..
382  COMPLEX*16 zlarnd
383  DOUBLE PRECISION dlarnd
384  LOGICAL lsame
385  EXTERNAL zlarnd, dlarnd, lsame
386 * ..
387 * .. External Subroutines ..
388  EXTERNAL dlatm7, dscal, xerbla, zlagge, zlaghe,
390 * ..
391 * .. Intrinsic Functions ..
392  INTRINSIC abs, cos, dble, dcmplx, dconjg, max, min, mod,
393  $ sin
394 * ..
395 * .. Executable Statements ..
396 *
397 * 1) Decode and Test the input parameters.
398 * Initialize flags & seed.
399 *
400  info = 0
401 *
402 * Quick return if possible
403 *
404  IF( m.EQ.0 .OR. n.EQ.0 )
405  $ RETURN
406 *
407 * Decode DIST
408 *
409  IF( lsame( dist, 'U' ) ) THEN
410  idist = 1
411  ELSE IF( lsame( dist, 'S' ) ) THEN
412  idist = 2
413  ELSE IF( lsame( dist, 'N' ) ) THEN
414  idist = 3
415  ELSE
416  idist = -1
417  END IF
418 *
419 * Decode SYM
420 *
421  IF( lsame( sym, 'N' ) ) THEN
422  isym = 1
423  irsign = 0
424  csym = .false.
425  ELSE IF( lsame( sym, 'P' ) ) THEN
426  isym = 2
427  irsign = 0
428  csym = .false.
429  ELSE IF( lsame( sym, 'S' ) ) THEN
430  isym = 2
431  irsign = 0
432  csym = .true.
433  ELSE IF( lsame( sym, 'H' ) ) THEN
434  isym = 2
435  irsign = 1
436  csym = .false.
437  ELSE
438  isym = -1
439  END IF
440 *
441 * Decode PACK
442 *
443  isympk = 0
444  IF( lsame( pack, 'N' ) ) THEN
445  ipack = 0
446  ELSE IF( lsame( pack, 'U' ) ) THEN
447  ipack = 1
448  isympk = 1
449  ELSE IF( lsame( pack, 'L' ) ) THEN
450  ipack = 2
451  isympk = 1
452  ELSE IF( lsame( pack, 'C' ) ) THEN
453  ipack = 3
454  isympk = 2
455  ELSE IF( lsame( pack, 'R' ) ) THEN
456  ipack = 4
457  isympk = 3
458  ELSE IF( lsame( pack, 'B' ) ) THEN
459  ipack = 5
460  isympk = 3
461  ELSE IF( lsame( pack, 'Q' ) ) THEN
462  ipack = 6
463  isympk = 2
464  ELSE IF( lsame( pack, 'Z' ) ) THEN
465  ipack = 7
466  ELSE
467  ipack = -1
468  END IF
469 *
470 * Set certain internal parameters
471 *
472  mnmin = min( m, n )
473  llb = min( kl, m-1 )
474  uub = min( ku, n-1 )
475  mr = min( m, n+llb )
476  nc = min( n, m+uub )
477 *
478  IF( ipack.EQ.5 .OR. ipack.EQ.6 ) THEN
479  minlda = uub + 1
480  ELSE IF( ipack.EQ.7 ) THEN
481  minlda = llb + uub + 1
482  ELSE
483  minlda = m
484  END IF
485 *
486 * Use Givens rotation method if bandwidth small enough,
487 * or if LDA is too small to store the matrix unpacked.
488 *
489  givens = .false.
490  IF( isym.EQ.1 ) THEN
491  IF( dble( llb+uub ).LT.0.3d0*dble( max( 1, mr+nc ) ) )
492  $ givens = .true.
493  ELSE
494  IF( 2*llb.LT.m )
495  $ givens = .true.
496  END IF
497  IF( lda.LT.m .AND. lda.GE.minlda )
498  $ givens = .true.
499 *
500 * Set INFO if an error
501 *
502  IF( m.LT.0 ) THEN
503  info = -1
504  ELSE IF( m.NE.n .AND. isym.NE.1 ) THEN
505  info = -1
506  ELSE IF( n.LT.0 ) THEN
507  info = -2
508  ELSE IF( idist.EQ.-1 ) THEN
509  info = -3
510  ELSE IF( isym.EQ.-1 ) THEN
511  info = -5
512  ELSE IF( abs( mode ).GT.6 ) THEN
513  info = -7
514  ELSE IF( ( mode.NE.0 .AND. abs( mode ).NE.6 ) .AND. cond.LT.one )
515  $ THEN
516  info = -8
517  ELSE IF( kl.LT.0 ) THEN
518  info = -10
519  ELSE IF( ku.LT.0 .OR. ( isym.NE.1 .AND. kl.NE.ku ) ) THEN
520  info = -11
521  ELSE IF( ipack.EQ.-1 .OR. ( isympk.EQ.1 .AND. isym.EQ.1 ) .OR.
522  $ ( isympk.EQ.2 .AND. isym.EQ.1 .AND. kl.GT.0 ) .OR.
523  $ ( isympk.EQ.3 .AND. isym.EQ.1 .AND. ku.GT.0 ) .OR.
524  $ ( isympk.NE.0 .AND. m.NE.n ) ) THEN
525  info = -12
526  ELSE IF( lda.LT.max( 1, minlda ) ) THEN
527  info = -14
528  END IF
529 *
530  IF( info.NE.0 ) THEN
531  CALL xerbla( 'ZLATMT', -info )
532  RETURN
533  END IF
534 *
535 * Initialize random number generator
536 *
537  DO 100 i = 1, 4
538  iseed( i ) = mod( abs( iseed( i ) ), 4096 )
539  100 CONTINUE
540 *
541  IF( mod( iseed( 4 ), 2 ).NE.1 )
542  $ iseed( 4 ) = iseed( 4 ) + 1
543 *
544 * 2) Set up D if indicated.
545 *
546 * Compute D according to COND and MODE
547 *
548  CALL dlatm7( mode, cond, irsign, idist, iseed, d, mnmin, rank,
549  $ iinfo )
550  IF( iinfo.NE.0 ) THEN
551  info = 1
552  RETURN
553  END IF
554 *
555 * Choose Top-Down if D is (apparently) increasing,
556 * Bottom-Up if D is (apparently) decreasing.
557 *
558  IF( abs( d( 1 ) ).LE.abs( d( rank ) ) ) THEN
559  topdwn = .true.
560  ELSE
561  topdwn = .false.
562  END IF
563 *
564  IF( mode.NE.0 .AND. abs( mode ).NE.6 ) THEN
565 *
566 * Scale by DMAX
567 *
568  temp = abs( d( 1 ) )
569  DO 110 i = 2, rank
570  temp = max( temp, abs( d( i ) ) )
571  110 CONTINUE
572 *
573  IF( temp.GT.zero ) THEN
574  alpha = dmax / temp
575  ELSE
576  info = 2
577  RETURN
578  END IF
579 *
580  CALL dscal( rank, alpha, d, 1 )
581 *
582  END IF
583 *
584  CALL zlaset( 'Full', lda, n, czero, czero, a, lda )
585 *
586 * 3) Generate Banded Matrix using Givens rotations.
587 * Also the special case of UUB=LLB=0
588 *
589 * Compute Addressing constants to cover all
590 * storage formats. Whether GE, HE, SY, GB, HB, or SB,
591 * upper or lower triangle or both,
592 * the (i,j)-th element is in
593 * A( i - ISKEW*j + IOFFST, j )
594 *
595  IF( ipack.GT.4 ) THEN
596  ilda = lda - 1
597  iskew = 1
598  IF( ipack.GT.5 ) THEN
599  ioffst = uub + 1
600  ELSE
601  ioffst = 1
602  END IF
603  ELSE
604  ilda = lda
605  iskew = 0
606  ioffst = 0
607  END IF
608 *
609 * IPACKG is the format that the matrix is generated in. If this is
610 * different from IPACK, then the matrix must be repacked at the
611 * end. It also signals how to compute the norm, for scaling.
612 *
613  ipackg = 0
614 *
615 * Diagonal Matrix -- We are done, unless it
616 * is to be stored HP/SP/PP/TP (PACK='R' or 'C')
617 *
618  IF( llb.EQ.0 .AND. uub.EQ.0 ) THEN
619  DO 120 j = 1, mnmin
620  a( ( 1-iskew )*j+ioffst, j ) = dcmplx( d( j ) )
621  120 CONTINUE
622 *
623  IF( ipack.LE.2 .OR. ipack.GE.5 )
624  $ ipackg = ipack
625 *
626  ELSE IF( givens ) THEN
627 *
628 * Check whether to use Givens rotations,
629 * Householder transformations, or nothing.
630 *
631  IF( isym.EQ.1 ) THEN
632 *
633 * Non-symmetric -- A = U D V
634 *
635  IF( ipack.GT.4 ) THEN
636  ipackg = ipack
637  ELSE
638  ipackg = 0
639  END IF
640 *
641  DO 130 j = 1, mnmin
642  a( ( 1-iskew )*j+ioffst, j ) = dcmplx( d( j ) )
643  130 CONTINUE
644 *
645  IF( topdwn ) THEN
646  jkl = 0
647  DO 160 jku = 1, uub
648 *
649 * Transform from bandwidth JKL, JKU-1 to JKL, JKU
650 *
651 * Last row actually rotated is M
652 * Last column actually rotated is MIN( M+JKU, N )
653 *
654  DO 150 jr = 1, min( m+jku, n ) + jkl - 1
655  extra = czero
656  angle = twopi*dlarnd( 1, iseed )
657  c = cos( angle )*zlarnd( 5, iseed )
658  s = sin( angle )*zlarnd( 5, iseed )
659  icol = max( 1, jr-jkl )
660  IF( jr.LT.m ) THEN
661  il = min( n, jr+jku ) + 1 - icol
662  CALL zlarot( .true., jr.GT.jkl, .false., il, c,
663  $ s, a( jr-iskew*icol+ioffst, icol ),
664  $ ilda, extra, dummy )
665  END IF
666 *
667 * Chase "EXTRA" back up
668 *
669  ir = jr
670  ic = icol
671  DO 140 jch = jr - jkl, 1, -jkl - jku
672  IF( ir.LT.m ) THEN
673  CALL zlartg( a( ir+1-iskew*( ic+1 )+ioffst,
674  $ ic+1 ), extra, realc, s, dummy )
675  dummy = dlarnd( 5, iseed )
676  c = dconjg( realc*dummy )
677  s = dconjg( -s*dummy )
678  END IF
679  irow = max( 1, jch-jku )
680  il = ir + 2 - irow
681  ztemp = czero
682  iltemp = jch.GT.jku
683  CALL zlarot( .false., iltemp, .true., il, c, s,
684  $ a( irow-iskew*ic+ioffst, ic ),
685  $ ilda, ztemp, extra )
686  IF( iltemp ) THEN
687  CALL zlartg( a( irow+1-iskew*( ic+1 )+ioffst,
688  $ ic+1 ), ztemp, realc, s, dummy )
689  dummy = zlarnd( 5, iseed )
690  c = dconjg( realc*dummy )
691  s = dconjg( -s*dummy )
692 *
693  icol = max( 1, jch-jku-jkl )
694  il = ic + 2 - icol
695  extra = czero
696  CALL zlarot( .true., jch.GT.jku+jkl, .true.,
697  $ il, c, s, a( irow-iskew*icol+
698  $ ioffst, icol ), ilda, extra,
699  $ ztemp )
700  ic = icol
701  ir = irow
702  END IF
703  140 CONTINUE
704  150 CONTINUE
705  160 CONTINUE
706 *
707  jku = uub
708  DO 190 jkl = 1, llb
709 *
710 * Transform from bandwidth JKL-1, JKU to JKL, JKU
711 *
712  DO 180 jc = 1, min( n+jkl, m ) + jku - 1
713  extra = czero
714  angle = twopi*dlarnd( 1, iseed )
715  c = cos( angle )*zlarnd( 5, iseed )
716  s = sin( angle )*zlarnd( 5, iseed )
717  irow = max( 1, jc-jku )
718  IF( jc.LT.n ) THEN
719  il = min( m, jc+jkl ) + 1 - irow
720  CALL zlarot( .false., jc.GT.jku, .false., il, c,
721  $ s, a( irow-iskew*jc+ioffst, jc ),
722  $ ilda, extra, dummy )
723  END IF
724 *
725 * Chase "EXTRA" back up
726 *
727  ic = jc
728  ir = irow
729  DO 170 jch = jc - jku, 1, -jkl - jku
730  IF( ic.LT.n ) THEN
731  CALL zlartg( a( ir+1-iskew*( ic+1 )+ioffst,
732  $ ic+1 ), extra, realc, s, dummy )
733  dummy = zlarnd( 5, iseed )
734  c = dconjg( realc*dummy )
735  s = dconjg( -s*dummy )
736  END IF
737  icol = max( 1, jch-jkl )
738  il = ic + 2 - icol
739  ztemp = czero
740  iltemp = jch.GT.jkl
741  CALL zlarot( .true., iltemp, .true., il, c, s,
742  $ a( ir-iskew*icol+ioffst, icol ),
743  $ ilda, ztemp, extra )
744  IF( iltemp ) THEN
745  CALL zlartg( a( ir+1-iskew*( icol+1 )+ioffst,
746  $ icol+1 ), ztemp, realc, s,
747  $ dummy )
748  dummy = zlarnd( 5, iseed )
749  c = dconjg( realc*dummy )
750  s = dconjg( -s*dummy )
751  irow = max( 1, jch-jkl-jku )
752  il = ir + 2 - irow
753  extra = czero
754  CALL zlarot( .false., jch.GT.jkl+jku, .true.,
755  $ il, c, s, a( irow-iskew*icol+
756  $ ioffst, icol ), ilda, extra,
757  $ ztemp )
758  ic = icol
759  ir = irow
760  END IF
761  170 CONTINUE
762  180 CONTINUE
763  190 CONTINUE
764 *
765  ELSE
766 *
767 * Bottom-Up -- Start at the bottom right.
768 *
769  jkl = 0
770  DO 220 jku = 1, uub
771 *
772 * Transform from bandwidth JKL, JKU-1 to JKL, JKU
773 *
774 * First row actually rotated is M
775 * First column actually rotated is MIN( M+JKU, N )
776 *
777  iendch = min( m, n+jkl ) - 1
778  DO 210 jc = min( m+jku, n ) - 1, 1 - jkl, -1
779  extra = czero
780  angle = twopi*dlarnd( 1, iseed )
781  c = cos( angle )*zlarnd( 5, iseed )
782  s = sin( angle )*zlarnd( 5, iseed )
783  irow = max( 1, jc-jku+1 )
784  IF( jc.GT.0 ) THEN
785  il = min( m, jc+jkl+1 ) + 1 - irow
786  CALL zlarot( .false., .false., jc+jkl.LT.m, il,
787  $ c, s, a( irow-iskew*jc+ioffst,
788  $ jc ), ilda, dummy, extra )
789  END IF
790 *
791 * Chase "EXTRA" back down
792 *
793  ic = jc
794  DO 200 jch = jc + jkl, iendch, jkl + jku
795  ilextr = ic.GT.0
796  IF( ilextr ) THEN
797  CALL zlartg( a( jch-iskew*ic+ioffst, ic ),
798  $ extra, realc, s, dummy )
799  dummy = zlarnd( 5, iseed )
800  c = realc*dummy
801  s = s*dummy
802  END IF
803  ic = max( 1, ic )
804  icol = min( n-1, jch+jku )
805  iltemp = jch + jku.LT.n
806  ztemp = czero
807  CALL zlarot( .true., ilextr, iltemp, icol+2-ic,
808  $ c, s, a( jch-iskew*ic+ioffst, ic ),
809  $ ilda, extra, ztemp )
810  IF( iltemp ) THEN
811  CALL zlartg( a( jch-iskew*icol+ioffst,
812  $ icol ), ztemp, realc, s, dummy )
813  dummy = zlarnd( 5, iseed )
814  c = realc*dummy
815  s = s*dummy
816  il = min( iendch, jch+jkl+jku ) + 2 - jch
817  extra = czero
818  CALL zlarot( .false., .true.,
819  $ jch+jkl+jku.LE.iendch, il, c, s,
820  $ a( jch-iskew*icol+ioffst,
821  $ icol ), ilda, ztemp, extra )
822  ic = icol
823  END IF
824  200 CONTINUE
825  210 CONTINUE
826  220 CONTINUE
827 *
828  jku = uub
829  DO 250 jkl = 1, llb
830 *
831 * Transform from bandwidth JKL-1, JKU to JKL, JKU
832 *
833 * First row actually rotated is MIN( N+JKL, M )
834 * First column actually rotated is N
835 *
836  iendch = min( n, m+jku ) - 1
837  DO 240 jr = min( n+jkl, m ) - 1, 1 - jku, -1
838  extra = czero
839  angle = twopi*dlarnd( 1, iseed )
840  c = cos( angle )*zlarnd( 5, iseed )
841  s = sin( angle )*zlarnd( 5, iseed )
842  icol = max( 1, jr-jkl+1 )
843  IF( jr.GT.0 ) THEN
844  il = min( n, jr+jku+1 ) + 1 - icol
845  CALL zlarot( .true., .false., jr+jku.LT.n, il,
846  $ c, s, a( jr-iskew*icol+ioffst,
847  $ icol ), ilda, dummy, extra )
848  END IF
849 *
850 * Chase "EXTRA" back down
851 *
852  ir = jr
853  DO 230 jch = jr + jku, iendch, jkl + jku
854  ilextr = ir.GT.0
855  IF( ilextr ) THEN
856  CALL zlartg( a( ir-iskew*jch+ioffst, jch ),
857  $ extra, realc, s, dummy )
858  dummy = zlarnd( 5, iseed )
859  c = realc*dummy
860  s = s*dummy
861  END IF
862  ir = max( 1, ir )
863  irow = min( m-1, jch+jkl )
864  iltemp = jch + jkl.LT.m
865  ztemp = czero
866  CALL zlarot( .false., ilextr, iltemp, irow+2-ir,
867  $ c, s, a( ir-iskew*jch+ioffst,
868  $ jch ), ilda, extra, ztemp )
869  IF( iltemp ) THEN
870  CALL zlartg( a( irow-iskew*jch+ioffst, jch ),
871  $ ztemp, realc, s, dummy )
872  dummy = zlarnd( 5, iseed )
873  c = realc*dummy
874  s = s*dummy
875  il = min( iendch, jch+jkl+jku ) + 2 - jch
876  extra = czero
877  CALL zlarot( .true., .true.,
878  $ jch+jkl+jku.LE.iendch, il, c, s,
879  $ a( irow-iskew*jch+ioffst, jch ),
880  $ ilda, ztemp, extra )
881  ir = irow
882  END IF
883  230 CONTINUE
884  240 CONTINUE
885  250 CONTINUE
886 *
887  END IF
888 *
889  ELSE
890 *
891 * Symmetric -- A = U D U'
892 * Hermitian -- A = U D U*
893 *
894  ipackg = ipack
895  ioffg = ioffst
896 *
897  IF( topdwn ) THEN
898 *
899 * Top-Down -- Generate Upper triangle only
900 *
901  IF( ipack.GE.5 ) THEN
902  ipackg = 6
903  ioffg = uub + 1
904  ELSE
905  ipackg = 1
906  END IF
907 *
908  DO 260 j = 1, mnmin
909  a( ( 1-iskew )*j+ioffg, j ) = dcmplx( d( j ) )
910  260 CONTINUE
911 *
912  DO 290 k = 1, uub
913  DO 280 jc = 1, n - 1
914  irow = max( 1, jc-k )
915  il = min( jc+1, k+2 )
916  extra = czero
917  ztemp = a( jc-iskew*( jc+1 )+ioffg, jc+1 )
918  angle = twopi*dlarnd( 1, iseed )
919  c = cos( angle )*zlarnd( 5, iseed )
920  s = sin( angle )*zlarnd( 5, iseed )
921  IF( csym ) THEN
922  ct = c
923  st = s
924  ELSE
925  ztemp = dconjg( ztemp )
926  ct = dconjg( c )
927  st = dconjg( s )
928  END IF
929  CALL zlarot( .false., jc.GT.k, .true., il, c, s,
930  $ a( irow-iskew*jc+ioffg, jc ), ilda,
931  $ extra, ztemp )
932  CALL zlarot( .true., .true., .false.,
933  $ min( k, n-jc )+1, ct, st,
934  $ a( ( 1-iskew )*jc+ioffg, jc ), ilda,
935  $ ztemp, dummy )
936 *
937 * Chase EXTRA back up the matrix
938 *
939  icol = jc
940  DO 270 jch = jc - k, 1, -k
941  CALL zlartg( a( jch+1-iskew*( icol+1 )+ioffg,
942  $ icol+1 ), extra, realc, s, dummy )
943  dummy = zlarnd( 5, iseed )
944  c = dconjg( realc*dummy )
945  s = dconjg( -s*dummy )
946  ztemp = a( jch-iskew*( jch+1 )+ioffg, jch+1 )
947  IF( csym ) THEN
948  ct = c
949  st = s
950  ELSE
951  ztemp = dconjg( ztemp )
952  ct = dconjg( c )
953  st = dconjg( s )
954  END IF
955  CALL zlarot( .true., .true., .true., k+2, c, s,
956  $ a( ( 1-iskew )*jch+ioffg, jch ),
957  $ ilda, ztemp, extra )
958  irow = max( 1, jch-k )
959  il = min( jch+1, k+2 )
960  extra = czero
961  CALL zlarot( .false., jch.GT.k, .true., il, ct,
962  $ st, a( irow-iskew*jch+ioffg, jch ),
963  $ ilda, extra, ztemp )
964  icol = jch
965  270 CONTINUE
966  280 CONTINUE
967  290 CONTINUE
968 *
969 * If we need lower triangle, copy from upper. Note that
970 * the order of copying is chosen to work for 'q' -> 'b'
971 *
972  IF( ipack.NE.ipackg .AND. ipack.NE.3 ) THEN
973  DO 320 jc = 1, n
974  irow = ioffst - iskew*jc
975  IF( csym ) THEN
976  DO 300 jr = jc, min( n, jc+uub )
977  a( jr+irow, jc ) = a( jc-iskew*jr+ioffg, jr )
978  300 CONTINUE
979  ELSE
980  DO 310 jr = jc, min( n, jc+uub )
981  a( jr+irow, jc ) = dconjg( a( jc-iskew*jr+
982  $ ioffg, jr ) )
983  310 CONTINUE
984  END IF
985  320 CONTINUE
986  IF( ipack.EQ.5 ) THEN
987  DO 340 jc = n - uub + 1, n
988  DO 330 jr = n + 2 - jc, uub + 1
989  a( jr, jc ) = czero
990  330 CONTINUE
991  340 CONTINUE
992  END IF
993  IF( ipackg.EQ.6 ) THEN
994  ipackg = ipack
995  ELSE
996  ipackg = 0
997  END IF
998  END IF
999  ELSE
1000 *
1001 * Bottom-Up -- Generate Lower triangle only
1002 *
1003  IF( ipack.GE.5 ) THEN
1004  ipackg = 5
1005  IF( ipack.EQ.6 )
1006  $ ioffg = 1
1007  ELSE
1008  ipackg = 2
1009  END IF
1010 *
1011  DO 350 j = 1, mnmin
1012  a( ( 1-iskew )*j+ioffg, j ) = dcmplx( d( j ) )
1013  350 CONTINUE
1014 *
1015  DO 380 k = 1, uub
1016  DO 370 jc = n - 1, 1, -1
1017  il = min( n+1-jc, k+2 )
1018  extra = czero
1019  ztemp = a( 1+( 1-iskew )*jc+ioffg, jc )
1020  angle = twopi*dlarnd( 1, iseed )
1021  c = cos( angle )*zlarnd( 5, iseed )
1022  s = sin( angle )*zlarnd( 5, iseed )
1023  IF( csym ) THEN
1024  ct = c
1025  st = s
1026  ELSE
1027  ztemp = dconjg( ztemp )
1028  ct = dconjg( c )
1029  st = dconjg( s )
1030  END IF
1031  CALL zlarot( .false., .true., n-jc.GT.k, il, c, s,
1032  $ a( ( 1-iskew )*jc+ioffg, jc ), ilda,
1033  $ ztemp, extra )
1034  icol = max( 1, jc-k+1 )
1035  CALL zlarot( .true., .false., .true., jc+2-icol,
1036  $ ct, st, a( jc-iskew*icol+ioffg,
1037  $ icol ), ilda, dummy, ztemp )
1038 *
1039 * Chase EXTRA back down the matrix
1040 *
1041  icol = jc
1042  DO 360 jch = jc + k, n - 1, k
1043  CALL zlartg( a( jch-iskew*icol+ioffg, icol ),
1044  $ extra, realc, s, dummy )
1045  dummy = zlarnd( 5, iseed )
1046  c = realc*dummy
1047  s = s*dummy
1048  ztemp = a( 1+( 1-iskew )*jch+ioffg, jch )
1049  IF( csym ) THEN
1050  ct = c
1051  st = s
1052  ELSE
1053  ztemp = dconjg( ztemp )
1054  ct = dconjg( c )
1055  st = dconjg( s )
1056  END IF
1057  CALL zlarot( .true., .true., .true., k+2, c, s,
1058  $ a( jch-iskew*icol+ioffg, icol ),
1059  $ ilda, extra, ztemp )
1060  il = min( n+1-jch, k+2 )
1061  extra = czero
1062  CALL zlarot( .false., .true., n-jch.GT.k, il,
1063  $ ct, st, a( ( 1-iskew )*jch+ioffg,
1064  $ jch ), ilda, ztemp, extra )
1065  icol = jch
1066  360 CONTINUE
1067  370 CONTINUE
1068  380 CONTINUE
1069 *
1070 * If we need upper triangle, copy from lower. Note that
1071 * the order of copying is chosen to work for 'b' -> 'q'
1072 *
1073  IF( ipack.NE.ipackg .AND. ipack.NE.4 ) THEN
1074  DO 410 jc = n, 1, -1
1075  irow = ioffst - iskew*jc
1076  IF( csym ) THEN
1077  DO 390 jr = jc, max( 1, jc-uub ), -1
1078  a( jr+irow, jc ) = a( jc-iskew*jr+ioffg, jr )
1079  390 CONTINUE
1080  ELSE
1081  DO 400 jr = jc, max( 1, jc-uub ), -1
1082  a( jr+irow, jc ) = dconjg( a( jc-iskew*jr+
1083  $ ioffg, jr ) )
1084  400 CONTINUE
1085  END IF
1086  410 CONTINUE
1087  IF( ipack.EQ.6 ) THEN
1088  DO 430 jc = 1, uub
1089  DO 420 jr = 1, uub + 1 - jc
1090  a( jr, jc ) = czero
1091  420 CONTINUE
1092  430 CONTINUE
1093  END IF
1094  IF( ipackg.EQ.5 ) THEN
1095  ipackg = ipack
1096  ELSE
1097  ipackg = 0
1098  END IF
1099  END IF
1100  END IF
1101 *
1102 * Ensure that the diagonal is real if Hermitian
1103 *
1104  IF( .NOT.csym ) THEN
1105  DO 440 jc = 1, n
1106  irow = ioffst + ( 1-iskew )*jc
1107  a( irow, jc ) = dcmplx( dble( a( irow, jc ) ) )
1108  440 CONTINUE
1109  END IF
1110 *
1111  END IF
1112 *
1113  ELSE
1114 *
1115 * 4) Generate Banded Matrix by first
1116 * Rotating by random Unitary matrices,
1117 * then reducing the bandwidth using Householder
1118 * transformations.
1119 *
1120 * Note: we should get here only if LDA .ge. N
1121 *
1122  IF( isym.EQ.1 ) THEN
1123 *
1124 * Non-symmetric -- A = U D V
1125 *
1126  CALL zlagge( mr, nc, llb, uub, d, a, lda, iseed, work,
1127  $ iinfo )
1128  ELSE
1129 *
1130 * Symmetric -- A = U D U' or
1131 * Hermitian -- A = U D U*
1132 *
1133  IF( csym ) THEN
1134  CALL zlagsy( m, llb, d, a, lda, iseed, work, iinfo )
1135  ELSE
1136  CALL zlaghe( m, llb, d, a, lda, iseed, work, iinfo )
1137  END IF
1138  END IF
1139 *
1140  IF( iinfo.NE.0 ) THEN
1141  info = 3
1142  RETURN
1143  END IF
1144  END IF
1145 *
1146 * 5) Pack the matrix
1147 *
1148  IF( ipack.NE.ipackg ) THEN
1149  IF( ipack.EQ.1 ) THEN
1150 *
1151 * 'U' -- Upper triangular, not packed
1152 *
1153  DO 460 j = 1, m
1154  DO 450 i = j + 1, m
1155  a( i, j ) = czero
1156  450 CONTINUE
1157  460 CONTINUE
1158 *
1159  ELSE IF( ipack.EQ.2 ) THEN
1160 *
1161 * 'L' -- Lower triangular, not packed
1162 *
1163  DO 480 j = 2, m
1164  DO 470 i = 1, j - 1
1165  a( i, j ) = czero
1166  470 CONTINUE
1167  480 CONTINUE
1168 *
1169  ELSE IF( ipack.EQ.3 ) THEN
1170 *
1171 * 'C' -- Upper triangle packed Columnwise.
1172 *
1173  icol = 1
1174  irow = 0
1175  DO 500 j = 1, m
1176  DO 490 i = 1, j
1177  irow = irow + 1
1178  IF( irow.GT.lda ) THEN
1179  irow = 1
1180  icol = icol + 1
1181  END IF
1182  a( irow, icol ) = a( i, j )
1183  490 CONTINUE
1184  500 CONTINUE
1185 *
1186  ELSE IF( ipack.EQ.4 ) THEN
1187 *
1188 * 'R' -- Lower triangle packed Columnwise.
1189 *
1190  icol = 1
1191  irow = 0
1192  DO 520 j = 1, m
1193  DO 510 i = j, m
1194  irow = irow + 1
1195  IF( irow.GT.lda ) THEN
1196  irow = 1
1197  icol = icol + 1
1198  END IF
1199  a( irow, icol ) = a( i, j )
1200  510 CONTINUE
1201  520 CONTINUE
1202 *
1203  ELSE IF( ipack.GE.5 ) THEN
1204 *
1205 * 'B' -- The lower triangle is packed as a band matrix.
1206 * 'Q' -- The upper triangle is packed as a band matrix.
1207 * 'Z' -- The whole matrix is packed as a band matrix.
1208 *
1209  IF( ipack.EQ.5 )
1210  $ uub = 0
1211  IF( ipack.EQ.6 )
1212  $ llb = 0
1213 *
1214  DO 540 j = 1, uub
1215  DO 530 i = min( j+llb, m ), 1, -1
1216  a( i-j+uub+1, j ) = a( i, j )
1217  530 CONTINUE
1218  540 CONTINUE
1219 *
1220  DO 560 j = uub + 2, n
1221  DO 550 i = j - uub, min( j+llb, m )
1222  a( i-j+uub+1, j ) = a( i, j )
1223  550 CONTINUE
1224  560 CONTINUE
1225  END IF
1226 *
1227 * If packed, zero out extraneous elements.
1228 *
1229 * Symmetric/Triangular Packed --
1230 * zero out everything after A(IROW,ICOL)
1231 *
1232  IF( ipack.EQ.3 .OR. ipack.EQ.4 ) THEN
1233  DO 580 jc = icol, m
1234  DO 570 jr = irow + 1, lda
1235  a( jr, jc ) = czero
1236  570 CONTINUE
1237  irow = 0
1238  580 CONTINUE
1239 *
1240  ELSE IF( ipack.GE.5 ) THEN
1241 *
1242 * Packed Band --
1243 * 1st row is now in A( UUB+2-j, j), zero above it
1244 * m-th row is now in A( M+UUB-j,j), zero below it
1245 * last non-zero diagonal is now in A( UUB+LLB+1,j ),
1246 * zero below it, too.
1247 *
1248  ir1 = uub + llb + 2
1249  ir2 = uub + m + 2
1250  DO 610 jc = 1, n
1251  DO 590 jr = 1, uub + 1 - jc
1252  a( jr, jc ) = czero
1253  590 CONTINUE
1254  DO 600 jr = max( 1, min( ir1, ir2-jc ) ), lda
1255  a( jr, jc ) = czero
1256  600 CONTINUE
1257  610 CONTINUE
1258  END IF
1259  END IF
1260 *
1261  RETURN
1262 *
1263 * End of ZLATMT
1264 *
subroutine zlaghe(N, K, D, A, LDA, ISEED, WORK, INFO)
ZLAGHE
Definition: zlaghe.f:104
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:108
subroutine zlagsy(N, K, D, A, LDA, ISEED, WORK, INFO)
ZLAGSY
Definition: zlagsy.f:104
subroutine xerbla(SRNAME, INFO)
XERBLA
Definition: xerbla.f:62
double precision function dlarnd(IDIST, ISEED)
DLARND
Definition: dlarnd.f:75
subroutine dscal(N, DA, DX, INCX)
DSCAL
Definition: dscal.f:55
subroutine dlatm7(MODE, COND, IRSIGN, IDIST, ISEED, D, N, RANK, INFO)
DLATM7
Definition: dlatm7.f:124
subroutine zlagge(M, N, KL, KU, D, A, LDA, ISEED, WORK, INFO)
ZLAGGE
Definition: zlagge.f:116
logical function lsame(CA, CB)
LSAME
Definition: lsame.f:55
complex *16 function zlarnd(IDIST, ISEED)
ZLARND
Definition: zlarnd.f:77
subroutine zlarot(LROWS, LLEFT, LRIGHT, NL, C, S, A, LDA, XLEFT, XRIGHT)
ZLAROT
Definition: zlarot.f:231
subroutine zlartg(F, G, CS, SN, R)
ZLARTG generates a plane rotation with real cosine and complex sine.
Definition: zlartg.f:105

Here is the call graph for this function:

Here is the caller graph for this function: