1       SUBROUTINE DLA_GBAMV( TRANS, M, N, KL, KU, ALPHA, AB, LDAB, X,
  2      $                      INCX, BETA, Y, INCY )
  3 *
  4 *     -- LAPACK routine (version 3.3.1)                                 --
  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, LDAB, M, N, KL, KU, TRANS
 17 *     ..
 18 *     .. Array Arguments ..
 19       DOUBLE PRECISION   AB( LDAB, * ), X( * ), Y( * )
 20 *     ..
 21 *
 22 *  Purpose
 23 *  =======
 24 *
 25 *  DLA_GBAMV  performs one of the matrix-vector operations
 26 *
 27 *          y := alpha*abs(A)*abs(x) + beta*abs(y),
 28 *     or   y := alpha*abs(A)**T*abs(x) + beta*abs(y),
 29 *
 30 *  where alpha and beta are scalars, x and y are vectors and A is an
 31 *  m by n matrix.
 32 *
 33 *  This function is primarily used in calculating error bounds.
 34 *  To protect against underflow during evaluation, components in
 35 *  the resulting vector are perturbed away from zero by (N+1)
 36 *  times the underflow threshold.  To prevent unnecessarily large
 37 *  errors for block-structure embedded in general matrices,
 38 *  "symbolically" zero components are not perturbed.  A zero
 39 *  entry is considered "symbolic" if all multiplications involved
 40 *  in computing that entry have at least one zero multiplicand.
 41 *
 42 *  Arguments
 43 *  ==========
 44 *
 45 *  TRANS   (input) INTEGER
 46 *           On entry, TRANS specifies the operation to be performed as
 47 *           follows:
 48 *
 49 *             BLAS_NO_TRANS      y := alpha*abs(A)*abs(x) + beta*abs(y)
 50 *             BLAS_TRANS         y := alpha*abs(A**T)*abs(x) + beta*abs(y)
 51 *             BLAS_CONJ_TRANS    y := alpha*abs(A**T)*abs(x) + beta*abs(y)
 52 *
 53 *           Unchanged on exit.
 54 *
 55 *  M        (input) INTEGER
 56 *           On entry, M specifies the number of rows of the matrix A.
 57 *           M must be at least zero.
 58 *           Unchanged on exit.
 59 *
 60 *  N        (input) INTEGER
 61 *           On entry, N specifies the number of columns of the matrix A.
 62 *           N must be at least zero.
 63 *           Unchanged on exit.
 64 *
 65 *  KL       (input) INTEGER
 66 *           The number of subdiagonals within the band of A.  KL >= 0.
 67 *
 68 *  KU       (input) INTEGER
 69 *           The number of superdiagonals within the band of A.  KU >= 0.
 70 *
 71 *  ALPHA    (input) DOUBLE PRECISION
 72 *           On entry, ALPHA specifies the scalar alpha.
 73 *           Unchanged on exit.
 74 *
 75 *  AB       (input) DOUBLE PRECISION array of DIMENSION ( LDAB, n )
 76 *           Before entry, the leading m by n part of the array AB must
 77 *           contain the matrix of coefficients.
 78 *           Unchanged on exit.
 79 *
 80 *  LDAB     (input) INTEGER
 81 *           On entry, LDA specifies the first dimension of AB as declared
 82 *           in the calling (sub) program. LDAB must be at least
 83 *           max( 1, m ).
 84 *           Unchanged on exit.
 85 *
 86 *  X        (input) DOUBLE PRECISION array, dimension
 87 *           ( 1 + ( n - 1 )*abs( INCX ) ) when TRANS = 'N' or 'n'
 88 *           and at least
 89 *           ( 1 + ( m - 1 )*abs( INCX ) ) otherwise.
 90 *           Before entry, the incremented array X must contain the
 91 *           vector x.
 92 *           Unchanged on exit.
 93 *
 94 *  INCX     (input) INTEGER
 95 *           On entry, INCX specifies the increment for the elements of
 96 *           X. INCX must not be zero.
 97 *           Unchanged on exit.
 98 *
 99 *  BETA     (input)  DOUBLE PRECISION
100 *           On entry, BETA specifies the scalar beta. When BETA is
101 *           supplied as zero then Y need not be set on input.
102 *           Unchanged on exit.
103 *
104 *  Y        (input/output) DOUBLE PRECISION  array, dimension
105 *           ( 1 + ( m - 1 )*abs( INCY ) ) when TRANS = 'N' or 'n'
106 *           and at least
107 *           ( 1 + ( n - 1 )*abs( INCY ) ) otherwise.
108 *           Before entry with BETA non-zero, the incremented array Y
109 *           must contain the vector y. On exit, Y is overwritten by the
110 *           updated vector y.
111 *
112 *  INCY     (input) INTEGER
113 *           On entry, INCY specifies the increment for the elements of
114 *           Y. INCY must not be zero.
115 *           Unchanged on exit.
116 *
117 *
118 *  Level 2 Blas routine.
119 *
120 *  =====================================================================
121 *
122 *     .. Parameters ..
123       DOUBLE PRECISION   ONE, ZERO
124       PARAMETER          ( ONE = 1.0D+0, ZERO = 0.0D+0 )
125 *     ..
126 *     .. Local Scalars ..
127       LOGICAL            SYMB_ZERO
128       DOUBLE PRECISION   TEMP, SAFE1
129       INTEGER            I, INFO, IY, J, JX, KX, KY, LENX, LENY, KD, KE
130 *     ..
131 *     .. External Subroutines ..
132       EXTERNAL           XERBLA, DLAMCH
133       DOUBLE PRECISION   DLAMCH
134 *     ..
135 *     .. External Functions ..
136       EXTERNAL           ILATRANS
137       INTEGER            ILATRANS
138 *     ..
139 *     .. Intrinsic Functions ..
140       INTRINSIC          MAXABSSIGN
141 *     ..
142 *     .. Executable Statements ..
143 *
144 *     Test the input parameters.
145 *
146       INFO = 0
147       IF     ( .NOT.( ( TRANS.EQ.ILATRANS( 'N' ) )
148      $           .OR. ( TRANS.EQ.ILATRANS( 'T' ) )
149      $           .OR. ( TRANS.EQ.ILATRANS( 'C' ) ) ) ) THEN
150          INFO = 1
151       ELSE IF( M.LT.0 )THEN
152          INFO = 2
153       ELSE IF( N.LT.0 )THEN
154          INFO = 3
155       ELSE IF( KL.LT.0 .OR. KL.GT.M-1 ) THEN
156          INFO = 4
157       ELSE IF( KU.LT.0 .OR. KU.GT.N-1 ) THEN
158          INFO = 5
159       ELSE IF( LDAB.LT.KL+KU+1 )THEN
160          INFO = 6
161       ELSE IF( INCX.EQ.0 )THEN
162          INFO = 8
163       ELSE IF( INCY.EQ.0 )THEN
164          INFO = 11
165       END IF
166       IF( INFO.NE.0 )THEN
167          CALL XERBLA( 'DLA_GBAMV ', INFO )
168          RETURN
169       END IF
170 *
171 *     Quick return if possible.
172 *
173       IF( ( M.EQ.0 ).OR.( N.EQ.0 ).OR.
174      $    ( ( ALPHA.EQ.ZERO ).AND.( BETA.EQ.ONE ) ) )
175      $   RETURN
176 *
177 *     Set  LENX  and  LENY, the lengths of the vectors x and y, and set
178 *     up the start points in  X  and  Y.
179 *
180       IF( TRANS.EQ.ILATRANS( 'N' ) )THEN
181          LENX = N
182          LENY = M
183       ELSE
184          LENX = M
185          LENY = N
186       END IF
187       IF( INCX.GT.0 )THEN
188          KX = 1
189       ELSE
190          KX = 1 - ( LENX - 1 )*INCX
191       END IF
192       IF( INCY.GT.0 )THEN
193          KY = 1
194       ELSE
195          KY = 1 - ( LENY - 1 )*INCY
196       END IF
197 *
198 *     Set SAFE1 essentially to be the underflow threshold times the
199 *     number of additions in each row.
200 *
201       SAFE1 = DLAMCH( 'Safe minimum' )
202       SAFE1 = (N+1)*SAFE1
203 *
204 *     Form  y := alpha*abs(A)*abs(x) + beta*abs(y).
205 *
206 *     The O(M*N) SYMB_ZERO tests could be replaced by O(N) queries to
207 *     the inexact flag.  Still doesn't help change the iteration order
208 *     to per-column.
209 *
210       KD = KU + 1
211       KE = KL + 1
212       IY = KY
213       IF ( INCX.EQ.1 ) THEN
214          IF( TRANS.EQ.ILATRANS( 'N' ) )THEN
215             DO I = 1, LENY
216                IF ( BETA .EQ. ZERO ) THEN
217                   SYMB_ZERO = .TRUE.
218                   Y( IY ) = 0.0D+0
219                ELSE IF ( Y( IY ) .EQ. ZERO ) THEN
220                   SYMB_ZERO = .TRUE.
221                ELSE
222                   SYMB_ZERO = .FALSE.
223                   Y( IY ) = BETA * ABS( Y( IY ) )
224                END IF
225                IF ( ALPHA .NE. ZERO ) THEN
226                   DO J = MAX( I-KL, 1 ), MIN( I+KU, LENX )
227                      TEMP = ABS( AB( KD+I-J, J ) )
228                      SYMB_ZERO = SYMB_ZERO .AND.
229      $                    ( X( J ) .EQ. ZERO .OR. TEMP .EQ. ZERO )
230 
231                      Y( IY ) = Y( IY ) + ALPHA*ABS( X( J ) )*TEMP
232                   END DO
233                END IF
234 
235                IF ( .NOT.SYMB_ZERO )
236      $              Y( IY ) = Y( IY ) + SIGN( SAFE1, Y( IY ) )
237                IY = IY + INCY
238             END DO
239          ELSE
240             DO I = 1, LENY
241                IF ( BETA .EQ. ZERO ) THEN
242                   SYMB_ZERO = .TRUE.
243                   Y( IY ) = 0.0D+0
244                ELSE IF ( Y( IY ) .EQ. ZERO ) THEN
245                   SYMB_ZERO = .TRUE.
246                ELSE
247                   SYMB_ZERO = .FALSE.
248                   Y( IY ) = BETA * ABS( Y( IY ) )
249                END IF
250                IF ( ALPHA .NE. ZERO ) THEN
251                   DO J = MAX( I-KL, 1 ), MIN( I+KU, LENX )
252                      TEMP = ABS( AB( KE-I+J, I ) )
253                      SYMB_ZERO = SYMB_ZERO .AND.
254      $                    ( X( J ) .EQ. ZERO .OR. TEMP .EQ. ZERO )
255 
256                      Y( IY ) = Y( IY ) + ALPHA*ABS( X( J ) )*TEMP
257                   END DO
258                END IF
259 
260                IF ( .NOT.SYMB_ZERO )
261      $              Y( IY ) = Y( IY ) + SIGN( SAFE1, Y( IY ) )
262                IY = IY + INCY
263             END DO
264          END IF
265       ELSE
266          IF( TRANS.EQ.ILATRANS( 'N' ) )THEN
267             DO I = 1, LENY
268                IF ( BETA .EQ. ZERO ) THEN
269                   SYMB_ZERO = .TRUE.
270                   Y( IY ) = 0.0D+0
271                ELSE IF ( Y( IY ) .EQ. ZERO ) THEN
272                   SYMB_ZERO = .TRUE.
273                ELSE
274                   SYMB_ZERO = .FALSE.
275                   Y( IY ) = BETA * ABS( Y( IY ) )
276                END IF
277                IF ( ALPHA .NE. ZERO ) THEN
278                   JX = KX
279                   DO J = MAX( I-KL, 1 ), MIN( I+KU, LENX )
280                      TEMP = ABS( AB( KD+I-J, J ) )
281                      SYMB_ZERO = SYMB_ZERO .AND.
282      $                    ( X( JX ) .EQ. ZERO .OR. TEMP .EQ. ZERO )
283 
284                      Y( IY ) = Y( IY ) + ALPHA*ABS( X( JX ) )*TEMP
285                      JX = JX + INCX
286                   END DO
287                END IF
288 
289                IF ( .NOT.SYMB_ZERO )
290      $              Y( IY ) = Y( IY ) + SIGN( SAFE1, Y( IY ) )
291 
292                IY = IY + INCY
293             END DO
294          ELSE
295             DO I = 1, LENY
296                IF ( BETA .EQ. ZERO ) THEN
297                   SYMB_ZERO = .TRUE.
298                   Y( IY ) = 0.0D+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                   JX = KX
307                   DO J = MAX( I-KL, 1 ), MIN( I+KU, LENX )
308                      TEMP = ABS( AB( KE-I+J, I ) )
309                      SYMB_ZERO = SYMB_ZERO .AND.
310      $                    ( X( JX ) .EQ. ZERO .OR. TEMP .EQ. ZERO )
311 
312                      Y( IY ) = Y( IY ) + ALPHA*ABS( X( JX ) )*TEMP
313                      JX = JX + INCX
314                   END DO
315                END IF
316 
317                IF ( .NOT.SYMB_ZERO )
318      $              Y( IY ) = Y( IY ) + SIGN( SAFE1, Y( IY ) )
319 
320                IY = IY + INCY
321             END DO
322          END IF
323 
324       END IF
325 *
326       RETURN
327 *
328 *     End of DLA_GBAMV
329 *
330       END