slattb.f

Go to the documentation of this file.

*> \brief \b SLATTB
*
*  =========== DOCUMENTATION ===========
*
* Online html documentation available at 
*            http://www.netlib.org/lapack/explore-html/ 
*
*  Definition:
*  ===========
*
*       SUBROUTINE SLATTB( IMAT, UPLO, TRANS, DIAG, ISEED, N, KD, AB,
*                          LDAB, B, WORK, INFO )
* 
*       .. Scalar Arguments ..
*       CHARACTER          DIAG, TRANS, UPLO
*       INTEGER            IMAT, INFO, KD, LDAB, N
*       ..
*       .. Array Arguments ..
*       INTEGER            ISEED( 4 )
*       REAL               AB( LDAB, * ), B( * ), WORK( * )
*       ..
*  
*
*> \par Purpose:
*  =============
*>
*> \verbatim
*>
*> SLATTB 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.
*> \endverbatim
*
*  Arguments:
*  ==========
*
*> \param[in] IMAT
*> \verbatim
*>          IMAT is INTEGER
*>          An integer key describing which matrix to generate for this
*>          path.
*> \endverbatim
*>
*> \param[in] UPLO
*> \verbatim
*>          UPLO is CHARACTER*1
*>          Specifies whether the matrix A will be upper or lower
*>          triangular.
*>          = 'U':  Upper triangular
*>          = 'L':  Lower triangular
*> \endverbatim
*>
*> \param[in] TRANS
*> \verbatim
*>          TRANS is CHARACTER*1
*>          Specifies whether the matrix or its transpose will be used.
*>          = 'N':  No transpose
*>          = 'T':  Transpose
*>          = 'C':  Conjugate transpose (= transpose)
*> \endverbatim
*>
*> \param[out] DIAG
*> \verbatim
*>          DIAG is CHARACTER*1
*>          Specifies whether or not the matrix A is unit triangular.
*>          = 'N':  Non-unit triangular
*>          = 'U':  Unit triangular
*> \endverbatim
*>
*> \param[in,out] ISEED
*> \verbatim
*>          ISEED is INTEGER array, dimension (4)
*>          The seed vector for the random number generator (used in
*>          SLATMS).  Modified on exit.
*> \endverbatim
*>
*> \param[in] N
*> \verbatim
*>          N is INTEGER
*>          The order of the matrix to be generated.
*> \endverbatim
*>
*> \param[in] KD
*> \verbatim
*>          KD is INTEGER
*>          The number of superdiagonals or subdiagonals of the banded
*>          triangular matrix A.  KD >= 0.
*> \endverbatim
*>
*> \param[out] AB
*> \verbatim
*>          AB is REAL 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).
*> \endverbatim
*>
*> \param[in] LDAB
*> \verbatim
*>          LDAB is INTEGER
*>          The leading dimension of the array AB.  LDAB >= KD+1.
*> \endverbatim
*>
*> \param[out] B
*> \verbatim
*>          B is REAL array, dimension (N)
*> \endverbatim
*>
*> \param[out] WORK
*> \verbatim
*>          WORK is REAL array, dimension (2*N)
*> \endverbatim
*>
*> \param[out] INFO
*> \verbatim
*>          INFO is INTEGER
*>          = 0:  successful exit
*>          < 0: if INFO = -k, the k-th argument had an illegal value
*> \endverbatim
*
*  Authors:
*  ========
*
*> \author Univ. of Tennessee 
*> \author Univ. of California Berkeley 
*> \author Univ. of Colorado Denver 
*> \author NAG Ltd. 
*
*> \date November 2011
*
*> \ingroup single_lin
*
*  =====================================================================
      SUBROUTINE slattb( IMAT, UPLO, TRANS, DIAG, ISEED, N, KD, AB,
     $                   ldab, b, work, info )
*
*  -- LAPACK test routine (version 3.4.0) --
*  -- LAPACK is a software package provided by Univ. of Tennessee,    --
*  -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
*     November 2011
*
*     .. Scalar Arguments ..
      CHARACTER          DIAG, TRANS, UPLO
      INTEGER            IMAT, INFO, KD, LDAB, N
*     ..
*     .. Array Arguments ..
      INTEGER            ISEED( 4 )
      REAL               AB( ldab, * ), B( * ), WORK( * )
*     ..
*
*  =====================================================================
*
*     .. 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, PLUS1,
     $                   plus2, rexp, sfac, smlnum, star1, texp, tleft,
     $                   tnorm, tscal, ulp, unfl
*     ..
*     .. External Functions ..
      LOGICAL            LSAME
      INTEGER            ISAMAX
      REAL               SLAMCH, SLARND
      EXTERNAL           lsame, isamax, slamch, slarnd
*     ..
*     .. External Subroutines ..
      EXTERNAL           scopy, slabad, slarnv, slatb4, slatms, sscal,
     $                   sswap
*     ..
*     .. Intrinsic Functions ..
      INTRINSIC          abs, max, min, REAL, SIGN, SQRT
*     ..
*     .. Executable Statements ..
*
      path( 1: 1 ) = 'Single 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 SLATB4 to set parameters for SLATMS.
*
      upper = lsame( uplo, 'U' )
      IF( upper ) THEN
         CALL slatb4( 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 slatb4( 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 slatms( n, n, dist, iseed, TYPE, B, 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 ) = sign( tnorm, slarnd( 2, iseed ) )
               lenj = ( n-3 ) / 2
               CALL slarnv( 2, iseed, lenj, work )
               DO 90 j = 1, lenj
                  ab( 1, 2*( j+1 ) ) = tnorm*work( j )
   90          CONTINUE
            ELSE
               ab( 2, 1 ) = sign( tnorm, slarnd( 2, iseed ) )
               lenj = ( n-3 ) / 2
               CALL slarnv( 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 = sign( tnorm, slarnd( 2, iseed ) )
            sfac = sqrt( tnorm )
            plus1 = sign( sfac, slarnd( 2, 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 )
                  ELSE
                     star1 = sfac**( one+rexp )
                  END IF
               END IF
  110       CONTINUE
*
*           Copy the tridiagonal T to AB.
*
            IF( upper ) THEN
               CALL scopy( n-1, work, 1, ab( kd, 2 ), ldab )
               CALL scopy( n-2, work( n+1 ), 1, ab( kd-1, 3 ), ldab )
            ELSE
               CALL scopy( n-1, work, 1, ab( 2, 1 ), ldab )
               CALL scopy( 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, kd+1 )
               CALL slarnv( 2, iseed, lenj, ab( kd+2-lenj, j ) )
               ab( kd+1, j ) = sign( two, ab( kd+1, j ) )
  120       CONTINUE
         ELSE
            DO 130 j = 1, n
               lenj = min( n-j+1, kd+1 )
               IF( lenj.GT.0 )
     $            CALL slarnv( 2, iseed, lenj, ab( 1, j ) )
               ab( 1, j ) = sign( two, ab( 1, j ) )
  130       CONTINUE
         END IF
*
*        Set the right hand side so that the largest value is BIGNUM.
*
         CALL slarnv( 2, iseed, n, b )
         iy = isamax( n, b, 1 )
         bnorm = abs( b( iy ) )
         bscal = bignum / max( one, bnorm )
         CALL sscal( n, bscal, b, 1 )
*
      ELSE IF( imat.EQ.11 ) THEN
*
*        Type 11:  Make the first diagonal element in the solve small to
*        cause immediate overflow when dividing by T(j,j).
*        In type 11, the offdiagonal elements are small (CNORM(j) < 1).
*
         CALL slarnv( 2, iseed, n, b )
         tscal = one / REAL( kd+1 )
         IF( upper ) THEN
            DO 140 j = 1, n
               lenj = min( j, kd+1 )
               CALL slarnv( 2, iseed, lenj, ab( kd+2-lenj, j ) )
               CALL sscal( lenj-1, tscal, ab( kd+2-lenj, j ), 1 )
               ab( kd+1, j ) = sign( one, ab( kd+1, j ) )
  140       CONTINUE
            ab( kd+1, n ) = smlnum*ab( kd+1, n )
         ELSE
            DO 150 j = 1, n
               lenj = min( n-j+1, kd+1 )
               CALL slarnv( 2, iseed, lenj, ab( 1, j ) )
               IF( lenj.GT.1 )
     $            CALL sscal( lenj-1, tscal, ab( 2, j ), 1 )
               ab( 1, j ) = sign( one, ab( 1, j ) )
  150       CONTINUE
            ab( 1, 1 ) = smlnum*ab( 1, 1 )
         END IF
*
      ELSE IF( imat.EQ.12 ) THEN
*
*        Type 12:  Make the first diagonal element in the solve small to
*        cause immediate overflow when dividing by T(j,j).
*        In type 12, the offdiagonal elements are O(1) (CNORM(j) > 1).
*
         CALL slarnv( 2, iseed, n, b )
         IF( upper ) THEN
            DO 160 j = 1, n
               lenj = min( j, kd+1 )
               CALL slarnv( 2, iseed, lenj, ab( kd+2-lenj, j ) )
               ab( kd+1, j ) = sign( one, ab( kd+1, j ) )
  160       CONTINUE
            ab( kd+1, n ) = smlnum*ab( kd+1, n )
         ELSE
            DO 170 j = 1, n
               lenj = min( n-j+1, kd+1 )
               CALL slarnv( 2, iseed, lenj, ab( 1, j ) )
               ab( 1, j ) = sign( one, ab( 1, j ) )
  170       CONTINUE
            ab( 1, 1 ) = smlnum*ab( 1, 1 )
         END IF
*
      ELSE IF( imat.EQ.13 ) THEN
*
*        Type 13:  T is diagonal with small numbers on the diagonal to
*        make the growth factor underflow, but a small right hand side
*        chosen so that the solution does not overflow.
*
         IF( upper ) THEN
            jcount = 1
            DO 190 j = n, 1, -1
               DO 180 i = max( 1, kd+1-( j-1 ) ), kd
                  ab( i, j ) = zero
  180          CONTINUE
               IF( jcount.LE.2 ) THEN
                  ab( kd+1, j ) = smlnum
               ELSE
                  ab( kd+1, j ) = one
               END IF
               jcount = jcount + 1
               IF( jcount.GT.4 )
     $            jcount = 1
  190       CONTINUE
         ELSE
            jcount = 1
            DO 210 j = 1, n
               DO 200 i = 2, min( n-j+1, kd+1 )
                  ab( i, j ) = zero
  200          CONTINUE
               IF( jcount.LE.2 ) THEN
                  ab( 1, j ) = smlnum
               ELSE
                  ab( 1, j ) = one
               END IF
               jcount = jcount + 1
               IF( jcount.GT.4 )
     $            jcount = 1
  210       CONTINUE
         END IF
*
*        Set the right hand side alternately zero and small.
*
         IF( upper ) THEN
            b( 1 ) = zero
            DO 220 i = n, 2, -2
               b( i ) = zero
               b( i-1 ) = smlnum
  220       CONTINUE
         ELSE
            b( n ) = zero
            DO 230 i = 1, n - 1, 2
               b( i ) = zero
               b( i+1 ) = smlnum
  230       CONTINUE
         END IF
*
      ELSE IF( imat.EQ.14 ) THEN
*
*        Type 14:  Make the diagonal elements small to cause gradual
*        overflow when dividing by T(j,j).  To control the amount of
*        scaling needed, the matrix is bidiagonal.
*
         texp = one / REAL( kd+1 )
         tscal = smlnum**texp
         CALL slarnv( 2, iseed, n, b )
         IF( upper ) THEN
            DO 250 j = 1, n
               DO 240 i = max( 1, kd+2-j ), kd
                  ab( i, j ) = zero
  240          CONTINUE
               IF( j.GT.1 .AND. kd.GT.0 )
     $            ab( kd, j ) = -one
               ab( kd+1, j ) = tscal
  250       CONTINUE
            b( n ) = one
         ELSE
            DO 270 j = 1, n
               DO 260 i = 3, min( n-j+1, kd+1 )
                  ab( i, j ) = zero
  260          CONTINUE
               IF( j.LT.n .AND. kd.GT.0 )
     $            ab( 2, j ) = -one
               ab( 1, j ) = tscal
  270       CONTINUE
            b( 1 ) = one
         END IF
*
      ELSE IF( imat.EQ.15 ) THEN
*
*        Type 15:  One zero diagonal element.
*
         iy = n / 2 + 1
         IF( upper ) THEN
            DO 280 j = 1, n
               lenj = min( j, kd+1 )
               CALL slarnv( 2, iseed, lenj, ab( kd+2-lenj, j ) )
               IF( j.NE.iy ) THEN
                  ab( kd+1, j ) = sign( two, ab( kd+1, j ) )
               ELSE
                  ab( kd+1, j ) = zero
               END IF
  280       CONTINUE
         ELSE
            DO 290 j = 1, n
               lenj = min( n-j+1, kd+1 )
               CALL slarnv( 2, iseed, lenj, ab( 1, j ) )
               IF( j.NE.iy ) THEN
                  ab( 1, j ) = sign( two, ab( 1, j ) )
               ELSE
                  ab( 1, j ) = zero
               END IF
  290       CONTINUE
         END IF
         CALL slarnv( 2, iseed, n, b )
         CALL sscal( n, two, b, 1 )
*
      ELSE IF( imat.EQ.16 ) THEN
*
*        Type 16:  Make the offdiagonal elements large to cause overflow
*        when adding a column of T.  In the non-transposed case, the
*        matrix is constructed to cause overflow when adding a column in
*        every other step.
*
         tscal = unfl / ulp
         tscal = ( one-ulp ) / tscal
         DO 310 j = 1, n
            DO 300 i = 1, kd + 1
               ab( i, j ) = zero
  300       CONTINUE
  310    CONTINUE
         texp = one
         IF( kd.GT.0 ) THEN
            IF( upper ) THEN
               DO 330 j = n, 1, -kd
                  DO 320 i = j, max( 1, j-kd+1 ), -2
                     ab( 1+( j-i ), i ) = -tscal / REAL( kd+2 )
                     ab( kd+1, i ) = one
                     b( i ) = texp*( one-ulp )
                     IF( i.GT.max( 1, j-kd+1 ) ) THEN
                        ab( 2+( j-i ), i-1 ) = -( tscal / REAL( KD+2 ) )
     $                                          / REAL( kd+3 )
                        ab( kd+1, i-1 ) = one
                        b( i-1 ) = texp*REAL( ( kd+1 )*( kd+1 )+kd )
                     END IF
                     texp = texp*two
  320             CONTINUE
                  b( max( 1, j-kd+1 ) ) = ( REAL( KD+2 ) /
     $                                    REAL( KD+3 ) )*tscal
  330          CONTINUE
            ELSE
               DO 350 j = 1, n, kd
                  texp = one
                  lenj = min( kd+1, n-j+1 )
                  DO 340 i = j, min( n, j+kd-1 ), 2
                     ab( lenj-( i-j ), j ) = -tscal / REAL( kd+2 )
                     ab( 1, j ) = one
                     b( j ) = texp*( one-ulp )
                     IF( i.LT.min( n, j+kd-1 ) ) THEN
                        ab( lenj-( i-j+1 ), i+1 ) = -( tscal /
     $                     REAL( KD+2 ) ) / REAL( KD+3 )
                        ab( 1, i+1 ) = one
                        b( i+1 ) = texp*REAL( ( kd+1 )*( kd+1 )+kd )
                     END IF
                     texp = texp*two
  340             CONTINUE
                  b( min( n, j+kd-1 ) ) = ( REAL( KD+2 ) /
     $                                    REAL( KD+3 ) )*tscal
  350          CONTINUE
            END IF
         ELSE
            DO 360 j = 1, n
               ab( 1, j ) = one
               b( j ) = REAL( j )
  360       CONTINUE
         END IF
*
      ELSE IF( imat.EQ.17 ) THEN
*
*        Type 17:  Generate a unit triangular matrix with elements
*        between -1 and 1, and make the right hand side large so that it
*        requires scaling.
*
         IF( upper ) THEN
            DO 370 j = 1, n
               lenj = min( j-1, kd )
               CALL slarnv( 2, iseed, lenj, ab( kd+1-lenj, j ) )
               ab( kd+1, j ) = REAL( j )
  370       CONTINUE
         ELSE
            DO 380 j = 1, n
               lenj = min( n-j, kd )
               IF( lenj.GT.0 )
     $            CALL slarnv( 2, iseed, lenj, ab( 2, j ) )
               ab( 1, j ) = REAL( j )
  380       CONTINUE
         END IF
*
*        Set the right hand side so that the largest value is BIGNUM.
*
         CALL slarnv( 2, iseed, n, b )
         iy = isamax( n, b, 1 )
         bnorm = abs( b( iy ) )
         bscal = bignum / max( one, bnorm )
         CALL sscal( n, bscal, b, 1 )
*
      ELSE IF( imat.EQ.18 ) THEN
*
*        Type 18:  Generate a triangular matrix with elements between
*        BIGNUM/KD and BIGNUM so that at least one of the column
*        norms will exceed BIGNUM.
*
         tleft = bignum / max( one, REAL( KD ) )
         tscal = bignum*( REAL( KD ) / REAL( KD+1 ) )
         IF( upper ) THEN
            DO 400 j = 1, n
               lenj = min( j, kd+1 )
               CALL slarnv( 2, iseed, lenj, ab( kd+2-lenj, j ) )
               DO 390 i = kd + 2 - lenj, kd + 1
                  ab( i, j ) = sign( tleft, ab( i, j ) ) +
     $                         tscal*ab( i, j )
  390          CONTINUE
  400       CONTINUE
         ELSE
            DO 420 j = 1, n
               lenj = min( n-j+1, kd+1 )
               CALL slarnv( 2, iseed, lenj, ab( 1, j ) )
               DO 410 i = 1, lenj
                  ab( i, j ) = sign( tleft, ab( i, j ) ) +
     $                         tscal*ab( i, j )
  410          CONTINUE
  420       CONTINUE
         END IF
         CALL slarnv( 2, iseed, n, b )
         CALL sscal( n, two, b, 1 )
      END IF
*
*     Flip the matrix if the transpose will be used.
*
      IF( .NOT.lsame( trans, 'N' ) ) THEN
         IF( upper ) THEN
            DO 430 j = 1, n / 2
               lenj = min( n-2*j+1, kd+1 )
               CALL sswap( lenj, ab( kd+1, j ), ldab-1,
     $                     ab( kd+2-lenj, n-j+1 ), -1 )
  430       CONTINUE
         ELSE
            DO 440 j = 1, n / 2
               lenj = min( n-2*j+1, kd+1 )
               CALL sswap( lenj, ab( 1, j ), 1, ab( lenj, n-j+2-lenj ),
     $                     -ldab+1 )
  440       CONTINUE
         END IF
      END IF
*
      RETURN
*
*     End of SLATTB
*
      END