001:       SUBROUTINE CLA_GBAMV( TRANS, M, N, KL, KU, ALPHA, AB, LDAB, X,
002:      $                      INCX, BETA, Y, INCY )
003: *
004: *     -- LAPACK routine (version 3.2)                                 --
005: *     -- Contributed by James Demmel, Deaglan Halligan, Yozo Hida and --
006: *     -- Jason Riedy of Univ. of California Berkeley.                 --
007: *     -- November 2008                                                --
008: *
009: *     -- LAPACK is a software package provided by Univ. of Tennessee, --
010: *     -- Univ. of California Berkeley and NAG Ltd.                    --
011: *
012:       IMPLICIT NONE
013: *     ..
014: *     .. Scalar Arguments ..
015:       REAL               ALPHA, BETA
016:       INTEGER            INCX, INCY, LDAB, M, N, KL, KU, TRANS
017: *     ..
018: *     .. Array Arguments ..
019:       COMPLEX            AB( LDAB, * ), X( * )
020:       REAL               Y( * )
021: *     ..
022: *
023: *  Purpose
024: *  =======
025: *
026: *  SLA_GEAMV  performs one of the matrix-vector operations
027: *
028: *          y := alpha*abs(A)*abs(x) + beta*abs(y),
029: *     or   y := alpha*abs(A)'*abs(x) + beta*abs(y),
030: *
031: *  where alpha and beta are scalars, x and y are vectors and A is an
032: *  m by n matrix.
033: *
034: *  This function is primarily used in calculating error bounds.
035: *  To protect against underflow during evaluation, components in
036: *  the resulting vector are perturbed away from zero by (N+1)
037: *  times the underflow threshold.  To prevent unnecessarily large
038: *  errors for block-structure embedded in general matrices,
039: *  "symbolically" zero components are not perturbed.  A zero
040: *  entry is considered "symbolic" if all multiplications involved
041: *  in computing that entry have at least one zero multiplicand.
042: *
043: *  Parameters
044: *  ==========
045: *
046: *  TRANS  - INTEGER
047: *           On entry, TRANS specifies the operation to be performed as
048: *           follows:
049: *
050: *             BLAS_NO_TRANS      y := alpha*abs(A)*abs(x) + beta*abs(y)
051: *             BLAS_TRANS         y := alpha*abs(A')*abs(x) + beta*abs(y)
052: *             BLAS_CONJ_TRANS    y := alpha*abs(A')*abs(x) + beta*abs(y)
053: *
054: *           Unchanged on exit.
055: *
056: *  M      - INTEGER
057: *           On entry, M specifies the number of rows of the matrix A.
058: *           M must be at least zero.
059: *           Unchanged on exit.
060: *
061: *  N      - INTEGER
062: *           On entry, N specifies the number of columns of the matrix A.
063: *           N must be at least zero.
064: *           Unchanged on exit.
065: *
066: *  KL     - INTEGER
067: *           The number of subdiagonals within the band of A.  KL >= 0.
068: *
069: *  KU     - INTEGER
070: *           The number of superdiagonals within the band of A.  KU >= 0.
071: *
072: *  ALPHA  - REAL
073: *           On entry, ALPHA specifies the scalar alpha.
074: *           Unchanged on exit.
075: *
076: *  A      - REAL             array of DIMENSION ( LDA, n )
077: *           Before entry, the leading m by n part of the array A must
078: *           contain the matrix of coefficients.
079: *           Unchanged on exit.
080: *
081: *  LDA    - INTEGER
082: *           On entry, LDA specifies the first dimension of A as declared
083: *           in the calling (sub) program. LDA must be at least
084: *           max( 1, m ).
085: *           Unchanged on exit.
086: *
087: *  X      - REAL             array of DIMENSION at least
088: *           ( 1 + ( n - 1 )*abs( INCX ) ) when TRANS = 'N' or 'n'
089: *           and at least
090: *           ( 1 + ( m - 1 )*abs( INCX ) ) otherwise.
091: *           Before entry, the incremented array X must contain the
092: *           vector x.
093: *           Unchanged on exit.
094: *
095: *  INCX   - INTEGER
096: *           On entry, INCX specifies the increment for the elements of
097: *           X. INCX must not be zero.
098: *           Unchanged on exit.
099: *
100: *  BETA   - REAL
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      - REAL             array of DIMENSION at least
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   - 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: *     .. Parameters ..
123:       COMPLEX            ONE, ZERO
124:       PARAMETER          ( ONE = 1.0E+0, ZERO = 0.0E+0 )
125: *     ..
126: *     .. Local Scalars ..
127:       LOGICAL            SYMB_ZERO
128:       REAL               TEMP, SAFE1
129:       INTEGER            I, INFO, IY, J, JX, KX, KY, LENX, LENY, KD
130:       COMPLEX            CDUM
131: *     ..
132: *     .. External Subroutines ..
133:       EXTERNAL           XERBLA, SLAMCH
134:       REAL               SLAMCH
135: *     ..
136: *     .. External Functions ..
137:       EXTERNAL           ILATRANS
138:       INTEGER            ILATRANS
139: *     ..
140: *     .. Intrinsic Functions ..
141:       INTRINSIC          MAX, ABS, REAL, AIMAG, SIGN
142: *     ..
143: *     .. Statement Functions
144:       REAL               CABS1
145: *     ..
146: *     .. Statement Function Definitions ..
147:       CABS1( CDUM ) = ABS( REAL( CDUM ) ) + ABS( AIMAG( CDUM ) )
148: *     ..
149: *     .. Executable Statements ..
150: *
151: *     Test the input parameters.
152: *
153:       INFO = 0
154:       IF     ( .NOT.( ( TRANS.EQ.ILATRANS( 'N' ) )
155:      $           .OR. ( TRANS.EQ.ILATRANS( 'T' ) )
156:      $           .OR. ( TRANS.EQ.ILATRANS( 'C' ) ) ) ) THEN
157:          INFO = 1
158:       ELSE IF( M.LT.0 )THEN
159:          INFO = 2
160:       ELSE IF( N.LT.0 )THEN
161:          INFO = 3
162:       ELSE IF( KL.LT.0 ) THEN
163:          INFO = 4
164:       ELSE IF( KU.LT.0 ) THEN
165:          INFO = 5
166:       ELSE IF( LDAB.LT.KL+KU+1 )THEN
167:          INFO = 6
168:       ELSE IF( INCX.EQ.0 )THEN
169:          INFO = 8
170:       ELSE IF( INCY.EQ.0 )THEN
171:          INFO = 11
172:       END IF
173:       IF( INFO.NE.0 )THEN
174:          CALL XERBLA( 'CLA_GBAMV ', INFO )
175:          RETURN
176:       END IF
177: *
178: *     Quick return if possible.
179: *
180:       IF( ( M.EQ.0 ).OR.( N.EQ.0 ).OR.
181:      $    ( ( ALPHA.EQ.ZERO ).AND.( BETA.EQ.ONE ) ) )
182:      $   RETURN
183: *
184: *     Set  LENX  and  LENY, the lengths of the vectors x and y, and set
185: *     up the start points in  X  and  Y.
186: *
187:       IF( TRANS.EQ.ILATRANS( 'N' ) )THEN
188:          LENX = N
189:          LENY = M
190:       ELSE
191:          LENX = M
192:          LENY = N
193:       END IF
194:       IF( INCX.GT.0 )THEN
195:          KX = 1
196:       ELSE
197:          KX = 1 - ( LENX - 1 )*INCX
198:       END IF
199:       IF( INCY.GT.0 )THEN
200:          KY = 1
201:       ELSE
202:          KY = 1 - ( LENY - 1 )*INCY
203:       END IF
204: *
205: *     Set SAFE1 essentially to be the underflow threshold times the
206: *     number of additions in each row.
207: *
208:       SAFE1 = SLAMCH( 'Safe minimum' )
209:       SAFE1 = (N+1)*SAFE1
210: *
211: *     Form  y := alpha*abs(A)*abs(x) + beta*abs(y).
212: *
213: *     The O(M*N) SYMB_ZERO tests could be replaced by O(N) queries to
214: *     the inexact flag.  Still doesn't help change the iteration order
215: *     to per-column.
216: *
217:       KD = KU + 1
218:       IY = KY
219:       IF ( INCX.EQ.1 ) THEN
220:          DO I = 1, LENY
221:             IF ( BETA .EQ. 0.0 ) THEN
222:                SYMB_ZERO = .TRUE.
223:                Y( IY ) = 0.0
224:             ELSE IF ( Y( IY ) .EQ. 0.0 ) THEN
225:                SYMB_ZERO = .TRUE.
226:             ELSE
227:                SYMB_ZERO = .FALSE.
228:                Y( IY ) = BETA * ABS( Y( IY ) )
229:             END IF
230:             IF ( ALPHA .NE. 0.0 ) THEN
231:                DO J = MAX( I-KU, 1 ), MIN( I+KL, LENX )
232:                   IF( TRANS.EQ.ILATRANS( 'N' ) )THEN
233:                      TEMP = CABS1( AB( KD+I-J, J ) )
234:                   ELSE
235:                      TEMP = CABS1( AB( J, KD+I-J ) )
236:                   END IF
237: 
238:                   SYMB_ZERO = SYMB_ZERO .AND.
239:      $                 ( X( J ) .EQ. ZERO .OR. TEMP .EQ. ZERO )
240: 
241:                   Y( IY ) = Y( IY ) + ALPHA*CABS1( X( J ) )*TEMP
242:                END DO
243:             END IF
244: 
245:             IF ( .NOT.SYMB_ZERO)
246:      $           Y( IY ) = Y( IY ) + SIGN( SAFE1, Y( IY ) )
247: 
248:             IY = IY + INCY
249:          END DO
250:       ELSE
251:          DO I = 1, LENY
252:             IF ( BETA .EQ. 0.0 ) THEN
253:                SYMB_ZERO = .TRUE.
254:                Y( IY ) = 0.0
255:             ELSE IF ( Y( IY ) .EQ. 0.0 ) THEN
256:                SYMB_ZERO = .TRUE.
257:             ELSE
258:                SYMB_ZERO = .FALSE.
259:                Y( IY ) = BETA * ABS( Y( IY ) )
260:             END IF
261:             IF ( ALPHA .NE. 0.0 ) THEN
262:                JX = KX
263:                DO J = MAX( I-KU, 1 ), MIN( I+KL, LENX )
264: 
265:                   IF( TRANS.EQ.ILATRANS( 'N' ) )THEN
266:                      TEMP = CABS1( AB( KD+I-J, J ) )
267:                   ELSE
268:                      TEMP = CABS1( AB( J, KD+I-J ) )
269:                   END IF
270: 
271:                   SYMB_ZERO = SYMB_ZERO .AND.
272:      $                 ( X( JX ) .EQ. ZERO .OR. TEMP .EQ. ZERO )
273: 
274:                   Y( IY ) = Y( IY ) + ALPHA*CABS1( X( JX ) )*TEMP
275:                   JX = JX + INCX
276:                END DO
277:             END IF
278: 
279:             IF ( .NOT.SYMB_ZERO )
280:      $           Y( IY ) = Y( IY ) + SIGN( SAFE1, Y( IY ) )
281: 
282:             IY = IY + INCY
283:          END DO
284:       END IF
285: *
286:       RETURN
287: *
288: *     End of CLA_GBAMV
289: *
290:       END
291: