SRC/sla_gbamv.f

   1 *> \brief \b SLA_GBAMV performs a matrix-vector operation 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 SLA_GBAMV + dependencies
  10 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/sla_gbamv.f">
  11 *> [TGZ]</a>
  12 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/sla_gbamv.f">
  13 *> [ZIP]</a>
  14 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/sla_gbamv.f">
  15 *> [TXT]</a>
  16 *> \endhtmlonly
  17 *
  18 *  Definition:
  19 *  ===========
  20 *
  21 *       SUBROUTINE SLA_GBAMV( TRANS, M, N, KL, KU, ALPHA, AB, LDAB, X,
  22 *                             INCX, BETA, Y, INCY )
  23 *
  24 *       .. Scalar Arguments ..
  25 *       REAL               ALPHA, BETA
  26 *       INTEGER            INCX, INCY, LDAB, M, N, KL, KU, TRANS
  27 *       ..
  28 *       .. Array Arguments ..
  29 *       REAL               AB( LDAB, * ), X( * ), Y( * )
  30 *       ..
  31 *
  32 *
  33 *> \par Purpose:
  34 *  =============
  35 *>
  36 *> \verbatim
  37 *>
  38 *> SLA_GBAMV  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] KL
  89 *> \verbatim
  90 *>          KL is INTEGER
  91 *>           The number of subdiagonals within the band of A.  KL >= 0.
  92 *> \endverbatim
  93 *>
  94 *> \param[in] KU
  95 *> \verbatim
  96 *>          KU is INTEGER
  97 *>           The number of superdiagonals within the band of A.  KU >= 0.
  98 *> \endverbatim
  99 *>
 100 *> \param[in] ALPHA
 101 *> \verbatim
 102 *>          ALPHA is REAL
 103 *>           On entry, ALPHA specifies the scalar alpha.
 104 *>           Unchanged on exit.
 105 *> \endverbatim
 106 *>
 107 *> \param[in] AB
 108 *> \verbatim
 109 *>          AB is REAL array of DIMENSION ( LDAB, n )
 110 *>           Before entry, the leading m by n part of the array AB must
 111 *>           contain the matrix of coefficients.
 112 *>           Unchanged on exit.
 113 *> \endverbatim
 114 *>
 115 *> \param[in] LDAB
 116 *> \verbatim
 117 *>          LDAB is INTEGER
 118 *>           On entry, LDA specifies the first dimension of AB as declared
 119 *>           in the calling (sub) program. LDAB must be at least
 120 *>           max( 1, m ).
 121 *>           Unchanged on exit.
 122 *> \endverbatim
 123 *>
 124 *> \param[in] X
 125 *> \verbatim
 126 *>          X is REAL array, dimension
 127 *>           ( 1 + ( n - 1 )*abs( INCX ) ) when TRANS = 'N' or 'n'
 128 *>           and at least
 129 *>           ( 1 + ( m - 1 )*abs( INCX ) ) otherwise.
 130 *>           Before entry, the incremented array X must contain the
 131 *>           vector x.
 132 *>           Unchanged on exit.
 133 *> \endverbatim
 134 *>
 135 *> \param[in] INCX
 136 *> \verbatim
 137 *>          INCX is INTEGER
 138 *>           On entry, INCX specifies the increment for the elements of
 139 *>           X. INCX must not be zero.
 140 *>           Unchanged on exit.
 141 *> \endverbatim
 142 *>
 143 *> \param[in] BETA
 144 *> \verbatim
 145 *>          BETA is REAL
 146 *>           On entry, BETA specifies the scalar beta. When BETA is
 147 *>           supplied as zero then Y need not be set on input.
 148 *>           Unchanged on exit.
 149 *> \endverbatim
 150 *>
 151 *> \param[in,out] Y
 152 *> \verbatim
 153 *>          Y is REAL array, dimension
 154 *>           ( 1 + ( m - 1 )*abs( INCY ) ) when TRANS = 'N' or 'n'
 155 *>           and at least
 156 *>           ( 1 + ( n - 1 )*abs( INCY ) ) otherwise.
 157 *>           Before entry with BETA non-zero, the incremented array Y
 158 *>           must contain the vector y. On exit, Y is overwritten by the
 159 *>           updated vector y.
 160 *> \endverbatim
 161 *>
 162 *> \param[in] INCY
 163 *> \verbatim
 164 *>          INCY is INTEGER
 165 *>           On entry, INCY specifies the increment for the elements of
 166 *>           Y. INCY must not be zero.
 167 *>           Unchanged on exit.
 168 *>
 169 *>  Level 2 Blas routine.
 170 *> \endverbatim
 171 *
 172 *  Authors:
 173 *  ========
 174 *
 175 *> \author Univ. of Tennessee
 176 *> \author Univ. of California Berkeley
 177 *> \author Univ. of Colorado Denver
 178 *> \author NAG Ltd.
 179 *
 180 *> \date September 2012
 181 *
 182 *> \ingroup realGBcomputational
 183 *
 184 *  =====================================================================
 185       SUBROUTINE SLA_GBAMV( TRANS, M, N, KL, KU, ALPHA, AB, LDAB, X,
 186      $                      INCX, BETA, Y, INCY )
 187 *
 188 *  -- LAPACK computational routine (version 3.4.2) --
 189 *  -- LAPACK is a software package provided by Univ. of Tennessee,    --
 190 *  -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
 191 *     September 2012
 192 *
 193 *     .. Scalar Arguments ..
 194       REAL               ALPHA, BETA
 195       INTEGER            INCX, INCY, LDAB, M, N, KL, KU, TRANS
 196 *     ..
 197 *     .. Array Arguments ..
 198       REAL               AB( LDAB, * ), X( * ), Y( * )
 199 *     ..
 200 *
 201 *  =====================================================================
 202 *     .. Parameters ..
 203       REAL               ONE, ZERO
 204       PARAMETER          ( ONE = 1.0E+0, ZERO = 0.0E+0 )
 205 *     ..
 206 *     .. Local Scalars ..
 207       LOGICAL            SYMB_ZERO
 208       REAL               TEMP, SAFE1
 209       INTEGER            I, INFO, IY, J, JX, KX, KY, LENX, LENY, KD, KE
 210 *     ..
 211 *     .. External Subroutines ..
 212       EXTERNAL           XERBLA, SLAMCH
 213       REAL               SLAMCH
 214 *     ..
 215 *     .. External Functions ..
 216       EXTERNAL           ILATRANS
 217       INTEGER            ILATRANS
 218 *     ..
 219 *     .. Intrinsic Functions ..
 220       INTRINSIC          MAX, ABS, SIGN
 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( KL.LT.0 .OR. KL.GT.M-1 ) THEN
 236          INFO = 4
 237       ELSE IF( KU.LT.0 .OR. KU.GT.N-1 ) THEN
 238          INFO = 5
 239       ELSE IF( LDAB.LT.KL+KU+1 )THEN
 240          INFO = 6
 241       ELSE IF( INCX.EQ.0 )THEN
 242          INFO = 8
 243       ELSE IF( INCY.EQ.0 )THEN
 244          INFO = 11
 245       END IF
 246       IF( INFO.NE.0 )THEN
 247          CALL XERBLA( 'SLA_GBAMV ', INFO )
 248          RETURN
 249       END IF
 250 *
 251 *     Quick return if possible.
 252 *
 253       IF( ( M.EQ.0 ).OR.( N.EQ.0 ).OR.
 254      $    ( ( ALPHA.EQ.ZERO ).AND.( BETA.EQ.ONE ) ) )
 255      $   RETURN
 256 *
 257 *     Set  LENX  and  LENY, the lengths of the vectors x and y, and set
 258 *     up the start points in  X  and  Y.
 259 *
 260       IF( TRANS.EQ.ILATRANS( 'N' ) )THEN
 261          LENX = N
 262          LENY = M
 263       ELSE
 264          LENX = M
 265          LENY = N
 266       END IF
 267       IF( INCX.GT.0 )THEN
 268          KX = 1
 269       ELSE
 270          KX = 1 - ( LENX - 1 )*INCX
 271       END IF
 272       IF( INCY.GT.0 )THEN
 273          KY = 1
 274       ELSE
 275          KY = 1 - ( LENY - 1 )*INCY
 276       END IF
 277 *
 278 *     Set SAFE1 essentially to be the underflow threshold times the
 279 *     number of additions in each row.
 280 *
 281       SAFE1 = SLAMCH( 'Safe minimum' )
 282       SAFE1 = (N+1)*SAFE1
 283 *
 284 *     Form  y := alpha*abs(A)*abs(x) + beta*abs(y).
 285 *
 286 *     The O(M*N) SYMB_ZERO tests could be replaced by O(N) queries to
 287 *     the inexact flag.  Still doesn't help change the iteration order
 288 *     to per-column.
 289 *
 290       KD = KU + 1
 291       KE = KL + 1
 292       IY = KY
 293       IF ( INCX.EQ.1 ) THEN
 294          IF( TRANS.EQ.ILATRANS( 'N' ) )THEN
 295             DO I = 1, LENY
 296                IF ( BETA .EQ. ZERO ) THEN
 297                   SYMB_ZERO = .TRUE.
 298                   Y( IY ) = 0.0
 299                ELSE IF ( Y( IY ) .EQ. ZERO ) THEN
 300                   SYMB_ZERO = .TRUE.
 301                ELSE
 302                   SYMB_ZERO = .FALSE.
 303                   Y( IY ) = BETA * ABS( Y( IY ) )
 304                END IF
 305                IF ( ALPHA .NE. ZERO ) THEN
 306                   DO J = MAX( I-KL, 1 ), MIN( I+KU, LENX )
 307                      TEMP = ABS( AB( KD+I-J, J ) )
 308                      SYMB_ZERO = SYMB_ZERO .AND.
 309      $                    ( X( J ) .EQ. ZERO .OR. TEMP .EQ. ZERO )
 310
 311                      Y( IY ) = Y( IY ) + ALPHA*ABS( X( J ) )*TEMP
 312                   END DO
 313                END IF
 314
 315                IF ( .NOT.SYMB_ZERO )
 316      $              Y( IY ) = Y( IY ) + SIGN( SAFE1, Y( IY ) )
 317                IY = IY + INCY
 318             END DO
 319          ELSE
 320             DO I = 1, LENY
 321                IF ( BETA .EQ. ZERO ) THEN
 322                   SYMB_ZERO = .TRUE.
 323                   Y( IY ) = 0.0
 324                ELSE IF ( Y( IY ) .EQ. ZERO ) THEN
 325                   SYMB_ZERO = .TRUE.
 326                ELSE
 327                   SYMB_ZERO = .FALSE.
 328                   Y( IY ) = BETA * ABS( Y( IY ) )
 329                END IF
 330                IF ( ALPHA .NE. ZERO ) THEN
 331                   DO J = MAX( I-KL, 1 ), MIN( I+KU, LENX )
 332                      TEMP = ABS( AB( KE-I+J, I ) )
 333                      SYMB_ZERO = SYMB_ZERO .AND.
 334      $                    ( X( J ) .EQ. ZERO .OR. TEMP .EQ. ZERO )
 335
 336                      Y( IY ) = Y( IY ) + ALPHA*ABS( X( J ) )*TEMP
 337                   END DO
 338                END IF
 339
 340                IF ( .NOT.SYMB_ZERO )
 341      $              Y( IY ) = Y( IY ) + SIGN( SAFE1, Y( IY ) )
 342                IY = IY + INCY
 343             END DO
 344          END IF
 345       ELSE
 346          IF( TRANS.EQ.ILATRANS( 'N' ) )THEN
 347             DO I = 1, LENY
 348                IF ( BETA .EQ. ZERO ) THEN
 349                   SYMB_ZERO = .TRUE.
 350                   Y( IY ) = 0.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                IF ( ALPHA .NE. ZERO ) THEN
 358                   JX = KX
 359                   DO J = MAX( I-KL, 1 ), MIN( I+KU, LENX )
 360                      TEMP = ABS( AB( KD+I-J, J ) )
 361                      SYMB_ZERO = SYMB_ZERO .AND.
 362      $                    ( X( JX ) .EQ. ZERO .OR. TEMP .EQ. ZERO )
 363
 364                      Y( IY ) = Y( IY ) + ALPHA*ABS( X( JX ) )*TEMP
 365                      JX = JX + INCX
 366                   END DO
 367                END IF
 368
 369                IF ( .NOT.SYMB_ZERO )
 370      $           Y( IY ) = Y( IY ) + SIGN( SAFE1, Y( IY ) )
 371
 372                IY = IY + INCY
 373             END DO
 374          ELSE
 375             DO I = 1, LENY
 376                IF ( BETA .EQ. ZERO ) THEN
 377                   SYMB_ZERO = .TRUE.
 378                   Y( IY ) = 0.0
 379                ELSE IF ( Y( IY ) .EQ. ZERO ) THEN
 380                   SYMB_ZERO = .TRUE.
 381                ELSE
 382                   SYMB_ZERO = .FALSE.
 383                   Y( IY ) = BETA * ABS( Y( IY ) )
 384                END IF
 385                IF ( ALPHA .NE. ZERO ) THEN
 386                   JX = KX
 387                   DO J = MAX( I-KL, 1 ), MIN( I+KU, LENX )
 388                      TEMP = ABS( AB( KE-I+J, I ) )
 389                      SYMB_ZERO = SYMB_ZERO .AND.
 390      $                    ( X( JX ) .EQ. ZERO .OR. TEMP .EQ. ZERO )
 391
 392                      Y( IY ) = Y( IY ) + ALPHA*ABS( X( JX ) )*TEMP
 393                      JX = JX + INCX
 394                   END DO
 395                END IF
 396
 397                IF ( .NOT.SYMB_ZERO )
 398      $           Y( IY ) = Y( IY ) + SIGN( SAFE1, Y( IY ) )
 399
 400                IY = IY + INCY
 401             END DO
 402          END IF
 403
 404       END IF
 405 *
 406       RETURN
 407 *
 408 *     End of SLA_GBAMV
 409 *
 410       END