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

◆ 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,
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 xerbla(srname, info)
Definition cblat2.f:3285
double precision function dlarnd(idist, iseed)
DLARND
Definition dlarnd.f:73
subroutine dlatm1(mode, cond, irsign, idist, iseed, d, n, info)
DLATM1
Definition dlatm1.f:135
subroutine zlartg(f, g, c, s, r)
ZLARTG generates a plane rotation with real cosine and complex sine.
Definition zlartg.f90:116
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
logical function lsame(ca, cb)
LSAME
Definition lsame.f:48
subroutine dscal(n, da, dx, incx)
DSCAL
Definition dscal.f:79
subroutine zlagge(m, n, kl, ku, d, a, lda, iseed, work, info)
ZLAGGE
Definition zlagge.f:114
subroutine zlaghe(n, k, d, a, lda, iseed, work, info)
ZLAGHE
Definition zlaghe.f:102
subroutine zlagsy(n, k, d, a, lda, iseed, work, info)
ZLAGSY
Definition zlagsy.f:102
complex *16 function zlarnd(idist, iseed)
ZLARND
Definition zlarnd.f:75
subroutine zlarot(lrows, lleft, lright, nl, c, s, a, lda, xleft, xright)
ZLAROT
Definition zlarot.f:229
Here is the call graph for this function: