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

◆ clatms()

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

CLATMS

Purpose:
!>
!>    CLATMS generates random matrices with specified singular values
!>    (or hermitian with specified eigenvalues)
!>    for testing LAPACK programs.
!>
!>    CLATMS operates by applying the following sequence of
!>    operations:
!>
!>      Set the diagonal to D, where D may be input or
!>         computed according to MODE, COND, DMAX, and SYM
!>         as described below.
!>
!>      Generate a matrix with the appropriate band structure, by one
!>         of two methods:
!>
!>      Method A:
!>          Generate a dense M x N matrix by multiplying D on the left
!>              and the right by random unitary matrices, then:
!>
!>          Reduce the bandwidth according to KL and KU, using
!>              Householder transformations.
!>
!>      Method B:
!>          Convert the bandwidth-0 (i.e., diagonal) matrix to a
!>              bandwidth-1 matrix using Givens rotations, 
!>              out-of-band elements back, much as in QR; then convert
!>              the bandwidth-1 to a bandwidth-2 matrix, etc.  Note
!>              that for reasonably small bandwidths (relative to M and
!>              N) this requires less storage, as a dense matrix is not
!>              generated.  Also, for hermitian or symmetric matrices,
!>              only one triangle is generated.
!>
!>      Method A is chosen if the bandwidth is a large fraction of the
!>          order of the matrix, and LDA is at least M (so a dense
!>          matrix can be stored.)  Method B is chosen if the bandwidth
!>          is small (< 1/2 N for hermitian or symmetric, < .3 N+M for
!>          non-symmetric), or LDA is less than M and not less than the
!>          bandwidth.
!>
!>      Pack the matrix if desired. Options specified by PACK are:
!>         no packing
!>         zero out upper half (if hermitian)
!>         zero out lower half (if hermitian)
!>         store the upper half columnwise (if hermitian or upper
!>               triangular)
!>         store the lower half columnwise (if hermitian or lower
!>               triangular)
!>         store the lower triangle in banded format (if hermitian or
!>               lower triangular)
!>         store the upper triangle in banded format (if hermitian or
!>               upper triangular)
!>         store the entire matrix in banded format
!>      If Method B is chosen, and band format is specified, then the
!>         matrix will be generated in the band format, so no repacking
!>         will be necessary.
!> 
Parameters
[in]M
!>          M is INTEGER
!>           The number of rows of A. Not modified.
!> 
[in]N
!>          N is INTEGER
!>           The number of columns of A. N must equal M if the matrix
!>           is symmetric or hermitian (i.e., if SYM is not 'N')
!>           Not modified.
!> 
[in]DIST
!>          DIST is CHARACTER*1
!>           On entry, DIST specifies the type of distribution to be used
!>           to generate the random eigen-/singular values.
!>           'U' => UNIFORM( 0, 1 )  ( 'U' for uniform )
!>           'S' => UNIFORM( -1, 1 ) ( 'S' for symmetric )
!>           'N' => NORMAL( 0, 1 )   ( 'N' for normal )
!>           Not modified.
!> 
[in,out]ISEED
!>          ISEED is INTEGER array, dimension ( 4 )
!>           On entry ISEED specifies the seed of the random number
!>           generator. They should lie between 0 and 4095 inclusive,
!>           and ISEED(4) should be odd. The random number generator
!>           uses a linear congruential sequence limited to small
!>           integers, and so should produce machine independent
!>           random numbers. The values of ISEED are changed on
!>           exit, and can be used in the next call to CLATMS
!>           to continue the same random number sequence.
!>           Changed on exit.
!> 
[in]SYM
!>          SYM is CHARACTER*1
!>           If SYM='H', the generated matrix is hermitian, with
!>             eigenvalues specified by D, COND, MODE, and DMAX; they
!>             may be positive, negative, or zero.
!>           If SYM='P', the generated matrix is hermitian, with
!>             eigenvalues (= singular values) specified by D, COND,
!>             MODE, and DMAX; they will not be negative.
!>           If SYM='N', the generated matrix is nonsymmetric, with
!>             singular values specified by D, COND, MODE, and DMAX;
!>             they will not be negative.
!>           If SYM='S', the generated matrix is (complex) symmetric,
!>             with singular values specified by D, COND, MODE, and
!>             DMAX; they will not be negative.
!>           Not modified.
!> 
[in,out]D
!>          D is REAL 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='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 REAL
!>           On entry, this is used as described under MODE above.
!>           If used, it must be >= 1. Not modified.
!> 
[in]DMAX
!>          DMAX is REAL
!>           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 or hermitian.
!>           Not modified.
!> 
[in]KU
!>          KU is INTEGER
!>           This specifies the upper bandwidth of the  matrix. For
!>           example, KU=0 implies lower triangular, KU=1 implies lower
!>           Hessenberg, and KU being at least N-1 means that the matrix
!>           has full upper bandwidth.  KL must equal KU if the matrix
!>           is symmetric or hermitian.
!>           Not modified.
!> 
[in]PACK
!>          PACK is CHARACTER*1
!>           This specifies packing of matrix as follows:
!>           'N' => no packing
!>           'U' => zero out all subdiagonal entries (if symmetric
!>                  or hermitian)
!>           'L' => zero out all superdiagonal entries (if symmetric
!>                  or hermitian)
!>           'C' => store the upper triangle columnwise (only if the
!>                  matrix is symmetric, hermitian, or upper triangular)
!>           'R' => store the lower triangle columnwise (only if the
!>                  matrix is symmetric, hermitian, or lower triangular)
!>           'B' => store the lower triangle in band storage scheme
!>                  (only if the matrix is symmetric, hermitian, or
!>                  lower triangular)
!>           'Q' => store the upper triangle in band storage scheme
!>                  (only if the matrix is symmetric, hermitian, or
!>                  upper triangular)
!>           'Z' => store the entire matrix in band storage scheme
!>                      (pivoting can be provided for by using this
!>                      option to store A in the trailing rows of
!>                      the allocated storage)
!>
!>           Using these options, the various LAPACK packed and banded
!>           storage schemes can be obtained:
!>           GB                    - use 'Z'
!>           PB, SB, HB, or TB     - use 'B' or 'Q'
!>           PP, SP, HB, or TP     - use 'C' or 'R'
!>
!>           If two calls to CLATMS differ only in the PACK parameter,
!>           they will generate mathematically equivalent matrices.
!>           Not modified.
!> 
[in,out]A
!>          A is COMPLEX array, dimension ( LDA, N )
!>           On exit A is the desired test matrix.  A is first generated
!>           in full (unpacked) form, and then packed, if so specified
!>           by PACK.  Thus, the first M elements of the first N
!>           columns will always be modified.  If PACK specifies a
!>           packed or banded storage scheme, all LDA elements of the
!>           first N columns will be modified; the elements of the
!>           array which do not correspond to elements of the generated
!>           matrix are set to zero.
!>           Modified.
!> 
[in]LDA
!>          LDA is INTEGER
!>           LDA specifies the first dimension of A as declared in the
!>           calling program.  If PACK='N', 'U', 'L', 'C', or 'R', then
!>           LDA must be at least M.  If PACK='B' or 'Q', then LDA must
!>           be at least MIN( KL, M-1) (which is equal to MIN(KU,N-1)).
!>           If PACK='Z', LDA must be large enough to hold the packed
!>           array: MIN( KU, N-1) + MIN( KL, M-1) + 1.
!>           Not modified.
!> 
[out]WORK
!>          WORK is COMPLEX array, dimension ( 3*MAX( N, M ) )
!>           Workspace.
!>           Modified.
!> 
[out]INFO
!>          INFO is INTEGER
!>           Error code.  On exit, INFO will be set to one of the
!>           following values:
!>             0 => normal return
!>            -1 => M negative or unequal to N and SYM='S', 'H', or 'P'
!>            -2 => N negative
!>            -3 => DIST illegal string
!>            -5 => SYM illegal string
!>            -7 => MODE not in range -6 to 6
!>            -8 => COND less than 1.0, and MODE neither -6, 0 nor 6
!>           -10 => KL negative
!>           -11 => KU negative, or SYM is not 'N' and KU is not equal to
!>                  KL
!>           -12 => PACK illegal string, or PACK='U' or 'L', and SYM='N';
!>                  or PACK='C' or 'Q' and SYM='N' and KL is not zero;
!>                  or PACK='R' or 'B' and SYM='N' and KU is not zero;
!>                  or PACK='U', 'L', 'C', 'R', 'B', or 'Q', and M is not
!>                  N.
!>           -14 => LDA is less than M, or PACK='Z' and LDA is less than
!>                  MIN(KU,N-1) + MIN(KL,M-1) + 1.
!>            1  => Error return from SLATM1
!>            2  => Cannot scale to DMAX (max. sing. value is 0)
!>            3  => Error return from CLAGGE, CLAGHE or CLAGSY
!> 
Author
Univ. of Tennessee
Univ. of California Berkeley
Univ. of Colorado Denver
NAG Ltd.

Definition at line 330 of file clatms.f.

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