summaryrefslogtreecommitdiff
path: root/TESTING/MATGEN/dlatmt.f
diff options
context:
space:
mode:
Diffstat (limited to 'TESTING/MATGEN/dlatmt.f')
-rw-r--r--TESTING/MATGEN/dlatmt.f1047
1 files changed, 1047 insertions, 0 deletions
diff --git a/TESTING/MATGEN/dlatmt.f b/TESTING/MATGEN/dlatmt.f
new file mode 100644
index 00000000..2ff32554
--- /dev/null
+++ b/TESTING/MATGEN/dlatmt.f
@@ -0,0 +1,1047 @@
+ SUBROUTINE DLATMT( M, N, DIST, ISEED, SYM, D, MODE, COND, DMAX,
+ $ RANK, KL, KU, PACK, A, LDA, WORK, INFO )
+*
+* -- LAPACK test routine (version 3.1) --
+* Craig Lucas, University of Manchester / NAG Ltd.
+* October, 2008
+*
+* .. Scalar Arguments ..
+ DOUBLE PRECISION COND, DMAX
+ INTEGER INFO, KL, KU, LDA, M, MODE, N, RANK
+ CHARACTER DIST, PACK, SYM
+* ..
+* .. Array Arguments ..
+ DOUBLE PRECISION A( LDA, * ), D( * ), WORK( * )
+ INTEGER ISEED( 4 )
+* ..
+*
+* 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.
+*
+* Arguments
+* =========
+*
+* M - INTEGER
+* The number of rows of A. Not modified.
+*
+* N - INTEGER
+* The number of columns of A. Not modified.
+*
+* DIST - 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.
+*
+* ISEED - 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.
+*
+* SYM - 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.
+*
+* D - 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.
+*
+* MODE - 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.
+*
+* COND - DOUBLE PRECISION
+* On entry, this is used as described under MODE above.
+* If used, it must be >= 1. Not modified.
+*
+* DMAX - 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.
+*
+* RANK - INTEGER
+* The rank of matrix to be generated for modes 1,2,3 only.
+* D( RANK+1:N ) = 0.
+* Not modified.
+*
+* KL - 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.
+*
+* KU - 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.
+*
+* PACK - 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.
+*
+* A - 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.
+*
+* LDA - 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.
+*
+* WORK - DOUBLE PRECISION array, dimension ( 3*MAX( N , M ) )
+* Workspace.
+* Modified.
+*
+* INFO - 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
+*
+* =====================================================================
+*
+* .. Parameters ..
+ DOUBLE PRECISION ZERO
+ PARAMETER ( ZERO = 0.0D0 )
+ DOUBLE PRECISION ONE
+ PARAMETER ( ONE = 1.0D0 )
+ DOUBLE PRECISION TWOPI
+ PARAMETER ( TWOPI = 6.2831853071795864769252867663D+0 )
+* ..
+* .. Local Scalars ..
+ DOUBLE PRECISION ALPHA, ANGLE, C, DUMMY, EXTRA, S, TEMP
+ INTEGER I, IC, ICOL, IDIST, IENDCH, IINFO, IL, ILDA,
+ $ IOFFG, IOFFST, IPACK, IPACKG, IR, IR1, IR2,
+ $ IROW, IRSIGN, ISKEW, ISYM, ISYMPK, J, JC, JCH,
+ $ JKL, JKU, JR, K, LLB, MINLDA, MNMIN, MR, NC,
+ $ UUB
+ LOGICAL GIVENS, ILEXTR, ILTEMP, TOPDWN
+* ..
+* .. External Functions ..
+ DOUBLE PRECISION DLARND
+ LOGICAL LSAME
+ EXTERNAL DLARND, LSAME
+* ..
+* .. External Subroutines ..
+ EXTERNAL DLATM7, DCOPY, DLAGGE, DLAGSY, DLAROT,
+ $ DLARTG, DLASET, DSCAL, XERBLA
+* ..
+* .. Intrinsic Functions ..
+ INTRINSIC ABS, COS, DBLE, MAX, MIN, MOD, SIN
+* ..
+* .. Executable Statements ..
+*
+* 1) Decode and Test the input parameters.
+* Initialize flags & seed.
+*
+ INFO = 0
+*
+* Quick return if possible
+*
+ IF( M.EQ.0 .OR. N.EQ.0 )
+ $ RETURN
+*
+* Decode DIST
+*
+ IF( LSAME( DIST, 'U' ) ) THEN
+ IDIST = 1
+ ELSE IF( LSAME( DIST, 'S' ) ) THEN
+ IDIST = 2
+ ELSE IF( LSAME( DIST, 'N' ) ) THEN
+ IDIST = 3
+ ELSE
+ IDIST = -1
+ END IF
+*
+* Decode SYM
+*
+ IF( LSAME( SYM, 'N' ) ) THEN
+ ISYM = 1
+ IRSIGN = 0
+ ELSE IF( LSAME( SYM, 'P' ) ) THEN
+ ISYM = 2
+ IRSIGN = 0
+ ELSE IF( LSAME( SYM, 'S' ) ) THEN
+ ISYM = 2
+ IRSIGN = 1
+ ELSE IF( LSAME( SYM, 'H' ) ) THEN
+ ISYM = 2
+ IRSIGN = 1
+ ELSE
+ ISYM = -1
+ END IF
+*
+* Decode PACK
+*
+ ISYMPK = 0
+ IF( LSAME( PACK, 'N' ) ) THEN
+ IPACK = 0
+ ELSE IF( LSAME( PACK, 'U' ) ) THEN
+ IPACK = 1
+ ISYMPK = 1
+ ELSE IF( LSAME( PACK, 'L' ) ) THEN
+ IPACK = 2
+ ISYMPK = 1
+ ELSE IF( LSAME( PACK, 'C' ) ) THEN
+ IPACK = 3
+ ISYMPK = 2
+ ELSE IF( LSAME( PACK, 'R' ) ) THEN
+ IPACK = 4
+ ISYMPK = 3
+ ELSE IF( LSAME( PACK, 'B' ) ) THEN
+ IPACK = 5
+ ISYMPK = 3
+ ELSE IF( LSAME( PACK, 'Q' ) ) THEN
+ IPACK = 6
+ ISYMPK = 2
+ ELSE IF( LSAME( PACK, 'Z' ) ) THEN
+ IPACK = 7
+ ELSE
+ IPACK = -1
+ END IF
+*
+* Set certain internal parameters
+*
+ MNMIN = MIN( M, N )
+ LLB = MIN( KL, M-1 )
+ UUB = MIN( KU, N-1 )
+ MR = MIN( M, N+LLB )
+ NC = MIN( N, M+UUB )
+*
+ IF( IPACK.EQ.5 .OR. IPACK.EQ.6 ) THEN
+ MINLDA = UUB + 1
+ ELSE IF( IPACK.EQ.7 ) THEN
+ MINLDA = LLB + UUB + 1
+ ELSE
+ MINLDA = M
+ END IF
+*
+* Use Givens rotation method if bandwidth small enough,
+* or if LDA is too small to store the matrix unpacked.
+*
+ GIVENS = .FALSE.
+ IF( ISYM.EQ.1 ) THEN
+ IF( DBLE( LLB+UUB ).LT.0.3D0*DBLE( MAX( 1, MR+NC ) ) )
+ $ GIVENS = .TRUE.
+ ELSE
+ IF( 2*LLB.LT.M )
+ $ GIVENS = .TRUE.
+ END IF
+ IF( LDA.LT.M .AND. LDA.GE.MINLDA )
+ $ GIVENS = .TRUE.
+*
+* Set INFO if an error
+*
+ IF( M.LT.0 ) THEN
+ INFO = -1
+ ELSE IF( M.NE.N .AND. ISYM.NE.1 ) THEN
+ INFO = -1
+ ELSE IF( N.LT.0 ) THEN
+ INFO = -2
+ ELSE IF( IDIST.EQ.-1 ) THEN
+ INFO = -3
+ ELSE IF( ISYM.EQ.-1 ) THEN
+ INFO = -5
+ ELSE IF( ABS( MODE ).GT.6 ) THEN
+ INFO = -7
+ ELSE IF( ( MODE.NE.0 .AND. ABS( MODE ).NE.6 ) .AND. COND.LT.ONE )
+ $ THEN
+ INFO = -8
+ ELSE IF( KL.LT.0 ) THEN
+ INFO = -10
+ ELSE IF( KU.LT.0 .OR. ( ISYM.NE.1 .AND. KL.NE.KU ) ) THEN
+ INFO = -11
+ ELSE IF( IPACK.EQ.-1 .OR. ( ISYMPK.EQ.1 .AND. ISYM.EQ.1 ) .OR.
+ $ ( ISYMPK.EQ.2 .AND. ISYM.EQ.1 .AND. KL.GT.0 ) .OR.
+ $ ( ISYMPK.EQ.3 .AND. ISYM.EQ.1 .AND. KU.GT.0 ) .OR.
+ $ ( ISYMPK.NE.0 .AND. M.NE.N ) ) THEN
+ INFO = -12
+ ELSE IF( LDA.LT.MAX( 1, MINLDA ) ) THEN
+ INFO = -14
+ END IF
+*
+ IF( INFO.NE.0 ) THEN
+ CALL XERBLA( 'DLATMT', -INFO )
+ RETURN
+ END IF
+*
+* Initialize random number generator
+*
+ DO 100 I = 1, 4
+ ISEED( I ) = MOD( ABS( ISEED( I ) ), 4096 )
+ 100 CONTINUE
+*
+ IF( MOD( ISEED( 4 ), 2 ).NE.1 )
+ $ ISEED( 4 ) = ISEED( 4 ) + 1
+*
+* 2) Set up D if indicated.
+*
+* Compute D according to COND and MODE
+*
+ CALL DLATM7( MODE, COND, IRSIGN, IDIST, ISEED, D, MNMIN, RANK,
+ $ IINFO )
+ IF( IINFO.NE.0 ) THEN
+ INFO = 1
+ RETURN
+ END IF
+*
+* Choose Top-Down if D is (apparently) increasing,
+* Bottom-Up if D is (apparently) decreasing.
+*
+ IF( ABS( D( 1 ) ).LE.ABS( D( RANK ) ) ) THEN
+ TOPDWN = .TRUE.
+ ELSE
+ TOPDWN = .FALSE.
+ END IF
+*
+ IF( MODE.NE.0 .AND. ABS( MODE ).NE.6 ) THEN
+*
+* Scale by DMAX
+*
+ TEMP = ABS( D( 1 ) )
+ DO 110 I = 2, RANK
+ TEMP = MAX( TEMP, ABS( D( I ) ) )
+ 110 CONTINUE
+*
+ IF( TEMP.GT.ZERO ) THEN
+ ALPHA = DMAX / TEMP
+ ELSE
+ INFO = 2
+ RETURN
+ END IF
+*
+ CALL DSCAL( RANK, ALPHA, D, 1 )
+*
+ END IF
+*
+* 3) Generate Banded Matrix using Givens rotations.
+* Also the special case of UUB=LLB=0
+*
+* Compute Addressing constants to cover all
+* storage formats. Whether GE, SY, GB, or SB,
+* upper or lower triangle or both,
+* the (i,j)-th element is in
+* A( i - ISKEW*j + IOFFST, j )
+*
+ IF( IPACK.GT.4 ) THEN
+ ILDA = LDA - 1
+ ISKEW = 1
+ IF( IPACK.GT.5 ) THEN
+ IOFFST = UUB + 1
+ ELSE
+ IOFFST = 1
+ END IF
+ ELSE
+ ILDA = LDA
+ ISKEW = 0
+ IOFFST = 0
+ END IF
+*
+* IPACKG is the format that the matrix is generated in. If this is
+* different from IPACK, then the matrix must be repacked at the
+* end. It also signals how to compute the norm, for scaling.
+*
+ IPACKG = 0
+ CALL DLASET( 'Full', LDA, N, ZERO, ZERO, A, LDA )
+*
+* Diagonal Matrix -- We are done, unless it
+* is to be stored SP/PP/TP (PACK='R' or 'C')
+*
+ IF( LLB.EQ.0 .AND. UUB.EQ.0 ) THEN
+ CALL DCOPY( MNMIN, D, 1, A( 1-ISKEW+IOFFST, 1 ), ILDA+1 )
+ IF( IPACK.LE.2 .OR. IPACK.GE.5 )
+ $ IPACKG = IPACK
+*
+ ELSE IF( GIVENS ) THEN
+*
+* Check whether to use Givens rotations,
+* Householder transformations, or nothing.
+*
+ IF( ISYM.EQ.1 ) THEN
+*
+* Non-symmetric -- A = U D V
+*
+ IF( IPACK.GT.4 ) THEN
+ IPACKG = IPACK
+ ELSE
+ IPACKG = 0
+ END IF
+*
+ CALL DCOPY( MNMIN, D, 1, A( 1-ISKEW+IOFFST, 1 ), ILDA+1 )
+*
+ IF( TOPDWN ) THEN
+ JKL = 0
+ DO 140 JKU = 1, UUB
+*
+* Transform from bandwidth JKL, JKU-1 to JKL, JKU
+*
+* Last row actually rotated is M
+* Last column actually rotated is MIN( M+JKU, N )
+*
+ DO 130 JR = 1, MIN( M+JKU, N ) + JKL - 1
+ EXTRA = ZERO
+ ANGLE = TWOPI*DLARND( 1, ISEED )
+ C = COS( ANGLE )
+ S = SIN( ANGLE )
+ ICOL = MAX( 1, JR-JKL )
+ IF( JR.LT.M ) THEN
+ IL = MIN( N, JR+JKU ) + 1 - ICOL
+ CALL DLAROT( .TRUE., JR.GT.JKL, .FALSE., IL, C,
+ $ S, A( JR-ISKEW*ICOL+IOFFST, ICOL ),
+ $ ILDA, EXTRA, DUMMY )
+ END IF
+*
+* Chase "EXTRA" back up
+*
+ IR = JR
+ IC = ICOL
+ DO 120 JCH = JR - JKL, 1, -JKL - JKU
+ IF( IR.LT.M ) THEN
+ CALL DLARTG( A( IR+1-ISKEW*( IC+1 )+IOFFST,
+ $ IC+1 ), EXTRA, C, S, DUMMY )
+ END IF
+ IROW = MAX( 1, JCH-JKU )
+ IL = IR + 2 - IROW
+ TEMP = ZERO
+ ILTEMP = JCH.GT.JKU
+ CALL DLAROT( .FALSE., ILTEMP, .TRUE., IL, C, -S,
+ $ A( IROW-ISKEW*IC+IOFFST, IC ),
+ $ ILDA, TEMP, EXTRA )
+ IF( ILTEMP ) THEN
+ CALL DLARTG( A( IROW+1-ISKEW*( IC+1 )+IOFFST,
+ $ IC+1 ), TEMP, C, S, DUMMY )
+ ICOL = MAX( 1, JCH-JKU-JKL )
+ IL = IC + 2 - ICOL
+ EXTRA = ZERO
+ CALL DLAROT( .TRUE., JCH.GT.JKU+JKL, .TRUE.,
+ $ IL, C, -S, A( IROW-ISKEW*ICOL+
+ $ IOFFST, ICOL ), ILDA, EXTRA,
+ $ TEMP )
+ IC = ICOL
+ IR = IROW
+ END IF
+ 120 CONTINUE
+ 130 CONTINUE
+ 140 CONTINUE
+*
+ JKU = UUB
+ DO 170 JKL = 1, LLB
+*
+* Transform from bandwidth JKL-1, JKU to JKL, JKU
+*
+ DO 160 JC = 1, MIN( N+JKL, M ) + JKU - 1
+ EXTRA = ZERO
+ ANGLE = TWOPI*DLARND( 1, ISEED )
+ C = COS( ANGLE )
+ S = SIN( ANGLE )
+ IROW = MAX( 1, JC-JKU )
+ IF( JC.LT.N ) THEN
+ IL = MIN( M, JC+JKL ) + 1 - IROW
+ CALL DLAROT( .FALSE., JC.GT.JKU, .FALSE., IL, C,
+ $ S, A( IROW-ISKEW*JC+IOFFST, JC ),
+ $ ILDA, EXTRA, DUMMY )
+ END IF
+*
+* Chase "EXTRA" back up
+*
+ IC = JC
+ IR = IROW
+ DO 150 JCH = JC - JKU, 1, -JKL - JKU
+ IF( IC.LT.N ) THEN
+ CALL DLARTG( A( IR+1-ISKEW*( IC+1 )+IOFFST,
+ $ IC+1 ), EXTRA, C, S, DUMMY )
+ END IF
+ ICOL = MAX( 1, JCH-JKL )
+ IL = IC + 2 - ICOL
+ TEMP = ZERO
+ ILTEMP = JCH.GT.JKL
+ CALL DLAROT( .TRUE., ILTEMP, .TRUE., IL, C, -S,
+ $ A( IR-ISKEW*ICOL+IOFFST, ICOL ),
+ $ ILDA, TEMP, EXTRA )
+ IF( ILTEMP ) THEN
+ CALL DLARTG( A( IR+1-ISKEW*( ICOL+1 )+IOFFST,
+ $ ICOL+1 ), TEMP, C, S, DUMMY )
+ IROW = MAX( 1, JCH-JKL-JKU )
+ IL = IR + 2 - IROW
+ EXTRA = ZERO
+ CALL DLAROT( .FALSE., JCH.GT.JKL+JKU, .TRUE.,
+ $ IL, C, -S, A( IROW-ISKEW*ICOL+
+ $ IOFFST, ICOL ), ILDA, EXTRA,
+ $ TEMP )
+ IC = ICOL
+ IR = IROW
+ END IF
+ 150 CONTINUE
+ 160 CONTINUE
+ 170 CONTINUE
+*
+ ELSE
+*
+* Bottom-Up -- Start at the bottom right.
+*
+ JKL = 0
+ DO 200 JKU = 1, UUB
+*
+* Transform from bandwidth JKL, JKU-1 to JKL, JKU
+*
+* First row actually rotated is M
+* First column actually rotated is MIN( M+JKU, N )
+*
+ IENDCH = MIN( M, N+JKL ) - 1
+ DO 190 JC = MIN( M+JKU, N ) - 1, 1 - JKL, -1
+ EXTRA = ZERO
+ ANGLE = TWOPI*DLARND( 1, ISEED )
+ C = COS( ANGLE )
+ S = SIN( ANGLE )
+ IROW = MAX( 1, JC-JKU+1 )
+ IF( JC.GT.0 ) THEN
+ IL = MIN( M, JC+JKL+1 ) + 1 - IROW
+ CALL DLAROT( .FALSE., .FALSE., JC+JKL.LT.M, IL,
+ $ C, S, A( IROW-ISKEW*JC+IOFFST,
+ $ JC ), ILDA, DUMMY, EXTRA )
+ END IF
+*
+* Chase "EXTRA" back down
+*
+ IC = JC
+ DO 180 JCH = JC + JKL, IENDCH, JKL + JKU
+ ILEXTR = IC.GT.0
+ IF( ILEXTR ) THEN
+ CALL DLARTG( A( JCH-ISKEW*IC+IOFFST, IC ),
+ $ EXTRA, C, S, DUMMY )
+ END IF
+ IC = MAX( 1, IC )
+ ICOL = MIN( N-1, JCH+JKU )
+ ILTEMP = JCH + JKU.LT.N
+ TEMP = ZERO
+ CALL DLAROT( .TRUE., ILEXTR, ILTEMP, ICOL+2-IC,
+ $ C, S, A( JCH-ISKEW*IC+IOFFST, IC ),
+ $ ILDA, EXTRA, TEMP )
+ IF( ILTEMP ) THEN
+ CALL DLARTG( A( JCH-ISKEW*ICOL+IOFFST,
+ $ ICOL ), TEMP, C, S, DUMMY )
+ IL = MIN( IENDCH, JCH+JKL+JKU ) + 2 - JCH
+ EXTRA = ZERO
+ CALL DLAROT( .FALSE., .TRUE.,
+ $ JCH+JKL+JKU.LE.IENDCH, IL, C, S,
+ $ A( JCH-ISKEW*ICOL+IOFFST,
+ $ ICOL ), ILDA, TEMP, EXTRA )
+ IC = ICOL
+ END IF
+ 180 CONTINUE
+ 190 CONTINUE
+ 200 CONTINUE
+*
+ JKU = UUB
+ DO 230 JKL = 1, LLB
+*
+* Transform from bandwidth JKL-1, JKU to JKL, JKU
+*
+* First row actually rotated is MIN( N+JKL, M )
+* First column actually rotated is N
+*
+ IENDCH = MIN( N, M+JKU ) - 1
+ DO 220 JR = MIN( N+JKL, M ) - 1, 1 - JKU, -1
+ EXTRA = ZERO
+ ANGLE = TWOPI*DLARND( 1, ISEED )
+ C = COS( ANGLE )
+ S = SIN( ANGLE )
+ ICOL = MAX( 1, JR-JKL+1 )
+ IF( JR.GT.0 ) THEN
+ IL = MIN( N, JR+JKU+1 ) + 1 - ICOL
+ CALL DLAROT( .TRUE., .FALSE., JR+JKU.LT.N, IL,
+ $ C, S, A( JR-ISKEW*ICOL+IOFFST,
+ $ ICOL ), ILDA, DUMMY, EXTRA )
+ END IF
+*
+* Chase "EXTRA" back down
+*
+ IR = JR
+ DO 210 JCH = JR + JKU, IENDCH, JKL + JKU
+ ILEXTR = IR.GT.0
+ IF( ILEXTR ) THEN
+ CALL DLARTG( A( IR-ISKEW*JCH+IOFFST, JCH ),
+ $ EXTRA, C, S, DUMMY )
+ END IF
+ IR = MAX( 1, IR )
+ IROW = MIN( M-1, JCH+JKL )
+ ILTEMP = JCH + JKL.LT.M
+ TEMP = ZERO
+ CALL DLAROT( .FALSE., ILEXTR, ILTEMP, IROW+2-IR,
+ $ C, S, A( IR-ISKEW*JCH+IOFFST,
+ $ JCH ), ILDA, EXTRA, TEMP )
+ IF( ILTEMP ) THEN
+ CALL DLARTG( A( IROW-ISKEW*JCH+IOFFST, JCH ),
+ $ TEMP, C, S, DUMMY )
+ IL = MIN( IENDCH, JCH+JKL+JKU ) + 2 - JCH
+ EXTRA = ZERO
+ CALL DLAROT( .TRUE., .TRUE.,
+ $ JCH+JKL+JKU.LE.IENDCH, IL, C, S,
+ $ A( IROW-ISKEW*JCH+IOFFST, JCH ),
+ $ ILDA, TEMP, EXTRA )
+ IR = IROW
+ END IF
+ 210 CONTINUE
+ 220 CONTINUE
+ 230 CONTINUE
+ END IF
+*
+ ELSE
+*
+* Symmetric -- A = U D U'
+*
+ IPACKG = IPACK
+ IOFFG = IOFFST
+*
+ IF( TOPDWN ) THEN
+*
+* Top-Down -- Generate Upper triangle only
+*
+ IF( IPACK.GE.5 ) THEN
+ IPACKG = 6
+ IOFFG = UUB + 1
+ ELSE
+ IPACKG = 1
+ END IF
+ CALL DCOPY( MNMIN, D, 1, A( 1-ISKEW+IOFFG, 1 ), ILDA+1 )
+*
+ DO 260 K = 1, UUB
+ DO 250 JC = 1, N - 1
+ IROW = MAX( 1, JC-K )
+ IL = MIN( JC+1, K+2 )
+ EXTRA = ZERO
+ TEMP = A( JC-ISKEW*( JC+1 )+IOFFG, JC+1 )
+ ANGLE = TWOPI*DLARND( 1, ISEED )
+ C = COS( ANGLE )
+ S = SIN( ANGLE )
+ CALL DLAROT( .FALSE., JC.GT.K, .TRUE., IL, C, S,
+ $ A( IROW-ISKEW*JC+IOFFG, JC ), ILDA,
+ $ EXTRA, TEMP )
+ CALL DLAROT( .TRUE., .TRUE., .FALSE.,
+ $ MIN( K, N-JC )+1, C, S,
+ $ A( ( 1-ISKEW )*JC+IOFFG, JC ), ILDA,
+ $ TEMP, DUMMY )
+*
+* Chase EXTRA back up the matrix
+*
+ ICOL = JC
+ DO 240 JCH = JC - K, 1, -K
+ CALL DLARTG( A( JCH+1-ISKEW*( ICOL+1 )+IOFFG,
+ $ ICOL+1 ), EXTRA, C, S, DUMMY )
+ TEMP = A( JCH-ISKEW*( JCH+1 )+IOFFG, JCH+1 )
+ CALL DLAROT( .TRUE., .TRUE., .TRUE., K+2, C, -S,
+ $ A( ( 1-ISKEW )*JCH+IOFFG, JCH ),
+ $ ILDA, TEMP, EXTRA )
+ IROW = MAX( 1, JCH-K )
+ IL = MIN( JCH+1, K+2 )
+ EXTRA = ZERO
+ CALL DLAROT( .FALSE., JCH.GT.K, .TRUE., IL, C,
+ $ -S, A( IROW-ISKEW*JCH+IOFFG, JCH ),
+ $ ILDA, EXTRA, TEMP )
+ ICOL = JCH
+ 240 CONTINUE
+ 250 CONTINUE
+ 260 CONTINUE
+*
+* If we need lower triangle, copy from upper. Note that
+* the order of copying is chosen to work for 'q' -> 'b'
+*
+ IF( IPACK.NE.IPACKG .AND. IPACK.NE.3 ) THEN
+ DO 280 JC = 1, N
+ IROW = IOFFST - ISKEW*JC
+ DO 270 JR = JC, MIN( N, JC+UUB )
+ A( JR+IROW, JC ) = A( JC-ISKEW*JR+IOFFG, JR )
+ 270 CONTINUE
+ 280 CONTINUE
+ IF( IPACK.EQ.5 ) THEN
+ DO 300 JC = N - UUB + 1, N
+ DO 290 JR = N + 2 - JC, UUB + 1
+ A( JR, JC ) = ZERO
+ 290 CONTINUE
+ 300 CONTINUE
+ END IF
+ IF( IPACKG.EQ.6 ) THEN
+ IPACKG = IPACK
+ ELSE
+ IPACKG = 0
+ END IF
+ END IF
+ ELSE
+*
+* Bottom-Up -- Generate Lower triangle only
+*
+ IF( IPACK.GE.5 ) THEN
+ IPACKG = 5
+ IF( IPACK.EQ.6 )
+ $ IOFFG = 1
+ ELSE
+ IPACKG = 2
+ END IF
+ CALL DCOPY( MNMIN, D, 1, A( 1-ISKEW+IOFFG, 1 ), ILDA+1 )
+*
+ DO 330 K = 1, UUB
+ DO 320 JC = N - 1, 1, -1
+ IL = MIN( N+1-JC, K+2 )
+ EXTRA = ZERO
+ TEMP = A( 1+( 1-ISKEW )*JC+IOFFG, JC )
+ ANGLE = TWOPI*DLARND( 1, ISEED )
+ C = COS( ANGLE )
+ S = -SIN( ANGLE )
+ CALL DLAROT( .FALSE., .TRUE., N-JC.GT.K, IL, C, S,
+ $ A( ( 1-ISKEW )*JC+IOFFG, JC ), ILDA,
+ $ TEMP, EXTRA )
+ ICOL = MAX( 1, JC-K+1 )
+ CALL DLAROT( .TRUE., .FALSE., .TRUE., JC+2-ICOL, C,
+ $ S, A( JC-ISKEW*ICOL+IOFFG, ICOL ),
+ $ ILDA, DUMMY, TEMP )
+*
+* Chase EXTRA back down the matrix
+*
+ ICOL = JC
+ DO 310 JCH = JC + K, N - 1, K
+ CALL DLARTG( A( JCH-ISKEW*ICOL+IOFFG, ICOL ),
+ $ EXTRA, C, S, DUMMY )
+ TEMP = A( 1+( 1-ISKEW )*JCH+IOFFG, JCH )
+ CALL DLAROT( .TRUE., .TRUE., .TRUE., K+2, C, S,
+ $ A( JCH-ISKEW*ICOL+IOFFG, ICOL ),
+ $ ILDA, EXTRA, TEMP )
+ IL = MIN( N+1-JCH, K+2 )
+ EXTRA = ZERO
+ CALL DLAROT( .FALSE., .TRUE., N-JCH.GT.K, IL, C,
+ $ S, A( ( 1-ISKEW )*JCH+IOFFG, JCH ),
+ $ ILDA, TEMP, EXTRA )
+ ICOL = JCH
+ 310 CONTINUE
+ 320 CONTINUE
+ 330 CONTINUE
+*
+* If we need upper triangle, copy from lower. Note that
+* the order of copying is chosen to work for 'b' -> 'q'
+*
+ IF( IPACK.NE.IPACKG .AND. IPACK.NE.4 ) THEN
+ DO 350 JC = N, 1, -1
+ IROW = IOFFST - ISKEW*JC
+ DO 340 JR = JC, MAX( 1, JC-UUB ), -1
+ A( JR+IROW, JC ) = A( JC-ISKEW*JR+IOFFG, JR )
+ 340 CONTINUE
+ 350 CONTINUE
+ IF( IPACK.EQ.6 ) THEN
+ DO 370 JC = 1, UUB
+ DO 360 JR = 1, UUB + 1 - JC
+ A( JR, JC ) = ZERO
+ 360 CONTINUE
+ 370 CONTINUE
+ END IF
+ IF( IPACKG.EQ.5 ) THEN
+ IPACKG = IPACK
+ ELSE
+ IPACKG = 0
+ END IF
+ END IF
+ END IF
+ END IF
+*
+ ELSE
+*
+* 4) Generate Banded Matrix by first
+* Rotating by random Unitary matrices,
+* then reducing the bandwidth using Householder
+* transformations.
+*
+* Note: we should get here only if LDA .ge. N
+*
+ IF( ISYM.EQ.1 ) THEN
+*
+* Non-symmetric -- A = U D V
+*
+ CALL DLAGGE( MR, NC, LLB, UUB, D, A, LDA, ISEED, WORK,
+ $ IINFO )
+ ELSE
+*
+* Symmetric -- A = U D U'
+*
+ CALL DLAGSY( M, LLB, D, A, LDA, ISEED, WORK, IINFO )
+*
+ END IF
+ IF( IINFO.NE.0 ) THEN
+ INFO = 3
+ RETURN
+ END IF
+ END IF
+*
+* 5) Pack the matrix
+*
+ IF( IPACK.NE.IPACKG ) THEN
+ IF( IPACK.EQ.1 ) THEN
+*
+* 'U' -- Upper triangular, not packed
+*
+ DO 390 J = 1, M
+ DO 380 I = J + 1, M
+ A( I, J ) = ZERO
+ 380 CONTINUE
+ 390 CONTINUE
+*
+ ELSE IF( IPACK.EQ.2 ) THEN
+*
+* 'L' -- Lower triangular, not packed
+*
+ DO 410 J = 2, M
+ DO 400 I = 1, J - 1
+ A( I, J ) = ZERO
+ 400 CONTINUE
+ 410 CONTINUE
+*
+ ELSE IF( IPACK.EQ.3 ) THEN
+*
+* 'C' -- Upper triangle packed Columnwise.
+*
+ ICOL = 1
+ IROW = 0
+ DO 430 J = 1, M
+ DO 420 I = 1, J
+ IROW = IROW + 1
+ IF( IROW.GT.LDA ) THEN
+ IROW = 1
+ ICOL = ICOL + 1
+ END IF
+ A( IROW, ICOL ) = A( I, J )
+ 420 CONTINUE
+ 430 CONTINUE
+*
+ ELSE IF( IPACK.EQ.4 ) THEN
+*
+* 'R' -- Lower triangle packed Columnwise.
+*
+ ICOL = 1
+ IROW = 0
+ DO 450 J = 1, M
+ DO 440 I = J, M
+ IROW = IROW + 1
+ IF( IROW.GT.LDA ) THEN
+ IROW = 1
+ ICOL = ICOL + 1
+ END IF
+ A( IROW, ICOL ) = A( I, J )
+ 440 CONTINUE
+ 450 CONTINUE
+*
+ ELSE IF( IPACK.GE.5 ) THEN
+*
+* 'B' -- The lower triangle is packed as a band matrix.
+* 'Q' -- The upper triangle is packed as a band matrix.
+* 'Z' -- The whole matrix is packed as a band matrix.
+*
+ IF( IPACK.EQ.5 )
+ $ UUB = 0
+ IF( IPACK.EQ.6 )
+ $ LLB = 0
+*
+ DO 470 J = 1, UUB
+ DO 460 I = MIN( J+LLB, M ), 1, -1
+ A( I-J+UUB+1, J ) = A( I, J )
+ 460 CONTINUE
+ 470 CONTINUE
+*
+ DO 490 J = UUB + 2, N
+ DO 480 I = J - UUB, MIN( J+LLB, M )
+ A( I-J+UUB+1, J ) = A( I, J )
+ 480 CONTINUE
+ 490 CONTINUE
+ END IF
+*
+* If packed, zero out extraneous elements.
+*
+* Symmetric/Triangular Packed --
+* zero out everything after A(IROW,ICOL)
+*
+ IF( IPACK.EQ.3 .OR. IPACK.EQ.4 ) THEN
+ DO 510 JC = ICOL, M
+ DO 500 JR = IROW + 1, LDA
+ A( JR, JC ) = ZERO
+ 500 CONTINUE
+ IROW = 0
+ 510 CONTINUE
+*
+ ELSE IF( IPACK.GE.5 ) THEN
+*
+* Packed Band --
+* 1st row is now in A( UUB+2-j, j), zero above it
+* m-th row is now in A( M+UUB-j,j), zero below it
+* last non-zero diagonal is now in A( UUB+LLB+1,j ),
+* zero below it, too.
+*
+ IR1 = UUB + LLB + 2
+ IR2 = UUB + M + 2
+ DO 540 JC = 1, N
+ DO 520 JR = 1, UUB + 1 - JC
+ A( JR, JC ) = ZERO
+ 520 CONTINUE
+ DO 530 JR = MAX( 1, MIN( IR1, IR2-JC ) ), LDA
+ A( JR, JC ) = ZERO
+ 530 CONTINUE
+ 540 CONTINUE
+ END IF
+ END IF
+*
+ RETURN
+*
+* End of DLATMT
+*
+ END
generated by cgit v1.2.3 (git 2.25.1) at 2024-11-22 22:04:26 +0000