LAPACK  3.8.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.
Date
December 2016

Definition at line 334 of file zlatms.f.

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