File:  [local] / rpl / lapack / blas / zgbmv.f
Revision 1.7: download - view: text, annotated - select for diffs - revision graph
Fri Jul 22 07:38:02 2011 UTC (12 years, 10 months ago) by bertrand
Branches: MAIN
CVS tags: rpl-4_1_3, rpl-4_1_2, rpl-4_1_1, HEAD
En route vers la 4.4.1.

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

CVSweb interface <joel.bertrand@systella.fr>