Annotation of rpl/lapack/lapack/dla_syamv.f, revision 1.17

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

CVSweb interface <joel.bertrand@systella.fr>