reference/zsbmvf.f

   1       SUBROUTINE ZSBMVF(UPLO, N, K, ALPHA, A, LDA, X, INCX, BETA, Y,
   2      $                  INCY )
   3 *
   4 *  -- LAPACK auxiliary routine (version 3.1) --
   5 *     Univ. of Tennessee, Univ. of California Berkeley and NAG Ltd..
   6 *     November 2006
   7 *
   8 *     .. Scalar Arguments ..
   9       CHARACTER          UPLO
  10       INTEGER            INCX, INCY, K, LDA, N
  11       COMPLEX*16         ALPHA, BETA
  12 *     ..
  13 *     .. Array Arguments ..
  14       COMPLEX*16         A( LDA, * ), X( * ), Y( * )
  15 *     ..
  16 *
  17 *  Purpose
  18 *  =======
  19 *
  20 *  ZSBMV  performs the matrix-vector  operation
  21 *
  22 *     y := alpha*A*x + beta*y,
  23 *
  24 *  where alpha and beta are scalars, x and y are n element vectors and
  25 *  A is an n by n symmetric band matrix, with k super-diagonals.
  26 *
  27 *  Arguments
  28 *  ==========
  29 *
  30 *  UPLO   - CHARACTER*1
  31 *           On entry, UPLO specifies whether the upper or lower
  32 *           triangular part of the band matrix A is being supplied as
  33 *           follows:
  34 *
  35 *              UPLO = 'U' or 'u'   The upper triangular part of A is
  36 *                                  being supplied.
  37 *
  38 *              UPLO = 'L' or 'l'   The lower triangular part of A is
  39 *                                  being supplied.
  40 *
  41 *           Unchanged on exit.
  42 *
  43 *  N      - INTEGER
  44 *           On entry, N specifies the order of the matrix A.
  45 *           N must be at least zero.
  46 *           Unchanged on exit.
  47 *
  48 *  K      - INTEGER
  49 *           On entry, K specifies the number of super-diagonals of the
  50 *           matrix A. K must satisfy  0 .le. K.
  51 *           Unchanged on exit.
  52 *
  53 *  ALPHA  - COMPLEX*16
  54 *           On entry, ALPHA specifies the scalar alpha.
  55 *           Unchanged on exit.
  56 *
  57 *  A      - COMPLEX*16 array, dimension( LDA, N )
  58 *           Before entry with UPLO = 'U' or 'u', the leading ( k + 1 )
  59 *           by n part of the array A must contain the upper triangular
  60 *           band part of the symmetric matrix, supplied column by
  61 *           column, with the leading diagonal of the matrix in row
  62 *           ( k + 1 ) of the array, the first super-diagonal starting at
  63 *           position 2 in row k, and so on. The top left k by k triangle
  64 *           of the array A is not referenced.
  65 *           The following program segment will transfer the upper
  66 *           triangular part of a symmetric band matrix from conventional
  67 *           full matrix storage to band storage:
  68 *
  69 *                 DO 20, J = 1, N
  70 *                    M = K + 1 - J
  71 *                    DO 10, I = MAX( 1, J - K ), J
  72 *                       A( M + I, J ) = matrix( I, J )
  73 *              10    CONTINUE
  74 *              20 CONTINUE
  75 *
  76 *           Before entry with UPLO = 'L' or 'l', the leading ( k + 1 )
  77 *           by n part of the array A must contain the lower triangular
  78 *           band part of the symmetric matrix, supplied column by
  79 *           column, with the leading diagonal of the matrix in row 1 of
  80 *           the array, the first sub-diagonal starting at position 1 in
  81 *           row 2, and so on. The bottom right k by k triangle of the
  82 *           array A is not referenced.
  83 *           The following program segment will transfer the lower
  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 = 1 - J
  89 *                    DO 10, I = J, MIN( N, J + K )
  90 *                       A( M + I, J ) = matrix( I, J )
  91 *              10    CONTINUE
  92 *              20 CONTINUE
  93 *
  94 *           Unchanged on exit.
  95 *
  96 *  LDA    - INTEGER
  97 *           On entry, LDA specifies the first dimension of A as declared
  98 *           in the calling (sub) program. LDA must be at least
  99 *           ( k + 1 ).
 100 *           Unchanged on exit.
 101 *
 102 *  X      - COMPLEX*16 array, dimension at least
 103 *           ( 1 + ( N - 1 )*abs( INCX ) ).
 104 *           Before entry, the incremented array X must contain the
 105 *           vector x.
 106 *           Unchanged on exit.
 107 *
 108 *  INCX   - INTEGER
 109 *           On entry, INCX specifies the increment for the elements of
 110 *           X. INCX must not be zero.
 111 *           Unchanged on exit.
 112 *
 113 *  BETA   - COMPLEX*16
 114 *           On entry, BETA specifies the scalar beta.
 115 *           Unchanged on exit.
 116 *
 117 *  Y      - COMPLEX*16 array, dimension at least
 118 *           ( 1 + ( N - 1 )*abs( INCY ) ).
 119 *           Before entry, the incremented array Y must contain the
 120 *           vector y. On exit, Y is overwritten by the updated vector y.
 121 *
 122 *  INCY   - INTEGER
 123 *           On entry, INCY specifies the increment for the elements of
 124 *           Y. INCY must not be zero.
 125 *           Unchanged on exit.
 126 *
 127 *  =====================================================================
 128 *
 129 *     .. Parameters ..
 130       COMPLEX*16         ONE
 131       PARAMETER          ( ONE = ( 1.0D+0, 0.0D+0 ) )
 132       COMPLEX*16         ZERO
 133       PARAMETER          ( ZERO = ( 0.0D+0, 0.0D+0 ) )
 134 *     ..
 135 *     .. Local Scalars ..
 136       INTEGER            I, INFO, IX, IY, J, JX, JY, KPLUS1, KX, KY, L
 137       COMPLEX*16         TEMP1, TEMP2
 138 *     ..
 139 *     .. External Functions ..
 140       LOGICAL            LSAME
 141       EXTERNAL           LSAME
 142 *     ..
 143 *     .. External Subroutines ..
 144       EXTERNAL           XERBLA
 145 *     ..
 146 *     .. Intrinsic Functions ..
 147       INTRINSIC          MAX, MIN
 148 *     ..
 149 *     .. Executable Statements ..
 150 *
 151 *     Test the input parameters.
 152 *
 153       INFO = 0
 154       IF( .NOT.LSAME( UPLO, 'U' ) .AND. .NOT.LSAME( UPLO, 'L' ) ) THEN
 155          INFO = 1
 156       ELSE IF( N.LT.0 ) THEN
 157          INFO = 2
 158       ELSE IF( K.LT.0 ) THEN
 159          INFO = 3
 160       ELSE IF( LDA.LT.( K+1 ) ) THEN
 161          INFO = 6
 162       ELSE IF( INCX.EQ.0 ) THEN
 163          INFO = 8
 164       ELSE IF( INCY.EQ.0 ) THEN
 165          INFO = 11
 166       END IF
 167       IF( INFO.NE.0 ) THEN
 168          CALL XERBLA( 'ZSBMV ', INFO )
 169          RETURN
 170       END IF
 171 *
 172 *     Quick return if possible.
 173 *
 174       IF( ( N.EQ.0 ) .OR. ( ( ALPHA.EQ.ZERO ) .AND. ( BETA.EQ.ONE ) ) )
 175      $   RETURN
 176 *
 177 *     Set up the start points in  X  and  Y.
 178 *
 179       IF( INCX.GT.0 ) THEN
 180          KX = 1
 181       ELSE
 182          KX = 1 - ( N-1 )*INCX
 183       END IF
 184       IF( INCY.GT.0 ) THEN
 185          KY = 1
 186       ELSE
 187          KY = 1 - ( N-1 )*INCY
 188       END IF
 189 *
 190 *     Start the operations. In this version the elements of the array A
 191 *     are accessed sequentially with one pass through A.
 192 *
 193 *     First form  y := beta*y.
 194 *
 195       IF( BETA.NE.ONE ) THEN
 196          IF( INCY.EQ.1 ) THEN
 197             IF( BETA.EQ.ZERO ) THEN
 198                DO 10 I = 1, N
 199                   Y( I ) = ZERO
 200    10          CONTINUE
 201             ELSE
 202                DO 20 I = 1, N
 203                   Y( I ) = BETA*Y( I )
 204    20          CONTINUE
 205             END IF
 206          ELSE
 207             IY = KY
 208             IF( BETA.EQ.ZERO ) THEN
 209                DO 30 I = 1, N
 210                   Y( IY ) = ZERO
 211                   IY = IY + INCY
 212    30          CONTINUE
 213             ELSE
 214                DO 40 I = 1, N
 215                   Y( IY ) = BETA*Y( IY )
 216                   IY = IY + INCY
 217    40          CONTINUE
 218             END IF
 219          END IF
 220       END IF
 221       IF( ALPHA.EQ.ZERO )
 222      $   RETURN
 223       IF( LSAME( UPLO, 'U' ) ) THEN
 224 *
 225 *        Form  y  when upper triangle of A is stored.
 226 *
 227          KPLUS1 = K + 1
 228          IF( ( INCX.EQ.1 ) .AND. ( INCY.EQ.1 ) ) THEN
 229             DO 60 J = 1, N
 230                TEMP1 = ALPHA*X( J )
 231                TEMP2 = ZERO
 232                L = KPLUS1 - J
 233                DO 50 I = MAX( 1, J-K ), J - 1
 234                   Y( I ) = Y( I ) + TEMP1*A( L+I, J )
 235                   TEMP2 = TEMP2 + A( L+I, J )*X( I )
 236    50          CONTINUE
 237                Y( J ) = Y( J ) + TEMP1*A( KPLUS1, J ) + ALPHA*TEMP2
 238    60       CONTINUE
 239          ELSE
 240             JX = KX
 241             JY = KY
 242             DO 80 J = 1, N
 243                TEMP1 = ALPHA*X( JX )
 244                TEMP2 = ZERO
 245                IX = KX
 246                IY = KY
 247                L = KPLUS1 - J
 248                DO 70 I = MAX( 1, J-K ), J - 1
 249                   Y( IY ) = Y( IY ) + TEMP1*A( L+I, J )
 250                   TEMP2 = TEMP2 + A( L+I, J )*X( IX )
 251                   IX = IX + INCX
 252                   IY = IY + INCY
 253    70          CONTINUE
 254                Y( JY ) = Y( JY ) + TEMP1*A( KPLUS1, J ) + ALPHA*TEMP2
 255                JX = JX + INCX
 256                JY = JY + INCY
 257                IF( J.GT.K ) THEN
 258                   KX = KX + INCX
 259                   KY = KY + INCY
 260                END IF
 261    80       CONTINUE
 262          END IF
 263       ELSE
 264 *
 265 *        Form  y  when lower triangle of A is stored.
 266 *
 267          IF( ( INCX.EQ.1 ) .AND. ( INCY.EQ.1 ) ) THEN
 268             DO 100 J = 1, N
 269                TEMP1 = ALPHA*X( J )
 270                TEMP2 = ZERO
 271                Y( J ) = Y( J ) + TEMP1*A( 1, J )
 272                L = 1 - J
 273                DO 90 I = J + 1, MIN( N, J+K )
 274                   Y( I ) = Y( I ) + TEMP1*A( L+I, J )
 275                   TEMP2 = TEMP2 + A( L+I, J )*X( I )
 276    90          CONTINUE
 277                Y( J ) = Y( J ) + ALPHA*TEMP2
 278   100       CONTINUE
 279          ELSE
 280             JX = KX
 281             JY = KY
 282             DO 120 J = 1, N
 283                TEMP1 = ALPHA*X( JX )
 284                TEMP2 = ZERO
 285                Y( JY ) = Y( JY ) + TEMP1*A( 1, J )
 286                L = 1 - J
 287                IX = JX
 288                IY = JY
 289                DO 110 I = J + 1, MIN( N, J+K )
 290                   IX = IX + INCX
 291                   IY = IY + INCY
 292                   Y( IY ) = Y( IY ) + TEMP1*A( L+I, J )
 293                   TEMP2 = TEMP2 + A( L+I, J )*X( IX )
 294   110          CONTINUE
 295                Y( JY ) = Y( JY ) + ALPHA*TEMP2
 296                JX = JX + INCX
 297                JY = JY + INCY
 298   120       CONTINUE
 299          END IF
 300       END IF
 301 *
 302       RETURN
 303 *
 304 *     End of ZSBMV
 305 *
 306       END