1*> \brief \b DLA_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 DLA_GBAMV + dependencies
10*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/dla_gbamv.f">
11*> [TGZ]</a>
12*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/dla_gbamv.f">
13*> [ZIP]</a>
14*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/dla_gbamv.f">
15*> [TXT]</a>
16*> \endhtmlonly
17*
18*  Definition:
19*  ===========
20*
21*       SUBROUTINE DLA_GBAMV( TRANS, M, N, KL, KU, ALPHA, AB, LDAB, X,
22*                             INCX, BETA, Y, INCY )
23*
24*       .. Scalar Arguments ..
25*       DOUBLE PRECISION   ALPHA, BETA
26*       INTEGER            INCX, INCY, LDAB, M, N, KL, KU, TRANS
27*       ..
28*       .. Array Arguments ..
29*       DOUBLE PRECISION   AB( LDAB, * ), X( * ), Y( * )
30*       ..
31*
32*
33*> \par Purpose:
34*  =============
35*>
36*> \verbatim
37*>
38*> DLA_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 DOUBLE PRECISION
103*>           On entry, ALPHA specifies the scalar alpha.
104*>           Unchanged on exit.
105*> \endverbatim
106*>
107*> \param[in] AB
108*> \verbatim
109*>          AB is DOUBLE PRECISION array, 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 DOUBLE PRECISION 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 DOUBLE PRECISION
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 DOUBLE PRECISION 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*> \ingroup doubleGBcomputational
181*
182*  =====================================================================
183      SUBROUTINE DLA_GBAMV( TRANS, M, N, KL, KU, ALPHA, AB, LDAB, X,
184     $                      INCX, BETA, Y, INCY )
185*
186*  -- LAPACK computational routine --
187*  -- LAPACK is a software package provided by Univ. of Tennessee,    --
188*  -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
189*
190*     .. Scalar Arguments ..
191      DOUBLE PRECISION   ALPHA, BETA
192      INTEGER            INCX, INCY, LDAB, M, N, KL, KU, TRANS
193*     ..
194*     .. Array Arguments ..
195      DOUBLE PRECISION   AB( LDAB, * ), X( * ), Y( * )
196*     ..
197*
198*  =====================================================================
199*
200*     .. Parameters ..
201      DOUBLE PRECISION   ONE, ZERO
202      PARAMETER          ( ONE = 1.0D+0, ZERO = 0.0D+0 )
203*     ..
204*     .. Local Scalars ..
205      LOGICAL            SYMB_ZERO
206      DOUBLE PRECISION   TEMP, SAFE1
207      INTEGER            I, INFO, IY, J, JX, KX, KY, LENX, LENY, KD, KE
208*     ..
209*     .. External Subroutines ..
210      EXTERNAL           XERBLA, DLAMCH
211      DOUBLE PRECISION   DLAMCH
212*     ..
213*     .. External Functions ..
214      EXTERNAL           ILATRANS
215      INTEGER            ILATRANS
216*     ..
217*     .. Intrinsic Functions ..
218      INTRINSIC          MAX, ABS, SIGN
219*     ..
220*     .. Executable Statements ..
221*
222*     Test the input parameters.
223*
224      INFO = 0
225      IF     ( .NOT.( ( TRANS.EQ.ILATRANS( 'N' ) )
226     $           .OR. ( TRANS.EQ.ILATRANS( 'T' ) )
227     $           .OR. ( TRANS.EQ.ILATRANS( 'C' ) ) ) ) THEN
228         INFO = 1
229      ELSE IF( M.LT.0 )THEN
230         INFO = 2
231      ELSE IF( N.LT.0 )THEN
232         INFO = 3
233      ELSE IF( KL.LT.0 .OR. KL.GT.M-1 ) THEN
234         INFO = 4
235      ELSE IF( KU.LT.0 .OR. KU.GT.N-1 ) THEN
236         INFO = 5
237      ELSE IF( LDAB.LT.KL+KU+1 )THEN
238         INFO = 6
239      ELSE IF( INCX.EQ.0 )THEN
240         INFO = 8
241      ELSE IF( INCY.EQ.0 )THEN
242         INFO = 11
243      END IF
244      IF( INFO.NE.0 )THEN
245         CALL XERBLA( 'DLA_GBAMV ', INFO )
246         RETURN
247      END IF
248*
249*     Quick return if possible.
250*
251      IF( ( M.EQ.0 ).OR.( N.EQ.0 ).OR.
252     $    ( ( ALPHA.EQ.ZERO ).AND.( BETA.EQ.ONE ) ) )
253     $   RETURN
254*
255*     Set  LENX  and  LENY, the lengths of the vectors x and y, and set
256*     up the start points in  X  and  Y.
257*
258      IF( TRANS.EQ.ILATRANS( 'N' ) )THEN
259         LENX = N
260         LENY = M
261      ELSE
262         LENX = M
263         LENY = N
264      END IF
265      IF( INCX.GT.0 )THEN
266         KX = 1
267      ELSE
268         KX = 1 - ( LENX - 1 )*INCX
269      END IF
270      IF( INCY.GT.0 )THEN
271         KY = 1
272      ELSE
273         KY = 1 - ( LENY - 1 )*INCY
274      END IF
275*
276*     Set SAFE1 essentially to be the underflow threshold times the
277*     number of additions in each row.
278*
279      SAFE1 = DLAMCH( 'Safe minimum' )
280      SAFE1 = (N+1)*SAFE1
281*
282*     Form  y := alpha*abs(A)*abs(x) + beta*abs(y).
283*
284*     The O(M*N) SYMB_ZERO tests could be replaced by O(N) queries to
285*     the inexact flag.  Still doesn't help change the iteration order
286*     to per-column.
287*
288      KD = KU + 1
289      KE = KL + 1
290      IY = KY
291      IF ( INCX.EQ.1 ) THEN
292         IF( TRANS.EQ.ILATRANS( 'N' ) )THEN
293            DO I = 1, LENY
294               IF ( BETA .EQ. ZERO ) THEN
295                  SYMB_ZERO = .TRUE.
296                  Y( IY ) = 0.0D+0
297               ELSE IF ( Y( IY ) .EQ. ZERO ) THEN
298                  SYMB_ZERO = .TRUE.
299               ELSE
300                  SYMB_ZERO = .FALSE.
301                  Y( IY ) = BETA * ABS( Y( IY ) )
302               END IF
303               IF ( ALPHA .NE. ZERO ) THEN
304                  DO J = MAX( I-KL, 1 ), MIN( I+KU, LENX )
305                     TEMP = ABS( AB( KD+I-J, J ) )
306                     SYMB_ZERO = SYMB_ZERO .AND.
307     $                    ( X( J ) .EQ. ZERO .OR. TEMP .EQ. ZERO )
308
309                     Y( IY ) = Y( IY ) + ALPHA*ABS( X( J ) )*TEMP
310                  END DO
311               END IF
312
313               IF ( .NOT.SYMB_ZERO )
314     $              Y( IY ) = Y( IY ) + SIGN( SAFE1, Y( IY ) )
315               IY = IY + INCY
316            END DO
317         ELSE
318            DO I = 1, LENY
319               IF ( BETA .EQ. ZERO ) THEN
320                  SYMB_ZERO = .TRUE.
321                  Y( IY ) = 0.0D+0
322               ELSE IF ( Y( IY ) .EQ. ZERO ) THEN
323                  SYMB_ZERO = .TRUE.
324               ELSE
325                  SYMB_ZERO = .FALSE.
326                  Y( IY ) = BETA * ABS( Y( IY ) )
327               END IF
328               IF ( ALPHA .NE. ZERO ) THEN
329                  DO J = MAX( I-KL, 1 ), MIN( I+KU, LENX )
330                     TEMP = ABS( AB( KE-I+J, I ) )
331                     SYMB_ZERO = SYMB_ZERO .AND.
332     $                    ( X( J ) .EQ. ZERO .OR. TEMP .EQ. ZERO )
333
334                     Y( IY ) = Y( IY ) + ALPHA*ABS( X( J ) )*TEMP
335                  END DO
336               END IF
337
338               IF ( .NOT.SYMB_ZERO )
339     $              Y( IY ) = Y( IY ) + SIGN( SAFE1, Y( IY ) )
340               IY = IY + INCY
341            END DO
342         END IF
343      ELSE
344         IF( TRANS.EQ.ILATRANS( 'N' ) )THEN
345            DO I = 1, LENY
346               IF ( BETA .EQ. ZERO ) THEN
347                  SYMB_ZERO = .TRUE.
348                  Y( IY ) = 0.0D+0
349               ELSE IF ( Y( IY ) .EQ. ZERO ) THEN
350                  SYMB_ZERO = .TRUE.
351               ELSE
352                  SYMB_ZERO = .FALSE.
353                  Y( IY ) = BETA * ABS( Y( IY ) )
354               END IF
355               IF ( ALPHA .NE. ZERO ) THEN
356                  JX = KX
357                  DO J = MAX( I-KL, 1 ), MIN( I+KU, LENX )
358                     TEMP = ABS( AB( KD+I-J, J ) )
359                     SYMB_ZERO = SYMB_ZERO .AND.
360     $                    ( X( JX ) .EQ. ZERO .OR. TEMP .EQ. ZERO )
361
362                     Y( IY ) = Y( IY ) + ALPHA*ABS( X( JX ) )*TEMP
363                     JX = JX + INCX
364                  END DO
365               END IF
366
367               IF ( .NOT.SYMB_ZERO )
368     $              Y( IY ) = Y( IY ) + SIGN( SAFE1, Y( IY ) )
369
370               IY = IY + INCY
371            END DO
372         ELSE
373            DO I = 1, LENY
374               IF ( BETA .EQ. ZERO ) THEN
375                  SYMB_ZERO = .TRUE.
376                  Y( IY ) = 0.0D+0
377               ELSE IF ( Y( IY ) .EQ. ZERO ) THEN
378                  SYMB_ZERO = .TRUE.
379               ELSE
380                  SYMB_ZERO = .FALSE.
381                  Y( IY ) = BETA * ABS( Y( IY ) )
382               END IF
383               IF ( ALPHA .NE. ZERO ) THEN
384                  JX = KX
385                  DO J = MAX( I-KL, 1 ), MIN( I+KU, LENX )
386                     TEMP = ABS( AB( KE-I+J, I ) )
387                     SYMB_ZERO = SYMB_ZERO .AND.
388     $                    ( X( JX ) .EQ. ZERO .OR. TEMP .EQ. ZERO )
389
390                     Y( IY ) = Y( IY ) + ALPHA*ABS( X( JX ) )*TEMP
391                     JX = JX + INCX
392                  END DO
393               END IF
394
395               IF ( .NOT.SYMB_ZERO )
396     $              Y( IY ) = Y( IY ) + SIGN( SAFE1, Y( IY ) )
397
398               IY = IY + INCY
399            END DO
400         END IF
401
402      END IF
403*
404      RETURN
405*
406*     End of DLA_GBAMV
407*
408      END
409