LAPACK: BLAS/SRC/dsbmv.f Source File

Go to the documentation of this file.
00001       SUBROUTINE DSBMV(UPLO,N,K,ALPHA,A,LDA,X,INCX,BETA,Y,INCY)
00002 *     .. Scalar Arguments ..
00003       DOUBLE PRECISION ALPHA,BETA
00004       INTEGER INCX,INCY,K,LDA,N
00005       CHARACTER UPLO
00006 *     ..
00007 *     .. Array Arguments ..
00008       DOUBLE PRECISION A(LDA,*),X(*),Y(*)
00009 *     ..
00010 *
00011 *  Purpose
00012 *  =======
00013 *
00014 *  DSBMV  performs the matrix-vector  operation
00015 *
00016 *     y := alpha*A*x + beta*y,
00017 *
00018 *  where alpha and beta are scalars, x and y are n element vectors and
00019 *  A is an n by n symmetric band matrix, with k super-diagonals.
00020 *
00021 *  Arguments
00022 *  ==========
00023 *
00024 *  UPLO   - CHARACTER*1.
00025 *           On entry, UPLO specifies whether the upper or lower
00026 *           triangular part of the band matrix A is being supplied as
00027 *           follows:
00028 *
00029 *              UPLO = 'U' or 'u'   The upper triangular part of A is
00030 *                                  being supplied.
00031 *
00032 *              UPLO = 'L' or 'l'   The lower triangular part of A is
00033 *                                  being supplied.
00034 *
00035 *           Unchanged on exit.
00036 *
00037 *  N      - INTEGER.
00038 *           On entry, N specifies the order of the matrix A.
00039 *           N must be at least zero.
00040 *           Unchanged on exit.
00041 *
00042 *  K      - INTEGER.
00043 *           On entry, K specifies the number of super-diagonals of the
00044 *           matrix A. K must satisfy  0 .le. K.
00045 *           Unchanged on exit.
00046 *
00047 *  ALPHA  - DOUBLE PRECISION.
00048 *           On entry, ALPHA specifies the scalar alpha.
00049 *           Unchanged on exit.
00050 *
00051 *  A      - DOUBLE PRECISION array of DIMENSION ( LDA, n ).
00052 *           Before entry with UPLO = 'U' or 'u', the leading ( k + 1 )
00053 *           by n part of the array A must contain the upper triangular
00054 *           band part of the symmetric matrix, supplied column by
00055 *           column, with the leading diagonal of the matrix in row
00056 *           ( k + 1 ) of the array, the first super-diagonal starting at
00057 *           position 2 in row k, and so on. The top left k by k triangle
00058 *           of the array A is not referenced.
00059 *           The following program segment will transfer the upper
00060 *           triangular part of a symmetric band matrix from conventional
00061 *           full matrix storage to band storage:
00062 *
00063 *                 DO 20, J = 1, N
00064 *                    M = K + 1 - J
00065 *                    DO 10, I = MAX( 1, J - K ), J
00066 *                       A( M + I, J ) = matrix( I, J )
00067 *              10    CONTINUE
00068 *              20 CONTINUE
00069 *
00070 *           Before entry with UPLO = 'L' or 'l', the leading ( k + 1 )
00071 *           by n part of the array A must contain the lower triangular
00072 *           band part of the symmetric matrix, supplied column by
00073 *           column, with the leading diagonal of the matrix in row 1 of
00074 *           the array, the first sub-diagonal starting at position 1 in
00075 *           row 2, and so on. The bottom right k by k triangle of the
00076 *           array A is not referenced.
00077 *           The following program segment will transfer the lower
00078 *           triangular part of a symmetric band matrix from conventional
00079 *           full matrix storage to band storage:
00080 *
00081 *                 DO 20, J = 1, N
00082 *                    M = 1 - J
00083 *                    DO 10, I = J, MIN( N, J + K )
00084 *                       A( M + I, J ) = matrix( I, J )
00085 *              10    CONTINUE
00086 *              20 CONTINUE
00087 *
00088 *           Unchanged on exit.
00089 *
00090 *  LDA    - INTEGER.
00091 *           On entry, LDA specifies the first dimension of A as declared
00092 *           in the calling (sub) program. LDA must be at least
00093 *           ( k + 1 ).
00094 *           Unchanged on exit.
00095 *
00096 *  X      - DOUBLE PRECISION array of DIMENSION at least
00097 *           ( 1 + ( n - 1 )*abs( INCX ) ).
00098 *           Before entry, the incremented array X must contain the
00099 *           vector x.
00100 *           Unchanged on exit.
00101 *
00102 *  INCX   - INTEGER.
00103 *           On entry, INCX specifies the increment for the elements of
00104 *           X. INCX must not be zero.
00105 *           Unchanged on exit.
00106 *
00107 *  BETA   - DOUBLE PRECISION.
00108 *           On entry, BETA specifies the scalar beta.
00109 *           Unchanged on exit.
00110 *
00111 *  Y      - DOUBLE PRECISION array of DIMENSION at least
00112 *           ( 1 + ( n - 1 )*abs( INCY ) ).
00113 *           Before entry, the incremented array Y must contain the
00114 *           vector y. On exit, Y is overwritten by the updated vector y.
00115 *
00116 *  INCY   - INTEGER.
00117 *           On entry, INCY specifies the increment for the elements of
00118 *           Y. INCY must not be zero.
00119 *           Unchanged on exit.
00120 *
00121 *
00122 *  Level 2 Blas routine.
00123 *
00124 *  -- Written on 22-October-1986.
00125 *     Jack Dongarra, Argonne National Lab.
00126 *     Jeremy Du Croz, Nag Central Office.
00127 *     Sven Hammarling, Nag Central Office.
00128 *     Richard Hanson, Sandia National Labs.
00129 *
00130 *  =====================================================================
00131 *
00132 *     .. Parameters ..
00133       DOUBLE PRECISION ONE,ZERO
00134       PARAMETER (ONE=1.0D+0,ZERO=0.0D+0)
00135 *     ..
00136 *     .. Local Scalars ..
00137       DOUBLE PRECISION TEMP1,TEMP2
00138       INTEGER I,INFO,IX,IY,J,JX,JY,KPLUS1,KX,KY,L
00139 *     ..
00140 *     .. External Functions ..
00141       LOGICAL LSAME
00142       EXTERNAL LSAME
00143 *     ..
00144 *     .. External Subroutines ..
00145       EXTERNAL XERBLA
00146 *     ..
00147 *     .. Intrinsic Functions ..
00148       INTRINSIC MAX,MIN
00149 *     ..
00150 *
00151 *     Test the input parameters.
00152 *
00153       INFO = 0
00154       IF (.NOT.LSAME(UPLO,'U') .AND. .NOT.LSAME(UPLO,'L')) THEN
00155           INFO = 1
00156       ELSE IF (N.LT.0) THEN
00157           INFO = 2
00158       ELSE IF (K.LT.0) THEN
00159           INFO = 3
00160       ELSE IF (LDA.LT. (K+1)) THEN
00161           INFO = 6
00162       ELSE IF (INCX.EQ.0) THEN
00163           INFO = 8
00164       ELSE IF (INCY.EQ.0) THEN
00165           INFO = 11
00166       END IF
00167       IF (INFO.NE.0) THEN
00168           CALL XERBLA('DSBMV ',INFO)
00169           RETURN
00170       END IF
00171 *
00172 *     Quick return if possible.
00173 *
00174       IF ((N.EQ.0) .OR. ((ALPHA.EQ.ZERO).AND. (BETA.EQ.ONE))) RETURN
00175 *
00176 *     Set up the start points in  X  and  Y.
00177 *
00178       IF (INCX.GT.0) THEN
00179           KX = 1
00180       ELSE
00181           KX = 1 - (N-1)*INCX
00182       END IF
00183       IF (INCY.GT.0) THEN
00184           KY = 1
00185       ELSE
00186           KY = 1 - (N-1)*INCY
00187       END IF
00188 *
00189 *     Start the operations. In this version the elements of the array A
00190 *     are accessed sequentially with one pass through A.
00191 *
00192 *     First form  y := beta*y.
00193 *
00194       IF (BETA.NE.ONE) THEN
00195           IF (INCY.EQ.1) THEN
00196               IF (BETA.EQ.ZERO) THEN
00197                   DO 10 I = 1,N
00198                       Y(I) = ZERO
00199    10             CONTINUE
00200               ELSE
00201                   DO 20 I = 1,N
00202                       Y(I) = BETA*Y(I)
00203    20             CONTINUE
00204               END IF
00205           ELSE
00206               IY = KY
00207               IF (BETA.EQ.ZERO) THEN
00208                   DO 30 I = 1,N
00209                       Y(IY) = ZERO
00210                       IY = IY + INCY
00211    30             CONTINUE
00212               ELSE
00213                   DO 40 I = 1,N
00214                       Y(IY) = BETA*Y(IY)
00215                       IY = IY + INCY
00216    40             CONTINUE
00217               END IF
00218           END IF
00219       END IF
00220       IF (ALPHA.EQ.ZERO) RETURN
00221       IF (LSAME(UPLO,'U')) THEN
00222 *
00223 *        Form  y  when upper triangle of A is stored.
00224 *
00225           KPLUS1 = K + 1
00226           IF ((INCX.EQ.1) .AND. (INCY.EQ.1)) THEN
00227               DO 60 J = 1,N
00228                   TEMP1 = ALPHA*X(J)
00229                   TEMP2 = ZERO
00230                   L = KPLUS1 - J
00231                   DO 50 I = MAX(1,J-K),J - 1
00232                       Y(I) = Y(I) + TEMP1*A(L+I,J)
00233                       TEMP2 = TEMP2 + A(L+I,J)*X(I)
00234    50             CONTINUE
00235                   Y(J) = Y(J) + TEMP1*A(KPLUS1,J) + ALPHA*TEMP2
00236    60         CONTINUE
00237           ELSE
00238               JX = KX
00239               JY = KY
00240               DO 80 J = 1,N
00241                   TEMP1 = ALPHA*X(JX)
00242                   TEMP2 = ZERO
00243                   IX = KX
00244                   IY = KY
00245                   L = KPLUS1 - J
00246                   DO 70 I = MAX(1,J-K),J - 1
00247                       Y(IY) = Y(IY) + TEMP1*A(L+I,J)
00248                       TEMP2 = TEMP2 + A(L+I,J)*X(IX)
00249                       IX = IX + INCX
00250                       IY = IY + INCY
00251    70             CONTINUE
00252                   Y(JY) = Y(JY) + TEMP1*A(KPLUS1,J) + ALPHA*TEMP2
00253                   JX = JX + INCX
00254                   JY = JY + INCY
00255                   IF (J.GT.K) THEN
00256                       KX = KX + INCX
00257                       KY = KY + INCY
00258                   END IF
00259    80         CONTINUE
00260           END IF
00261       ELSE
00262 *
00263 *        Form  y  when lower triangle of A is stored.
00264 *
00265           IF ((INCX.EQ.1) .AND. (INCY.EQ.1)) THEN
00266               DO 100 J = 1,N
00267                   TEMP1 = ALPHA*X(J)
00268                   TEMP2 = ZERO
00269                   Y(J) = Y(J) + TEMP1*A(1,J)
00270                   L = 1 - J
00271                   DO 90 I = J + 1,MIN(N,J+K)
00272                       Y(I) = Y(I) + TEMP1*A(L+I,J)
00273                       TEMP2 = TEMP2 + A(L+I,J)*X(I)
00274    90             CONTINUE
00275                   Y(J) = Y(J) + ALPHA*TEMP2
00276   100         CONTINUE
00277           ELSE
00278               JX = KX
00279               JY = KY
00280               DO 120 J = 1,N
00281                   TEMP1 = ALPHA*X(JX)
00282                   TEMP2 = ZERO
00283                   Y(JY) = Y(JY) + TEMP1*A(1,J)
00284                   L = 1 - J
00285                   IX = JX
00286                   IY = JY
00287                   DO 110 I = J + 1,MIN(N,J+K)
00288                       IX = IX + INCX
00289                       IY = IY + INCY
00290                       Y(IY) = Y(IY) + TEMP1*A(L+I,J)
00291                       TEMP2 = TEMP2 + A(L+I,J)*X(IX)
00292   110             CONTINUE
00293                   Y(JY) = Y(JY) + ALPHA*TEMP2
00294                   JX = JX + INCX
00295                   JY = JY + INCY
00296   120         CONTINUE
00297           END IF
00298       END IF
00299 *
00300       RETURN
00301 *
00302 *     End of DSBMV .
00303 *
00304       END