share/lapack/blas/fortran/ztbsv.f

   1       SUBROUTINE ZTBSV ( UPLO, TRANS, DIAG, N, K, A, LDA, X, INCX )
   2 *     .. Scalar Arguments ..
   3       INTEGER            INCX, K, LDA, N
   4       CHARACTER*1        DIAG, TRANS, UPLO
   5 *     .. Array Arguments ..
   6       COMPLEX*16         A( LDA, * ), X( * )
   7 *     ..
   8 *
   9 *  Purpose
  10 *  =======
  11 *
  12 *  ZTBSV  solves one of the systems of equations
  13 *
  14 *     A*x = b,   or   A'*x = b,   or   conjg( A' )*x = b,
  15 *
  16 *  where b and x are n element vectors and A is an n by n unit, or
  17 *  non-unit, upper or lower triangular band matrix, with ( k + 1 )
  18 *  diagonals.
  19 *
  20 *  No test for singularity or near-singularity is included in this
  21 *  routine. Such tests must be performed before calling this routine.
  22 *
  23 *  Parameters
  24 *  ==========
  25 *
  26 *  UPLO   - CHARACTER*1.
  27 *           On entry, UPLO specifies whether the matrix is an upper or
  28 *           lower triangular matrix as follows:
  29 *
  30 *              UPLO = 'U' or 'u'   A is an upper triangular matrix.
  31 *
  32 *              UPLO = 'L' or 'l'   A is a lower triangular matrix.
  33 *
  34 *           Unchanged on exit.
  35 *
  36 *  TRANS  - CHARACTER*1.
  37 *           On entry, TRANS specifies the equations to be solved as
  38 *           follows:
  39 *
  40 *              TRANS = 'N' or 'n'   A*x = b.
  41 *
  42 *              TRANS = 'T' or 't'   A'*x = b.
  43 *
  44 *              TRANS = 'C' or 'c'   conjg( A' )*x = b.
  45 *
  46 *           Unchanged on exit.
  47 *
  48 *  DIAG   - CHARACTER*1.
  49 *           On entry, DIAG specifies whether or not A is unit
  50 *           triangular as follows:
  51 *
  52 *              DIAG = 'U' or 'u'   A is assumed to be unit triangular.
  53 *
  54 *              DIAG = 'N' or 'n'   A is not assumed to be unit
  55 *                                  triangular.
  56 *
  57 *           Unchanged on exit.
  58 *
  59 *  N      - INTEGER.
  60 *           On entry, N specifies the order of the matrix A.
  61 *           N must be at least zero.
  62 *           Unchanged on exit.
  63 *
  64 *  K      - INTEGER.
  65 *           On entry with UPLO = 'U' or 'u', K specifies the number of
  66 *           super-diagonals of the matrix A.
  67 *           On entry with UPLO = 'L' or 'l', K specifies the number of
  68 *           sub-diagonals of the matrix A.
  69 *           K must satisfy  0 .le. K.
  70 *           Unchanged on exit.
  71 *
  72 *  A      - COMPLEX*16       array of DIMENSION ( LDA, n ).
  73 *           Before entry with UPLO = 'U' or 'u', the leading ( k + 1 )
  74 *           by n part of the array A must contain the upper triangular
  75 *           band part of the matrix of coefficients, supplied column by
  76 *           column, with the leading diagonal of the matrix in row
  77 *           ( k + 1 ) of the array, the first super-diagonal starting at
  78 *           position 2 in row k, and so on. The top left k by k triangle
  79 *           of the array A is not referenced.
  80 *           The following program segment will transfer an upper
  81 *           triangular band matrix from conventional full matrix storage
  82 *           to band storage:
  83 *
  84 *                 DO 20, J = 1, N
  85 *                    M = K + 1 - J
  86 *                    DO 10, I = MAX( 1, J - K ), J
  87 *                       A( M + I, J ) = matrix( I, J )
  88 *              10    CONTINUE
  89 *              20 CONTINUE
  90 *
  91 *           Before entry with UPLO = 'L' or 'l', the leading ( k + 1 )
  92 *           by n part of the array A must contain the lower triangular
  93 *           band part of the matrix of coefficients, supplied column by
  94 *           column, with the leading diagonal of the matrix in row 1 of
  95 *           the array, the first sub-diagonal starting at position 1 in
  96 *           row 2, and so on. The bottom right k by k triangle of the
  97 *           array A is not referenced.
  98 *           The following program segment will transfer a lower
  99 *           triangular band matrix from conventional full matrix storage
 100 *           to band storage:
 101 *
 102 *                 DO 20, J = 1, N
 103 *                    M = 1 - J
 104 *                    DO 10, I = J, MIN( N, J + K )
 105 *                       A( M + I, J ) = matrix( I, J )
 106 *              10    CONTINUE
 107 *              20 CONTINUE
 108 *
 109 *           Note that when DIAG = 'U' or 'u' the elements of the array A
 110 *           corresponding to the diagonal elements of the matrix are not
 111 *           referenced, but are assumed to be unity.
 112 *           Unchanged on exit.
 113 *
 114 *  LDA    - INTEGER.
 115 *           On entry, LDA specifies the first dimension of A as declared
 116 *           in the calling (sub) program. LDA must be at least
 117 *           ( k + 1 ).
 118 *           Unchanged on exit.
 119 *
 120 *  X      - COMPLEX*16       array of dimension at least
 121 *           ( 1 + ( n - 1 )*abs( INCX ) ).
 122 *           Before entry, the incremented array X must contain the n
 123 *           element right-hand side vector b. On exit, X is overwritten
 124 *           with the solution vector x.
 125 *
 126 *  INCX   - INTEGER.
 127 *           On entry, INCX specifies the increment for the elements of
 128 *           X. INCX must not be zero.
 129 *           Unchanged on exit.
 130 *
 131 *
 132 *  Level 2 Blas routine.
 133 *
 134 *  -- Written on 22-October-1986.
 135 *     Jack Dongarra, Argonne National Lab.
 136 *     Jeremy Du Croz, Nag Central Office.
 137 *     Sven Hammarling, Nag Central Office.
 138 *     Richard Hanson, Sandia National Labs.
 139 *
 140 *
 141 *     .. Parameters ..
 142       COMPLEX*16         ZERO
 143       PARAMETER        ( ZERO = ( 0.0D+0, 0.0D+0 ) )
 144 *     .. Local Scalars ..
 145       COMPLEX*16         TEMP
 146       INTEGER            I, INFO, IX, J, JX, KPLUS1, KX, L
 147       LOGICAL            NOCONJ, NOUNIT
 148 *     .. External Functions ..
 149       LOGICAL            LSAME
 150       EXTERNAL           LSAME
 151 *     .. External Subroutines ..
 152       EXTERNAL           XERBLA
 153 *     .. Intrinsic Functions ..
 154       INTRINSIC          DCONJG, MAX, MIN
 155 *     ..
 156 *     .. Executable Statements ..
 157 *
 158 *     Test the input parameters.
 159 *
 160       INFO = 0
 161       IF     ( .NOT.LSAME( UPLO , 'U' ).AND.
 162      $         .NOT.LSAME( UPLO , 'L' )      )THEN
 163          INFO = 1
 164       ELSE IF( .NOT.LSAME( TRANS, 'N' ).AND.
 165      $         .NOT.LSAME( TRANS, 'T' ).AND.
 166      $         .NOT.LSAME( TRANS, 'C' )      )THEN
 167          INFO = 2
 168       ELSE IF( .NOT.LSAME( DIAG , 'U' ).AND.
 169      $         .NOT.LSAME( DIAG , 'N' )      )THEN
 170          INFO = 3
 171       ELSE IF( N.LT.0 )THEN
 172          INFO = 4
 173       ELSE IF( K.LT.0 )THEN
 174          INFO = 5
 175       ELSE IF( LDA.LT.( K + 1 ) )THEN
 176          INFO = 7
 177       ELSE IF( INCX.EQ.0 )THEN
 178          INFO = 9
 179       END IF
 180       IF( INFO.NE.0 )THEN
 181          CALL XERBLA( 'ZTBSV ', INFO )
 182          RETURN
 183       END IF
 184 *
 185 *     Quick return if possible.
 186 *
 187       IF( N.EQ.0 )
 188      $   RETURN
 189 *
 190       NOCONJ = LSAME( TRANS, 'T' )
 191       NOUNIT = LSAME( DIAG , 'N' )
 192 *
 193 *     Set up the start point in X if the increment is not unity. This
 194 *     will be  ( N - 1 )*INCX  too small for descending loops.
 195 *
 196       IF( INCX.LE.0 )THEN
 197          KX = 1 - ( N - 1 )*INCX
 198       ELSE IF( INCX.NE.1 )THEN
 199          KX = 1
 200       END IF
 201 *
 202 *     Start the operations. In this version the elements of A are
 203 *     accessed by sequentially with one pass through A.
 204 *
 205       IF( LSAME( TRANS, 'N' ) )THEN
 206 *
 207 *        Form  x := inv( A )*x.
 208 *
 209          IF( LSAME( UPLO, 'U' ) )THEN
 210             KPLUS1 = K + 1
 211             IF( INCX.EQ.1 )THEN
 212                DO 20, J = N, 1, -1
 213                   IF( X( J ).NE.ZERO )THEN
 214                      L = KPLUS1 - J
 215                      IF( NOUNIT )
 216      $                  X( J ) = X( J )/A( KPLUS1, J )
 217                      TEMP = X( J )
 218                      DO 10, I = J - 1, MAX( 1, J - K ), -1
 219                         X( I ) = X( I ) - TEMP*A( L + I, J )
 220    10                CONTINUE
 221                   END IF
 222    20          CONTINUE
 223             ELSE
 224                KX = KX + ( N - 1 )*INCX
 225                JX = KX
 226                DO 40, J = N, 1, -1
 227                   KX = KX - INCX
 228                   IF( X( JX ).NE.ZERO )THEN
 229                      IX = KX
 230                      L  = KPLUS1 - J
 231                      IF( NOUNIT )
 232      $                  X( JX ) = X( JX )/A( KPLUS1, J )
 233                      TEMP = X( JX )
 234                      DO 30, I = J - 1, MAX( 1, J - K ), -1
 235                         X( IX ) = X( IX ) - TEMP*A( L + I, J )
 236                         IX      = IX      - INCX
 237    30                CONTINUE
 238                   END IF
 239                   JX = JX - INCX
 240    40          CONTINUE
 241             END IF
 242          ELSE
 243             IF( INCX.EQ.1 )THEN
 244                DO 60, J = 1, N
 245                   IF( X( J ).NE.ZERO )THEN
 246                      L = 1 - J
 247                      IF( NOUNIT )
 248      $                  X( J ) = X( J )/A( 1, J )
 249                      TEMP = X( J )
 250                      DO 50, I = J + 1, MIN( N, J + K )
 251                         X( I ) = X( I ) - TEMP*A( L + I, J )
 252    50                CONTINUE
 253                   END IF
 254    60          CONTINUE
 255             ELSE
 256                JX = KX
 257                DO 80, J = 1, N
 258                   KX = KX + INCX
 259                   IF( X( JX ).NE.ZERO )THEN
 260                      IX = KX
 261                      L  = 1  - J
 262                      IF( NOUNIT )
 263      $                  X( JX ) = X( JX )/A( 1, J )
 264                      TEMP = X( JX )
 265                      DO 70, I = J + 1, MIN( N, J + K )
 266                         X( IX ) = X( IX ) - TEMP*A( L + I, J )
 267                         IX      = IX      + INCX
 268    70                CONTINUE
 269                   END IF
 270                   JX = JX + INCX
 271    80          CONTINUE
 272             END IF
 273          END IF
 274       ELSE
 275 *
 276 *        Form  x := inv( A' )*x  or  x := inv( conjg( A') )*x.
 277 *
 278          IF( LSAME( UPLO, 'U' ) )THEN
 279             KPLUS1 = K + 1
 280             IF( INCX.EQ.1 )THEN
 281                DO 110, J = 1, N
 282                   TEMP = X( J )
 283                   L    = KPLUS1 - J
 284                   IF( NOCONJ )THEN
 285                      DO 90, I = MAX( 1, J - K ), J - 1
 286                         TEMP = TEMP - A( L + I, J )*X( I )
 287    90                CONTINUE
 288                      IF( NOUNIT )
 289      $                  TEMP = TEMP/A( KPLUS1, J )
 290                   ELSE
 291                      DO 100, I = MAX( 1, J - K ), J - 1
 292                         TEMP = TEMP - DCONJG( A( L + I, J ) )*X( I )
 293   100                CONTINUE
 294                      IF( NOUNIT )
 295      $                  TEMP = TEMP/DCONJG( A( KPLUS1, J ) )
 296                   END IF
 297                   X( J ) = TEMP
 298   110          CONTINUE
 299             ELSE
 300                JX = KX
 301                DO 140, J = 1, N
 302                   TEMP = X( JX )
 303                   IX   = KX
 304                   L    = KPLUS1  - J
 305                   IF( NOCONJ )THEN
 306                      DO 120, I = MAX( 1, J - K ), J - 1
 307                         TEMP = TEMP - A( L + I, J )*X( IX )
 308                         IX   = IX   + INCX
 309   120                CONTINUE
 310                      IF( NOUNIT )
 311      $                  TEMP = TEMP/A( KPLUS1, J )
 312                   ELSE
 313                      DO 130, I = MAX( 1, J - K ), J - 1
 314                         TEMP = TEMP - DCONJG( A( L + I, J ) )*X( IX )
 315                         IX   = IX   + INCX
 316   130                CONTINUE
 317                      IF( NOUNIT )
 318      $                  TEMP = TEMP/DCONJG( A( KPLUS1, J ) )
 319                   END IF
 320                   X( JX ) = TEMP
 321                   JX      = JX   + INCX
 322                   IF( J.GT.K )
 323      $               KX = KX + INCX
 324   140          CONTINUE
 325             END IF
 326          ELSE
 327             IF( INCX.EQ.1 )THEN
 328                DO 170, J = N, 1, -1
 329                   TEMP = X( J )
 330                   L    = 1      - J
 331                   IF( NOCONJ )THEN
 332                      DO 150, I = MIN( N, J + K ), J + 1, -1
 333                         TEMP = TEMP - A( L + I, J )*X( I )
 334   150                CONTINUE
 335                      IF( NOUNIT )
 336      $                  TEMP = TEMP/A( 1, J )
 337                   ELSE
 338                      DO 160, I = MIN( N, J + K ), J + 1, -1
 339                         TEMP = TEMP - DCONJG( A( L + I, J ) )*X( I )
 340   160                CONTINUE
 341                      IF( NOUNIT )
 342      $                  TEMP = TEMP/DCONJG( A( 1, J ) )
 343                   END IF
 344                   X( J ) = TEMP
 345   170          CONTINUE
 346             ELSE
 347                KX = KX + ( N - 1 )*INCX
 348                JX = KX
 349                DO 200, J = N, 1, -1
 350                   TEMP = X( JX )
 351                   IX   = KX
 352                   L    = 1       - J
 353                   IF( NOCONJ )THEN
 354                      DO 180, I = MIN( N, J + K ), J + 1, -1
 355                         TEMP = TEMP - A( L + I, J )*X( IX )
 356                         IX   = IX   - INCX
 357   180                CONTINUE
 358                      IF( NOUNIT )
 359      $                  TEMP = TEMP/A( 1, J )
 360                   ELSE
 361                      DO 190, I = MIN( N, J + K ), J + 1, -1
 362                         TEMP = TEMP - DCONJG( A( L + I, J ) )*X( IX )
 363                         IX   = IX   - INCX
 364   190                CONTINUE
 365                      IF( NOUNIT )
 366      $                  TEMP = TEMP/DCONJG( A( 1, J ) )
 367                   END IF
 368                   X( JX ) = TEMP
 369                   JX      = JX   - INCX
 370                   IF( ( N - J ).GE.K )
 371      $               KX = KX - INCX
 372   200          CONTINUE
 373             END IF
 374          END IF
 375       END IF
 376 *
 377       RETURN
 378 *
 379 *     End of ZTBSV .
 380 *
 381       END