zlattr.f

来自「famous linear algebra library (LAPACK) p」· F 代码 · 共 659 行 · 第 1/2 页

      SUBROUTINE ZLATTR( IMAT, UPLO, TRANS, DIAG, ISEED, N, A, LDA, B,
     $                   WORK, RWORK, INFO )
*
*  -- LAPACK test routine (version 3.1) --
*     Univ. of Tennessee, Univ. of California Berkeley and NAG Ltd..
*     November 2006
*
*     .. Scalar Arguments ..
      CHARACTER          DIAG, TRANS, UPLO
      INTEGER            IMAT, INFO, LDA, N
*     ..
*     .. Array Arguments ..
      INTEGER            ISEED( 4 )
      DOUBLE PRECISION   RWORK( * )
      COMPLEX*16         A( LDA, * ), B( * ), WORK( * )
*     ..
*
*  Purpose
*  =======
*
*  ZLATTR generates a triangular test matrix in 2-dimensional storage.
*  IMAT and UPLO uniquely specify the properties of the test matrix,
*  which is returned in the array A.
*
*  Arguments
*  =========
*
*  IMAT    (input) INTEGER
*          An integer key describing which matrix to generate for this
*          path.
*
*  UPLO    (input) CHARACTER*1
*          Specifies whether the matrix A will be upper or lower
*          triangular.
*          = 'U':  Upper triangular
*          = 'L':  Lower triangular
*
*  TRANS   (input) CHARACTER*1
*          Specifies whether the matrix or its transpose will be used.
*          = 'N':  No transpose
*          = 'T':  Transpose
*          = 'C':  Conjugate transpose
*
*  DIAG    (output) CHARACTER*1
*          Specifies whether or not the matrix A is unit triangular.
*          = 'N':  Non-unit triangular
*          = 'U':  Unit triangular
*
*  ISEED   (input/output) INTEGER array, dimension (4)
*          The seed vector for the random number generator (used in
*          ZLATMS).  Modified on exit.
*
*  N       (input) INTEGER
*          The order of the matrix to be generated.
*
*  A       (output) COMPLEX*16 array, dimension (LDA,N)
*          The triangular matrix A.  If UPLO = 'U', the leading N x N
*          upper triangular part of the array A contains the upper
*          triangular matrix, and the strictly lower triangular part of
*          A is not referenced.  If UPLO = 'L', the leading N x N lower
*          triangular part of the array A contains the lower triangular
*          matrix and the strictly upper triangular part of A is not
*          referenced.
*
*  LDA     (input) INTEGER
*          The leading dimension of the array A.  LDA >= max(1,N).
*
*  B       (output) COMPLEX*16 array, dimension (N)
*          The right hand side vector, if IMAT > 10.
*
*  WORK    (workspace) COMPLEX*16 array, dimension (2*N)
*
*  RWORK   (workspace) DOUBLE PRECISION array, dimension (N)
*
*  INFO    (output) INTEGER
*          = 0:  successful exit
*          < 0:  if INFO = -i, the i-th argument had an illegal value
*
*  =====================================================================
*
*     .. Parameters ..
      DOUBLE PRECISION   ONE, TWO, ZERO
      PARAMETER          ( ONE = 1.0D+0, TWO = 2.0D+0, ZERO = 0.0D+0 )
*     ..
*     .. Local Scalars ..
      LOGICAL            UPPER
      CHARACTER          DIST, TYPE
      CHARACTER*3        PATH
      INTEGER            I, IY, J, JCOUNT, KL, KU, MODE
      DOUBLE PRECISION   ANORM, BIGNUM, BNORM, BSCAL, C, CNDNUM, REXP,
     $                   SFAC, SMLNUM, TEXP, TLEFT, TSCAL, ULP, UNFL, X,
     $                   Y, Z
      COMPLEX*16         PLUS1, PLUS2, RA, RB, S, STAR1
*     ..
*     .. External Functions ..
      LOGICAL            LSAME
      INTEGER            IZAMAX
      DOUBLE PRECISION   DLAMCH, DLARND
      COMPLEX*16         ZLARND
      EXTERNAL           LSAME, IZAMAX, DLAMCH, DLARND, ZLARND
*     ..
*     .. External Subroutines ..
      EXTERNAL           DLABAD, DLARNV, ZCOPY, ZDSCAL, ZLARNV, ZLATB4,
     $                   ZLATMS, ZROT, ZROTG, ZSWAP
*     ..
*     .. Intrinsic Functions ..
      INTRINSIC          ABS, DBLE, DCMPLX, DCONJG, MAX, SQRT
*     ..
*     .. Executable Statements ..
*
      PATH( 1: 1 ) = 'Zomplex precision'
      PATH( 2: 3 ) = 'TR'
      UNFL = DLAMCH( 'Safe minimum' )
      ULP = DLAMCH( 'Epsilon' )*DLAMCH( 'Base' )
      SMLNUM = UNFL
      BIGNUM = ( ONE-ULP ) / SMLNUM
      CALL DLABAD( SMLNUM, BIGNUM )
      IF( ( IMAT.GE.7 .AND. IMAT.LE.10 ) .OR. IMAT.EQ.18 ) THEN
         DIAG = 'U'
      ELSE
         DIAG = 'N'
      END IF
      INFO = 0
*
*     Quick return if N.LE.0.
*
      IF( N.LE.0 )
     $   RETURN
*
*     Call ZLATB4 to set parameters for CLATMS.
*
      UPPER = LSAME( UPLO, 'U' )
      IF( UPPER ) THEN
         CALL ZLATB4( PATH, IMAT, N, N, TYPE, KL, KU, ANORM, MODE,
     $                CNDNUM, DIST )
      ELSE
         CALL ZLATB4( PATH, -IMAT, N, N, TYPE, KL, KU, ANORM, MODE,
     $                CNDNUM, DIST )
      END IF
*
*     IMAT <= 6:  Non-unit triangular matrix
*
      IF( IMAT.LE.6 ) THEN
         CALL ZLATMS( N, N, DIST, ISEED, TYPE, RWORK, MODE, CNDNUM,
     $                ANORM, KL, KU, 'No packing', A, LDA, WORK, INFO )
*
*     IMAT > 6:  Unit triangular matrix
*     The diagonal is deliberately set to something other than 1.
*
*     IMAT = 7:  Matrix is the identity
*
      ELSE IF( IMAT.EQ.7 ) THEN
         IF( UPPER ) THEN
            DO 20 J = 1, N
               DO 10 I = 1, J - 1
                  A( I, J ) = ZERO
   10          CONTINUE
               A( J, J ) = J
   20       CONTINUE
         ELSE
            DO 40 J = 1, N
               A( J, J ) = J
               DO 30 I = J + 1, N
                  A( I, J ) = ZERO
   30          CONTINUE
   40       CONTINUE
         END IF
*
*     IMAT > 7:  Non-trivial unit triangular matrix
*
*     Generate a unit triangular matrix T with condition CNDNUM by
*     forming a triangular matrix with known singular values and
*     filling in the zero entries with Givens rotations.
*
      ELSE IF( IMAT.LE.10 ) THEN
         IF( UPPER ) THEN
            DO 60 J = 1, N
               DO 50 I = 1, J - 1
                  A( I, J ) = ZERO
   50          CONTINUE
               A( J, J ) = J
   60       CONTINUE
         ELSE
            DO 80 J = 1, N
               A( J, J ) = J
               DO 70 I = J + 1, N
                  A( I, J ) = ZERO
   70          CONTINUE
   80       CONTINUE
         END IF
*
*        Since the trace of a unit triangular matrix is 1, the product
*        of its singular values must be 1.  Let s = sqrt(CNDNUM),
*        x = sqrt(s) - 1/sqrt(s), y = sqrt(2/(n-2))*x, and z = x**2.
*        The following triangular matrix has singular values s, 1, 1,
*        ..., 1, 1/s:
*
*        1  y  y  y  ...  y  y  z
*           1  0  0  ...  0  0  y
*              1  0  ...  0  0  y
*                 .  ...  .  .  .
*                     .   .  .  .
*                         1  0  y
*                            1  y
*                               1
*
*        To fill in the zeros, we first multiply by a matrix with small
*        condition number of the form
*
*        1  0  0  0  0  ...
*           1  +  *  0  0  ...
*              1  +  0  0  0
*                 1  +  *  0  0
*                    1  +  0  0
*                       ...
*                          1  +  0
*                             1  0
*                                1
*
*        Each element marked with a '*' is formed by taking the product
*        of the adjacent elements marked with '+'.  The '*'s can be
*        chosen freely, and the '+'s are chosen so that the inverse of
*        T will have elements of the same magnitude as T.  If the *'s in
*        both T and inv(T) have small magnitude, T is well conditioned.
*        The two offdiagonals of T are stored in WORK.
*
*        The product of these two matrices has the form
*
*        1  y  y  y  y  y  .  y  y  z
*           1  +  *  0  0  .  0  0  y
*              1  +  0  0  .  0  0  y
*                 1  +  *  .  .  .  .
*                    1  +  .  .  .  .
*                       .  .  .  .  .
*                          .  .  .  .
*                             1  +  y
*                                1  y
*                                   1
*
*        Now we multiply by Givens rotations, using the fact that
*
*              [  c   s ] [  1   w ] [ -c  -s ] =  [  1  -w ]
*              [ -s   c ] [  0   1 ] [  s  -c ]    [  0   1 ]
*        and
*              [ -c  -s ] [  1   0 ] [  c   s ] =  [  1   0 ]
*              [  s  -c ] [  w   1 ] [ -s   c ]    [ -w   1 ]
*
*        where c = w / sqrt(w**2+4) and s = 2 / sqrt(w**2+4).
*
         STAR1 = 0.25D0*ZLARND( 5, ISEED )
         SFAC = 0.5D0
         PLUS1 = SFAC*ZLARND( 5, ISEED )
         DO 90 J = 1, N, 2
            PLUS2 = STAR1 / PLUS1
            WORK( J ) = PLUS1
            WORK( N+J ) = STAR1
            IF( J+1.LE.N ) THEN
               WORK( J+1 ) = PLUS2
               WORK( N+J+1 ) = ZERO
               PLUS1 = STAR1 / PLUS2
               REXP = DLARND( 2, ISEED )
               IF( REXP.LT.ZERO ) THEN
                  STAR1 = -SFAC**( ONE-REXP )*ZLARND( 5, ISEED )
               ELSE
                  STAR1 = SFAC**( ONE+REXP )*ZLARND( 5, ISEED )
               END IF
            END IF
   90    CONTINUE
*
         X = SQRT( CNDNUM ) - 1 / SQRT( CNDNUM )
         IF( N.GT.2 ) THEN
            Y = SQRT( 2.D0 / ( N-2 ) )*X
         ELSE
            Y = ZERO
         END IF
         Z = X*X
*
         IF( UPPER ) THEN
            IF( N.GT.3 ) THEN
               CALL ZCOPY( N-3, WORK, 1, A( 2, 3 ), LDA+1 )
               IF( N.GT.4 )
     $            CALL ZCOPY( N-4, WORK( N+1 ), 1, A( 2, 4 ), LDA+1 )
            END IF
            DO 100 J = 2, N - 1
               A( 1, J ) = Y
               A( J, N ) = Y
  100       CONTINUE
            A( 1, N ) = Z
         ELSE
            IF( N.GT.3 ) THEN
               CALL ZCOPY( N-3, WORK, 1, A( 3, 2 ), LDA+1 )
               IF( N.GT.4 )
     $            CALL ZCOPY( N-4, WORK( N+1 ), 1, A( 4, 2 ), LDA+1 )
            END IF
            DO 110 J = 2, N - 1
               A( J, 1 ) = Y
               A( N, J ) = Y
  110       CONTINUE
            A( N, 1 ) = Z
         END IF
*
*        Fill in the zeros using Givens rotations.
*
         IF( UPPER ) THEN
            DO 120 J = 1, N - 1
               RA = A( J, J+1 )
               RB = 2.0D0
               CALL ZROTG( RA, RB, C, S )
*
*              Multiply by [ c  s; -conjg(s)  c] on the left.
*
               IF( N.GT.J+1 )
     $            CALL ZROT( N-J-1, A( J, J+2 ), LDA, A( J+1, J+2 ),
     $                       LDA, C, S )
*
*              Multiply by [-c -s;  conjg(s) -c] on the right.
*
               IF( J.GT.1 )
     $            CALL ZROT( J-1, A( 1, J+1 ), 1, A( 1, J ), 1, -C, -S )
*
*              Negate A(J,J+1).
*
               A( J, J+1 ) = -A( J, J+1 )
  120       CONTINUE
         ELSE
            DO 130 J = 1, N - 1
               RA = A( J+1, J )
               RB = 2.0D0
               CALL ZROTG( RA, RB, C, S )
               S = DCONJG( S )