LAPACK: SRC/sla_syamv.f Source File

Go to the documentation of this file.
00001       SUBROUTINE SLA_SYAMV( UPLO, N, ALPHA, A, LDA, X, INCX, BETA, Y,
00002      $                      INCY )
00003 *
00004 *     -- LAPACK routine (version 3.2.2)                                 --
00005 *     -- Contributed by James Demmel, Deaglan Halligan, Yozo Hida and --
00006 *     -- Jason Riedy of Univ. of California Berkeley.                 --
00007 *     -- June 2010                                                    --
00008 *
00009 *     -- LAPACK is a software package provided by Univ. of Tennessee, --
00010 *     -- Univ. of California Berkeley and NAG Ltd.                    --
00011 *
00012       IMPLICIT NONE
00013 *     ..
00014 *     .. Scalar Arguments ..
00015       REAL               ALPHA, BETA
00016       INTEGER            INCX, INCY, LDA, N, UPLO
00017 *     ..
00018 *     .. Array Arguments ..
00019       REAL               A( LDA, * ), X( * ), Y( * )
00020 *     ..
00021 *
00022 *  Purpose
00023 *  =======
00024 *
00025 *  SLA_SYAMV  performs the matrix-vector operation
00026 *
00027 *          y := alpha*abs(A)*abs(x) + beta*abs(y),
00028 *
00029 *  where alpha and beta are scalars, x and y are vectors and A is an
00030 *  n by n symmetric matrix.
00031 *
00032 *  This function is primarily used in calculating error bounds.
00033 *  To protect against underflow during evaluation, components in
00034 *  the resulting vector are perturbed away from zero by (N+1)
00035 *  times the underflow threshold.  To prevent unnecessarily large
00036 *  errors for block-structure embedded in general matrices,
00037 *  "symbolically" zero components are not perturbed.  A zero
00038 *  entry is considered "symbolic" if all multiplications involved
00039 *  in computing that entry have at least one zero multiplicand.
00040 *
00041 *  Arguments
00042 *  ==========
00043 *
00044 *  UPLO    (input) INTEGER
00045 *           On entry, UPLO specifies whether the upper or lower
00046 *           triangular part of the array A is to be referenced as
00047 *           follows:
00048 *
00049 *              UPLO = BLAS_UPPER   Only the upper triangular part of A
00050 *                                  is to be referenced.
00051 *
00052 *              UPLO = BLAS_LOWER   Only the lower triangular part of A
00053 *                                  is to be referenced.
00054 *
00055 *           Unchanged on exit.
00056 *
00057 *  N       (input) INTEGER
00058 *           On entry, N specifies the number of columns of the matrix A.
00059 *           N must be at least zero.
00060 *           Unchanged on exit.
00061 *
00062 *  ALPHA   (input) REAL            .
00063 *           On entry, ALPHA specifies the scalar alpha.
00064 *           Unchanged on exit.
00065 *
00066 *  A      - REAL             array of DIMENSION ( LDA, n ).
00067 *           Before entry, the leading m by n part of the array A must
00068 *           contain the matrix of coefficients.
00069 *           Unchanged on exit.
00070 *
00071 *  LDA     (input) INTEGER
00072 *           On entry, LDA specifies the first dimension of A as declared
00073 *           in the calling (sub) program. LDA must be at least
00074 *           max( 1, n ).
00075 *           Unchanged on exit.
00076 *
00077 *  X       (input) REAL array, dimension
00078 *           ( 1 + ( n - 1 )*abs( INCX ) )
00079 *           Before entry, the incremented array X must contain the
00080 *           vector x.
00081 *           Unchanged on exit.
00082 *
00083 *  INCX    (input) INTEGER
00084 *           On entry, INCX specifies the increment for the elements of
00085 *           X. INCX must not be zero.
00086 *           Unchanged on exit.
00087 *
00088 *  BETA    (input) REAL            .
00089 *           On entry, BETA specifies the scalar beta. When BETA is
00090 *           supplied as zero then Y need not be set on input.
00091 *           Unchanged on exit.
00092 *
00093 *  Y       (input/output) REAL array, dimension
00094 *           ( 1 + ( n - 1 )*abs( INCY ) )
00095 *           Before entry with BETA non-zero, the incremented array Y
00096 *           must contain the vector y. On exit, Y is overwritten by the
00097 *           updated vector y.
00098 *
00099 *  INCY    (input) INTEGER
00100 *           On entry, INCY specifies the increment for the elements of
00101 *           Y. INCY must not be zero.
00102 *           Unchanged on exit.
00103 *
00104 *  Further Details
00105 *  ===============
00106 *
00107 *  Level 2 Blas routine.
00108 *
00109 *  -- Written on 22-October-1986.
00110 *     Jack Dongarra, Argonne National Lab.
00111 *     Jeremy Du Croz, Nag Central Office.
00112 *     Sven Hammarling, Nag Central Office.
00113 *     Richard Hanson, Sandia National Labs.
00114 *  -- Modified for the absolute-value product, April 2006
00115 *     Jason Riedy, UC Berkeley
00116 *
00117 *  =====================================================================
00118 *
00119 *     .. Parameters ..
00120       REAL               ONE, ZERO
00121       PARAMETER          ( ONE = 1.0E+0, ZERO = 0.0E+0 )
00122 *     ..
00123 *     .. Local Scalars ..
00124       LOGICAL            SYMB_ZERO
00125       REAL               TEMP, SAFE1
00126       INTEGER            I, INFO, IY, J, JX, KX, KY
00127 *     ..
00128 *     .. External Subroutines ..
00129       EXTERNAL           XERBLA, SLAMCH
00130       REAL               SLAMCH
00131 *     ..
00132 *     .. External Functions ..
00133       EXTERNAL           ILAUPLO
00134       INTEGER            ILAUPLO
00135 *     ..
00136 *     .. Intrinsic Functions ..
00137       INTRINSIC          MAX, ABS, SIGN
00138 *     ..
00139 *     .. Executable Statements ..
00140 *
00141 *     Test the input parameters.
00142 *
00143       INFO = 0
00144       IF     ( UPLO.NE.ILAUPLO( 'U' ) .AND.
00145      $         UPLO.NE.ILAUPLO( 'L' ) ) THEN
00146          INFO = 1
00147       ELSE IF( N.LT.0 )THEN
00148          INFO = 2
00149       ELSE IF( LDA.LT.MAX( 1, N ) )THEN
00150          INFO = 5
00151       ELSE IF( INCX.EQ.0 )THEN
00152          INFO = 7
00153       ELSE IF( INCY.EQ.0 )THEN
00154          INFO = 10
00155       END IF
00156       IF( INFO.NE.0 )THEN
00157          CALL XERBLA( 'SSYMV ', INFO )
00158          RETURN
00159       END IF
00160 *
00161 *     Quick return if possible.
00162 *
00163       IF( ( N.EQ.0 ).OR.( ( ALPHA.EQ.ZERO ).AND.( BETA.EQ.ONE ) ) )
00164      $   RETURN
00165 *
00166 *     Set up the start points in  X  and  Y.
00167 *
00168       IF( INCX.GT.0 )THEN
00169          KX = 1
00170       ELSE
00171          KX = 1 - ( N - 1 )*INCX
00172       END IF
00173       IF( INCY.GT.0 )THEN
00174          KY = 1
00175       ELSE
00176          KY = 1 - ( N - 1 )*INCY
00177       END IF
00178 *
00179 *     Set SAFE1 essentially to be the underflow threshold times the
00180 *     number of additions in each row.
00181 *
00182       SAFE1 = SLAMCH( 'Safe minimum' )
00183       SAFE1 = (N+1)*SAFE1
00184 *
00185 *     Form  y := alpha*abs(A)*abs(x) + beta*abs(y).
00186 *
00187 *     The O(N^2) SYMB_ZERO tests could be replaced by O(N) queries to
00188 *     the inexact flag.  Still doesn't help change the iteration order
00189 *     to per-column.
00190 *
00191       IY = KY
00192       IF ( INCX.EQ.1 ) THEN
00193          IF ( UPLO .EQ. ILAUPLO( 'U' ) ) THEN
00194             DO I = 1, N
00195                IF ( BETA .EQ. ZERO ) THEN
00196                   SYMB_ZERO = .TRUE.
00197                   Y( IY ) = 0.0
00198                ELSE IF ( Y( IY ) .EQ. ZERO ) THEN
00199                   SYMB_ZERO = .TRUE.
00200                ELSE
00201                   SYMB_ZERO = .FALSE.
00202                   Y( IY ) = BETA * ABS( Y( IY ) )
00203                END IF
00204                IF ( ALPHA .NE. ZERO ) THEN
00205                   DO J = 1, I
00206                      TEMP = ABS( A( J, I ) )
00207                      SYMB_ZERO = SYMB_ZERO .AND.
00208      $                    ( X( J ) .EQ. ZERO .OR. TEMP .EQ. ZERO )
00209 
00210                      Y( IY ) = Y( IY ) + ALPHA*ABS( X( J ) )*TEMP
00211                   END DO
00212                   DO J = I+1, N
00213                      TEMP = ABS( A( I, J ) )
00214                      SYMB_ZERO = SYMB_ZERO .AND.
00215      $                    ( X( J ) .EQ. ZERO .OR. TEMP .EQ. ZERO )
00216 
00217                      Y( IY ) = Y( IY ) + ALPHA*ABS( X( J ) )*TEMP
00218                   END DO
00219                END IF
00220 
00221                IF ( .NOT.SYMB_ZERO )
00222      $              Y( IY ) = Y( IY ) + SIGN( SAFE1, Y( IY ) )
00223 
00224                IY = IY + INCY
00225             END DO
00226          ELSE
00227             DO I = 1, N
00228                IF ( BETA .EQ. ZERO ) THEN
00229                   SYMB_ZERO = .TRUE.
00230                   Y( IY ) = 0.0
00231                ELSE IF ( Y( IY ) .EQ. ZERO ) THEN
00232                   SYMB_ZERO = .TRUE.
00233                ELSE
00234                   SYMB_ZERO = .FALSE.
00235                   Y( IY ) = BETA * ABS( Y( IY ) )
00236                END IF
00237                IF ( ALPHA .NE. ZERO ) THEN
00238                   DO J = 1, I
00239                      TEMP = ABS( A( I, J ) )
00240                      SYMB_ZERO = SYMB_ZERO .AND.
00241      $                    ( X( J ) .EQ. ZERO .OR. TEMP .EQ. ZERO )
00242 
00243                      Y( IY ) = Y( IY ) + ALPHA*ABS( X( J ) )*TEMP
00244                   END DO
00245                   DO J = I+1, N
00246                      TEMP = ABS( A( J, I ) )
00247                      SYMB_ZERO = SYMB_ZERO .AND.
00248      $                    ( X( J ) .EQ. ZERO .OR. TEMP .EQ. ZERO )
00249 
00250                      Y( IY ) = Y( IY ) + ALPHA*ABS( X( J ) )*TEMP
00251                   END DO
00252                END IF
00253 
00254                IF ( .NOT.SYMB_ZERO )
00255      $              Y( IY ) = Y( IY ) + SIGN( SAFE1, Y( IY ) )
00256 
00257                IY = IY + INCY
00258             END DO
00259          END IF
00260       ELSE
00261          IF ( UPLO .EQ. ILAUPLO( 'U' ) ) THEN
00262             DO I = 1, N
00263                IF ( BETA .EQ. ZERO ) THEN
00264                   SYMB_ZERO = .TRUE.
00265                   Y( IY ) = 0.0
00266                ELSE IF ( Y( IY ) .EQ. ZERO ) THEN
00267                   SYMB_ZERO = .TRUE.
00268                ELSE
00269                   SYMB_ZERO = .FALSE.
00270                   Y( IY ) = BETA * ABS( Y( IY ) )
00271                END IF
00272                JX = KX
00273                IF ( ALPHA .NE. ZERO ) THEN
00274                   DO J = 1, I
00275                      TEMP = ABS( A( J, I ) )
00276                      SYMB_ZERO = SYMB_ZERO .AND.
00277      $                    ( X( J ) .EQ. ZERO .OR. TEMP .EQ. ZERO )
00278 
00279                      Y( IY ) = Y( IY ) + ALPHA*ABS( X( JX ) )*TEMP
00280                      JX = JX + INCX
00281                   END DO
00282                   DO J = I+1, N
00283                      TEMP = ABS( A( I, J ) )
00284                      SYMB_ZERO = SYMB_ZERO .AND.
00285      $                    ( X( J ) .EQ. ZERO .OR. TEMP .EQ. ZERO )
00286 
00287                      Y( IY ) = Y( IY ) + ALPHA*ABS( X( JX ) )*TEMP
00288                      JX = JX + INCX
00289                   END DO
00290                END IF
00291 
00292                IF ( .NOT.SYMB_ZERO )
00293      $              Y( IY ) = Y( IY ) + SIGN( SAFE1, Y( IY ) )
00294 
00295                IY = IY + INCY
00296             END DO
00297          ELSE
00298             DO I = 1, N
00299                IF ( BETA .EQ. ZERO ) THEN
00300                   SYMB_ZERO = .TRUE.
00301                   Y( IY ) = 0.0
00302                ELSE IF ( Y( IY ) .EQ. ZERO ) THEN
00303                   SYMB_ZERO = .TRUE.
00304                ELSE
00305                   SYMB_ZERO = .FALSE.
00306                   Y( IY ) = BETA * ABS( Y( IY ) )
00307                END IF
00308                JX = KX
00309                IF ( ALPHA .NE. ZERO ) THEN
00310                   DO J = 1, I
00311                      TEMP = ABS( A( I, J ) )
00312                      SYMB_ZERO = SYMB_ZERO .AND.
00313      $                    ( X( J ) .EQ. ZERO .OR. TEMP .EQ. ZERO )
00314 
00315                      Y( IY ) = Y( IY ) + ALPHA*ABS( X( JX ) )*TEMP
00316                      JX = JX + INCX
00317                   END DO
00318                   DO J = I+1, N
00319                      TEMP = ABS( A( J, I ) )
00320                      SYMB_ZERO = SYMB_ZERO .AND.
00321      $                    ( X( J ) .EQ. ZERO .OR. TEMP .EQ. ZERO )
00322 
00323                      Y( IY ) = Y( IY ) + ALPHA*ABS( X( JX ) )*TEMP
00324                      JX = JX + INCX
00325                   END DO
00326                END IF
00327 
00328                IF ( .NOT.SYMB_ZERO )
00329      $              Y( IY ) = Y( IY ) + SIGN( SAFE1, Y( IY ) )
00330 
00331                IY = IY + INCY
00332             END DO
00333          END IF
00334 
00335       END IF
00336 *
00337       RETURN
00338 *
00339 *     End of SLA_SYAMV
00340 *
00341       END