File:  [local] / rpl / lapack / blas / zhbmv.f
Revision 1.8: download - view: text, annotated - select for diffs - revision graph
Mon Nov 21 20:37:09 2011 UTC (12 years, 6 months ago) by bertrand
Branches: MAIN
CVS tags: rpl-4_1_8, rpl-4_1_7, rpl-4_1_6, rpl-4_1_5, rpl-4_1_4, HEAD
Mise à jour de 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 of 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 of 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 of 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: *> \date November 2011
  169: *
  170: *> \ingroup complex16_blas_level2
  171: *
  172: *> \par Further Details:
  173: *  =====================
  174: *>
  175: *> \verbatim
  176: *>
  177: *>  Level 2 Blas routine.
  178: *>  The vector and matrix arguments are not referenced when N = 0, or M = 0
  179: *>
  180: *>  -- Written on 22-October-1986.
  181: *>     Jack Dongarra, Argonne National Lab.
  182: *>     Jeremy Du Croz, Nag Central Office.
  183: *>     Sven Hammarling, Nag Central Office.
  184: *>     Richard Hanson, Sandia National Labs.
  185: *> \endverbatim
  186: *>
  187: *  =====================================================================
  188:       SUBROUTINE ZHBMV(UPLO,N,K,ALPHA,A,LDA,X,INCX,BETA,Y,INCY)
  189: *
  190: *  -- Reference BLAS level2 routine (version 3.4.0) --
  191: *  -- Reference BLAS is a software package provided by Univ. of Tennessee,    --
  192: *  -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
  193: *     November 2011
  194: *
  195: *     .. Scalar Arguments ..
  196:       COMPLEX*16 ALPHA,BETA
  197:       INTEGER INCX,INCY,K,LDA,N
  198:       CHARACTER UPLO
  199: *     ..
  200: *     .. Array Arguments ..
  201:       COMPLEX*16 A(LDA,*),X(*),Y(*)
  202: *     ..
  203: *
  204: *  =====================================================================
  205: *
  206: *     .. Parameters ..
  207:       COMPLEX*16 ONE
  208:       PARAMETER (ONE= (1.0D+0,0.0D+0))
  209:       COMPLEX*16 ZERO
  210:       PARAMETER (ZERO= (0.0D+0,0.0D+0))
  211: *     ..
  212: *     .. Local Scalars ..
  213:       COMPLEX*16 TEMP1,TEMP2
  214:       INTEGER I,INFO,IX,IY,J,JX,JY,KPLUS1,KX,KY,L
  215: *     ..
  216: *     .. External Functions ..
  217:       LOGICAL LSAME
  218:       EXTERNAL LSAME
  219: *     ..
  220: *     .. External Subroutines ..
  221:       EXTERNAL XERBLA
  222: *     ..
  223: *     .. Intrinsic Functions ..
  224:       INTRINSIC DBLE,DCONJG,MAX,MIN
  225: *     ..
  226: *
  227: *     Test the input parameters.
  228: *
  229:       INFO = 0
  230:       IF (.NOT.LSAME(UPLO,'U') .AND. .NOT.LSAME(UPLO,'L')) THEN
  231:           INFO = 1
  232:       ELSE IF (N.LT.0) THEN
  233:           INFO = 2
  234:       ELSE IF (K.LT.0) THEN
  235:           INFO = 3
  236:       ELSE IF (LDA.LT. (K+1)) THEN
  237:           INFO = 6
  238:       ELSE IF (INCX.EQ.0) THEN
  239:           INFO = 8
  240:       ELSE IF (INCY.EQ.0) THEN
  241:           INFO = 11
  242:       END IF
  243:       IF (INFO.NE.0) THEN
  244:           CALL XERBLA('ZHBMV ',INFO)
  245:           RETURN
  246:       END IF
  247: *
  248: *     Quick return if possible.
  249: *
  250:       IF ((N.EQ.0) .OR. ((ALPHA.EQ.ZERO).AND. (BETA.EQ.ONE))) RETURN
  251: *
  252: *     Set up the start points in  X  and  Y.
  253: *
  254:       IF (INCX.GT.0) THEN
  255:           KX = 1
  256:       ELSE
  257:           KX = 1 - (N-1)*INCX
  258:       END IF
  259:       IF (INCY.GT.0) THEN
  260:           KY = 1
  261:       ELSE
  262:           KY = 1 - (N-1)*INCY
  263:       END IF
  264: *
  265: *     Start the operations. In this version the elements of the array A
  266: *     are accessed sequentially with one pass through A.
  267: *
  268: *     First form  y := beta*y.
  269: *
  270:       IF (BETA.NE.ONE) THEN
  271:           IF (INCY.EQ.1) THEN
  272:               IF (BETA.EQ.ZERO) THEN
  273:                   DO 10 I = 1,N
  274:                       Y(I) = ZERO
  275:    10             CONTINUE
  276:               ELSE
  277:                   DO 20 I = 1,N
  278:                       Y(I) = BETA*Y(I)
  279:    20             CONTINUE
  280:               END IF
  281:           ELSE
  282:               IY = KY
  283:               IF (BETA.EQ.ZERO) THEN
  284:                   DO 30 I = 1,N
  285:                       Y(IY) = ZERO
  286:                       IY = IY + INCY
  287:    30             CONTINUE
  288:               ELSE
  289:                   DO 40 I = 1,N
  290:                       Y(IY) = BETA*Y(IY)
  291:                       IY = IY + INCY
  292:    40             CONTINUE
  293:               END IF
  294:           END IF
  295:       END IF
  296:       IF (ALPHA.EQ.ZERO) RETURN
  297:       IF (LSAME(UPLO,'U')) THEN
  298: *
  299: *        Form  y  when upper triangle of A is stored.
  300: *
  301:           KPLUS1 = K + 1
  302:           IF ((INCX.EQ.1) .AND. (INCY.EQ.1)) THEN
  303:               DO 60 J = 1,N
  304:                   TEMP1 = ALPHA*X(J)
  305:                   TEMP2 = ZERO
  306:                   L = KPLUS1 - J
  307:                   DO 50 I = MAX(1,J-K),J - 1
  308:                       Y(I) = Y(I) + TEMP1*A(L+I,J)
  309:                       TEMP2 = TEMP2 + DCONJG(A(L+I,J))*X(I)
  310:    50             CONTINUE
  311:                   Y(J) = Y(J) + TEMP1*DBLE(A(KPLUS1,J)) + ALPHA*TEMP2
  312:    60         CONTINUE
  313:           ELSE
  314:               JX = KX
  315:               JY = KY
  316:               DO 80 J = 1,N
  317:                   TEMP1 = ALPHA*X(JX)
  318:                   TEMP2 = ZERO
  319:                   IX = KX
  320:                   IY = KY
  321:                   L = KPLUS1 - J
  322:                   DO 70 I = MAX(1,J-K),J - 1
  323:                       Y(IY) = Y(IY) + TEMP1*A(L+I,J)
  324:                       TEMP2 = TEMP2 + DCONJG(A(L+I,J))*X(IX)
  325:                       IX = IX + INCX
  326:                       IY = IY + INCY
  327:    70             CONTINUE
  328:                   Y(JY) = Y(JY) + TEMP1*DBLE(A(KPLUS1,J)) + ALPHA*TEMP2
  329:                   JX = JX + INCX
  330:                   JY = JY + INCY
  331:                   IF (J.GT.K) THEN
  332:                       KX = KX + INCX
  333:                       KY = KY + INCY
  334:                   END IF
  335:    80         CONTINUE
  336:           END IF
  337:       ELSE
  338: *
  339: *        Form  y  when lower triangle of A is stored.
  340: *
  341:           IF ((INCX.EQ.1) .AND. (INCY.EQ.1)) THEN
  342:               DO 100 J = 1,N
  343:                   TEMP1 = ALPHA*X(J)
  344:                   TEMP2 = ZERO
  345:                   Y(J) = Y(J) + TEMP1*DBLE(A(1,J))
  346:                   L = 1 - J
  347:                   DO 90 I = J + 1,MIN(N,J+K)
  348:                       Y(I) = Y(I) + TEMP1*A(L+I,J)
  349:                       TEMP2 = TEMP2 + DCONJG(A(L+I,J))*X(I)
  350:    90             CONTINUE
  351:                   Y(J) = Y(J) + ALPHA*TEMP2
  352:   100         CONTINUE
  353:           ELSE
  354:               JX = KX
  355:               JY = KY
  356:               DO 120 J = 1,N
  357:                   TEMP1 = ALPHA*X(JX)
  358:                   TEMP2 = ZERO
  359:                   Y(JY) = Y(JY) + TEMP1*DBLE(A(1,J))
  360:                   L = 1 - J
  361:                   IX = JX
  362:                   IY = JY
  363:                   DO 110 I = J + 1,MIN(N,J+K)
  364:                       IX = IX + INCX
  365:                       IY = IY + INCY
  366:                       Y(IY) = Y(IY) + TEMP1*A(L+I,J)
  367:                       TEMP2 = TEMP2 + DCONJG(A(L+I,J))*X(IX)
  368:   110             CONTINUE
  369:                   Y(JY) = Y(JY) + ALPHA*TEMP2
  370:                   JX = JX + INCX
  371:                   JY = JY + INCY
  372:   120         CONTINUE
  373:           END IF
  374:       END IF
  375: *
  376:       RETURN
  377: *
  378: *     End of ZHBMV .
  379: *
  380:       END

CVSweb interface <joel.bertrand@systella.fr>