BLAS/SRC/ctpsv.f

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