File:  [local] / rpl / lapack / blas / ztpmv.f
Revision 1.7: download - view: text, annotated - select for diffs - revision graph
Fri Jul 22 07:38:03 2011 UTC (12 years, 9 months ago) by bertrand
Branches: MAIN
CVS tags: rpl-4_1_3, rpl-4_1_2, rpl-4_1_1, HEAD
En route vers la 4.4.1.

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

CVSweb interface <joel.bertrand@systella.fr>