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

CVSweb interface <joel.bertrand@systella.fr>