1       SUBROUTINE DLA_GEAMV ( TRANS, M, N, ALPHA, A, LDA, X, INCX, BETA,
  2      $                       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, LDA, M, N, TRANS
 17 *     ..
 18 *     .. Array Arguments ..
 19       DOUBLE PRECISION   A( LDA, * ), X( * ), Y( * )
 20 *     ..
 21 *
 22 *  Purpose
 23 *  =======
 24 *
 25 *  DLA_GEAMV  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 *  ALPHA    (input) DOUBLE PRECISION
 66 *           On entry, ALPHA specifies the scalar alpha.
 67 *           Unchanged on exit.
 68 *
 69 *  A        (input) DOUBLE PRECISION array of DIMENSION ( LDA, n )
 70 *           Before entry, the leading m by n part of the array A must
 71 *           contain the matrix of coefficients.
 72 *           Unchanged on exit.
 73 *
 74 *  LDA      (input) INTEGER
 75 *           On entry, LDA specifies the first dimension of A as declared
 76 *           in the calling (sub) program. LDA must be at least
 77 *           max( 1, m ).
 78 *           Unchanged on exit.
 79 *
 80 *  X        (input) DOUBLE PRECISION array, dimension
 81 *           ( 1 + ( n - 1 )*abs( INCX ) ) when TRANS = 'N' or 'n'
 82 *           and at least
 83 *           ( 1 + ( m - 1 )*abs( INCX ) ) otherwise.
 84 *           Before entry, the incremented array X must contain the
 85 *           vector x.
 86 *           Unchanged on exit.
 87 *
 88 *  INCX     (input) INTEGER
 89 *           On entry, INCX specifies the increment for the elements of
 90 *           X. INCX must not be zero.
 91 *           Unchanged on exit.
 92 *
 93 *  BETA     (input) DOUBLE PRECISION
 94 *           On entry, BETA specifies the scalar beta. When BETA is
 95 *           supplied as zero then Y need not be set on input.
 96 *           Unchanged on exit.
 97 *
 98 *  Y        (input/output) DOUBLE PRECISION
 99 *           Array of DIMENSION at least
100 *           ( 1 + ( m - 1 )*abs( INCY ) ) when TRANS = 'N' or 'n'
101 *           and at least
102 *           ( 1 + ( n - 1 )*abs( INCY ) ) otherwise.
103 *           Before entry with BETA non-zero, the incremented array Y
104 *           must contain the vector y. On exit, Y is overwritten by the
105 *           updated vector y.
106 *
107 *  INCY     (input) INTEGER
108 *           On entry, INCY specifies the increment for the elements of
109 *           Y. INCY must not be zero.
110 *           Unchanged on exit.
111 *
112 *  Level 2 Blas routine.
113 *
114 *  =====================================================================
115 *
116 *     .. Parameters ..
117       DOUBLE PRECISION   ONE, ZERO
118       PARAMETER          ( ONE = 1.0D+0, ZERO = 0.0D+0 )
119 *     ..
120 *     .. Local Scalars ..
121       LOGICAL            SYMB_ZERO
122       DOUBLE PRECISION   TEMP, SAFE1
123       INTEGER            I, INFO, IY, J, JX, KX, KY, LENX, LENY
124 *     ..
125 *     .. External Subroutines ..
126       EXTERNAL           XERBLA, DLAMCH
127       DOUBLE PRECISION   DLAMCH
128 *     ..
129 *     .. External Functions ..
130       EXTERNAL           ILATRANS
131       INTEGER            ILATRANS
132 *     ..
133 *     .. Intrinsic Functions ..
134       INTRINSIC          MAXABSSIGN
135 *     ..
136 *     .. Executable Statements ..
137 *
138 *     Test the input parameters.
139 *
140       INFO = 0
141       IF     ( .NOT.( ( TRANS.EQ.ILATRANS( 'N' ) )
142      $           .OR. ( TRANS.EQ.ILATRANS( 'T' ) )
143      $           .OR. ( TRANS.EQ.ILATRANS( 'C' )) ) ) THEN
144          INFO = 1
145       ELSE IF( M.LT.0 )THEN
146          INFO = 2
147       ELSE IF( N.LT.0 )THEN
148          INFO = 3
149       ELSE IF( LDA.LT.MAX1, M ) )THEN
150          INFO = 6
151       ELSE IF( INCX.EQ.0 )THEN
152          INFO = 8
153       ELSE IF( INCY.EQ.0 )THEN
154          INFO = 11
155       END IF
156       IF( INFO.NE.0 )THEN
157          CALL XERBLA( 'DLA_GEAMV ', INFO )
158          RETURN
159       END IF
160 *
161 *     Quick return if possible.
162 *
163       IF( ( M.EQ.0 ).OR.( N.EQ.0 ).OR.
164      $    ( ( ALPHA.EQ.ZERO ).AND.( BETA.EQ.ONE ) ) )
165      $   RETURN
166 *
167 *     Set  LENX  and  LENY, the lengths of the vectors x and y, and set
168 *     up the start points in  X  and  Y.
169 *
170       IF( TRANS.EQ.ILATRANS( 'N' ) )THEN
171          LENX = N
172          LENY = M
173       ELSE
174          LENX = M
175          LENY = N
176       END IF
177       IF( INCX.GT.0 )THEN
178          KX = 1
179       ELSE
180          KX = 1 - ( LENX - 1 )*INCX
181       END IF
182       IF( INCY.GT.0 )THEN
183          KY = 1
184       ELSE
185          KY = 1 - ( LENY - 1 )*INCY
186       END IF
187 *
188 *     Set SAFE1 essentially to be the underflow threshold times the
189 *     number of additions in each row.
190 *
191       SAFE1 = DLAMCH( 'Safe minimum' )
192       SAFE1 = (N+1)*SAFE1
193 *
194 *     Form  y := alpha*abs(A)*abs(x) + beta*abs(y).
195 *
196 *     The O(M*N) SYMB_ZERO tests could be replaced by O(N) queries to
197 *     the inexact flag.  Still doesn't help change the iteration order
198 *     to per-column.
199 *
200       IY = KY
201       IF ( INCX.EQ.1 ) THEN
202          IF( TRANS.EQ.ILATRANS( 'N' ) )THEN
203             DO I = 1, LENY
204                IF ( BETA .EQ. ZERO ) THEN
205                   SYMB_ZERO = .TRUE.
206                   Y( IY ) = 0.0D+0
207                ELSE IF ( Y( IY ) .EQ. ZERO ) THEN
208                   SYMB_ZERO = .TRUE.
209                ELSE
210                   SYMB_ZERO = .FALSE.
211                   Y( IY ) = BETA * ABS( Y( IY ) )
212                END IF
213                IF ( ALPHA .NE. ZERO ) THEN
214                   DO J = 1, LENX
215                      TEMP = ABS( A( I, J ) )
216                      SYMB_ZERO = SYMB_ZERO .AND.
217      $                    ( X( J ) .EQ. ZERO .OR. TEMP .EQ. ZERO )
218 
219                      Y( IY ) = Y( IY ) + ALPHA*ABS( X( J ) )*TEMP
220                   END DO
221                END IF
222 
223                IF ( .NOT.SYMB_ZERO )
224      $              Y( IY ) = Y( IY ) + SIGN( SAFE1, Y( IY ) )
225 
226                IY = IY + INCY
227             END DO
228          ELSE
229             DO I = 1, LENY
230                IF ( BETA .EQ. ZERO ) THEN
231                   SYMB_ZERO = .TRUE.
232                   Y( IY ) = 0.0D+0
233                ELSE IF ( Y( IY ) .EQ. ZERO ) THEN
234                   SYMB_ZERO = .TRUE.
235                ELSE
236                   SYMB_ZERO = .FALSE.
237                   Y( IY ) = BETA * ABS( Y( IY ) )
238                END IF
239                IF ( ALPHA .NE. ZERO ) THEN
240                   DO J = 1, LENX
241                      TEMP = ABS( A( J, I ) )
242                      SYMB_ZERO = SYMB_ZERO .AND.
243      $                    ( X( J ) .EQ. ZERO .OR. TEMP .EQ. ZERO )
244 
245                      Y( IY ) = Y( IY ) + ALPHA*ABS( X( J ) )*TEMP
246                   END DO
247                END IF
248 
249                IF ( .NOT.SYMB_ZERO )
250      $              Y( IY ) = Y( IY ) + SIGN( SAFE1, Y( IY ) )
251 
252                IY = IY + INCY
253             END DO
254          END IF
255       ELSE
256          IF( TRANS.EQ.ILATRANS( 'N' ) )THEN
257             DO I = 1, LENY
258                IF ( BETA .EQ. ZERO ) THEN
259                   SYMB_ZERO = .TRUE.
260                   Y( IY ) = 0.0D+0
261                ELSE IF ( Y( IY ) .EQ. ZERO ) THEN
262                   SYMB_ZERO = .TRUE.
263                ELSE
264                   SYMB_ZERO = .FALSE.
265                   Y( IY ) = BETA * ABS( Y( IY ) )
266                END IF
267                IF ( ALPHA .NE. ZERO ) THEN
268                   JX = KX
269                   DO J = 1, LENX
270                      TEMP = ABS( A( I, J ) )
271                      SYMB_ZERO = SYMB_ZERO .AND.
272      $                    ( X( JX ) .EQ. ZERO .OR. TEMP .EQ. ZERO )
273 
274                      Y( IY ) = Y( IY ) + ALPHA*ABS( 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          ELSE
285             DO I = 1, LENY
286                IF ( BETA .EQ. ZERO ) THEN
287                   SYMB_ZERO = .TRUE.
288                   Y( IY ) = 0.0D+0
289                ELSE IF ( Y( IY ) .EQ. ZERO ) THEN
290                   SYMB_ZERO = .TRUE.
291                ELSE
292                   SYMB_ZERO = .FALSE.
293                   Y( IY ) = BETA * ABS( Y( IY ) )
294                END IF
295                IF ( ALPHA .NE. ZERO ) THEN
296                   JX = KX
297                   DO J = 1, LENX
298                      TEMP = ABS( A( J, I ) )
299                      SYMB_ZERO = SYMB_ZERO .AND.
300      $                    ( X( JX ) .EQ. ZERO .OR. TEMP .EQ. ZERO )
301 
302                      Y( IY ) = Y( IY ) + ALPHA*ABS( X( JX ) )*TEMP
303                      JX = JX + INCX
304                   END DO
305                END IF
306 
307                IF (.NOT.SYMB_ZERO)
308      $              Y( IY ) = Y( IY ) + SIGN( SAFE1, Y( IY ) )
309 
310                IY = IY + INCY
311             END DO
312          END IF
313 
314       END IF
315 *
316       RETURN
317 *
318 *     End of DLA_GEAMV
319 *
320       END