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

CVSweb interface <joel.bertrand@systella.fr>