SRC/dla_geamv.f

   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