Return to dla_geamv.f CVS log

Up to [local] / rpl / lapack / lapack

Cohérence.

    1: *> \brief \b DLA_GEAMV computes a matrix-vector product using a general matrix 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_GEAMV + dependencies 
   10: *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/dla_geamv.f"> 
   11: *> [TGZ]</a> 
   12: *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/dla_geamv.f"> 
   13: *> [ZIP]</a> 
   14: *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/dla_geamv.f"> 
   15: *> [TXT]</a>
   16: *> \endhtmlonly 
   17: *
   18: *  Definition:
   19: *  ===========
   20: *
   21: *       SUBROUTINE DLA_GEAMV ( TRANS, M, N, ALPHA, A, LDA, X, INCX, BETA,
   22: *                              Y, INCY )
   23: * 
   24: *       .. Scalar Arguments ..
   25: *       DOUBLE PRECISION   ALPHA, BETA
   26: *       INTEGER            INCX, INCY, LDA, M, N, TRANS
   27: *       ..
   28: *       .. Array Arguments ..
   29: *       DOUBLE PRECISION   A( LDA, * ), X( * ), Y( * )
   30: *       ..
   31: *  
   32: *
   33: *> \par Purpose:
   34: *  =============
   35: *>
   36: *> \verbatim
   37: *>
   38: *> DLA_GEAMV  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] ALPHA
   89: *> \verbatim
   90: *>          ALPHA is DOUBLE PRECISION
   91: *>           On entry, ALPHA specifies the scalar alpha.
   92: *>           Unchanged on exit.
   93: *> \endverbatim
   94: *>
   95: *> \param[in] A
   96: *> \verbatim
   97: *>          A is DOUBLE PRECISION array of DIMENSION ( LDA, n )
   98: *>           Before entry, the leading m by n part of the array A must
   99: *>           contain the matrix of coefficients.
  100: *>           Unchanged on exit.
  101: *> \endverbatim
  102: *>
  103: *> \param[in] LDA
  104: *> \verbatim
  105: *>          LDA is INTEGER
  106: *>           On entry, LDA specifies the first dimension of A as declared
  107: *>           in the calling (sub) program. LDA must be at least
  108: *>           max( 1, m ).
  109: *>           Unchanged on exit.
  110: *> \endverbatim
  111: *>
  112: *> \param[in] X
  113: *> \verbatim
  114: *>          X is DOUBLE PRECISION array, dimension
  115: *>           ( 1 + ( n - 1 )*abs( INCX ) ) when TRANS = 'N' or 'n'
  116: *>           and at least
  117: *>           ( 1 + ( m - 1 )*abs( INCX ) ) otherwise.
  118: *>           Before entry, the incremented array X must contain the
  119: *>           vector x.
  120: *>           Unchanged on exit.
  121: *> \endverbatim
  122: *>
  123: *> \param[in] INCX
  124: *> \verbatim
  125: *>          INCX is INTEGER
  126: *>           On entry, INCX specifies the increment for the elements of
  127: *>           X. INCX must not be zero.
  128: *>           Unchanged on exit.
  129: *> \endverbatim
  130: *>
  131: *> \param[in] BETA
  132: *> \verbatim
  133: *>          BETA is DOUBLE PRECISION
  134: *>           On entry, BETA specifies the scalar beta. When BETA is
  135: *>           supplied as zero then Y need not be set on input.
  136: *>           Unchanged on exit.
  137: *> \endverbatim
  138: *>
  139: *> \param[in,out] Y
  140: *> \verbatim
  141: *>          Y is DOUBLE PRECISION
  142: *>           Array of DIMENSION at least
  143: *>           ( 1 + ( m - 1 )*abs( INCY ) ) when TRANS = 'N' or 'n'
  144: *>           and at least
  145: *>           ( 1 + ( n - 1 )*abs( INCY ) ) otherwise.
  146: *>           Before entry with BETA non-zero, the incremented array Y
  147: *>           must contain the vector y. On exit, Y is overwritten by the
  148: *>           updated vector y.
  149: *> \endverbatim
  150: *>
  151: *> \param[in] INCY
  152: *> \verbatim
  153: *>          INCY is INTEGER
  154: *>           On entry, INCY specifies the increment for the elements of
  155: *>           Y. INCY must not be zero.
  156: *>           Unchanged on exit.
  157: *>
  158: *>  Level 2 Blas routine.
  159: *> \endverbatim
  160: *
  161: *  Authors:
  162: *  ========
  163: *
  164: *> \author Univ. of Tennessee 
  165: *> \author Univ. of California Berkeley 
  166: *> \author Univ. of Colorado Denver 
  167: *> \author NAG Ltd. 
  168: *
  169: *> \date September 2012
  170: *
  171: *> \ingroup doubleGEcomputational
  172: *
  173: *  =====================================================================
  174:       SUBROUTINE DLA_GEAMV ( TRANS, M, N, ALPHA, A, LDA, X, INCX, BETA,
  175:      $                       Y, INCY )
  176: *
  177: *  -- LAPACK computational routine (version 3.4.2) --
  178: *  -- LAPACK is a software package provided by Univ. of Tennessee,    --
  179: *  -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
  180: *     September 2012
  181: *
  182: *     .. Scalar Arguments ..
  183:       DOUBLE PRECISION   ALPHA, BETA
  184:       INTEGER            INCX, INCY, LDA, M, N, TRANS
  185: *     ..
  186: *     .. Array Arguments ..
  187:       DOUBLE PRECISION   A( LDA, * ), X( * ), Y( * )
  188: *     ..
  189: *
  190: *  =====================================================================
  191: *
  192: *     .. Parameters ..
  193:       DOUBLE PRECISION   ONE, ZERO
  194:       PARAMETER          ( ONE = 1.0D+0, ZERO = 0.0D+0 )
  195: *     ..
  196: *     .. Local Scalars ..
  197:       LOGICAL            SYMB_ZERO
  198:       DOUBLE PRECISION   TEMP, SAFE1
  199:       INTEGER            I, INFO, IY, J, JX, KX, KY, LENX, LENY
  200: *     ..
  201: *     .. External Subroutines ..
  202:       EXTERNAL           XERBLA, DLAMCH
  203:       DOUBLE PRECISION   DLAMCH
  204: *     ..
  205: *     .. External Functions ..
  206:       EXTERNAL           ILATRANS
  207:       INTEGER            ILATRANS
  208: *     ..
  209: *     .. Intrinsic Functions ..
  210:       INTRINSIC          MAX, ABS, SIGN
  211: *     ..
  212: *     .. Executable Statements ..
  213: *
  214: *     Test the input parameters.
  215: *
  216:       INFO = 0
  217:       IF     ( .NOT.( ( TRANS.EQ.ILATRANS( 'N' ) )
  218:      $           .OR. ( TRANS.EQ.ILATRANS( 'T' ) )
  219:      $           .OR. ( TRANS.EQ.ILATRANS( 'C' )) ) ) THEN
  220:          INFO = 1
  221:       ELSE IF( M.LT.0 )THEN
  222:          INFO = 2
  223:       ELSE IF( N.LT.0 )THEN
  224:          INFO = 3
  225:       ELSE IF( LDA.LT.MAX( 1, M ) )THEN
  226:          INFO = 6
  227:       ELSE IF( INCX.EQ.0 )THEN
  228:          INFO = 8
  229:       ELSE IF( INCY.EQ.0 )THEN
  230:          INFO = 11
  231:       END IF
  232:       IF( INFO.NE.0 )THEN
  233:          CALL XERBLA( 'DLA_GEAMV ', INFO )
  234:          RETURN
  235:       END IF
  236: *
  237: *     Quick return if possible.
  238: *
  239:       IF( ( M.EQ.0 ).OR.( N.EQ.0 ).OR.
  240:      $    ( ( ALPHA.EQ.ZERO ).AND.( BETA.EQ.ONE ) ) )
  241:      $   RETURN
  242: *
  243: *     Set  LENX  and  LENY, the lengths of the vectors x and y, and set
  244: *     up the start points in  X  and  Y.
  245: *
  246:       IF( TRANS.EQ.ILATRANS( 'N' ) )THEN
  247:          LENX = N
  248:          LENY = M
  249:       ELSE
  250:          LENX = M
  251:          LENY = N
  252:       END IF
  253:       IF( INCX.GT.0 )THEN
  254:          KX = 1
  255:       ELSE
  256:          KX = 1 - ( LENX - 1 )*INCX
  257:       END IF
  258:       IF( INCY.GT.0 )THEN
  259:          KY = 1
  260:       ELSE
  261:          KY = 1 - ( LENY - 1 )*INCY
  262:       END IF
  263: *
  264: *     Set SAFE1 essentially to be the underflow threshold times the
  265: *     number of additions in each row.
  266: *
  267:       SAFE1 = DLAMCH( 'Safe minimum' )
  268:       SAFE1 = (N+1)*SAFE1
  269: *
  270: *     Form  y := alpha*abs(A)*abs(x) + beta*abs(y).
  271: *
  272: *     The O(M*N) SYMB_ZERO tests could be replaced by O(N) queries to
  273: *     the inexact flag.  Still doesn't help change the iteration order
  274: *     to per-column.
  275: *
  276:       IY = KY
  277:       IF ( INCX.EQ.1 ) THEN
  278:          IF( TRANS.EQ.ILATRANS( 'N' ) )THEN
  279:             DO I = 1, LENY
  280:                IF ( BETA .EQ. ZERO ) THEN
  281:                   SYMB_ZERO = .TRUE.
  282:                   Y( IY ) = 0.0D+0
  283:                ELSE IF ( Y( IY ) .EQ. ZERO ) THEN
  284:                   SYMB_ZERO = .TRUE.
  285:                ELSE
  286:                   SYMB_ZERO = .FALSE.
  287:                   Y( IY ) = BETA * ABS( Y( IY ) )
  288:                END IF
  289:                IF ( ALPHA .NE. ZERO ) THEN
  290:                   DO J = 1, LENX
  291:                      TEMP = ABS( A( I, J ) )
  292:                      SYMB_ZERO = SYMB_ZERO .AND.
  293:      $                    ( X( J ) .EQ. ZERO .OR. TEMP .EQ. ZERO )
  294: 
  295:                      Y( IY ) = Y( IY ) + ALPHA*ABS( X( J ) )*TEMP
  296:                   END DO
  297:                END IF
  298: 
  299:                IF ( .NOT.SYMB_ZERO )
  300:      $              Y( IY ) = Y( IY ) + SIGN( SAFE1, Y( IY ) )
  301: 
  302:                IY = IY + INCY
  303:             END DO
  304:          ELSE
  305:             DO I = 1, LENY
  306:                IF ( BETA .EQ. ZERO ) THEN
  307:                   SYMB_ZERO = .TRUE.
  308:                   Y( IY ) = 0.0D+0
  309:                ELSE IF ( Y( IY ) .EQ. ZERO ) THEN
  310:                   SYMB_ZERO = .TRUE.
  311:                ELSE
  312:                   SYMB_ZERO = .FALSE.
  313:                   Y( IY ) = BETA * ABS( Y( IY ) )
  314:                END IF
  315:                IF ( ALPHA .NE. ZERO ) THEN
  316:                   DO J = 1, LENX
  317:                      TEMP = ABS( A( J, I ) )
  318:                      SYMB_ZERO = SYMB_ZERO .AND.
  319:      $                    ( X( J ) .EQ. ZERO .OR. TEMP .EQ. ZERO )
  320: 
  321:                      Y( IY ) = Y( IY ) + ALPHA*ABS( X( J ) )*TEMP
  322:                   END DO
  323:                END IF
  324: 
  325:                IF ( .NOT.SYMB_ZERO )
  326:      $              Y( IY ) = Y( IY ) + SIGN( SAFE1, Y( IY ) )
  327: 
  328:                IY = IY + INCY
  329:             END DO
  330:          END IF
  331:       ELSE
  332:          IF( TRANS.EQ.ILATRANS( 'N' ) )THEN
  333:             DO I = 1, LENY
  334:                IF ( BETA .EQ. ZERO ) THEN
  335:                   SYMB_ZERO = .TRUE.
  336:                   Y( IY ) = 0.0D+0
  337:                ELSE IF ( Y( IY ) .EQ. ZERO ) THEN
  338:                   SYMB_ZERO = .TRUE.
  339:                ELSE
  340:                   SYMB_ZERO = .FALSE.
  341:                   Y( IY ) = BETA * ABS( Y( IY ) )
  342:                END IF
  343:                IF ( ALPHA .NE. ZERO ) THEN
  344:                   JX = KX
  345:                   DO J = 1, LENX
  346:                      TEMP = ABS( A( I, J ) )
  347:                      SYMB_ZERO = SYMB_ZERO .AND.
  348:      $                    ( X( JX ) .EQ. ZERO .OR. TEMP .EQ. ZERO )
  349: 
  350:                      Y( IY ) = Y( IY ) + ALPHA*ABS( X( JX ) )*TEMP
  351:                      JX = JX + INCX
  352:                   END DO
  353:                END IF
  354: 
  355:                IF (.NOT.SYMB_ZERO)
  356:      $              Y( IY ) = Y( IY ) + SIGN( SAFE1, Y( IY ) )
  357: 
  358:                IY = IY + INCY
  359:             END DO
  360:          ELSE
  361:             DO I = 1, LENY
  362:                IF ( BETA .EQ. ZERO ) THEN
  363:                   SYMB_ZERO = .TRUE.
  364:                   Y( IY ) = 0.0D+0
  365:                ELSE IF ( Y( IY ) .EQ. ZERO ) THEN
  366:                   SYMB_ZERO = .TRUE.
  367:                ELSE
  368:                   SYMB_ZERO = .FALSE.
  369:                   Y( IY ) = BETA * ABS( Y( IY ) )
  370:                END IF
  371:                IF ( ALPHA .NE. ZERO ) THEN
  372:                   JX = KX
  373:                   DO J = 1, LENX
  374:                      TEMP = ABS( A( J, I ) )
  375:                      SYMB_ZERO = SYMB_ZERO .AND.
  376:      $                    ( X( JX ) .EQ. ZERO .OR. TEMP .EQ. ZERO )
  377: 
  378:                      Y( IY ) = Y( IY ) + ALPHA*ABS( X( JX ) )*TEMP
  379:                      JX = JX + INCX
  380:                   END DO
  381:                END IF
  382: 
  383:                IF (.NOT.SYMB_ZERO)
  384:      $              Y( IY ) = Y( IY ) + SIGN( SAFE1, Y( IY ) )
  385: 
  386:                IY = IY + INCY
  387:             END DO
  388:          END IF
  389: 
  390:       END IF
  391: *
  392:       RETURN
  393: *
  394: *     End of DLA_GEAMV
  395: *
  396:       END