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

◆ zlatmt()

subroutine zlatmt ( integer  m,
integer  n,
character  dist,
integer, dimension( 4 )  iseed,
character  sym,
double precision, dimension( * )  d,
integer  mode,
double precision  cond,
double precision  dmax,
integer  rank,
integer  kl,
integer  ku,
character  pack,
complex*16, dimension( lda, * )  a,
integer  lda,
complex*16, dimension( * )  work,
integer  info 
)

ZLATMT

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

    ZLATMT operates by applying the following sequence of
    operations:

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

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

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

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

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

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

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

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

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

Definition at line 338 of file zlatmt.f.

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