SRC/cla_geamv.f

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