LAPACK  3.8.0
LAPACK: Linear Algebra PACKage

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

Definition at line 333 of file dlatmt.f.

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