LAPACK  3.10.0
LAPACK: Linear Algebra PACKage

◆ zlatmt()

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.

Definition at line 338 of file zlatmt.f.

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