File:  [local] / rpl / lapack / blas / ztbmv.f
Revision 1.16: download - view: text, annotated - select for diffs - revision graph
Mon Aug 7 08:38:46 2023 UTC (8 months, 3 weeks ago) by bertrand
Branches: MAIN
CVS tags: rpl-4_1_35, rpl-4_1_34, HEAD
Première mise à jour de lapack et blas.

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

CVSweb interface <joel.bertrand@systella.fr>