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

CVSweb interface <joel.bertrand@systella.fr>