share/lapack/blas/fortran/zgbmv.f

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