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