LAPACK: SRC/zla_heamv.f Source File

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