LAPACK  3.4.2
LAPACK: Linear Algebra PACKage
 All Files Functions Groups
sla_syamv.f
Go to the documentation of this file.
1 *> \brief \b SLA_SYAMV computes a matrix-vector product using a symmetric indefinite matrix 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 SLA_SYAMV + dependencies
10 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/sla_syamv.f">
11 *> [TGZ]</a>
12 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/sla_syamv.f">
13 *> [ZIP]</a>
14 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/sla_syamv.f">
15 *> [TXT]</a>
16 *> \endhtmlonly
17 *
18 * Definition:
19 * ===========
20 *
21 * SUBROUTINE SLA_SYAMV( UPLO, N, ALPHA, A, LDA, X, INCX, BETA, Y,
22 * INCY )
23 *
24 * .. Scalar Arguments ..
25 * REAL ALPHA, BETA
26 * INTEGER INCX, INCY, LDA, N, UPLO
27 * ..
28 * .. Array Arguments ..
29 * REAL A( LDA, * ), X( * ), Y( * )
30 * ..
31 *
32 *
33 *> \par Purpose:
34 * =============
35 *>
36 *> \verbatim
37 *>
38 *> SLA_SYAMV performs the matrix-vector operation
39 *>
40 *> y := alpha*abs(A)*abs(x) + beta*abs(y),
41 *>
42 *> where alpha and beta are scalars, x and y are vectors and A is an
43 *> n by n symmetric matrix.
44 *>
45 *> This function is primarily used in calculating error bounds.
46 *> To protect against underflow during evaluation, components in
47 *> the resulting vector are perturbed away from zero by (N+1)
48 *> times the underflow threshold. To prevent unnecessarily large
49 *> errors for block-structure embedded in general matrices,
50 *> "symbolically" zero components are not perturbed. A zero
51 *> entry is considered "symbolic" if all multiplications involved
52 *> in computing that entry have at least one zero multiplicand.
53 *> \endverbatim
54 *
55 * Arguments:
56 * ==========
57 *
58 *> \param[in] UPLO
59 *> \verbatim
60 *> UPLO is INTEGER
61 *> On entry, UPLO specifies whether the upper or lower
62 *> triangular part of the array A is to be referenced as
63 *> follows:
64 *>
65 *> UPLO = BLAS_UPPER Only the upper triangular part of A
66 *> is to be referenced.
67 *>
68 *> UPLO = BLAS_LOWER Only the lower triangular part of A
69 *> is to be referenced.
70 *>
71 *> Unchanged on exit.
72 *> \endverbatim
73 *>
74 *> \param[in] N
75 *> \verbatim
76 *> N is INTEGER
77 *> On entry, N specifies the number of columns of the matrix A.
78 *> N must be at least zero.
79 *> Unchanged on exit.
80 *> \endverbatim
81 *>
82 *> \param[in] ALPHA
83 *> \verbatim
84 *> ALPHA is REAL .
85 *> On entry, ALPHA specifies the scalar alpha.
86 *> Unchanged on exit.
87 *> \endverbatim
88 *>
89 *> \param[in] A
90 *> \verbatim
91 *> A is REAL array of DIMENSION ( LDA, n ).
92 *> Before entry, the leading m by n part of the array A must
93 *> contain the matrix of coefficients.
94 *> Unchanged on exit.
95 *> \endverbatim
96 *>
97 *> \param[in] LDA
98 *> \verbatim
99 *> LDA is INTEGER
100 *> On entry, LDA specifies the first dimension of A as declared
101 *> in the calling (sub) program. LDA must be at least
102 *> max( 1, n ).
103 *> Unchanged on exit.
104 *> \endverbatim
105 *>
106 *> \param[in] X
107 *> \verbatim
108 *> X is REAL array, dimension
109 *> ( 1 + ( n - 1 )*abs( INCX ) )
110 *> Before entry, the incremented array X must contain the
111 *> vector x.
112 *> Unchanged on exit.
113 *> \endverbatim
114 *>
115 *> \param[in] INCX
116 *> \verbatim
117 *> INCX is INTEGER
118 *> On entry, INCX specifies the increment for the elements of
119 *> X. INCX must not be zero.
120 *> Unchanged on exit.
121 *> \endverbatim
122 *>
123 *> \param[in] BETA
124 *> \verbatim
125 *> BETA is REAL .
126 *> On entry, BETA specifies the scalar beta. When BETA is
127 *> supplied as zero then Y need not be set on input.
128 *> Unchanged on exit.
129 *> \endverbatim
130 *>
131 *> \param[in,out] Y
132 *> \verbatim
133 *> Y is REAL array, dimension
134 *> ( 1 + ( n - 1 )*abs( INCY ) )
135 *> Before entry with BETA non-zero, the incremented array Y
136 *> must contain the vector y. On exit, Y is overwritten by the
137 *> updated vector y.
138 *> \endverbatim
139 *>
140 *> \param[in] INCY
141 *> \verbatim
142 *> INCY is INTEGER
143 *> On entry, INCY specifies the increment for the elements of
144 *> Y. INCY must not be zero.
145 *> Unchanged on exit.
146 *> \endverbatim
147 *
148 * Authors:
149 * ========
150 *
151 *> \author Univ. of Tennessee
152 *> \author Univ. of California Berkeley
153 *> \author Univ. of Colorado Denver
154 *> \author NAG Ltd.
155 *
156 *> \date September 2012
157 *
158 *> \ingroup realSYcomputational
159 *
160 *> \par Further Details:
161 * =====================
162 *>
163 *> \verbatim
164 *>
165 *> Level 2 Blas routine.
166 *>
167 *> -- Written on 22-October-1986.
168 *> Jack Dongarra, Argonne National Lab.
169 *> Jeremy Du Croz, Nag Central Office.
170 *> Sven Hammarling, Nag Central Office.
171 *> Richard Hanson, Sandia National Labs.
172 *> -- Modified for the absolute-value product, April 2006
173 *> Jason Riedy, UC Berkeley
174 *> \endverbatim
175 *>
176 * =====================================================================
177  SUBROUTINE sla_syamv( UPLO, N, ALPHA, A, LDA, X, INCX, BETA, Y,
178  $ incy )
179 *
180 * -- LAPACK computational routine (version 3.4.2) --
181 * -- LAPACK is a software package provided by Univ. of Tennessee, --
182 * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
183 * September 2012
184 *
185 * .. Scalar Arguments ..
186  REAL alpha, beta
187  INTEGER incx, incy, lda, n, uplo
188 * ..
189 * .. Array Arguments ..
190  REAL a( lda, * ), x( * ), y( * )
191 * ..
192 *
193 * =====================================================================
194 *
195 * .. Parameters ..
196  REAL one, zero
197  parameter( one = 1.0e+0, zero = 0.0e+0 )
198 * ..
199 * .. Local Scalars ..
200  LOGICAL symb_zero
201  REAL temp, safe1
202  INTEGER i, info, iy, j, jx, kx, ky
203 * ..
204 * .. External Subroutines ..
205  EXTERNAL xerbla, slamch
206  REAL slamch
207 * ..
208 * .. External Functions ..
209  EXTERNAL ilauplo
210  INTEGER ilauplo
211 * ..
212 * .. Intrinsic Functions ..
213  INTRINSIC max, abs, sign
214 * ..
215 * .. Executable Statements ..
216 *
217 * Test the input parameters.
218 *
219  info = 0
220  IF ( uplo.NE.ilauplo( 'U' ) .AND.
221  $ uplo.NE.ilauplo( 'L' ) ) THEN
222  info = 1
223  ELSE IF( n.LT.0 )THEN
224  info = 2
225  ELSE IF( lda.LT.max( 1, n ) )THEN
226  info = 5
227  ELSE IF( incx.EQ.0 )THEN
228  info = 7
229  ELSE IF( incy.EQ.0 )THEN
230  info = 10
231  END IF
232  IF( info.NE.0 )THEN
233  CALL xerbla( 'SSYMV ', info )
234  return
235  END IF
236 *
237 * Quick return if possible.
238 *
239  IF( ( n.EQ.0 ).OR.( ( alpha.EQ.zero ).AND.( beta.EQ.one ) ) )
240  $ return
241 *
242 * Set up the start points in X and Y.
243 *
244  IF( incx.GT.0 )THEN
245  kx = 1
246  ELSE
247  kx = 1 - ( n - 1 )*incx
248  END IF
249  IF( incy.GT.0 )THEN
250  ky = 1
251  ELSE
252  ky = 1 - ( n - 1 )*incy
253  END IF
254 *
255 * Set SAFE1 essentially to be the underflow threshold times the
256 * number of additions in each row.
257 *
258  safe1 = slamch( 'Safe minimum' )
259  safe1 = (n+1)*safe1
260 *
261 * Form y := alpha*abs(A)*abs(x) + beta*abs(y).
262 *
263 * The O(N^2) SYMB_ZERO tests could be replaced by O(N) queries to
264 * the inexact flag. Still doesn't help change the iteration order
265 * to per-column.
266 *
267  iy = ky
268  IF ( incx.EQ.1 ) THEN
269  IF ( uplo .EQ. ilauplo( 'U' ) ) THEN
270  DO i = 1, n
271  IF ( beta .EQ. zero ) THEN
272  symb_zero = .true.
273  y( iy ) = 0.0
274  ELSE IF ( y( iy ) .EQ. zero ) THEN
275  symb_zero = .true.
276  ELSE
277  symb_zero = .false.
278  y( iy ) = beta * abs( y( iy ) )
279  END IF
280  IF ( alpha .NE. zero ) THEN
281  DO j = 1, i
282  temp = abs( a( j, i ) )
283  symb_zero = symb_zero .AND.
284  $ ( x( j ) .EQ. zero .OR. temp .EQ. zero )
285 
286  y( iy ) = y( iy ) + alpha*abs( x( j ) )*temp
287  END DO
288  DO j = i+1, n
289  temp = abs( a( i, j ) )
290  symb_zero = symb_zero .AND.
291  $ ( x( j ) .EQ. zero .OR. temp .EQ. zero )
292 
293  y( iy ) = y( iy ) + alpha*abs( x( j ) )*temp
294  END DO
295  END IF
296 
297  IF ( .NOT.symb_zero )
298  $ y( iy ) = y( iy ) + sign( safe1, y( iy ) )
299 
300  iy = iy + incy
301  END DO
302  ELSE
303  DO i = 1, n
304  IF ( beta .EQ. zero ) THEN
305  symb_zero = .true.
306  y( iy ) = 0.0
307  ELSE IF ( y( iy ) .EQ. zero ) THEN
308  symb_zero = .true.
309  ELSE
310  symb_zero = .false.
311  y( iy ) = beta * abs( y( iy ) )
312  END IF
313  IF ( alpha .NE. zero ) THEN
314  DO j = 1, i
315  temp = abs( a( i, j ) )
316  symb_zero = symb_zero .AND.
317  $ ( x( j ) .EQ. zero .OR. temp .EQ. zero )
318 
319  y( iy ) = y( iy ) + alpha*abs( x( j ) )*temp
320  END DO
321  DO j = i+1, n
322  temp = abs( a( j, i ) )
323  symb_zero = symb_zero .AND.
324  $ ( x( j ) .EQ. zero .OR. temp .EQ. zero )
325 
326  y( iy ) = y( iy ) + alpha*abs( x( j ) )*temp
327  END DO
328  END IF
329 
330  IF ( .NOT.symb_zero )
331  $ y( iy ) = y( iy ) + sign( safe1, y( iy ) )
332 
333  iy = iy + incy
334  END DO
335  END IF
336  ELSE
337  IF ( uplo .EQ. ilauplo( 'U' ) ) THEN
338  DO i = 1, n
339  IF ( beta .EQ. zero ) THEN
340  symb_zero = .true.
341  y( iy ) = 0.0
342  ELSE IF ( y( iy ) .EQ. zero ) THEN
343  symb_zero = .true.
344  ELSE
345  symb_zero = .false.
346  y( iy ) = beta * abs( y( iy ) )
347  END IF
348  jx = kx
349  IF ( alpha .NE. zero ) THEN
350  DO j = 1, i
351  temp = abs( a( j, i ) )
352  symb_zero = symb_zero .AND.
353  $ ( x( j ) .EQ. zero .OR. temp .EQ. zero )
354 
355  y( iy ) = y( iy ) + alpha*abs( x( jx ) )*temp
356  jx = jx + incx
357  END DO
358  DO j = i+1, n
359  temp = abs( a( i, j ) )
360  symb_zero = symb_zero .AND.
361  $ ( x( j ) .EQ. zero .OR. temp .EQ. zero )
362 
363  y( iy ) = y( iy ) + alpha*abs( x( jx ) )*temp
364  jx = jx + incx
365  END DO
366  END IF
367 
368  IF ( .NOT.symb_zero )
369  $ y( iy ) = y( iy ) + sign( safe1, y( iy ) )
370 
371  iy = iy + incy
372  END DO
373  ELSE
374  DO i = 1, n
375  IF ( beta .EQ. zero ) THEN
376  symb_zero = .true.
377  y( iy ) = 0.0
378  ELSE IF ( y( iy ) .EQ. zero ) THEN
379  symb_zero = .true.
380  ELSE
381  symb_zero = .false.
382  y( iy ) = beta * abs( y( iy ) )
383  END IF
384  jx = kx
385  IF ( alpha .NE. zero ) THEN
386  DO j = 1, i
387  temp = abs( a( i, j ) )
388  symb_zero = symb_zero .AND.
389  $ ( x( j ) .EQ. zero .OR. temp .EQ. zero )
390 
391  y( iy ) = y( iy ) + alpha*abs( x( jx ) )*temp
392  jx = jx + incx
393  END DO
394  DO j = i+1, n
395  temp = abs( a( j, i ) )
396  symb_zero = symb_zero .AND.
397  $ ( x( j ) .EQ. zero .OR. temp .EQ. zero )
398 
399  y( iy ) = y( iy ) + alpha*abs( x( jx ) )*temp
400  jx = jx + incx
401  END DO
402  END IF
403 
404  IF ( .NOT.symb_zero )
405  $ y( iy ) = y( iy ) + sign( safe1, y( iy ) )
406 
407  iy = iy + incy
408  END DO
409  END IF
410 
411  END IF
412 *
413  return
414 *
415 * End of SLA_SYAMV
416 *
417  END