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

◆ dlatms()

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

DLATMS

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

    DLATMS 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 DLATMS
           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:N)=1.0/COND
           MODE = 2 sets D(1:N-1)=1 and D(N)=1.0/COND
           MODE = 3 sets D(I)=COND**(-(I-1)/(N-1))
           MODE = 4 sets D(i)=1 - (i-1)/(N-1)*(1 - 1/COND)
           MODE = 5 sets D to random numbers in the range
                    ( 1/COND , 1 ) such that their logarithms
                    are uniformly distributed.
           MODE = 6 set D to random numbers from same distribution
                    as the rest of the matrix.
           MODE < 0 has the same meaning as ABS(MODE), except that
              the order of the elements of D is reversed.
           Thus if MODE is positive, D has entries ranging from
              1 to 1/COND, if negative, from 1/COND to 1,
           If SYM='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]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 DLATMS 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 DLATM1
            2  => Cannot scale to DMAX (max. sing. value is 0)
            3  => Error return from DLAGGE or SLAGSY
Author
Univ. of Tennessee
Univ. of California Berkeley
Univ. of Colorado Denver
NAG Ltd.

Definition at line 319 of file dlatms.f.

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