File:  [local] / rpl / lapack / blas / dsbmv.f
Revision 1.9: download - view: text, annotated - select for diffs - revision graph
Wed Aug 22 09:36:40 2012 UTC (11 years, 9 months ago) by bertrand
Branches: MAIN
CVS tags: rpl-4_1_9, rpl-4_1_10, HEAD
Cohérence

    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 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 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 of 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 of 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: *> \date November 2011
  166: *
  167: *> \ingroup double_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 DSBMV(UPLO,N,K,ALPHA,A,LDA,X,INCX,BETA,Y,INCY)
  186: *
  187: *  -- Reference BLAS level2 routine (version 3.4.0) --
  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: *     November 2011
  191: *
  192: *     .. Scalar Arguments ..
  193:       DOUBLE PRECISION ALPHA,BETA
  194:       INTEGER INCX,INCY,K,LDA,N
  195:       CHARACTER UPLO
  196: *     ..
  197: *     .. Array Arguments ..
  198:       DOUBLE PRECISION A(LDA,*),X(*),Y(*)
  199: *     ..
  200: *
  201: *  =====================================================================
  202: *
  203: *     .. Parameters ..
  204:       DOUBLE PRECISION ONE,ZERO
  205:       PARAMETER (ONE=1.0D+0,ZERO=0.0D+0)
  206: *     ..
  207: *     .. Local Scalars ..
  208:       DOUBLE PRECISION TEMP1,TEMP2
  209:       INTEGER I,INFO,IX,IY,J,JX,JY,KPLUS1,KX,KY,L
  210: *     ..
  211: *     .. External Functions ..
  212:       LOGICAL LSAME
  213:       EXTERNAL LSAME
  214: *     ..
  215: *     .. External Subroutines ..
  216:       EXTERNAL XERBLA
  217: *     ..
  218: *     .. Intrinsic Functions ..
  219:       INTRINSIC MAX,MIN
  220: *     ..
  221: *
  222: *     Test the input parameters.
  223: *
  224:       INFO = 0
  225:       IF (.NOT.LSAME(UPLO,'U') .AND. .NOT.LSAME(UPLO,'L')) THEN
  226:           INFO = 1
  227:       ELSE IF (N.LT.0) THEN
  228:           INFO = 2
  229:       ELSE IF (K.LT.0) THEN
  230:           INFO = 3
  231:       ELSE IF (LDA.LT. (K+1)) THEN
  232:           INFO = 6
  233:       ELSE IF (INCX.EQ.0) THEN
  234:           INFO = 8
  235:       ELSE IF (INCY.EQ.0) THEN
  236:           INFO = 11
  237:       END IF
  238:       IF (INFO.NE.0) THEN
  239:           CALL XERBLA('DSBMV ',INFO)
  240:           RETURN
  241:       END IF
  242: *
  243: *     Quick return if possible.
  244: *
  245:       IF ((N.EQ.0) .OR. ((ALPHA.EQ.ZERO).AND. (BETA.EQ.ONE))) RETURN
  246: *
  247: *     Set up the start points in  X  and  Y.
  248: *
  249:       IF (INCX.GT.0) THEN
  250:           KX = 1
  251:       ELSE
  252:           KX = 1 - (N-1)*INCX
  253:       END IF
  254:       IF (INCY.GT.0) THEN
  255:           KY = 1
  256:       ELSE
  257:           KY = 1 - (N-1)*INCY
  258:       END IF
  259: *
  260: *     Start the operations. In this version the elements of the array A
  261: *     are accessed sequentially with one pass through A.
  262: *
  263: *     First form  y := beta*y.
  264: *
  265:       IF (BETA.NE.ONE) THEN
  266:           IF (INCY.EQ.1) THEN
  267:               IF (BETA.EQ.ZERO) THEN
  268:                   DO 10 I = 1,N
  269:                       Y(I) = ZERO
  270:    10             CONTINUE
  271:               ELSE
  272:                   DO 20 I = 1,N
  273:                       Y(I) = BETA*Y(I)
  274:    20             CONTINUE
  275:               END IF
  276:           ELSE
  277:               IY = KY
  278:               IF (BETA.EQ.ZERO) THEN
  279:                   DO 30 I = 1,N
  280:                       Y(IY) = ZERO
  281:                       IY = IY + INCY
  282:    30             CONTINUE
  283:               ELSE
  284:                   DO 40 I = 1,N
  285:                       Y(IY) = BETA*Y(IY)
  286:                       IY = IY + INCY
  287:    40             CONTINUE
  288:               END IF
  289:           END IF
  290:       END IF
  291:       IF (ALPHA.EQ.ZERO) RETURN
  292:       IF (LSAME(UPLO,'U')) THEN
  293: *
  294: *        Form  y  when upper triangle of A is stored.
  295: *
  296:           KPLUS1 = K + 1
  297:           IF ((INCX.EQ.1) .AND. (INCY.EQ.1)) THEN
  298:               DO 60 J = 1,N
  299:                   TEMP1 = ALPHA*X(J)
  300:                   TEMP2 = ZERO
  301:                   L = KPLUS1 - J
  302:                   DO 50 I = MAX(1,J-K),J - 1
  303:                       Y(I) = Y(I) + TEMP1*A(L+I,J)
  304:                       TEMP2 = TEMP2 + A(L+I,J)*X(I)
  305:    50             CONTINUE
  306:                   Y(J) = Y(J) + TEMP1*A(KPLUS1,J) + ALPHA*TEMP2
  307:    60         CONTINUE
  308:           ELSE
  309:               JX = KX
  310:               JY = KY
  311:               DO 80 J = 1,N
  312:                   TEMP1 = ALPHA*X(JX)
  313:                   TEMP2 = ZERO
  314:                   IX = KX
  315:                   IY = KY
  316:                   L = KPLUS1 - J
  317:                   DO 70 I = MAX(1,J-K),J - 1
  318:                       Y(IY) = Y(IY) + TEMP1*A(L+I,J)
  319:                       TEMP2 = TEMP2 + A(L+I,J)*X(IX)
  320:                       IX = IX + INCX
  321:                       IY = IY + INCY
  322:    70             CONTINUE
  323:                   Y(JY) = Y(JY) + TEMP1*A(KPLUS1,J) + ALPHA*TEMP2
  324:                   JX = JX + INCX
  325:                   JY = JY + INCY
  326:                   IF (J.GT.K) THEN
  327:                       KX = KX + INCX
  328:                       KY = KY + INCY
  329:                   END IF
  330:    80         CONTINUE
  331:           END IF
  332:       ELSE
  333: *
  334: *        Form  y  when lower triangle of A is stored.
  335: *
  336:           IF ((INCX.EQ.1) .AND. (INCY.EQ.1)) THEN
  337:               DO 100 J = 1,N
  338:                   TEMP1 = ALPHA*X(J)
  339:                   TEMP2 = ZERO
  340:                   Y(J) = Y(J) + TEMP1*A(1,J)
  341:                   L = 1 - J
  342:                   DO 90 I = J + 1,MIN(N,J+K)
  343:                       Y(I) = Y(I) + TEMP1*A(L+I,J)
  344:                       TEMP2 = TEMP2 + A(L+I,J)*X(I)
  345:    90             CONTINUE
  346:                   Y(J) = Y(J) + ALPHA*TEMP2
  347:   100         CONTINUE
  348:           ELSE
  349:               JX = KX
  350:               JY = KY
  351:               DO 120 J = 1,N
  352:                   TEMP1 = ALPHA*X(JX)
  353:                   TEMP2 = ZERO
  354:                   Y(JY) = Y(JY) + TEMP1*A(1,J)
  355:                   L = 1 - J
  356:                   IX = JX
  357:                   IY = JY
  358:                   DO 110 I = J + 1,MIN(N,J+K)
  359:                       IX = IX + INCX
  360:                       IY = IY + INCY
  361:                       Y(IY) = Y(IY) + TEMP1*A(L+I,J)
  362:                       TEMP2 = TEMP2 + A(L+I,J)*X(IX)
  363:   110             CONTINUE
  364:                   Y(JY) = Y(JY) + ALPHA*TEMP2
  365:                   JX = JX + INCX
  366:                   JY = JY + INCY
  367:   120         CONTINUE
  368:           END IF
  369:       END IF
  370: *
  371:       RETURN
  372: *
  373: *     End of DSBMV .
  374: *
  375:       END

CVSweb interface <joel.bertrand@systella.fr>