LAPACK  3.10.0
LAPACK: Linear Algebra PACKage

◆ zlatms()

subroutine zlatms ( 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  KL,
integer  KU,
character  PACK,
complex*16, dimension( lda, * )  A,
integer  LDA,
complex*16, dimension( * )  WORK,
integer  INFO 
)

ZLATMS

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

    ZLATMS 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 ZLATMS
           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:N)=1.0/COND
           MODE = 2 sets D(1:N-1)=1 and D(N)=1.0/COND
           MODE = 3 sets D(I)=COND**(-(I-1)/(N-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]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 ZLATMS 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 DLATM1
            2  => Cannot scale to DMAX (max. sing. value is 0)
            3  => Error return from ZLAGGE, CLAGHE or CLAGSY
Author
Univ. of Tennessee
Univ. of California Berkeley
Univ. of Colorado Denver
NAG Ltd.

Definition at line 330 of file zlatms.f.

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