LAPACK  3.8.0
LAPACK: Linear Algebra PACKage

◆ dlatme()

subroutine dlatme ( integer  N,
character  DIST,
integer, dimension( 4 )  ISEED,
double precision, dimension( * )  D,
integer  MODE,
double precision  COND,
double precision  DMAX,
character, dimension( * )  EI,
character  RSIGN,
character  UPPER,
character  SIM,
double precision, dimension( * )  DS,
integer  MODES,
double precision  CONDS,
integer  KL,
integer  KU,
double precision  ANORM,
double precision, dimension( lda, * )  A,
integer  LDA,
double precision, dimension( * )  WORK,
integer  INFO 
)

DLATME

Purpose:
    DLATME generates random non-symmetric square matrices with
    specified eigenvalues for testing LAPACK programs.

    DLATME operates by applying the following sequence of
    operations:

    1. Set the diagonal to D, where D may be input or
         computed according to MODE, COND, DMAX, and RSIGN
         as described below.

    2. If complex conjugate pairs are desired (MODE=0 and EI(1)='R',
         or MODE=5), certain pairs of adjacent elements of D are
         interpreted as the real and complex parts of a complex
         conjugate pair; A thus becomes block diagonal, with 1x1
         and 2x2 blocks.

    3. If UPPER='T', the upper triangle of A is set to random values
         out of distribution DIST.

    4. If SIM='T', A is multiplied on the left by a random matrix
         X, whose singular values are specified by DS, MODES, and
         CONDS, and on the right by X inverse.

    5. If KL < N-1, the lower bandwidth is reduced to KL using
         Householder transformations.  If KU < N-1, the upper
         bandwidth is reduced to KU.

    6. If ANORM is not negative, the matrix is scaled to have
         maximum-element-norm ANORM.

    (Note: since the matrix cannot be reduced beyond Hessenberg form,
     no packing options are available.)
Parameters
[in]N
          N is INTEGER
           The number of columns (or rows) 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, and for the
           upper triangle (see UPPER).
           '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 DLATME
           to continue the same random number sequence.
           Changed on exit.
[in,out]D
          D is DOUBLE PRECISION array, dimension ( N )
           This array is used to specify the eigenvalues of A.  If
           MODE=0, then D is assumed to contain the eigenvalues (but
           see the description of EI), otherwise they will be
           computed according to MODE, COND, DMAX, and RSIGN and
           placed in D.
           Modified if MODE is nonzero.
[in]MODE
          MODE is INTEGER
           On entry this describes how the eigenvalues are to
           be specified:
           MODE = 0 means use D (with EI) 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.  Each odd-even pair
                    of elements will be either used as two real
                    eigenvalues or as the real and imaginary part
                    of a complex conjugate pair of eigenvalues;
                    the choice of which is done is random, with
                    50-50 probability, for each pair.
           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 between 1 and 4, D has entries ranging
              from 1 to 1/COND, if between -1 and -4, D has entries
              ranging from 1/COND to 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))).  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]EI
          EI is CHARACTER*1 array, dimension ( N )
           If MODE is 0, and EI(1) is not ' ' (space character),
           this array specifies which elements of D (on input) are
           real eigenvalues and which are the real and imaginary parts
           of a complex conjugate pair of eigenvalues.  The elements
           of EI may then only have the values 'R' and 'I'.  If
           EI(j)='R' and EI(j+1)='I', then the j-th eigenvalue is
           CMPLX( D(j) , D(j+1) ), and the (j+1)-th is the complex
           conjugate thereof.  If EI(j)=EI(j+1)='R', then the j-th
           eigenvalue is D(j) (i.e., real).  EI(1) may not be 'I',
           nor may two adjacent elements of EI both have the value 'I'.
           If MODE is not 0, then EI is ignored.  If MODE is 0 and
           EI(1)=' ', then the eigenvalues will all be real.
           Not modified.
[in]RSIGN
          RSIGN is CHARACTER*1
           If MODE is not 0, 6, or -6, and RSIGN='T', then the
           elements of D, as computed according to MODE and COND, will
           be multiplied by a random sign (+1 or -1).  If RSIGN='F',
           they will not be.  RSIGN may only have the values 'T' or
           'F'.
           Not modified.
[in]UPPER
          UPPER is CHARACTER*1
           If UPPER='T', then the elements of A above the diagonal
           (and above the 2x2 diagonal blocks, if A has complex
           eigenvalues) will be set to random numbers out of DIST.
           If UPPER='F', they will not.  UPPER may only have the
           values 'T' or 'F'.
           Not modified.
[in]SIM
          SIM is CHARACTER*1
           If SIM='T', then A will be operated on by a "similarity
           transform", i.e., multiplied on the left by a matrix X and
           on the right by X inverse.  X = U S V, where U and V are
           random unitary matrices and S is a (diagonal) matrix of
           singular values specified by DS, MODES, and CONDS.  If
           SIM='F', then A will not be transformed.
           Not modified.
[in,out]DS
          DS is DOUBLE PRECISION array, dimension ( N )
           This array is used to specify the singular values of X,
           in the same way that D specifies the eigenvalues of A.
           If MODE=0, the DS contains the singular values, which
           may not be zero.
           Modified if MODE is nonzero.
[in]MODES
          MODES is INTEGER
[in]CONDS
          CONDS is DOUBLE PRECISION
           Same as MODE and COND, but for specifying the diagonal
           of S.  MODES=-6 and +6 are not allowed (since they would
           result in randomly ill-conditioned eigenvalues.)
[in]KL
          KL is INTEGER
           This specifies the lower bandwidth of the  matrix.  KL=1
           specifies upper Hessenberg form.  If KL is at least N-1,
           then A will have full lower bandwidth.  KL must be at
           least 1.
           Not modified.
[in]KU
          KU is INTEGER
           This specifies the upper bandwidth of the  matrix.  KU=1
           specifies lower Hessenberg form.  If KU is at least N-1,
           then A will have full upper bandwidth; if KU and KL
           are both at least N-1, then A will be dense.  Only one of
           KU and KL may be less than N-1.  KU must be at least 1.
           Not modified.
[in]ANORM
          ANORM is DOUBLE PRECISION
           If ANORM is not negative, then A will be scaled by a non-
           negative real number to make the maximum-element-norm of A
           to be ANORM.
           Not modified.
[out]A
          A is DOUBLE PRECISION array, dimension ( LDA, N )
           On exit A is the desired test matrix.
           Modified.
[in]LDA
          LDA is INTEGER
           LDA specifies the first dimension of A as declared in the
           calling program.  LDA must be at least N.
           Not modified.
[out]WORK
          WORK is DOUBLE PRECISION array, dimension ( 3*N )
           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 => N negative
            -2 => DIST illegal string
            -5 => MODE not in range -6 to 6
            -6 => COND less than 1.0, and MODE neither -6, 0 nor 6
            -8 => EI(1) is not ' ' or 'R', EI(j) is not 'R' or 'I', or
                  two adjacent elements of EI are 'I'.
            -9 => RSIGN is not 'T' or 'F'
           -10 => UPPER is not 'T' or 'F'
           -11 => SIM   is not 'T' or 'F'
           -12 => MODES=0 and DS has a zero singular value.
           -13 => MODES is not in the range -5 to 5.
           -14 => MODES is nonzero and CONDS is less than 1.
           -15 => KL is less than 1.
           -16 => KU is less than 1, or KL and KU are both less than
                  N-1.
           -19 => LDA is less than N.
            1  => Error return from DLATM1 (computing D)
            2  => Cannot scale to DMAX (max. eigenvalue is 0)
            3  => Error return from DLATM1 (computing DS)
            4  => Error return from DLARGE
            5  => Zero singular value from DLATM1.
Author
Univ. of Tennessee
Univ. of California Berkeley
Univ. of Colorado Denver
NAG Ltd.
Date
December 2016

Definition at line 334 of file dlatme.f.

334 *
335 * -- LAPACK computational routine (version 3.7.0) --
336 * -- LAPACK is a software package provided by Univ. of Tennessee, --
337 * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
338 * December 2016
339 *
340 * .. Scalar Arguments ..
341  CHARACTER dist, rsign, sim, upper
342  INTEGER info, kl, ku, lda, mode, modes, n
343  DOUBLE PRECISION anorm, cond, conds, dmax
344 * ..
345 * .. Array Arguments ..
346  CHARACTER ei( * )
347  INTEGER iseed( 4 )
348  DOUBLE PRECISION a( lda, * ), d( * ), ds( * ), work( * )
349 * ..
350 *
351 * =====================================================================
352 *
353 * .. Parameters ..
354  DOUBLE PRECISION zero
355  parameter( zero = 0.0d0 )
356  DOUBLE PRECISION one
357  parameter( one = 1.0d0 )
358  DOUBLE PRECISION half
359  parameter( half = 1.0d0 / 2.0d0 )
360 * ..
361 * .. Local Scalars ..
362  LOGICAL badei, bads, useei
363  INTEGER i, ic, icols, idist, iinfo, ir, irows, irsign,
364  $ isim, iupper, j, jc, jcr, jr
365  DOUBLE PRECISION alpha, tau, temp, xnorms
366 * ..
367 * .. Local Arrays ..
368  DOUBLE PRECISION tempa( 1 )
369 * ..
370 * .. External Functions ..
371  LOGICAL lsame
372  DOUBLE PRECISION dlange, dlaran
373  EXTERNAL lsame, dlange, dlaran
374 * ..
375 * .. External Subroutines ..
376  EXTERNAL dcopy, dgemv, dger, dlarfg, dlarge, dlarnv,
378 * ..
379 * .. Intrinsic Functions ..
380  INTRINSIC abs, max, mod
381 * ..
382 * .. Executable Statements ..
383 *
384 * 1) Decode and Test the input parameters.
385 * Initialize flags & seed.
386 *
387  info = 0
388 *
389 * Quick return if possible
390 *
391  IF( n.EQ.0 )
392  $ RETURN
393 *
394 * Decode DIST
395 *
396  IF( lsame( dist, 'U' ) ) THEN
397  idist = 1
398  ELSE IF( lsame( dist, 'S' ) ) THEN
399  idist = 2
400  ELSE IF( lsame( dist, 'N' ) ) THEN
401  idist = 3
402  ELSE
403  idist = -1
404  END IF
405 *
406 * Check EI
407 *
408  useei = .true.
409  badei = .false.
410  IF( lsame( ei( 1 ), ' ' ) .OR. mode.NE.0 ) THEN
411  useei = .false.
412  ELSE
413  IF( lsame( ei( 1 ), 'R' ) ) THEN
414  DO 10 j = 2, n
415  IF( lsame( ei( j ), 'I' ) ) THEN
416  IF( lsame( ei( j-1 ), 'I' ) )
417  $ badei = .true.
418  ELSE
419  IF( .NOT.lsame( ei( j ), 'R' ) )
420  $ badei = .true.
421  END IF
422  10 CONTINUE
423  ELSE
424  badei = .true.
425  END IF
426  END IF
427 *
428 * Decode RSIGN
429 *
430  IF( lsame( rsign, 'T' ) ) THEN
431  irsign = 1
432  ELSE IF( lsame( rsign, 'F' ) ) THEN
433  irsign = 0
434  ELSE
435  irsign = -1
436  END IF
437 *
438 * Decode UPPER
439 *
440  IF( lsame( upper, 'T' ) ) THEN
441  iupper = 1
442  ELSE IF( lsame( upper, 'F' ) ) THEN
443  iupper = 0
444  ELSE
445  iupper = -1
446  END IF
447 *
448 * Decode SIM
449 *
450  IF( lsame( sim, 'T' ) ) THEN
451  isim = 1
452  ELSE IF( lsame( sim, 'F' ) ) THEN
453  isim = 0
454  ELSE
455  isim = -1
456  END IF
457 *
458 * Check DS, if MODES=0 and ISIM=1
459 *
460  bads = .false.
461  IF( modes.EQ.0 .AND. isim.EQ.1 ) THEN
462  DO 20 j = 1, n
463  IF( ds( j ).EQ.zero )
464  $ bads = .true.
465  20 CONTINUE
466  END IF
467 *
468 * Set INFO if an error
469 *
470  IF( n.LT.0 ) THEN
471  info = -1
472  ELSE IF( idist.EQ.-1 ) THEN
473  info = -2
474  ELSE IF( abs( mode ).GT.6 ) THEN
475  info = -5
476  ELSE IF( ( mode.NE.0 .AND. abs( mode ).NE.6 ) .AND. cond.LT.one )
477  $ THEN
478  info = -6
479  ELSE IF( badei ) THEN
480  info = -8
481  ELSE IF( irsign.EQ.-1 ) THEN
482  info = -9
483  ELSE IF( iupper.EQ.-1 ) THEN
484  info = -10
485  ELSE IF( isim.EQ.-1 ) THEN
486  info = -11
487  ELSE IF( bads ) THEN
488  info = -12
489  ELSE IF( isim.EQ.1 .AND. abs( modes ).GT.5 ) THEN
490  info = -13
491  ELSE IF( isim.EQ.1 .AND. modes.NE.0 .AND. conds.LT.one ) THEN
492  info = -14
493  ELSE IF( kl.LT.1 ) THEN
494  info = -15
495  ELSE IF( ku.LT.1 .OR. ( ku.LT.n-1 .AND. kl.LT.n-1 ) ) THEN
496  info = -16
497  ELSE IF( lda.LT.max( 1, n ) ) THEN
498  info = -19
499  END IF
500 *
501  IF( info.NE.0 ) THEN
502  CALL xerbla( 'DLATME', -info )
503  RETURN
504  END IF
505 *
506 * Initialize random number generator
507 *
508  DO 30 i = 1, 4
509  iseed( i ) = mod( abs( iseed( i ) ), 4096 )
510  30 CONTINUE
511 *
512  IF( mod( iseed( 4 ), 2 ).NE.1 )
513  $ iseed( 4 ) = iseed( 4 ) + 1
514 *
515 * 2) Set up diagonal of A
516 *
517 * Compute D according to COND and MODE
518 *
519  CALL dlatm1( mode, cond, irsign, idist, iseed, d, n, iinfo )
520  IF( iinfo.NE.0 ) THEN
521  info = 1
522  RETURN
523  END IF
524  IF( mode.NE.0 .AND. abs( mode ).NE.6 ) THEN
525 *
526 * Scale by DMAX
527 *
528  temp = abs( d( 1 ) )
529  DO 40 i = 2, n
530  temp = max( temp, abs( d( i ) ) )
531  40 CONTINUE
532 *
533  IF( temp.GT.zero ) THEN
534  alpha = dmax / temp
535  ELSE IF( dmax.NE.zero ) THEN
536  info = 2
537  RETURN
538  ELSE
539  alpha = zero
540  END IF
541 *
542  CALL dscal( n, alpha, d, 1 )
543 *
544  END IF
545 *
546  CALL dlaset( 'Full', n, n, zero, zero, a, lda )
547  CALL dcopy( n, d, 1, a, lda+1 )
548 *
549 * Set up complex conjugate pairs
550 *
551  IF( mode.EQ.0 ) THEN
552  IF( useei ) THEN
553  DO 50 j = 2, n
554  IF( lsame( ei( j ), 'I' ) ) THEN
555  a( j-1, j ) = a( j, j )
556  a( j, j-1 ) = -a( j, j )
557  a( j, j ) = a( j-1, j-1 )
558  END IF
559  50 CONTINUE
560  END IF
561 *
562  ELSE IF( abs( mode ).EQ.5 ) THEN
563 *
564  DO 60 j = 2, n, 2
565  IF( dlaran( iseed ).GT.half ) THEN
566  a( j-1, j ) = a( j, j )
567  a( j, j-1 ) = -a( j, j )
568  a( j, j ) = a( j-1, j-1 )
569  END IF
570  60 CONTINUE
571  END IF
572 *
573 * 3) If UPPER='T', set upper triangle of A to random numbers.
574 * (but don't modify the corners of 2x2 blocks.)
575 *
576  IF( iupper.NE.0 ) THEN
577  DO 70 jc = 2, n
578  IF( a( jc-1, jc ).NE.zero ) THEN
579  jr = jc - 2
580  ELSE
581  jr = jc - 1
582  END IF
583  CALL dlarnv( idist, iseed, jr, a( 1, jc ) )
584  70 CONTINUE
585  END IF
586 *
587 * 4) If SIM='T', apply similarity transformation.
588 *
589 * -1
590 * Transform is X A X , where X = U S V, thus
591 *
592 * it is U S V A V' (1/S) U'
593 *
594  IF( isim.NE.0 ) THEN
595 *
596 * Compute S (singular values of the eigenvector matrix)
597 * according to CONDS and MODES
598 *
599  CALL dlatm1( modes, conds, 0, 0, iseed, ds, n, iinfo )
600  IF( iinfo.NE.0 ) THEN
601  info = 3
602  RETURN
603  END IF
604 *
605 * Multiply by V and V'
606 *
607  CALL dlarge( n, a, lda, iseed, work, iinfo )
608  IF( iinfo.NE.0 ) THEN
609  info = 4
610  RETURN
611  END IF
612 *
613 * Multiply by S and (1/S)
614 *
615  DO 80 j = 1, n
616  CALL dscal( n, ds( j ), a( j, 1 ), lda )
617  IF( ds( j ).NE.zero ) THEN
618  CALL dscal( n, one / ds( j ), a( 1, j ), 1 )
619  ELSE
620  info = 5
621  RETURN
622  END IF
623  80 CONTINUE
624 *
625 * Multiply by U and U'
626 *
627  CALL dlarge( n, a, lda, iseed, work, iinfo )
628  IF( iinfo.NE.0 ) THEN
629  info = 4
630  RETURN
631  END IF
632  END IF
633 *
634 * 5) Reduce the bandwidth.
635 *
636  IF( kl.LT.n-1 ) THEN
637 *
638 * Reduce bandwidth -- kill column
639 *
640  DO 90 jcr = kl + 1, n - 1
641  ic = jcr - kl
642  irows = n + 1 - jcr
643  icols = n + kl - jcr
644 *
645  CALL dcopy( irows, a( jcr, ic ), 1, work, 1 )
646  xnorms = work( 1 )
647  CALL dlarfg( irows, xnorms, work( 2 ), 1, tau )
648  work( 1 ) = one
649 *
650  CALL dgemv( 'T', irows, icols, one, a( jcr, ic+1 ), lda,
651  $ work, 1, zero, work( irows+1 ), 1 )
652  CALL dger( irows, icols, -tau, work, 1, work( irows+1 ), 1,
653  $ a( jcr, ic+1 ), lda )
654 *
655  CALL dgemv( 'N', n, irows, one, a( 1, jcr ), lda, work, 1,
656  $ zero, work( irows+1 ), 1 )
657  CALL dger( n, irows, -tau, work( irows+1 ), 1, work, 1,
658  $ a( 1, jcr ), lda )
659 *
660  a( jcr, ic ) = xnorms
661  CALL dlaset( 'Full', irows-1, 1, zero, zero, a( jcr+1, ic ),
662  $ lda )
663  90 CONTINUE
664  ELSE IF( ku.LT.n-1 ) THEN
665 *
666 * Reduce upper bandwidth -- kill a row at a time.
667 *
668  DO 100 jcr = ku + 1, n - 1
669  ir = jcr - ku
670  irows = n + ku - jcr
671  icols = n + 1 - jcr
672 *
673  CALL dcopy( icols, a( ir, jcr ), lda, work, 1 )
674  xnorms = work( 1 )
675  CALL dlarfg( icols, xnorms, work( 2 ), 1, tau )
676  work( 1 ) = one
677 *
678  CALL dgemv( 'N', irows, icols, one, a( ir+1, jcr ), lda,
679  $ work, 1, zero, work( icols+1 ), 1 )
680  CALL dger( irows, icols, -tau, work( icols+1 ), 1, work, 1,
681  $ a( ir+1, jcr ), lda )
682 *
683  CALL dgemv( 'C', icols, n, one, a( jcr, 1 ), lda, work, 1,
684  $ zero, work( icols+1 ), 1 )
685  CALL dger( icols, n, -tau, work, 1, work( icols+1 ), 1,
686  $ a( jcr, 1 ), lda )
687 *
688  a( ir, jcr ) = xnorms
689  CALL dlaset( 'Full', 1, icols-1, zero, zero, a( ir, jcr+1 ),
690  $ lda )
691  100 CONTINUE
692  END IF
693 *
694 * Scale the matrix to have norm ANORM
695 *
696  IF( anorm.GE.zero ) THEN
697  temp = dlange( 'M', n, n, a, lda, tempa )
698  IF( temp.GT.zero ) THEN
699  alpha = anorm / temp
700  DO 110 j = 1, n
701  CALL dscal( n, alpha, a( 1, j ), 1 )
702  110 CONTINUE
703  END IF
704  END IF
705 *
706  RETURN
707 *
708 * End of DLATME
709 *
subroutine dcopy(N, DX, INCX, DY, INCY)
DCOPY
Definition: dcopy.f:84
subroutine dlatm1(MODE, COND, IRSIGN, IDIST, ISEED, D, N, INFO)
DLATM1
Definition: dlatm1.f:137
subroutine dgemv(TRANS, M, N, ALPHA, A, LDA, X, INCX, BETA, Y, INCY)
DGEMV
Definition: dgemv.f:158
subroutine dlarnv(IDIST, ISEED, N, X)
DLARNV returns a vector of random numbers from a uniform or normal distribution.
Definition: dlarnv.f:99
subroutine dlaset(UPLO, M, N, ALPHA, BETA, A, LDA)
DLASET initializes the off-diagonal elements and the diagonal elements of a matrix to given values...
Definition: dlaset.f:112
double precision function dlange(NORM, M, N, A, LDA, WORK)
DLANGE returns the value of the 1-norm, Frobenius norm, infinity-norm, or the largest absolute value ...
Definition: dlange.f:116
subroutine dger(M, N, ALPHA, X, INCX, Y, INCY, A, LDA)
DGER
Definition: dger.f:132
subroutine xerbla(SRNAME, INFO)
XERBLA
Definition: xerbla.f:62
logical function lsame(CA, CB)
LSAME
Definition: lsame.f:55
double precision function dlaran(ISEED)
DLARAN
Definition: dlaran.f:69
subroutine dscal(N, DA, DX, INCX)
DSCAL
Definition: dscal.f:81
subroutine dlarfg(N, ALPHA, X, INCX, TAU)
DLARFG generates an elementary reflector (Householder matrix).
Definition: dlarfg.f:108
subroutine dlarge(N, A, LDA, ISEED, WORK, INFO)
DLARGE
Definition: dlarge.f:89
Here is the call graph for this function:
Here is the caller graph for this function: