SRC/zla_heamv.f

   1 *> \brief \b ZLA_HEAMV computes a matrix-vector product using a Hermitian indefinite 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 ZLA_HEAMV + dependencies
  10 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/zla_heamv.f">
  11 *> [TGZ]</a>
  12 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/zla_heamv.f">
  13 *> [ZIP]</a>
  14 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/zla_heamv.f">
  15 *> [TXT]</a>
  16 *> \endhtmlonly
  17 *
  18 *  Definition:
  19 *  ===========
  20 *
  21 *       SUBROUTINE ZLA_HEAMV( UPLO, N, ALPHA, A, LDA, X, INCX, BETA, Y,
  22 *                             INCY )
  23 *
  24 *       .. Scalar Arguments ..
  25 *       DOUBLE PRECISION   ALPHA, BETA
  26 *       INTEGER            INCX, INCY, LDA, N, UPLO
  27 *       ..
  28 *       .. Array Arguments ..
  29 *       COMPLEX*16         A( LDA, * ), X( * )
  30 *       DOUBLE PRECISION   Y( * )
  31 *       ..
  32 *
  33 *
  34 *> \par Purpose:
  35 *  =============
  36 *>
  37 *> \verbatim
  38 *>
  39 *> ZLA_SYAMV  performs the matrix-vector operation
  40 *>
  41 *>         y := alpha*abs(A)*abs(x) + beta*abs(y),
  42 *>
  43 *> where alpha and beta are scalars, x and y are vectors and A is an
  44 *> n by n symmetric 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] UPLO
  60 *> \verbatim
  61 *>          UPLO is INTEGER
  62 *>           On entry, UPLO specifies whether the upper or lower
  63 *>           triangular part of the array A is to be referenced as
  64 *>           follows:
  65 *>
  66 *>              UPLO = BLAS_UPPER   Only the upper triangular part of A
  67 *>                                  is to be referenced.
  68 *>
  69 *>              UPLO = BLAS_LOWER   Only the lower triangular part of A
  70 *>                                  is to be referenced.
  71 *>
  72 *>           Unchanged on exit.
  73 *> \endverbatim
  74 *>
  75 *> \param[in] N
  76 *> \verbatim
  77 *>          N is INTEGER
  78 *>           On entry, N specifies the number of columns of the matrix A.
  79 *>           N must be at least zero.
  80 *>           Unchanged on exit.
  81 *> \endverbatim
  82 *>
  83 *> \param[in] ALPHA
  84 *> \verbatim
  85 *>          ALPHA is DOUBLE PRECISION .
  86 *>           On entry, ALPHA specifies the scalar alpha.
  87 *>           Unchanged on exit.
  88 *> \endverbatim
  89 *>
  90 *> \param[in] A
  91 *> \verbatim
  92 *>          A is COMPLEX*16 array, DIMENSION ( LDA, n ).
  93 *>           Before entry, the leading m by n part of the array A must
  94 *>           contain the matrix of coefficients.
  95 *>           Unchanged on exit.
  96 *> \endverbatim
  97 *>
  98 *> \param[in] LDA
  99 *> \verbatim
 100 *>          LDA is INTEGER
 101 *>           On entry, LDA specifies the first dimension of A as declared
 102 *>           in the calling (sub) program. LDA must be at least
 103 *>           max( 1, n ).
 104 *>           Unchanged on exit.
 105 *> \endverbatim
 106 *>
 107 *> \param[in] X
 108 *> \verbatim
 109 *>          X is COMPLEX*16 array, DIMENSION at least
 110 *>           ( 1 + ( n - 1 )*abs( INCX ) )
 111 *>           Before entry, the incremented array X must contain the
 112 *>           vector x.
 113 *>           Unchanged on exit.
 114 *> \endverbatim
 115 *>
 116 *> \param[in] INCX
 117 *> \verbatim
 118 *>          INCX is INTEGER
 119 *>           On entry, INCX specifies the increment for the elements of
 120 *>           X. INCX must not be zero.
 121 *>           Unchanged on exit.
 122 *> \endverbatim
 123 *>
 124 *> \param[in] BETA
 125 *> \verbatim
 126 *>          BETA is DOUBLE PRECISION .
 127 *>           On entry, BETA specifies the scalar beta. When BETA is
 128 *>           supplied as zero then Y need not be set on input.
 129 *>           Unchanged on exit.
 130 *> \endverbatim
 131 *>
 132 *> \param[in,out] Y
 133 *> \verbatim
 134 *>          Y is DOUBLE PRECISION array, dimension
 135 *>           ( 1 + ( n - 1 )*abs( INCY ) )
 136 *>           Before entry with BETA non-zero, the incremented array Y
 137 *>           must contain the vector y. On exit, Y is overwritten by the
 138 *>           updated vector y.
 139 *> \endverbatim
 140 *>
 141 *> \param[in] INCY
 142 *> \verbatim
 143 *>          INCY is INTEGER
 144 *>           On entry, INCY specifies the increment for the elements of
 145 *>           Y. INCY must not be zero.
 146 *>           Unchanged on exit.
 147 *> \endverbatim
 148 *
 149 *  Authors:
 150 *  ========
 151 *
 152 *> \author Univ. of Tennessee
 153 *> \author Univ. of California Berkeley
 154 *> \author Univ. of Colorado Denver
 155 *> \author NAG Ltd.
 156 *
 157 *> \date September 2012
 158 *
 159 *> \ingroup complex16HEcomputational
 160 *
 161 *> \par Further Details:
 162 *  =====================
 163 *>
 164 *> \verbatim
 165 *>
 166 *>  Level 2 Blas routine.
 167 *>
 168 *>  -- Written on 22-October-1986.
 169 *>     Jack Dongarra, Argonne National Lab.
 170 *>     Jeremy Du Croz, Nag Central Office.
 171 *>     Sven Hammarling, Nag Central Office.
 172 *>     Richard Hanson, Sandia National Labs.
 173 *>  -- Modified for the absolute-value product, April 2006
 174 *>     Jason Riedy, UC Berkeley
 175 *> \endverbatim
 176 *>
 177 *  =====================================================================
 178       SUBROUTINE ZLA_HEAMV( UPLO, N, ALPHA, A, LDA, X, INCX, BETA, Y,
 179      $                      INCY )
 180 *
 181 *  -- LAPACK computational routine (version 3.4.2) --
 182 *  -- LAPACK is a software package provided by Univ. of Tennessee,    --
 183 *  -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
 184 *     September 2012
 185 *
 186 *     .. Scalar Arguments ..
 187       DOUBLE PRECISION   ALPHA, BETA
 188       INTEGER            INCX, INCY, LDA, N, UPLO
 189 *     ..
 190 *     .. Array Arguments ..
 191       COMPLEX*16         A( LDA, * ), X( * )
 192       DOUBLE PRECISION   Y( * )
 193 *     ..
 194 *
 195 *  =====================================================================
 196 *
 197 *     .. Parameters ..
 198       DOUBLE PRECISION   ONE, ZERO
 199       PARAMETER          ( ONE = 1.0D+0, ZERO = 0.0D+0 )
 200 *     ..
 201 *     .. Local Scalars ..
 202       LOGICAL            SYMB_ZERO
 203       DOUBLE PRECISION   TEMP, SAFE1
 204       INTEGER            I, INFO, IY, J, JX, KX, KY
 205       COMPLEX*16         ZDUM
 206 *     ..
 207 *     .. External Subroutines ..
 208       EXTERNAL           XERBLA, DLAMCH
 209       DOUBLE PRECISION   DLAMCH
 210 *     ..
 211 *     .. External Functions ..
 212       EXTERNAL           ILAUPLO
 213       INTEGER            ILAUPLO
 214 *     ..
 215 *     .. Intrinsic Functions ..
 216       INTRINSIC          MAX, ABS, SIGN, REAL, DIMAG
 217 *     ..
 218 *     .. Statement Functions ..
 219       DOUBLE PRECISION   CABS1
 220 *     ..
 221 *     .. Statement Function Definitions ..
 222       CABS1( ZDUM ) = ABS( DBLE ( ZDUM ) ) + ABS( DIMAG ( ZDUM ) )
 223 *     ..
 224 *     .. Executable Statements ..
 225 *
 226 *     Test the input parameters.
 227 *
 228       INFO = 0
 229       IF     ( UPLO.NE.ILAUPLO( 'U' ) .AND.
 230      $         UPLO.NE.ILAUPLO( 'L' ) )THEN
 231          INFO = 1
 232       ELSE IF( N.LT.0 )THEN
 233          INFO = 2
 234       ELSE IF( LDA.LT.MAX( 1, N ) )THEN
 235          INFO = 5
 236       ELSE IF( INCX.EQ.0 )THEN
 237          INFO = 7
 238       ELSE IF( INCY.EQ.0 )THEN
 239          INFO = 10
 240       END IF
 241       IF( INFO.NE.0 )THEN
 242          CALL XERBLA( 'ZHEMV ', INFO )
 243          RETURN
 244       END IF
 245 *
 246 *     Quick return if possible.
 247 *
 248       IF( ( N.EQ.0 ).OR.( ( ALPHA.EQ.ZERO ).AND.( BETA.EQ.ONE ) ) )
 249      $   RETURN
 250 *
 251 *     Set up the start points in  X  and  Y.
 252 *
 253       IF( INCX.GT.0 )THEN
 254          KX = 1
 255       ELSE
 256          KX = 1 - ( N - 1 )*INCX
 257       END IF
 258       IF( INCY.GT.0 )THEN
 259          KY = 1
 260       ELSE
 261          KY = 1 - ( N - 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(N^2) 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 ( UPLO .EQ. ILAUPLO( 'U' ) ) THEN
 279             DO I = 1, N
 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, I
 291                      TEMP = CABS1( A( J, I ) )
 292                      SYMB_ZERO = SYMB_ZERO .AND.
 293      $                    ( X( J ) .EQ. ZERO .OR. TEMP .EQ. ZERO )
 294
 295                      Y( IY ) = Y( IY ) + ALPHA*CABS1( X( J ) )*TEMP
 296                   END DO
 297                   DO J = I+1, N
 298                      TEMP = CABS1( A( I, J ) )
 299                      SYMB_ZERO = SYMB_ZERO .AND.
 300      $                    ( X( J ) .EQ. ZERO .OR. TEMP .EQ. ZERO )
 301
 302                      Y( IY ) = Y( IY ) + ALPHA*CABS1( X( J ) )*TEMP
 303                   END DO
 304                END IF
 305
 306                IF (.NOT.SYMB_ZERO)
 307      $              Y( IY ) = Y( IY ) + SIGN( SAFE1, Y( IY ) )
 308
 309                IY = IY + INCY
 310             END DO
 311          ELSE
 312             DO I = 1, N
 313                IF ( BETA .EQ. ZERO ) THEN
 314                   SYMB_ZERO = .TRUE.
 315                   Y( IY ) = 0.0D+0
 316                ELSE IF ( Y( IY ) .EQ. ZERO ) THEN
 317                   SYMB_ZERO = .TRUE.
 318                ELSE
 319                   SYMB_ZERO = .FALSE.
 320                   Y( IY ) = BETA * ABS( Y( IY ) )
 321                END IF
 322                IF ( ALPHA .NE. ZERO ) THEN
 323                   DO J = 1, I
 324                      TEMP = CABS1( A( I, J ) )
 325                      SYMB_ZERO = SYMB_ZERO .AND.
 326      $                    ( X( J ) .EQ. ZERO .OR. TEMP .EQ. ZERO )
 327
 328                      Y( IY ) = Y( IY ) + ALPHA*CABS1( X( J ) )*TEMP
 329                   END DO
 330                   DO J = I+1, N
 331                      TEMP = CABS1( A( J, I ) )
 332                      SYMB_ZERO = SYMB_ZERO .AND.
 333      $                    ( X( J ) .EQ. ZERO .OR. TEMP .EQ. ZERO )
 334
 335                      Y( IY ) = Y( IY ) + ALPHA*CABS1( X( J ) )*TEMP
 336                   END DO
 337                END IF
 338
 339                IF (.NOT.SYMB_ZERO)
 340      $              Y( IY ) = Y( IY ) + SIGN( SAFE1, Y( IY ) )
 341
 342                IY = IY + INCY
 343             END DO
 344          END IF
 345       ELSE
 346          IF ( UPLO .EQ. ILAUPLO( 'U' ) ) THEN
 347             DO I = 1, N
 348                IF ( BETA .EQ. ZERO ) THEN
 349                   SYMB_ZERO = .TRUE.
 350                   Y( IY ) = 0.0D+0
 351                ELSE IF ( Y( IY ) .EQ. ZERO ) THEN
 352                   SYMB_ZERO = .TRUE.
 353                ELSE
 354                   SYMB_ZERO = .FALSE.
 355                   Y( IY ) = BETA * ABS( Y( IY ) )
 356                END IF
 357                JX = KX
 358                IF ( ALPHA .NE. ZERO ) THEN
 359                   DO J = 1, I
 360                      TEMP = CABS1( A( J, I ) )
 361                      SYMB_ZERO = SYMB_ZERO .AND.
 362      $                    ( X( J ) .EQ. ZERO .OR. TEMP .EQ. ZERO )
 363
 364                      Y( IY ) = Y( IY ) + ALPHA*CABS1( X( JX ) )*TEMP
 365                      JX = JX + INCX
 366                   END DO
 367                   DO J = I+1, N
 368                      TEMP = CABS1( A( I, J ) )
 369                      SYMB_ZERO = SYMB_ZERO .AND.
 370      $                    ( X( J ) .EQ. ZERO .OR. TEMP .EQ. ZERO )
 371
 372                      Y( IY ) = Y( IY ) + ALPHA*CABS1( X( JX ) )*TEMP
 373                      JX = JX + INCX
 374                   END DO
 375                END IF
 376
 377                IF ( .NOT.SYMB_ZERO )
 378      $              Y( IY ) = Y( IY ) + SIGN( SAFE1, Y( IY ) )
 379
 380                IY = IY + INCY
 381             END DO
 382          ELSE
 383             DO I = 1, N
 384                IF ( BETA .EQ. ZERO ) THEN
 385                   SYMB_ZERO = .TRUE.
 386                   Y( IY ) = 0.0D+0
 387                ELSE IF ( Y( IY ) .EQ. ZERO ) THEN
 388                   SYMB_ZERO = .TRUE.
 389                ELSE
 390                   SYMB_ZERO = .FALSE.
 391                   Y( IY ) = BETA * ABS( Y( IY ) )
 392                END IF
 393                JX = KX
 394                IF ( ALPHA .NE. ZERO ) THEN
 395                   DO J = 1, I
 396                      TEMP = CABS1( A( I, J ) )
 397                      SYMB_ZERO = SYMB_ZERO .AND.
 398      $                    ( X( J ) .EQ. ZERO .OR. TEMP .EQ. ZERO )
 399
 400                      Y( IY ) = Y( IY ) + ALPHA*CABS1( X( JX ) )*TEMP
 401                      JX = JX + INCX
 402                   END DO
 403                   DO J = I+1, N
 404                      TEMP = CABS1( A( J, I ) )
 405                      SYMB_ZERO = SYMB_ZERO .AND.
 406      $                    ( X( J ) .EQ. ZERO .OR. TEMP .EQ. ZERO )
 407
 408                      Y( IY ) = Y( IY ) + ALPHA*CABS1( X( JX ) )*TEMP
 409                      JX = JX + INCX
 410                   END DO
 411                END IF
 412
 413                IF ( .NOT.SYMB_ZERO )
 414      $              Y( IY ) = Y( IY ) + SIGN( SAFE1, Y( IY ) )
 415
 416                IY = IY + INCY
 417             END DO
 418          END IF
 419
 420       END IF
 421 *
 422       RETURN
 423 *
 424 *     End of ZLA_HEAMV
 425 *
 426       END