BLAS/SRC/ctbsv.f

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