File:  [local] / rpl / lapack / lapack / zla_syamv.f
Revision 1.3: download - view: text, annotated - select for diffs - revision graph
Fri Aug 13 21:04:07 2010 UTC (13 years, 9 months ago) by bertrand
Branches: MAIN
CVS tags: rpl-4_0_19, rpl-4_0_18, HEAD
Patches pour OS/2

    1:       SUBROUTINE ZLA_SYAMV( UPLO, N, ALPHA, A, LDA, X, INCX, BETA, Y,
    2:      $                      INCY )
    3: *
    4: *     -- LAPACK routine (version 3.2.2)                                 --
    5: *     -- Contributed by James Demmel, Deaglan Halligan, Yozo Hida and --
    6: *     -- Jason Riedy of Univ. of California Berkeley.                 --
    7: *     -- June 2010                                                    --
    8: *
    9: *     -- LAPACK is a software package provided by Univ. of Tennessee, --
   10: *     -- Univ. of California Berkeley and NAG Ltd.                    --
   11: *
   12:       IMPLICIT NONE
   13: *     ..
   14: *     .. Scalar Arguments ..
   15:       DOUBLE PRECISION   ALPHA, BETA
   16:       INTEGER            INCX, INCY, LDA, N
   17:       INTEGER            UPLO
   18: *     ..
   19: *     .. Array Arguments ..
   20:       COMPLEX*16         A( LDA, * ), X( * )
   21:       DOUBLE PRECISION   Y( * )
   22: *     ..
   23: *
   24: *  Purpose
   25: *  =======
   26: *
   27: *  ZLA_SYAMV  performs the matrix-vector operation
   28: *
   29: *          y := alpha*abs(A)*abs(x) + beta*abs(y),
   30: *
   31: *  where alpha and beta are scalars, x and y are vectors and A is an
   32: *  n by n symmetric matrix.
   33: *
   34: *  This function is primarily used in calculating error bounds.
   35: *  To protect against underflow during evaluation, components in
   36: *  the resulting vector are perturbed away from zero by (N+1)
   37: *  times the underflow threshold.  To prevent unnecessarily large
   38: *  errors for block-structure embedded in general matrices,
   39: *  "symbolically" zero components are not perturbed.  A zero
   40: *  entry is considered "symbolic" if all multiplications involved
   41: *  in computing that entry have at least one zero multiplicand.
   42: *
   43: *  Arguments
   44: *  ==========
   45: *
   46: *  UPLO    (input) INTEGER
   47: *           On entry, UPLO specifies whether the upper or lower
   48: *           triangular part of the array A is to be referenced as
   49: *           follows:
   50: *
   51: *              UPLO = BLAS_UPPER   Only the upper triangular part of A
   52: *                                  is to be referenced.
   53: *
   54: *              UPLO = BLAS_LOWER   Only the lower triangular part of A
   55: *                                  is to be referenced.
   56: *
   57: *           Unchanged on exit.
   58: *
   59: *  N       (input) INTEGER
   60: *           On entry, N specifies the number of columns of the matrix A.
   61: *           N must be at least zero.
   62: *           Unchanged on exit.
   63: *
   64: *  ALPHA  - DOUBLE PRECISION   .
   65: *           On entry, ALPHA specifies the scalar alpha.
   66: *           Unchanged on exit.
   67: *
   68: *  A      - COMPLEX*16         array of DIMENSION ( LDA, n ).
   69: *           Before entry, the leading m by n part of the array A must
   70: *           contain the matrix of coefficients.
   71: *           Unchanged on exit.
   72: *
   73: *  LDA     (input) INTEGER
   74: *           On entry, LDA specifies the first dimension of A as declared
   75: *           in the calling (sub) program. LDA must be at least
   76: *           max( 1, n ).
   77: *           Unchanged on exit.
   78: *
   79: *  X      - COMPLEX*16         array of DIMENSION at least
   80: *           ( 1 + ( n - 1 )*abs( INCX ) )
   81: *           Before entry, the incremented array X must contain the
   82: *           vector x.
   83: *           Unchanged on exit.
   84: *
   85: *  INCX    (input) INTEGER
   86: *           On entry, INCX specifies the increment for the elements of
   87: *           X. INCX must not be zero.
   88: *           Unchanged on exit.
   89: *
   90: *  BETA   - DOUBLE PRECISION   .
   91: *           On entry, BETA specifies the scalar beta. When BETA is
   92: *           supplied as zero then Y need not be set on input.
   93: *           Unchanged on exit.
   94: *
   95: *  Y       (input/output) DOUBLE PRECISION  array, dimension
   96: *           ( 1 + ( n - 1 )*abs( INCY ) )
   97: *           Before entry with BETA non-zero, the incremented array Y
   98: *           must contain the vector y. On exit, Y is overwritten by the
   99: *           updated vector y.
  100: *
  101: *  INCY    (input) INTEGER
  102: *           On entry, INCY specifies the increment for the elements of
  103: *           Y. INCY must not be zero.
  104: *           Unchanged on exit.
  105: *
  106: *  Further Details
  107: *  ===============
  108: *
  109: *  Level 2 Blas routine.
  110: *
  111: *  -- Written on 22-October-1986.
  112: *     Jack Dongarra, Argonne National Lab.
  113: *     Jeremy Du Croz, Nag Central Office.
  114: *     Sven Hammarling, Nag Central Office.
  115: *     Richard Hanson, Sandia National Labs.
  116: *  -- Modified for the absolute-value product, April 2006
  117: *     Jason Riedy, UC Berkeley
  118: *
  119: *  =====================================================================
  120: *
  121: *     .. Parameters ..
  122:       DOUBLE PRECISION   ONE, ZERO
  123:       PARAMETER          ( ONE = 1.0D+0, ZERO = 0.0D+0 )
  124: *     ..
  125: *     .. Local Scalars ..
  126:       LOGICAL            SYMB_ZERO
  127:       DOUBLE PRECISION   TEMP, SAFE1
  128:       INTEGER            I, INFO, IY, J, JX, KX, KY
  129:       COMPLEX*16         ZDUM
  130: *     ..
  131: *     .. External Subroutines ..
  132:       EXTERNAL           XERBLA, DLAMCH
  133:       DOUBLE PRECISION   DLAMCH
  134: *     ..
  135: *     .. External Functions ..
  136:       EXTERNAL           ILAUPLO
  137:       INTEGER            ILAUPLO
  138: *     ..
  139: *     .. Intrinsic Functions ..
  140:       INTRINSIC          MAX, ABS, SIGN, REAL, DIMAG
  141: *     ..
  142: *     .. Statement Functions ..
  143:       DOUBLE PRECISION   CABS1
  144: *     ..
  145: *     .. Statement Function Definitions ..
  146:       CABS1( ZDUM ) = ABS( DBLE ( ZDUM ) ) + ABS( DIMAG ( ZDUM ) )
  147: *     ..
  148: *     .. Executable Statements ..
  149: *
  150: *     Test the input parameters.
  151: *
  152:       INFO = 0
  153:       IF     ( UPLO.NE.ILAUPLO( 'U' ) .AND.
  154:      $         UPLO.NE.ILAUPLO( 'L' ) )THEN
  155:          INFO = 1
  156:       ELSE IF( N.LT.0 )THEN
  157:          INFO = 2
  158:       ELSE IF( LDA.LT.MAX( 1, N ) )THEN
  159:          INFO = 5
  160:       ELSE IF( INCX.EQ.0 )THEN
  161:          INFO = 7
  162:       ELSE IF( INCY.EQ.0 )THEN
  163:          INFO = 10
  164:       END IF
  165:       IF( INFO.NE.0 )THEN
  166:          CALL XERBLA( 'DSYMV ', INFO )
  167:          RETURN
  168:       END IF
  169: *
  170: *     Quick return if possible.
  171: *
  172:       IF( ( N.EQ.0 ).OR.( ( ALPHA.EQ.ZERO ).AND.( BETA.EQ.ONE ) ) )
  173:      $   RETURN
  174: *
  175: *     Set up the start points in  X  and  Y.
  176: *
  177:       IF( INCX.GT.0 )THEN
  178:          KX = 1
  179:       ELSE
  180:          KX = 1 - ( N - 1 )*INCX
  181:       END IF
  182:       IF( INCY.GT.0 )THEN
  183:          KY = 1
  184:       ELSE
  185:          KY = 1 - ( N - 1 )*INCY
  186:       END IF
  187: *
  188: *     Set SAFE1 essentially to be the underflow threshold times the
  189: *     number of additions in each row.
  190: *
  191:       SAFE1 = DLAMCH( 'Safe minimum' )
  192:       SAFE1 = (N+1)*SAFE1
  193: *
  194: *     Form  y := alpha*abs(A)*abs(x) + beta*abs(y).
  195: *
  196: *     The O(N^2) SYMB_ZERO tests could be replaced by O(N) queries to
  197: *     the inexact flag.  Still doesn't help change the iteration order
  198: *     to per-column.
  199: *
  200:       IY = KY
  201:       IF ( INCX.EQ.1 ) THEN
  202:          IF ( UPLO .EQ. ILAUPLO( 'U' ) ) THEN
  203:             DO I = 1, N
  204:                IF ( BETA .EQ. ZERO ) THEN
  205:                   SYMB_ZERO = .TRUE.
  206:                   Y( IY ) = 0.0D+0
  207:                ELSE IF ( Y( IY ) .EQ. ZERO ) THEN
  208:                   SYMB_ZERO = .TRUE.
  209:                ELSE
  210:                   SYMB_ZERO = .FALSE.
  211:                   Y( IY ) = BETA * ABS( Y( IY ) )
  212:                END IF
  213:                IF ( ALPHA .NE. ZERO ) THEN
  214:                   DO J = 1, I
  215:                      TEMP = CABS1( A( J, I ) )
  216:                      SYMB_ZERO = SYMB_ZERO .AND.
  217:      $                    ( X( J ) .EQ. ZERO .OR. TEMP .EQ. ZERO )
  218: 
  219:                      Y( IY ) = Y( IY ) + ALPHA*CABS1( X( J ) )*TEMP
  220:                   END DO
  221:                   DO J = I+1, N
  222:                      TEMP = CABS1( A( I, J ) )
  223:                      SYMB_ZERO = SYMB_ZERO .AND.
  224:      $                    ( X( J ) .EQ. ZERO .OR. TEMP .EQ. ZERO )
  225: 
  226:                      Y( IY ) = Y( IY ) + ALPHA*CABS1( X( J ) )*TEMP
  227:                   END DO
  228:                END IF
  229: 
  230:                IF ( .NOT.SYMB_ZERO )
  231:      $              Y( IY ) = Y( IY ) + SIGN( SAFE1, Y( IY ) )
  232: 
  233:                IY = IY + INCY
  234:             END DO
  235:          ELSE
  236:             DO I = 1, N
  237:                IF ( BETA .EQ. ZERO ) THEN
  238:                   SYMB_ZERO = .TRUE.
  239:                   Y( IY ) = 0.0D+0
  240:                ELSE IF ( Y( IY ) .EQ. ZERO ) THEN
  241:                   SYMB_ZERO = .TRUE.
  242:                ELSE
  243:                   SYMB_ZERO = .FALSE.
  244:                   Y( IY ) = BETA * ABS( Y( IY ) )
  245:                END IF
  246:                IF ( ALPHA .NE. ZERO ) THEN
  247:                   DO J = 1, I
  248:                      TEMP = CABS1( A( I, J ) )
  249:                      SYMB_ZERO = SYMB_ZERO .AND.
  250:      $                    ( X( J ) .EQ. ZERO .OR. TEMP .EQ. ZERO )
  251: 
  252:                      Y( IY ) = Y( IY ) + ALPHA*CABS1( X( J ) )*TEMP
  253:                   END DO
  254:                   DO J = I+1, N
  255:                      TEMP = CABS1( A( J, I ) )
  256:                      SYMB_ZERO = SYMB_ZERO .AND.
  257:      $                    ( X( J ) .EQ. ZERO .OR. TEMP .EQ. ZERO )
  258: 
  259:                      Y( IY ) = Y( IY ) + ALPHA*CABS1( X( J ) )*TEMP
  260:                   END DO
  261:                END IF
  262: 
  263:                IF ( .NOT.SYMB_ZERO )
  264:      $              Y( IY ) = Y( IY ) + SIGN( SAFE1, Y( IY ) )
  265: 
  266:                IY = IY + INCY
  267:             END DO
  268:          END IF
  269:       ELSE
  270:          IF ( UPLO .EQ. ILAUPLO( 'U' ) ) THEN
  271:             DO I = 1, N
  272:                IF ( BETA .EQ. ZERO ) THEN
  273:                   SYMB_ZERO = .TRUE.
  274:                   Y( IY ) = 0.0D+0
  275:                ELSE IF ( Y( IY ) .EQ. ZERO ) THEN
  276:                   SYMB_ZERO = .TRUE.
  277:                ELSE
  278:                   SYMB_ZERO = .FALSE.
  279:                   Y( IY ) = BETA * ABS( Y( IY ) )
  280:                END IF
  281:                JX = KX
  282:                IF ( ALPHA .NE. ZERO ) THEN
  283:                   DO J = 1, I
  284:                      TEMP = CABS1( A( J, I ) )
  285:                      SYMB_ZERO = SYMB_ZERO .AND.
  286:      $                    ( X( J ) .EQ. ZERO .OR. TEMP .EQ. ZERO )
  287: 
  288:                      Y( IY ) = Y( IY ) + ALPHA*CABS1( X( JX ) )*TEMP
  289:                      JX = JX + INCX
  290:                   END DO
  291:                   DO J = I+1, N
  292:                      TEMP = CABS1( A( I, J ) )
  293:                      SYMB_ZERO = SYMB_ZERO .AND.
  294:      $                    ( X( J ) .EQ. ZERO .OR. TEMP .EQ. ZERO )
  295: 
  296:                      Y( IY ) = Y( IY ) + ALPHA*CABS1( X( JX ) )*TEMP
  297:                      JX = JX + INCX
  298:                   END DO
  299:                END IF
  300: 
  301:                IF ( .NOT.SYMB_ZERO )
  302:      $              Y( IY ) = Y( IY ) + SIGN( SAFE1, Y( IY ) )
  303: 
  304:                IY = IY + INCY
  305:             END DO
  306:          ELSE
  307:             DO I = 1, N
  308:                IF ( BETA .EQ. ZERO ) THEN
  309:                   SYMB_ZERO = .TRUE.
  310:                   Y( IY ) = 0.0D+0
  311:                ELSE IF ( Y( IY ) .EQ. ZERO ) THEN
  312:                   SYMB_ZERO = .TRUE.
  313:                ELSE
  314:                   SYMB_ZERO = .FALSE.
  315:                   Y( IY ) = BETA * ABS( Y( IY ) )
  316:                END IF
  317:                JX = KX
  318:                IF ( ALPHA .NE. ZERO ) THEN
  319:                   DO J = 1, I
  320:                      TEMP = CABS1( A( I, J ) )
  321:                      SYMB_ZERO = SYMB_ZERO .AND.
  322:      $                    ( X( J ) .EQ. ZERO .OR. TEMP .EQ. ZERO )
  323: 
  324:                      Y( IY ) = Y( IY ) + ALPHA*CABS1( X( JX ) )*TEMP
  325:                      JX = JX + INCX
  326:                   END DO
  327:                   DO J = I+1, N
  328:                      TEMP = CABS1( A( J, I ) )
  329:                      SYMB_ZERO = SYMB_ZERO .AND.
  330:      $                    ( X( J ) .EQ. ZERO .OR. TEMP .EQ. ZERO )
  331: 
  332:                      Y( IY ) = Y( IY ) + ALPHA*CABS1( X( JX ) )*TEMP
  333:                      JX = JX + INCX
  334:                   END DO
  335:                END IF
  336: 
  337:                IF ( .NOT.SYMB_ZERO )
  338:      $              Y( IY ) = Y( IY ) + SIGN( SAFE1, Y( IY ) )
  339: 
  340:                IY = IY + INCY
  341:             END DO
  342:          END IF
  343: 
  344:       END IF
  345: *
  346:       RETURN
  347: *
  348: *     End of ZLA_SYAMV
  349: *
  350:       END

CVSweb interface <joel.bertrand@systella.fr>