LAPACK: SRC/cla_syamv.f Source File

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