clattb.f

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

      SUBROUTINE CLATTB( IMAT, UPLO, TRANS, DIAG, ISEED, N, KD, AB,
     $                   LDAB, 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, KD, LDAB, N
*     ..
*     .. Array Arguments ..
      INTEGER            ISEED( 4 )
      REAL               RWORK( * )
      COMPLEX            AB( LDAB, * ), B( * ), WORK( * )
*     ..
*
*  Purpose
*  =======
*
*  CLATTB 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 (= 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
*          CLATMS).  Modified on exit.
*
*  N       (input) INTEGER
*          The order of the matrix to be generated.
*
*  KD      (input) INTEGER
*          The number of superdiagonals or subdiagonals of the banded
*          triangular matrix A.  KD >= 0.
*
*  AB      (output) COMPLEX array, dimension (LDAB,N)
*          The upper or lower triangular banded matrix A, stored in the
*          first KD+1 rows of AB.  Let j be a column of A, 1<=j<=n.
*          If UPLO = 'U', AB(kd+1+i-j,j) = A(i,j) for max(1,j-kd)<=i<=j.
*          If UPLO = 'L', AB(1+i-j,j)    = A(i,j) for j<=i<=min(n,j+kd).
*
*  LDAB    (input) INTEGER
*          The leading dimension of the array AB.  LDAB >= KD+1.
*
*  B       (workspace) COMPLEX array, dimension (N)
*
*  WORK    (workspace) COMPLEX array, dimension (2*N)
*
*  RWORK   (workspace) REAL array, dimension (N)
*
*  INFO    (output) INTEGER
*          = 0:  successful exit
*          < 0:  if INFO = -i, the i-th argument had an illegal value
*
*  =====================================================================
*
*     .. Parameters ..
      REAL               ONE, TWO, ZERO
      PARAMETER          ( ONE = 1.0E+0, TWO = 2.0E+0, ZERO = 0.0E+0 )
*     ..
*     .. Local Scalars ..
      LOGICAL            UPPER
      CHARACTER          DIST, PACKIT, TYPE
      CHARACTER*3        PATH
      INTEGER            I, IOFF, IY, J, JCOUNT, KL, KU, LENJ, MODE
      REAL               ANORM, BIGNUM, BNORM, BSCAL, CNDNUM, REXP,
     $                   SFAC, SMLNUM, TEXP, TLEFT, TNORM, TSCAL, ULP,
     $                   UNFL
      COMPLEX            PLUS1, PLUS2, STAR1
*     ..
*     .. External Functions ..
      LOGICAL            LSAME
      INTEGER            ICAMAX
      REAL               SLAMCH, SLARND
      COMPLEX            CLARND
      EXTERNAL           LSAME, ICAMAX, SLAMCH, SLARND, CLARND
*     ..
*     .. External Subroutines ..
      EXTERNAL           CCOPY, CLARNV, CLATB4, CLATMS, CSSCAL, CSWAP,
     $                   SLABAD, SLARNV
*     ..
*     .. Intrinsic Functions ..
      INTRINSIC          ABS, CMPLX, MAX, MIN, REAL, SQRT
*     ..
*     .. Executable Statements ..
*
      PATH( 1: 1 ) = 'Complex precision'
      PATH( 2: 3 ) = 'TB'
      UNFL = SLAMCH( 'Safe minimum' )
      ULP = SLAMCH( 'Epsilon' )*SLAMCH( 'Base' )
      SMLNUM = UNFL
      BIGNUM = ( ONE-ULP ) / SMLNUM
      CALL SLABAD( SMLNUM, BIGNUM )
      IF( ( IMAT.GE.6 .AND. IMAT.LE.9 ) .OR. IMAT.EQ.17 ) THEN
         DIAG = 'U'
      ELSE
         DIAG = 'N'
      END IF
      INFO = 0
*
*     Quick return if N.LE.0.
*
      IF( N.LE.0 )
     $   RETURN
*
*     Call CLATB4 to set parameters for CLATMS.
*
      UPPER = LSAME( UPLO, 'U' )
      IF( UPPER ) THEN
         CALL CLATB4( PATH, IMAT, N, N, TYPE, KL, KU, ANORM, MODE,
     $                CNDNUM, DIST )
         KU = KD
         IOFF = 1 + MAX( 0, KD-N+1 )
         KL = 0
         PACKIT = 'Q'
      ELSE
         CALL CLATB4( PATH, -IMAT, N, N, TYPE, KL, KU, ANORM, MODE,
     $                CNDNUM, DIST )
         KL = KD
         IOFF = 1
         KU = 0
         PACKIT = 'B'
      END IF
*
*     IMAT <= 5:  Non-unit triangular matrix
*
      IF( IMAT.LE.5 ) THEN
         CALL CLATMS( N, N, DIST, ISEED, TYPE, RWORK, MODE, CNDNUM,
     $                ANORM, KL, KU, PACKIT, AB( IOFF, 1 ), LDAB, WORK,
     $                INFO )
*
*     IMAT > 5:  Unit triangular matrix
*     The diagonal is deliberately set to something other than 1.
*
*     IMAT = 6:  Matrix is the identity
*
      ELSE IF( IMAT.EQ.6 ) THEN
         IF( UPPER ) THEN
            DO 20 J = 1, N
               DO 10 I = MAX( 1, KD+2-J ), KD
                  AB( I, J ) = ZERO
   10          CONTINUE
               AB( KD+1, J ) = J
   20       CONTINUE
         ELSE
            DO 40 J = 1, N
               AB( 1, J ) = J
               DO 30 I = 2, MIN( KD+1, N-J+1 )
                  AB( I, J ) = ZERO
   30          CONTINUE
   40       CONTINUE
         END IF
*
*     IMAT > 6:  Non-trivial unit triangular matrix
*
*     A unit triangular matrix T with condition CNDNUM is formed.
*     In this version, T only has bandwidth 2, the rest of it is zero.
*
      ELSE IF( IMAT.LE.9 ) THEN
         TNORM = SQRT( CNDNUM )
*
*        Initialize AB to zero.
*
         IF( UPPER ) THEN
            DO 60 J = 1, N
               DO 50 I = MAX( 1, KD+2-J ), KD
                  AB( I, J ) = ZERO
   50          CONTINUE
               AB( KD+1, J ) = REAL( J )
   60       CONTINUE
         ELSE
            DO 80 J = 1, N
               DO 70 I = 2, MIN( KD+1, N-J+1 )
                  AB( I, J ) = ZERO
   70          CONTINUE
               AB( 1, J ) = REAL( J )
   80       CONTINUE
         END IF
*
*        Special case:  T is tridiagonal.  Set every other offdiagonal
*        so that the matrix has norm TNORM+1.
*
         IF( KD.EQ.1 ) THEN
            IF( UPPER ) THEN
               AB( 1, 2 ) = TNORM*CLARND( 5, ISEED )
               LENJ = ( N-3 ) / 2
               CALL CLARNV( 2, ISEED, LENJ, WORK )
               DO 90 J = 1, LENJ
                  AB( 1, 2*( J+1 ) ) = TNORM*WORK( J )
   90          CONTINUE
            ELSE
               AB( 2, 1 ) = TNORM*CLARND( 5, ISEED )
               LENJ = ( N-3 ) / 2
               CALL CLARNV( 2, ISEED, LENJ, WORK )
               DO 100 J = 1, LENJ
                  AB( 2, 2*J+1 ) = TNORM*WORK( J )
  100          CONTINUE
            END IF
         ELSE IF( KD.GT.1 ) THEN
*
*           Form a unit triangular matrix T with condition CNDNUM.  T is
*           given by
*                   | 1   +   *                      |
*                   |     1   +                      |
*               T = |         1   +   *              |
*                   |             1   +              |
*                   |                 1   +   *      |
*                   |                     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.
*
*        The two offdiagonals of T are stored in WORK.
*
            STAR1 = TNORM*CLARND( 5, ISEED )
            SFAC = SQRT( TNORM )
            PLUS1 = SFAC*CLARND( 5, ISEED )
            DO 110 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
*
*                 Generate a new *-value with norm between sqrt(TNORM)
*                 and TNORM.
*
                  REXP = SLARND( 2, ISEED )
                  IF( REXP.LT.ZERO ) THEN
                     STAR1 = -SFAC**( ONE-REXP )*CLARND( 5, ISEED )
                  ELSE
                     STAR1 = SFAC**( ONE+REXP )*CLARND( 5, ISEED )
                  END IF
               END IF
  110       CONTINUE
*
*           Copy the tridiagonal T to AB.
*
            IF( UPPER ) THEN
               CALL CCOPY( N-1, WORK, 1, AB( KD, 2 ), LDAB )
               CALL CCOPY( N-2, WORK( N+1 ), 1, AB( KD-1, 3 ), LDAB )
            ELSE
               CALL CCOPY( N-1, WORK, 1, AB( 2, 1 ), LDAB )
               CALL CCOPY( N-2, WORK( N+1 ), 1, AB( 3, 1 ), LDAB )
            END IF
         END IF
*
*     IMAT > 9:  Pathological test cases.  These triangular matrices
*     are badly scaled or badly conditioned, so when used in solving a
*     triangular system they may cause overflow in the solution vector.
*
      ELSE IF( IMAT.EQ.10 ) THEN
*
*        Type 10:  Generate a triangular matrix with elements between
*        -1 and 1. Give the diagonal norm 2 to make it well-conditioned.
*        Make the right hand side large so that it requires scaling.
*
         IF( UPPER ) THEN
            DO 120 J = 1, N
               LENJ = MIN( J-1, KD )
               CALL CLARNV( 4, ISEED, LENJ, AB( KD+1-LENJ, J ) )
               AB( KD+1, J ) = CLARND( 5, ISEED )*TWO
  120       CONTINUE
         ELSE
            DO 130 J = 1, N
               LENJ = MIN( N-J, KD )
               IF( LENJ.GT.0 )
     $            CALL CLARNV( 4, ISEED, LENJ, AB( 2, J ) )
               AB( 1, J ) = CLARND( 5, ISEED )*TWO
  130       CONTINUE
         END IF
*
*        Set the right hand side so that the largest value is BIGNUM.
*
         CALL CLARNV( 2, ISEED, N, B )
         IY = ICAMAX( N, B, 1 )
         BNORM = ABS( B( IY ) )
         BSCAL = BIGNUM / MAX( ONE, BNORM )
         CALL CSSCAL( N, BSCAL, B, 1 )
*
      ELSE IF( IMAT.EQ.11 ) THEN