Return to dla_gbamv.f CVS log

Up to [local] / rpl / lapack / lapack

Cohérence.

    1: *> \brief \b DLA_GBAMV performs a matrix-vector operation to calculate error bounds.
    2: *
    3: *  =========== DOCUMENTATION ===========
    4: *
    5: * Online html documentation available at 
    6: *            http://www.netlib.org/lapack/explore-html/ 
    7: *
    8: *> \htmlonly
    9: *> Download DLA_GBAMV + dependencies 
   10: *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/dla_gbamv.f"> 
   11: *> [TGZ]</a> 
   12: *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/dla_gbamv.f"> 
   13: *> [ZIP]</a> 
   14: *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/dla_gbamv.f"> 
   15: *> [TXT]</a>
   16: *> \endhtmlonly 
   17: *
   18: *  Definition:
   19: *  ===========
   20: *
   21: *       SUBROUTINE DLA_GBAMV( TRANS, M, N, KL, KU, ALPHA, AB, LDAB, X,
   22: *                             INCX, BETA, Y, INCY )
   23: * 
   24: *       .. Scalar Arguments ..
   25: *       DOUBLE PRECISION   ALPHA, BETA
   26: *       INTEGER            INCX, INCY, LDAB, M, N, KL, KU, TRANS
   27: *       ..
   28: *       .. Array Arguments ..
   29: *       DOUBLE PRECISION   AB( LDAB, * ), X( * ), Y( * )
   30: *       ..
   31: *  
   32: *
   33: *> \par Purpose:
   34: *  =============
   35: *>
   36: *> \verbatim
   37: *>
   38: *> DLA_GBAMV  performs one of the matrix-vector operations
   39: *>
   40: *>         y := alpha*abs(A)*abs(x) + beta*abs(y),
   41: *>    or   y := alpha*abs(A)**T*abs(x) + beta*abs(y),
   42: *>
   43: *> where alpha and beta are scalars, x and y are vectors and A is an
   44: *> m by n matrix.
   45: *>
   46: *> This function is primarily used in calculating error bounds.
   47: *> To protect against underflow during evaluation, components in
   48: *> the resulting vector are perturbed away from zero by (N+1)
   49: *> times the underflow threshold.  To prevent unnecessarily large
   50: *> errors for block-structure embedded in general matrices,
   51: *> "symbolically" zero components are not perturbed.  A zero
   52: *> entry is considered "symbolic" if all multiplications involved
   53: *> in computing that entry have at least one zero multiplicand.
   54: *> \endverbatim
   55: *
   56: *  Arguments:
   57: *  ==========
   58: *
   59: *> \param[in] TRANS
   60: *> \verbatim
   61: *>          TRANS is INTEGER
   62: *>           On entry, TRANS specifies the operation to be performed as
   63: *>           follows:
   64: *>
   65: *>             BLAS_NO_TRANS      y := alpha*abs(A)*abs(x) + beta*abs(y)
   66: *>             BLAS_TRANS         y := alpha*abs(A**T)*abs(x) + beta*abs(y)
   67: *>             BLAS_CONJ_TRANS    y := alpha*abs(A**T)*abs(x) + beta*abs(y)
   68: *>
   69: *>           Unchanged on exit.
   70: *> \endverbatim
   71: *>
   72: *> \param[in] M
   73: *> \verbatim
   74: *>          M is INTEGER
   75: *>           On entry, M specifies the number of rows of the matrix A.
   76: *>           M must be at least zero.
   77: *>           Unchanged on exit.
   78: *> \endverbatim
   79: *>
   80: *> \param[in] N
   81: *> \verbatim
   82: *>          N is INTEGER
   83: *>           On entry, N specifies the number of columns of the matrix A.
   84: *>           N must be at least zero.
   85: *>           Unchanged on exit.
   86: *> \endverbatim
   87: *>
   88: *> \param[in] KL
   89: *> \verbatim
   90: *>          KL is INTEGER
   91: *>           The number of subdiagonals within the band of A.  KL >= 0.
   92: *> \endverbatim
   93: *>
   94: *> \param[in] KU
   95: *> \verbatim
   96: *>          KU is INTEGER
   97: *>           The number of superdiagonals within the band of A.  KU >= 0.
   98: *> \endverbatim
   99: *>
  100: *> \param[in] ALPHA
  101: *> \verbatim
  102: *>          ALPHA is DOUBLE PRECISION
  103: *>           On entry, ALPHA specifies the scalar alpha.
  104: *>           Unchanged on exit.
  105: *> \endverbatim
  106: *>
  107: *> \param[in] AB
  108: *> \verbatim
  109: *>          AB is DOUBLE PRECISION array of DIMENSION ( LDAB, n )
  110: *>           Before entry, the leading m by n part of the array AB must
  111: *>           contain the matrix of coefficients.
  112: *>           Unchanged on exit.
  113: *> \endverbatim
  114: *>
  115: *> \param[in] LDAB
  116: *> \verbatim
  117: *>          LDAB is INTEGER
  118: *>           On entry, LDA specifies the first dimension of AB as declared
  119: *>           in the calling (sub) program. LDAB must be at least
  120: *>           max( 1, m ).
  121: *>           Unchanged on exit.
  122: *> \endverbatim
  123: *>
  124: *> \param[in] X
  125: *> \verbatim
  126: *>          X is DOUBLE PRECISION array, dimension
  127: *>           ( 1 + ( n - 1 )*abs( INCX ) ) when TRANS = 'N' or 'n'
  128: *>           and at least
  129: *>           ( 1 + ( m - 1 )*abs( INCX ) ) otherwise.
  130: *>           Before entry, the incremented array X must contain the
  131: *>           vector x.
  132: *>           Unchanged on exit.
  133: *> \endverbatim
  134: *>
  135: *> \param[in] INCX
  136: *> \verbatim
  137: *>          INCX is INTEGER
  138: *>           On entry, INCX specifies the increment for the elements of
  139: *>           X. INCX must not be zero.
  140: *>           Unchanged on exit.
  141: *> \endverbatim
  142: *>
  143: *> \param[in] BETA
  144: *> \verbatim
  145: *>          BETA is DOUBLE PRECISION
  146: *>           On entry, BETA specifies the scalar beta. When BETA is
  147: *>           supplied as zero then Y need not be set on input.
  148: *>           Unchanged on exit.
  149: *> \endverbatim
  150: *>
  151: *> \param[in,out] Y
  152: *> \verbatim
  153: *>          Y is DOUBLE PRECISION array, dimension
  154: *>           ( 1 + ( m - 1 )*abs( INCY ) ) when TRANS = 'N' or 'n'
  155: *>           and at least
  156: *>           ( 1 + ( n - 1 )*abs( INCY ) ) otherwise.
  157: *>           Before entry with BETA non-zero, the incremented array Y
  158: *>           must contain the vector y. On exit, Y is overwritten by the
  159: *>           updated vector y.
  160: *> \endverbatim
  161: *>
  162: *> \param[in] INCY
  163: *> \verbatim
  164: *>          INCY is INTEGER
  165: *>           On entry, INCY specifies the increment for the elements of
  166: *>           Y. INCY must not be zero.
  167: *>           Unchanged on exit.
  168: *>
  169: *>  Level 2 Blas routine.
  170: *> \endverbatim
  171: *
  172: *  Authors:
  173: *  ========
  174: *
  175: *> \author Univ. of Tennessee 
  176: *> \author Univ. of California Berkeley 
  177: *> \author Univ. of Colorado Denver 
  178: *> \author NAG Ltd. 
  179: *
  180: *> \date September 2012
  181: *
  182: *> \ingroup doubleGBcomputational
  183: *
  184: *  =====================================================================
  185:       SUBROUTINE DLA_GBAMV( TRANS, M, N, KL, KU, ALPHA, AB, LDAB, X,
  186:      $                      INCX, BETA, Y, INCY )
  187: *
  188: *  -- LAPACK computational routine (version 3.4.2) --
  189: *  -- LAPACK is a software package provided by Univ. of Tennessee,    --
  190: *  -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
  191: *     September 2012
  192: *
  193: *     .. Scalar Arguments ..
  194:       DOUBLE PRECISION   ALPHA, BETA
  195:       INTEGER            INCX, INCY, LDAB, M, N, KL, KU, TRANS
  196: *     ..
  197: *     .. Array Arguments ..
  198:       DOUBLE PRECISION   AB( LDAB, * ), X( * ), Y( * )
  199: *     ..
  200: *
  201: *  =====================================================================
  202: *
  203: *     .. Parameters ..
  204:       DOUBLE PRECISION   ONE, ZERO
  205:       PARAMETER          ( ONE = 1.0D+0, ZERO = 0.0D+0 )
  206: *     ..
  207: *     .. Local Scalars ..
  208:       LOGICAL            SYMB_ZERO
  209:       DOUBLE PRECISION   TEMP, SAFE1
  210:       INTEGER            I, INFO, IY, J, JX, KX, KY, LENX, LENY, KD, KE
  211: *     ..
  212: *     .. External Subroutines ..
  213:       EXTERNAL           XERBLA, DLAMCH
  214:       DOUBLE PRECISION   DLAMCH
  215: *     ..
  216: *     .. External Functions ..
  217:       EXTERNAL           ILATRANS
  218:       INTEGER            ILATRANS
  219: *     ..
  220: *     .. Intrinsic Functions ..
  221:       INTRINSIC          MAX, ABS, SIGN
  222: *     ..
  223: *     .. Executable Statements ..
  224: *
  225: *     Test the input parameters.
  226: *
  227:       INFO = 0
  228:       IF     ( .NOT.( ( TRANS.EQ.ILATRANS( 'N' ) )
  229:      $           .OR. ( TRANS.EQ.ILATRANS( 'T' ) )
  230:      $           .OR. ( TRANS.EQ.ILATRANS( 'C' ) ) ) ) THEN
  231:          INFO = 1
  232:       ELSE IF( M.LT.0 )THEN
  233:          INFO = 2
  234:       ELSE IF( N.LT.0 )THEN
  235:          INFO = 3
  236:       ELSE IF( KL.LT.0 .OR. KL.GT.M-1 ) THEN
  237:          INFO = 4
  238:       ELSE IF( KU.LT.0 .OR. KU.GT.N-1 ) THEN
  239:          INFO = 5
  240:       ELSE IF( LDAB.LT.KL+KU+1 )THEN
  241:          INFO = 6
  242:       ELSE IF( INCX.EQ.0 )THEN
  243:          INFO = 8
  244:       ELSE IF( INCY.EQ.0 )THEN
  245:          INFO = 11
  246:       END IF
  247:       IF( INFO.NE.0 )THEN
  248:          CALL XERBLA( 'DLA_GBAMV ', INFO )
  249:          RETURN
  250:       END IF
  251: *
  252: *     Quick return if possible.
  253: *
  254:       IF( ( M.EQ.0 ).OR.( N.EQ.0 ).OR.
  255:      $    ( ( ALPHA.EQ.ZERO ).AND.( BETA.EQ.ONE ) ) )
  256:      $   RETURN
  257: *
  258: *     Set  LENX  and  LENY, the lengths of the vectors x and y, and set
  259: *     up the start points in  X  and  Y.
  260: *
  261:       IF( TRANS.EQ.ILATRANS( 'N' ) )THEN
  262:          LENX = N
  263:          LENY = M
  264:       ELSE
  265:          LENX = M
  266:          LENY = N
  267:       END IF
  268:       IF( INCX.GT.0 )THEN
  269:          KX = 1
  270:       ELSE
  271:          KX = 1 - ( LENX - 1 )*INCX
  272:       END IF
  273:       IF( INCY.GT.0 )THEN
  274:          KY = 1
  275:       ELSE
  276:          KY = 1 - ( LENY - 1 )*INCY
  277:       END IF
  278: *
  279: *     Set SAFE1 essentially to be the underflow threshold times the
  280: *     number of additions in each row.
  281: *
  282:       SAFE1 = DLAMCH( 'Safe minimum' )
  283:       SAFE1 = (N+1)*SAFE1
  284: *
  285: *     Form  y := alpha*abs(A)*abs(x) + beta*abs(y).
  286: *
  287: *     The O(M*N) SYMB_ZERO tests could be replaced by O(N) queries to
  288: *     the inexact flag.  Still doesn't help change the iteration order
  289: *     to per-column.
  290: *
  291:       KD = KU + 1
  292:       KE = KL + 1
  293:       IY = KY
  294:       IF ( INCX.EQ.1 ) THEN
  295:          IF( TRANS.EQ.ILATRANS( 'N' ) )THEN
  296:             DO I = 1, LENY
  297:                IF ( BETA .EQ. ZERO ) THEN
  298:                   SYMB_ZERO = .TRUE.
  299:                   Y( IY ) = 0.0D+0
  300:                ELSE IF ( Y( IY ) .EQ. ZERO ) THEN
  301:                   SYMB_ZERO = .TRUE.
  302:                ELSE
  303:                   SYMB_ZERO = .FALSE.
  304:                   Y( IY ) = BETA * ABS( Y( IY ) )
  305:                END IF
  306:                IF ( ALPHA .NE. ZERO ) THEN
  307:                   DO J = MAX( I-KL, 1 ), MIN( I+KU, LENX )
  308:                      TEMP = ABS( AB( KD+I-J, J ) )
  309:                      SYMB_ZERO = SYMB_ZERO .AND.
  310:      $                    ( X( J ) .EQ. ZERO .OR. TEMP .EQ. ZERO )
  311: 
  312:                      Y( IY ) = Y( IY ) + ALPHA*ABS( X( J ) )*TEMP
  313:                   END DO
  314:                END IF
  315: 
  316:                IF ( .NOT.SYMB_ZERO )
  317:      $              Y( IY ) = Y( IY ) + SIGN( SAFE1, Y( IY ) )
  318:                IY = IY + INCY
  319:             END DO
  320:          ELSE
  321:             DO I = 1, LENY
  322:                IF ( BETA .EQ. ZERO ) THEN
  323:                   SYMB_ZERO = .TRUE.
  324:                   Y( IY ) = 0.0D+0
  325:                ELSE IF ( Y( IY ) .EQ. ZERO ) THEN
  326:                   SYMB_ZERO = .TRUE.
  327:                ELSE
  328:                   SYMB_ZERO = .FALSE.
  329:                   Y( IY ) = BETA * ABS( Y( IY ) )
  330:                END IF
  331:                IF ( ALPHA .NE. ZERO ) THEN
  332:                   DO J = MAX( I-KL, 1 ), MIN( I+KU, LENX )
  333:                      TEMP = ABS( AB( KE-I+J, I ) )
  334:                      SYMB_ZERO = SYMB_ZERO .AND.
  335:      $                    ( X( J ) .EQ. ZERO .OR. TEMP .EQ. ZERO )
  336: 
  337:                      Y( IY ) = Y( IY ) + ALPHA*ABS( X( J ) )*TEMP
  338:                   END DO
  339:                END IF
  340: 
  341:                IF ( .NOT.SYMB_ZERO )
  342:      $              Y( IY ) = Y( IY ) + SIGN( SAFE1, Y( IY ) )
  343:                IY = IY + INCY
  344:             END DO
  345:          END IF
  346:       ELSE
  347:          IF( TRANS.EQ.ILATRANS( 'N' ) )THEN
  348:             DO I = 1, LENY
  349:                IF ( BETA .EQ. ZERO ) THEN
  350:                   SYMB_ZERO = .TRUE.
  351:                   Y( IY ) = 0.0D+0
  352:                ELSE IF ( Y( IY ) .EQ. ZERO ) THEN
  353:                   SYMB_ZERO = .TRUE.
  354:                ELSE
  355:                   SYMB_ZERO = .FALSE.
  356:                   Y( IY ) = BETA * ABS( Y( IY ) )
  357:                END IF
  358:                IF ( ALPHA .NE. ZERO ) THEN
  359:                   JX = KX
  360:                   DO J = MAX( I-KL, 1 ), MIN( I+KU, LENX )
  361:                      TEMP = ABS( AB( KD+I-J, J ) )
  362:                      SYMB_ZERO = SYMB_ZERO .AND.
  363:      $                    ( X( JX ) .EQ. ZERO .OR. TEMP .EQ. ZERO )
  364: 
  365:                      Y( IY ) = Y( IY ) + ALPHA*ABS( X( JX ) )*TEMP
  366:                      JX = JX + INCX
  367:                   END DO
  368:                END IF
  369: 
  370:                IF ( .NOT.SYMB_ZERO )
  371:      $              Y( IY ) = Y( IY ) + SIGN( SAFE1, Y( IY ) )
  372: 
  373:                IY = IY + INCY
  374:             END DO
  375:          ELSE
  376:             DO I = 1, LENY
  377:                IF ( BETA .EQ. ZERO ) THEN
  378:                   SYMB_ZERO = .TRUE.
  379:                   Y( IY ) = 0.0D+0
  380:                ELSE IF ( Y( IY ) .EQ. ZERO ) THEN
  381:                   SYMB_ZERO = .TRUE.
  382:                ELSE
  383:                   SYMB_ZERO = .FALSE.
  384:                   Y( IY ) = BETA * ABS( Y( IY ) )
  385:                END IF
  386:                IF ( ALPHA .NE. ZERO ) THEN
  387:                   JX = KX
  388:                   DO J = MAX( I-KL, 1 ), MIN( I+KU, LENX )
  389:                      TEMP = ABS( AB( KE-I+J, I ) )
  390:                      SYMB_ZERO = SYMB_ZERO .AND.
  391:      $                    ( X( JX ) .EQ. ZERO .OR. TEMP .EQ. ZERO )
  392: 
  393:                      Y( IY ) = Y( IY ) + ALPHA*ABS( X( JX ) )*TEMP
  394:                      JX = JX + INCX
  395:                   END DO
  396:                END IF
  397: 
  398:                IF ( .NOT.SYMB_ZERO )
  399:      $              Y( IY ) = Y( IY ) + SIGN( SAFE1, Y( IY ) )
  400: 
  401:                IY = IY + INCY
  402:             END DO
  403:          END IF
  404: 
  405:       END IF
  406: *
  407:       RETURN
  408: *
  409: *     End of DLA_GBAMV
  410: *
  411:       END