BLAS/SRC/ssbmv.f

   1 *> \brief \b SSBMV
   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 SSBMV(UPLO,N,K,ALPHA,A,LDA,X,INCX,BETA,Y,INCY)
  12 *
  13 *       .. Scalar Arguments ..
  14 *       REAL ALPHA,BETA
  15 *       INTEGER INCX,INCY,K,LDA,N
  16 *       CHARACTER UPLO
  17 *       ..
  18 *       .. Array Arguments ..
  19 *       REAL A(LDA,*),X(*),Y(*)
  20 *       ..
  21 *
  22 *
  23 *> \par Purpose:
  24 *  =============
  25 *>
  26 *> \verbatim
  27 *>
  28 *> SSBMV  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 REAL
  70 *>           On entry, ALPHA specifies the scalar alpha.
  71 *> \endverbatim
  72 *>
  73 *> \param[in] A
  74 *> \verbatim
  75 *>          A is REAL 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 REAL 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 REAL
 139 *>           On entry, BETA specifies the scalar beta.
 140 *> \endverbatim
 141 *>
 142 *> \param[in,out] Y
 143 *> \verbatim
 144 *>          Y is REAL 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 single_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 SSBMV(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       REAL ALPHA,BETA
 194       INTEGER INCX,INCY,K,LDA,N
 195       CHARACTER UPLO
 196 *     ..
 197 *     .. Array Arguments ..
 198       REAL A(LDA,*),X(*),Y(*)
 199 *     ..
 200 *
 201 *  =====================================================================
 202 *
 203 *     .. Parameters ..
 204       REAL ONE,ZERO
 205       PARAMETER (ONE=1.0E+0,ZERO=0.0E+0)
 206 *     ..
 207 *     .. Local Scalars ..
 208       REAL 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('SSBMV ',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 SSBMV .
 374 *
 375       END