LAPACK 3.12.1
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, 
!>              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,
660 $ c,
661 $ s, a( jr-iskew*icol+ioffst, icol ),
662 $ ilda, extra, dummy )
663 END IF
664*
665* Chase "EXTRA" back up
666*
667 ir = jr
668 ic = icol
669 DO 140 jch = jr - jkl, 1, -jkl - jku
670 IF( ir.LT.m ) THEN
671 CALL zlartg( a( ir+1-iskew*( ic+1 )+ioffst,
672 $ ic+1 ), extra, realc, s, dummy )
673 dummy = dlarnd( 5, iseed )
674 c = dconjg( realc*dummy )
675 s = dconjg( -s*dummy )
676 END IF
677 irow = max( 1, jch-jku )
678 il = ir + 2 - irow
679 ztemp = czero
680 iltemp = jch.GT.jku
681 CALL zlarot( .false., iltemp, .true., il, c,
682 $ s,
683 $ a( irow-iskew*ic+ioffst, ic ),
684 $ ilda, ztemp, extra )
685 IF( iltemp ) THEN
686 CALL zlartg( a( irow+1-iskew*( ic+1 )+ioffst,
687 $ ic+1 ), ztemp, realc, s, dummy )
688 dummy = zlarnd( 5, iseed )
689 c = dconjg( realc*dummy )
690 s = dconjg( -s*dummy )
691*
692 icol = max( 1, jch-jku-jkl )
693 il = ic + 2 - icol
694 extra = czero
695 CALL zlarot( .true., jch.GT.jku+jkl,
696 $ .true.,
697 $ il, c, s, a( irow-iskew*icol+
698 $ ioffst, icol ), ilda, extra,
699 $ ztemp )
700 ic = icol
701 ir = irow
702 END IF
703 140 CONTINUE
704 150 CONTINUE
705 160 CONTINUE
706*
707 jku = uub
708 DO 190 jkl = 1, llb
709*
710* Transform from bandwidth JKL-1, JKU to JKL, JKU
711*
712 DO 180 jc = 1, min( n+jkl, m ) + jku - 1
713 extra = czero
714 angle = twopi*dlarnd( 1, iseed )
715 c = cos( angle )*zlarnd( 5, iseed )
716 s = sin( angle )*zlarnd( 5, iseed )
717 irow = max( 1, jc-jku )
718 IF( jc.LT.n ) THEN
719 il = min( m, jc+jkl ) + 1 - irow
720 CALL zlarot( .false., jc.GT.jku, .false., il,
721 $ c,
722 $ s, a( irow-iskew*jc+ioffst, jc ),
723 $ ilda, extra, dummy )
724 END IF
725*
726* Chase "EXTRA" back up
727*
728 ic = jc
729 ir = irow
730 DO 170 jch = jc - jku, 1, -jkl - jku
731 IF( ic.LT.n ) THEN
732 CALL zlartg( a( ir+1-iskew*( ic+1 )+ioffst,
733 $ ic+1 ), extra, realc, s, dummy )
734 dummy = zlarnd( 5, iseed )
735 c = dconjg( realc*dummy )
736 s = dconjg( -s*dummy )
737 END IF
738 icol = max( 1, jch-jkl )
739 il = ic + 2 - icol
740 ztemp = czero
741 iltemp = jch.GT.jkl
742 CALL zlarot( .true., iltemp, .true., il, c,
743 $ s,
744 $ a( ir-iskew*icol+ioffst, icol ),
745 $ ilda, ztemp, extra )
746 IF( iltemp ) THEN
747 CALL zlartg( a( ir+1-iskew*( icol+1 )+ioffst,
748 $ icol+1 ), ztemp, realc, s,
749 $ dummy )
750 dummy = zlarnd( 5, iseed )
751 c = dconjg( realc*dummy )
752 s = dconjg( -s*dummy )
753 irow = max( 1, jch-jkl-jku )
754 il = ir + 2 - irow
755 extra = czero
756 CALL zlarot( .false., jch.GT.jkl+jku,
757 $ .true.,
758 $ il, c, s, a( irow-iskew*icol+
759 $ ioffst, icol ), ilda, extra,
760 $ ztemp )
761 ic = icol
762 ir = irow
763 END IF
764 170 CONTINUE
765 180 CONTINUE
766 190 CONTINUE
767*
768 ELSE
769*
770* Bottom-Up -- Start at the bottom right.
771*
772 jkl = 0
773 DO 220 jku = 1, uub
774*
775* Transform from bandwidth JKL, JKU-1 to JKL, JKU
776*
777* First row actually rotated is M
778* First column actually rotated is MIN( M+JKU, N )
779*
780 iendch = min( m, n+jkl ) - 1
781 DO 210 jc = min( m+jku, n ) - 1, 1 - jkl, -1
782 extra = czero
783 angle = twopi*dlarnd( 1, iseed )
784 c = cos( angle )*zlarnd( 5, iseed )
785 s = sin( angle )*zlarnd( 5, iseed )
786 irow = max( 1, jc-jku+1 )
787 IF( jc.GT.0 ) THEN
788 il = min( m, jc+jkl+1 ) + 1 - irow
789 CALL zlarot( .false., .false., jc+jkl.LT.m,
790 $ il,
791 $ c, s, a( irow-iskew*jc+ioffst,
792 $ jc ), ilda, dummy, extra )
793 END IF
794*
795* Chase "EXTRA" back down
796*
797 ic = jc
798 DO 200 jch = jc + jkl, iendch, jkl + jku
799 ilextr = ic.GT.0
800 IF( ilextr ) THEN
801 CALL zlartg( a( jch-iskew*ic+ioffst, ic ),
802 $ extra, realc, s, dummy )
803 dummy = zlarnd( 5, iseed )
804 c = realc*dummy
805 s = s*dummy
806 END IF
807 ic = max( 1, ic )
808 icol = min( n-1, jch+jku )
809 iltemp = jch + jku.LT.n
810 ztemp = czero
811 CALL zlarot( .true., ilextr, iltemp,
812 $ icol+2-ic,
813 $ c, s, a( jch-iskew*ic+ioffst, ic ),
814 $ ilda, extra, ztemp )
815 IF( iltemp ) THEN
816 CALL zlartg( a( jch-iskew*icol+ioffst,
817 $ icol ), ztemp, realc, s, dummy )
818 dummy = zlarnd( 5, iseed )
819 c = realc*dummy
820 s = s*dummy
821 il = min( iendch, jch+jkl+jku ) + 2 - jch
822 extra = czero
823 CALL zlarot( .false., .true.,
824 $ jch+jkl+jku.LE.iendch, il, c, s,
825 $ a( jch-iskew*icol+ioffst,
826 $ icol ), ilda, ztemp, extra )
827 ic = icol
828 END IF
829 200 CONTINUE
830 210 CONTINUE
831 220 CONTINUE
832*
833 jku = uub
834 DO 250 jkl = 1, llb
835*
836* Transform from bandwidth JKL-1, JKU to JKL, JKU
837*
838* First row actually rotated is MIN( N+JKL, M )
839* First column actually rotated is N
840*
841 iendch = min( n, m+jku ) - 1
842 DO 240 jr = min( n+jkl, m ) - 1, 1 - jku, -1
843 extra = czero
844 angle = twopi*dlarnd( 1, iseed )
845 c = cos( angle )*zlarnd( 5, iseed )
846 s = sin( angle )*zlarnd( 5, iseed )
847 icol = max( 1, jr-jkl+1 )
848 IF( jr.GT.0 ) THEN
849 il = min( n, jr+jku+1 ) + 1 - icol
850 CALL zlarot( .true., .false., jr+jku.LT.n,
851 $ il,
852 $ c, s, a( jr-iskew*icol+ioffst,
853 $ icol ), ilda, dummy, extra )
854 END IF
855*
856* Chase "EXTRA" back down
857*
858 ir = jr
859 DO 230 jch = jr + jku, iendch, jkl + jku
860 ilextr = ir.GT.0
861 IF( ilextr ) THEN
862 CALL zlartg( a( ir-iskew*jch+ioffst, jch ),
863 $ extra, realc, s, dummy )
864 dummy = zlarnd( 5, iseed )
865 c = realc*dummy
866 s = s*dummy
867 END IF
868 ir = max( 1, ir )
869 irow = min( m-1, jch+jkl )
870 iltemp = jch + jkl.LT.m
871 ztemp = czero
872 CALL zlarot( .false., ilextr, iltemp,
873 $ irow+2-ir,
874 $ c, s, a( ir-iskew*jch+ioffst,
875 $ jch ), ilda, extra, ztemp )
876 IF( iltemp ) THEN
877 CALL zlartg( a( irow-iskew*jch+ioffst, jch ),
878 $ ztemp, realc, s, dummy )
879 dummy = zlarnd( 5, iseed )
880 c = realc*dummy
881 s = s*dummy
882 il = min( iendch, jch+jkl+jku ) + 2 - jch
883 extra = czero
884 CALL zlarot( .true., .true.,
885 $ jch+jkl+jku.LE.iendch, il, c, s,
886 $ a( irow-iskew*jch+ioffst, jch ),
887 $ ilda, ztemp, extra )
888 ir = irow
889 END IF
890 230 CONTINUE
891 240 CONTINUE
892 250 CONTINUE
893*
894 END IF
895*
896 ELSE
897*
898* Symmetric -- A = U D U'
899* Hermitian -- A = U D U*
900*
901 ipackg = ipack
902 ioffg = ioffst
903*
904 IF( topdwn ) THEN
905*
906* Top-Down -- Generate Upper triangle only
907*
908 IF( ipack.GE.5 ) THEN
909 ipackg = 6
910 ioffg = uub + 1
911 ELSE
912 ipackg = 1
913 END IF
914*
915 DO 260 j = 1, mnmin
916 a( ( 1-iskew )*j+ioffg, j ) = dcmplx( d( j ) )
917 260 CONTINUE
918*
919 DO 290 k = 1, uub
920 DO 280 jc = 1, n - 1
921 irow = max( 1, jc-k )
922 il = min( jc+1, k+2 )
923 extra = czero
924 ztemp = a( jc-iskew*( jc+1 )+ioffg, jc+1 )
925 angle = twopi*dlarnd( 1, iseed )
926 c = cos( angle )*zlarnd( 5, iseed )
927 s = sin( angle )*zlarnd( 5, iseed )
928 IF( csym ) THEN
929 ct = c
930 st = s
931 ELSE
932 ztemp = dconjg( ztemp )
933 ct = dconjg( c )
934 st = dconjg( s )
935 END IF
936 CALL zlarot( .false., jc.GT.k, .true., il, c, s,
937 $ a( irow-iskew*jc+ioffg, jc ), ilda,
938 $ extra, ztemp )
939 CALL zlarot( .true., .true., .false.,
940 $ min( k, n-jc )+1, ct, st,
941 $ a( ( 1-iskew )*jc+ioffg, jc ), ilda,
942 $ ztemp, dummy )
943*
944* Chase EXTRA back up the matrix
945*
946 icol = jc
947 DO 270 jch = jc - k, 1, -k
948 CALL zlartg( a( jch+1-iskew*( icol+1 )+ioffg,
949 $ icol+1 ), extra, realc, s, dummy )
950 dummy = zlarnd( 5, iseed )
951 c = dconjg( realc*dummy )
952 s = dconjg( -s*dummy )
953 ztemp = a( jch-iskew*( jch+1 )+ioffg, jch+1 )
954 IF( csym ) THEN
955 ct = c
956 st = s
957 ELSE
958 ztemp = dconjg( ztemp )
959 ct = dconjg( c )
960 st = dconjg( s )
961 END IF
962 CALL zlarot( .true., .true., .true., k+2, c,
963 $ s,
964 $ a( ( 1-iskew )*jch+ioffg, jch ),
965 $ ilda, ztemp, extra )
966 irow = max( 1, jch-k )
967 il = min( jch+1, k+2 )
968 extra = czero
969 CALL zlarot( .false., jch.GT.k, .true., il,
970 $ ct,
971 $ st, a( irow-iskew*jch+ioffg, jch ),
972 $ ilda, extra, ztemp )
973 icol = jch
974 270 CONTINUE
975 280 CONTINUE
976 290 CONTINUE
977*
978* If we need lower triangle, copy from upper. Note that
979* the order of copying is chosen to work for 'q' -> 'b'
980*
981 IF( ipack.NE.ipackg .AND. ipack.NE.3 ) THEN
982 DO 320 jc = 1, n
983 irow = ioffst - iskew*jc
984 IF( csym ) THEN
985 DO 300 jr = jc, min( n, jc+uub )
986 a( jr+irow, jc ) = a( jc-iskew*jr+ioffg, jr )
987 300 CONTINUE
988 ELSE
989 DO 310 jr = jc, min( n, jc+uub )
990 a( jr+irow, jc ) = dconjg( a( jc-iskew*jr+
991 $ ioffg, jr ) )
992 310 CONTINUE
993 END IF
994 320 CONTINUE
995 IF( ipack.EQ.5 ) THEN
996 DO 340 jc = n - uub + 1, n
997 DO 330 jr = n + 2 - jc, uub + 1
998 a( jr, jc ) = czero
999 330 CONTINUE
1000 340 CONTINUE
1001 END IF
1002 IF( ipackg.EQ.6 ) THEN
1003 ipackg = ipack
1004 ELSE
1005 ipackg = 0
1006 END IF
1007 END IF
1008 ELSE
1009*
1010* Bottom-Up -- Generate Lower triangle only
1011*
1012 IF( ipack.GE.5 ) THEN
1013 ipackg = 5
1014 IF( ipack.EQ.6 )
1015 $ ioffg = 1
1016 ELSE
1017 ipackg = 2
1018 END IF
1019*
1020 DO 350 j = 1, mnmin
1021 a( ( 1-iskew )*j+ioffg, j ) = dcmplx( d( j ) )
1022 350 CONTINUE
1023*
1024 DO 380 k = 1, uub
1025 DO 370 jc = n - 1, 1, -1
1026 il = min( n+1-jc, k+2 )
1027 extra = czero
1028 ztemp = a( 1+( 1-iskew )*jc+ioffg, jc )
1029 angle = twopi*dlarnd( 1, iseed )
1030 c = cos( angle )*zlarnd( 5, iseed )
1031 s = sin( angle )*zlarnd( 5, iseed )
1032 IF( csym ) THEN
1033 ct = c
1034 st = s
1035 ELSE
1036 ztemp = dconjg( ztemp )
1037 ct = dconjg( c )
1038 st = dconjg( s )
1039 END IF
1040 CALL zlarot( .false., .true., n-jc.GT.k, il, c,
1041 $ s,
1042 $ a( ( 1-iskew )*jc+ioffg, jc ), ilda,
1043 $ ztemp, extra )
1044 icol = max( 1, jc-k+1 )
1045 CALL zlarot( .true., .false., .true., jc+2-icol,
1046 $ ct, st, a( jc-iskew*icol+ioffg,
1047 $ icol ), ilda, dummy, ztemp )
1048*
1049* Chase EXTRA back down the matrix
1050*
1051 icol = jc
1052 DO 360 jch = jc + k, n - 1, k
1053 CALL zlartg( a( jch-iskew*icol+ioffg, icol ),
1054 $ extra, realc, s, dummy )
1055 dummy = zlarnd( 5, iseed )
1056 c = realc*dummy
1057 s = s*dummy
1058 ztemp = a( 1+( 1-iskew )*jch+ioffg, jch )
1059 IF( csym ) THEN
1060 ct = c
1061 st = s
1062 ELSE
1063 ztemp = dconjg( ztemp )
1064 ct = dconjg( c )
1065 st = dconjg( s )
1066 END IF
1067 CALL zlarot( .true., .true., .true., k+2, c,
1068 $ s,
1069 $ a( jch-iskew*icol+ioffg, icol ),
1070 $ ilda, extra, ztemp )
1071 il = min( n+1-jch, k+2 )
1072 extra = czero
1073 CALL zlarot( .false., .true., n-jch.GT.k, il,
1074 $ ct, st, a( ( 1-iskew )*jch+ioffg,
1075 $ jch ), ilda, ztemp, extra )
1076 icol = jch
1077 360 CONTINUE
1078 370 CONTINUE
1079 380 CONTINUE
1080*
1081* If we need upper triangle, copy from lower. Note that
1082* the order of copying is chosen to work for 'b' -> 'q'
1083*
1084 IF( ipack.NE.ipackg .AND. ipack.NE.4 ) THEN
1085 DO 410 jc = n, 1, -1
1086 irow = ioffst - iskew*jc
1087 IF( csym ) THEN
1088 DO 390 jr = jc, max( 1, jc-uub ), -1
1089 a( jr+irow, jc ) = a( jc-iskew*jr+ioffg, jr )
1090 390 CONTINUE
1091 ELSE
1092 DO 400 jr = jc, max( 1, jc-uub ), -1
1093 a( jr+irow, jc ) = dconjg( a( jc-iskew*jr+
1094 $ ioffg, jr ) )
1095 400 CONTINUE
1096 END IF
1097 410 CONTINUE
1098 IF( ipack.EQ.6 ) THEN
1099 DO 430 jc = 1, uub
1100 DO 420 jr = 1, uub + 1 - jc
1101 a( jr, jc ) = czero
1102 420 CONTINUE
1103 430 CONTINUE
1104 END IF
1105 IF( ipackg.EQ.5 ) THEN
1106 ipackg = ipack
1107 ELSE
1108 ipackg = 0
1109 END IF
1110 END IF
1111 END IF
1112*
1113* Ensure that the diagonal is real if Hermitian
1114*
1115 IF( .NOT.csym ) THEN
1116 DO 440 jc = 1, n
1117 irow = ioffst + ( 1-iskew )*jc
1118 a( irow, jc ) = dcmplx( dble( a( irow, jc ) ) )
1119 440 CONTINUE
1120 END IF
1121*
1122 END IF
1123*
1124 ELSE
1125*
1126* 4) Generate Banded Matrix by first
1127* Rotating by random Unitary matrices,
1128* then reducing the bandwidth using Householder
1129* transformations.
1130*
1131* Note: we should get here only if LDA .ge. N
1132*
1133 IF( isym.EQ.1 ) THEN
1134*
1135* Non-symmetric -- A = U D V
1136*
1137 CALL zlagge( mr, nc, llb, uub, d, a, lda, iseed, work,
1138 $ iinfo )
1139 ELSE
1140*
1141* Symmetric -- A = U D U' or
1142* Hermitian -- A = U D U*
1143*
1144 IF( csym ) THEN
1145 CALL zlagsy( m, llb, d, a, lda, iseed, work, iinfo )
1146 ELSE
1147 CALL zlaghe( m, llb, d, a, lda, iseed, work, iinfo )
1148 END IF
1149 END IF
1150*
1151 IF( iinfo.NE.0 ) THEN
1152 info = 3
1153 RETURN
1154 END IF
1155 END IF
1156*
1157* 5) Pack the matrix
1158*
1159 IF( ipack.NE.ipackg ) THEN
1160 IF( ipack.EQ.1 ) THEN
1161*
1162* 'U' -- Upper triangular, not packed
1163*
1164 DO 460 j = 1, m
1165 DO 450 i = j + 1, m
1166 a( i, j ) = czero
1167 450 CONTINUE
1168 460 CONTINUE
1169*
1170 ELSE IF( ipack.EQ.2 ) THEN
1171*
1172* 'L' -- Lower triangular, not packed
1173*
1174 DO 480 j = 2, m
1175 DO 470 i = 1, j - 1
1176 a( i, j ) = czero
1177 470 CONTINUE
1178 480 CONTINUE
1179*
1180 ELSE IF( ipack.EQ.3 ) THEN
1181*
1182* 'C' -- Upper triangle packed Columnwise.
1183*
1184 icol = 1
1185 irow = 0
1186 DO 500 j = 1, m
1187 DO 490 i = 1, j
1188 irow = irow + 1
1189 IF( irow.GT.lda ) THEN
1190 irow = 1
1191 icol = icol + 1
1192 END IF
1193 a( irow, icol ) = a( i, j )
1194 490 CONTINUE
1195 500 CONTINUE
1196*
1197 ELSE IF( ipack.EQ.4 ) THEN
1198*
1199* 'R' -- Lower triangle packed Columnwise.
1200*
1201 icol = 1
1202 irow = 0
1203 DO 520 j = 1, m
1204 DO 510 i = j, m
1205 irow = irow + 1
1206 IF( irow.GT.lda ) THEN
1207 irow = 1
1208 icol = icol + 1
1209 END IF
1210 a( irow, icol ) = a( i, j )
1211 510 CONTINUE
1212 520 CONTINUE
1213*
1214 ELSE IF( ipack.GE.5 ) THEN
1215*
1216* 'B' -- The lower triangle is packed as a band matrix.
1217* 'Q' -- The upper triangle is packed as a band matrix.
1218* 'Z' -- The whole matrix is packed as a band matrix.
1219*
1220 IF( ipack.EQ.5 )
1221 $ uub = 0
1222 IF( ipack.EQ.6 )
1223 $ llb = 0
1224*
1225 DO 540 j = 1, uub
1226 DO 530 i = min( j+llb, m ), 1, -1
1227 a( i-j+uub+1, j ) = a( i, j )
1228 530 CONTINUE
1229 540 CONTINUE
1230*
1231 DO 560 j = uub + 2, n
1232 DO 550 i = j - uub, min( j+llb, m )
1233 a( i-j+uub+1, j ) = a( i, j )
1234 550 CONTINUE
1235 560 CONTINUE
1236 END IF
1237*
1238* If packed, zero out extraneous elements.
1239*
1240* Symmetric/Triangular Packed --
1241* zero out everything after A(IROW,ICOL)
1242*
1243 IF( ipack.EQ.3 .OR. ipack.EQ.4 ) THEN
1244 DO 580 jc = icol, m
1245 DO 570 jr = irow + 1, lda
1246 a( jr, jc ) = czero
1247 570 CONTINUE
1248 irow = 0
1249 580 CONTINUE
1250*
1251 ELSE IF( ipack.GE.5 ) THEN
1252*
1253* Packed Band --
1254* 1st row is now in A( UUB+2-j, j), zero above it
1255* m-th row is now in A( M+UUB-j,j), zero below it
1256* last non-zero diagonal is now in A( UUB+LLB+1,j ),
1257* zero below it, too.
1258*
1259 ir1 = uub + llb + 2
1260 ir2 = uub + m + 2
1261 DO 610 jc = 1, n
1262 DO 590 jr = 1, uub + 1 - jc
1263 a( jr, jc ) = czero
1264 590 CONTINUE
1265 DO 600 jr = max( 1, min( ir1, ir2-jc ) ), lda
1266 a( jr, jc ) = czero
1267 600 CONTINUE
1268 610 CONTINUE
1269 END IF
1270 END IF
1271*
1272 RETURN
1273*
1274* End of ZLATMT
1275*
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:104
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:230
Here is the call graph for this function:
Here is the caller graph for this function: