File:  [local] / rpl / lapack / blas / zgbmv.f
Revision 1.1: download - view: text, annotated - select for diffs - revision graph
Tue Jan 26 15:22:45 2010 UTC (14 years, 3 months ago) by bertrand
Branches: MAIN
CVS tags: HEAD
Initial revision

    1:       SUBROUTINE ZGBMV(TRANS,M,N,KL,KU,ALPHA,A,LDA,X,INCX,BETA,Y,INCY)
    2: *     .. Scalar Arguments ..
    3:       DOUBLE COMPLEX ALPHA,BETA
    4:       INTEGER INCX,INCY,KL,KU,LDA,M,N
    5:       CHARACTER TRANS
    6: *     ..
    7: *     .. Array Arguments ..
    8:       DOUBLE COMPLEX 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: *  Arguments
   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: *  Further Details
  122: *  ===============
  123: *
  124: *  Level 2 Blas routine.
  125: *
  126: *  -- Written on 22-October-1986.
  127: *     Jack Dongarra, Argonne National Lab.
  128: *     Jeremy Du Croz, Nag Central Office.
  129: *     Sven Hammarling, Nag Central Office.
  130: *     Richard Hanson, Sandia National Labs.
  131: *
  132: *  =====================================================================
  133: *
  134: *     .. Parameters ..
  135:       DOUBLE COMPLEX ONE
  136:       PARAMETER (ONE= (1.0D+0,0.0D+0))
  137:       DOUBLE COMPLEX ZERO
  138:       PARAMETER (ZERO= (0.0D+0,0.0D+0))
  139: *     ..
  140: *     .. Local Scalars ..
  141:       DOUBLE COMPLEX TEMP
  142:       INTEGER I,INFO,IX,IY,J,JX,JY,K,KUP1,KX,KY,LENX,LENY
  143:       LOGICAL NOCONJ
  144: *     ..
  145: *     .. External Functions ..
  146:       LOGICAL LSAME
  147:       EXTERNAL LSAME
  148: *     ..
  149: *     .. External Subroutines ..
  150:       EXTERNAL XERBLA
  151: *     ..
  152: *     .. Intrinsic Functions ..
  153:       INTRINSIC DCONJG,MAX,MIN
  154: *     ..
  155: *
  156: *     Test the input parameters.
  157: *
  158:       INFO = 0
  159:       IF (.NOT.LSAME(TRANS,'N') .AND. .NOT.LSAME(TRANS,'T') .AND.
  160:      +    .NOT.LSAME(TRANS,'C')) THEN
  161:           INFO = 1
  162:       ELSE IF (M.LT.0) THEN
  163:           INFO = 2
  164:       ELSE IF (N.LT.0) THEN
  165:           INFO = 3
  166:       ELSE IF (KL.LT.0) THEN
  167:           INFO = 4
  168:       ELSE IF (KU.LT.0) THEN
  169:           INFO = 5
  170:       ELSE IF (LDA.LT. (KL+KU+1)) THEN
  171:           INFO = 8
  172:       ELSE IF (INCX.EQ.0) THEN
  173:           INFO = 10
  174:       ELSE IF (INCY.EQ.0) THEN
  175:           INFO = 13
  176:       END IF
  177:       IF (INFO.NE.0) THEN
  178:           CALL XERBLA('ZGBMV ',INFO)
  179:           RETURN
  180:       END IF
  181: *
  182: *     Quick return if possible.
  183: *
  184:       IF ((M.EQ.0) .OR. (N.EQ.0) .OR.
  185:      +    ((ALPHA.EQ.ZERO).AND. (BETA.EQ.ONE))) RETURN
  186: *
  187:       NOCONJ = LSAME(TRANS,'T')
  188: *
  189: *     Set  LENX  and  LENY, the lengths of the vectors x and y, and set
  190: *     up the start points in  X  and  Y.
  191: *
  192:       IF (LSAME(TRANS,'N')) THEN
  193:           LENX = N
  194:           LENY = M
  195:       ELSE
  196:           LENX = M
  197:           LENY = N
  198:       END IF
  199:       IF (INCX.GT.0) THEN
  200:           KX = 1
  201:       ELSE
  202:           KX = 1 - (LENX-1)*INCX
  203:       END IF
  204:       IF (INCY.GT.0) THEN
  205:           KY = 1
  206:       ELSE
  207:           KY = 1 - (LENY-1)*INCY
  208:       END IF
  209: *
  210: *     Start the operations. In this version the elements of A are
  211: *     accessed sequentially with one pass through the band part of A.
  212: *
  213: *     First form  y := beta*y.
  214: *
  215:       IF (BETA.NE.ONE) THEN
  216:           IF (INCY.EQ.1) THEN
  217:               IF (BETA.EQ.ZERO) THEN
  218:                   DO 10 I = 1,LENY
  219:                       Y(I) = ZERO
  220:    10             CONTINUE
  221:               ELSE
  222:                   DO 20 I = 1,LENY
  223:                       Y(I) = BETA*Y(I)
  224:    20             CONTINUE
  225:               END IF
  226:           ELSE
  227:               IY = KY
  228:               IF (BETA.EQ.ZERO) THEN
  229:                   DO 30 I = 1,LENY
  230:                       Y(IY) = ZERO
  231:                       IY = IY + INCY
  232:    30             CONTINUE
  233:               ELSE
  234:                   DO 40 I = 1,LENY
  235:                       Y(IY) = BETA*Y(IY)
  236:                       IY = IY + INCY
  237:    40             CONTINUE
  238:               END IF
  239:           END IF
  240:       END IF
  241:       IF (ALPHA.EQ.ZERO) RETURN
  242:       KUP1 = KU + 1
  243:       IF (LSAME(TRANS,'N')) THEN
  244: *
  245: *        Form  y := alpha*A*x + y.
  246: *
  247:           JX = KX
  248:           IF (INCY.EQ.1) THEN
  249:               DO 60 J = 1,N
  250:                   IF (X(JX).NE.ZERO) THEN
  251:                       TEMP = ALPHA*X(JX)
  252:                       K = KUP1 - J
  253:                       DO 50 I = MAX(1,J-KU),MIN(M,J+KL)
  254:                           Y(I) = Y(I) + TEMP*A(K+I,J)
  255:    50                 CONTINUE
  256:                   END IF
  257:                   JX = JX + INCX
  258:    60         CONTINUE
  259:           ELSE
  260:               DO 80 J = 1,N
  261:                   IF (X(JX).NE.ZERO) THEN
  262:                       TEMP = ALPHA*X(JX)
  263:                       IY = KY
  264:                       K = KUP1 - J
  265:                       DO 70 I = MAX(1,J-KU),MIN(M,J+KL)
  266:                           Y(IY) = Y(IY) + TEMP*A(K+I,J)
  267:                           IY = IY + INCY
  268:    70                 CONTINUE
  269:                   END IF
  270:                   JX = JX + INCX
  271:                   IF (J.GT.KU) KY = KY + INCY
  272:    80         CONTINUE
  273:           END IF
  274:       ELSE
  275: *
  276: *        Form  y := alpha*A'*x + y  or  y := alpha*conjg( A' )*x + y.
  277: *
  278:           JY = KY
  279:           IF (INCX.EQ.1) THEN
  280:               DO 110 J = 1,N
  281:                   TEMP = ZERO
  282:                   K = KUP1 - J
  283:                   IF (NOCONJ) THEN
  284:                       DO 90 I = MAX(1,J-KU),MIN(M,J+KL)
  285:                           TEMP = TEMP + A(K+I,J)*X(I)
  286:    90                 CONTINUE
  287:                   ELSE
  288:                       DO 100 I = MAX(1,J-KU),MIN(M,J+KL)
  289:                           TEMP = TEMP + DCONJG(A(K+I,J))*X(I)
  290:   100                 CONTINUE
  291:                   END IF
  292:                   Y(JY) = Y(JY) + ALPHA*TEMP
  293:                   JY = JY + INCY
  294:   110         CONTINUE
  295:           ELSE
  296:               DO 140 J = 1,N
  297:                   TEMP = ZERO
  298:                   IX = KX
  299:                   K = KUP1 - J
  300:                   IF (NOCONJ) THEN
  301:                       DO 120 I = MAX(1,J-KU),MIN(M,J+KL)
  302:                           TEMP = TEMP + A(K+I,J)*X(IX)
  303:                           IX = IX + INCX
  304:   120                 CONTINUE
  305:                   ELSE
  306:                       DO 130 I = MAX(1,J-KU),MIN(M,J+KL)
  307:                           TEMP = TEMP + DCONJG(A(K+I,J))*X(IX)
  308:                           IX = IX + INCX
  309:   130                 CONTINUE
  310:                   END IF
  311:                   Y(JY) = Y(JY) + ALPHA*TEMP
  312:                   JY = JY + INCY
  313:                   IF (J.GT.KU) KX = KX + INCX
  314:   140         CONTINUE
  315:           END IF
  316:       END IF
  317: *
  318:       RETURN
  319: *
  320: *     End of ZGBMV .
  321: *
  322:       END

CVSweb interface <joel.bertrand@systella.fr>