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

◆ dlatmt()

subroutine dlatmt ( 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,
double precision, dimension( lda, * )  A,
integer  LDA,
double precision, dimension( * )  WORK,
integer  INFO 
)

DLATMT

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

    DLATMT 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 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 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 symmetric)
         zero out lower half (if symmetric)
         store the upper half columnwise (if symmetric or upper
               triangular)
         store the lower half columnwise (if symmetric or lower
               triangular)
         store the lower triangle in banded format (if symmetric
               or lower triangular)
         store the upper triangle in banded format (if symmetric
               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. 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 DLATMT
           to continue the same random number sequence.
           Changed on exit.
[in]SYM
          SYM is CHARACTER*1
           If SYM='S' or 'H', the generated matrix is symmetric, with
             eigenvalues specified by D, COND, MODE, and DMAX; they
             may be positive, negative, or zero.
           If SYM='P', the generated matrix is symmetric, 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.
           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='S' or '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.
           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.
           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)
           'L' => zero out all superdiagonal entries (if symmetric)
           'C' => store the upper triangle columnwise
                  (only if the matrix is symmetric or upper triangular)
           'R' => store the lower triangle columnwise
                  (only if the matrix is symmetric or lower triangular)
           'B' => store the lower triangle in band storage scheme
                  (only if matrix symmetric or lower triangular)
           'Q' => store the upper triangle in band storage scheme
                  (only if matrix symmetric 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 or TB     - use 'B' or 'Q'
           PP, SP or TP     - use 'C' or 'R'

           If two calls to DLATMT differ only in the PACK parameter,
           they will generate mathematically equivalent matrices.
           Not modified.
[in,out]A
          A is DOUBLE PRECISION 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 DOUBLE PRECISION 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='S' or 'H' and KU 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 DLAGGE or DLAGSY
Author
Univ. of Tennessee
Univ. of California Berkeley
Univ. of Colorado Denver
NAG Ltd.

Definition at line 329 of file dlatmt.f.

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