LAPACK: SRC/dla_geamv.f Source File

Go to the documentation of this file.
00001       SUBROUTINE DLA_GEAMV ( TRANS, M, N, ALPHA, A, LDA, X, INCX, BETA,
00002      $                       Y, 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, M, N, TRANS
00017 *     ..
00018 *     .. Array Arguments ..
00019       DOUBLE PRECISION   A( LDA, * ), X( * ), Y( * )
00020 *     ..
00021 *
00022 *  Purpose
00023 *  =======
00024 *
00025 *  DLA_GEAMV  performs one of the matrix-vector operations
00026 *
00027 *          y := alpha*abs(A)*abs(x) + beta*abs(y),
00028 *     or   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 *  m by n 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 *  TRANS   (input) INTEGER
00046 *           On entry, TRANS specifies the operation to be performed as
00047 *           follows:
00048 *
00049 *             BLAS_NO_TRANS      y := alpha*abs(A)*abs(x) + beta*abs(y)
00050 *             BLAS_TRANS         y := alpha*abs(A')*abs(x) + beta*abs(y)
00051 *             BLAS_CONJ_TRANS    y := alpha*abs(A')*abs(x) + beta*abs(y)
00052 *
00053 *           Unchanged on exit.
00054 *
00055 *  M       (input) INTEGER
00056 *           On entry, M specifies the number of rows of the matrix A.
00057 *           M must be at least zero.
00058 *           Unchanged on exit.
00059 *
00060 *  N       (input) INTEGER
00061 *           On entry, N specifies the number of columns of the matrix A.
00062 *           N must be at least zero.
00063 *           Unchanged on exit.
00064 *
00065 *  ALPHA  - DOUBLE PRECISION
00066 *           On entry, ALPHA specifies the scalar alpha.
00067 *           Unchanged on exit.
00068 *
00069 *  A      - DOUBLE PRECISION   array of DIMENSION ( LDA, n )
00070 *           Before entry, the leading m by n part of the array A must
00071 *           contain the matrix of coefficients.
00072 *           Unchanged on exit.
00073 *
00074 *  LDA     (input) INTEGER
00075 *           On entry, LDA specifies the first dimension of A as declared
00076 *           in the calling (sub) program. LDA must be at least
00077 *           max( 1, m ).
00078 *           Unchanged on exit.
00079 *
00080 *  X       (input) DOUBLE PRECISION array, dimension
00081 *           ( 1 + ( n - 1 )*abs( INCX ) ) when TRANS = 'N' or 'n'
00082 *           and at least
00083 *           ( 1 + ( m - 1 )*abs( INCX ) ) otherwise.
00084 *           Before entry, the incremented array X must contain the
00085 *           vector x.
00086 *           Unchanged on exit.
00087 *
00088 *  INCX    (input) INTEGER
00089 *           On entry, INCX specifies the increment for the elements of
00090 *           X. INCX must not be zero.
00091 *           Unchanged on exit.
00092 *
00093 *  BETA   - DOUBLE PRECISION
00094 *           On entry, BETA specifies the scalar beta. When BETA is
00095 *           supplied as zero then Y need not be set on input.
00096 *           Unchanged on exit.
00097 *
00098 *  Y      - DOUBLE PRECISION
00099 *           Array of DIMENSION at least
00100 *           ( 1 + ( m - 1 )*abs( INCY ) ) when TRANS = 'N' or 'n'
00101 *           and at least
00102 *           ( 1 + ( n - 1 )*abs( INCY ) ) otherwise.
00103 *           Before entry with BETA non-zero, the incremented array Y
00104 *           must contain the vector y. On exit, Y is overwritten by the
00105 *           updated vector y.
00106 *
00107 *  INCY    (input) INTEGER
00108 *           On entry, INCY specifies the increment for the elements of
00109 *           Y. INCY must not be zero.
00110 *           Unchanged on exit.
00111 *
00112 *  Level 2 Blas routine.
00113 *
00114 *  =====================================================================
00115 *
00116 *     .. Parameters ..
00117       DOUBLE PRECISION   ONE, ZERO
00118       PARAMETER          ( ONE = 1.0D+0, ZERO = 0.0D+0 )
00119 *     ..
00120 *     .. Local Scalars ..
00121       LOGICAL            SYMB_ZERO
00122       DOUBLE PRECISION   TEMP, SAFE1
00123       INTEGER            I, INFO, IY, J, JX, KX, KY, LENX, LENY
00124 *     ..
00125 *     .. External Subroutines ..
00126       EXTERNAL           XERBLA, DLAMCH
00127       DOUBLE PRECISION   DLAMCH
00128 *     ..
00129 *     .. External Functions ..
00130       EXTERNAL           ILATRANS
00131       INTEGER            ILATRANS
00132 *     ..
00133 *     .. Intrinsic Functions ..
00134       INTRINSIC          MAX, ABS, SIGN
00135 *     ..
00136 *     .. Executable Statements ..
00137 *
00138 *     Test the input parameters.
00139 *
00140       INFO = 0
00141       IF     ( .NOT.( ( TRANS.EQ.ILATRANS( 'N' ) )
00142      $           .OR. ( TRANS.EQ.ILATRANS( 'T' ) )
00143      $           .OR. ( TRANS.EQ.ILATRANS( 'C' )) ) ) THEN
00144          INFO = 1
00145       ELSE IF( M.LT.0 )THEN
00146          INFO = 2
00147       ELSE IF( N.LT.0 )THEN
00148          INFO = 3
00149       ELSE IF( LDA.LT.MAX( 1, M ) )THEN
00150          INFO = 6
00151       ELSE IF( INCX.EQ.0 )THEN
00152          INFO = 8
00153       ELSE IF( INCY.EQ.0 )THEN
00154          INFO = 11
00155       END IF
00156       IF( INFO.NE.0 )THEN
00157          CALL XERBLA( 'DLA_GEAMV ', INFO )
00158          RETURN
00159       END IF
00160 *
00161 *     Quick return if possible.
00162 *
00163       IF( ( M.EQ.0 ).OR.( N.EQ.0 ).OR.
00164      $    ( ( ALPHA.EQ.ZERO ).AND.( BETA.EQ.ONE ) ) )
00165      $   RETURN
00166 *
00167 *     Set  LENX  and  LENY, the lengths of the vectors x and y, and set
00168 *     up the start points in  X  and  Y.
00169 *
00170       IF( TRANS.EQ.ILATRANS( 'N' ) )THEN
00171          LENX = N
00172          LENY = M
00173       ELSE
00174          LENX = M
00175          LENY = N
00176       END IF
00177       IF( INCX.GT.0 )THEN
00178          KX = 1
00179       ELSE
00180          KX = 1 - ( LENX - 1 )*INCX
00181       END IF
00182       IF( INCY.GT.0 )THEN
00183          KY = 1
00184       ELSE
00185          KY = 1 - ( LENY - 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 = DLAMCH( 'Safe minimum' )
00192       SAFE1 = (N+1)*SAFE1
00193 *
00194 *     Form  y := alpha*abs(A)*abs(x) + beta*abs(y).
00195 *
00196 *     The O(M*N) 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( TRANS.EQ.ILATRANS( 'N' ) )THEN
00203             DO I = 1, LENY
00204                IF ( BETA .EQ. ZERO ) THEN
00205                   SYMB_ZERO = .TRUE.
00206                   Y( IY ) = 0.0D+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, LENX
00215                      TEMP = ABS( A( I, J ) )
00216                      SYMB_ZERO = SYMB_ZERO .AND.
00217      $                    ( X( J ) .EQ. ZERO .OR. TEMP .EQ. ZERO )
00218 
00219                      Y( IY ) = Y( IY ) + ALPHA*ABS( X( J ) )*TEMP
00220                   END DO
00221                END IF
00222 
00223                IF ( .NOT.SYMB_ZERO )
00224      $              Y( IY ) = Y( IY ) + SIGN( SAFE1, Y( IY ) )
00225 
00226                IY = IY + INCY
00227             END DO
00228          ELSE
00229             DO I = 1, LENY
00230                IF ( BETA .EQ. ZERO ) THEN
00231                   SYMB_ZERO = .TRUE.
00232                   Y( IY ) = 0.0D+0
00233                ELSE IF ( Y( IY ) .EQ. ZERO ) THEN
00234                   SYMB_ZERO = .TRUE.
00235                ELSE
00236                   SYMB_ZERO = .FALSE.
00237                   Y( IY ) = BETA * ABS( Y( IY ) )
00238                END IF
00239                IF ( ALPHA .NE. ZERO ) THEN
00240                   DO J = 1, LENX
00241                      TEMP = ABS( A( J, I ) )
00242                      SYMB_ZERO = SYMB_ZERO .AND.
00243      $                    ( X( J ) .EQ. ZERO .OR. TEMP .EQ. ZERO )
00244 
00245                      Y( IY ) = Y( IY ) + ALPHA*ABS( X( J ) )*TEMP
00246                   END DO
00247                END IF
00248 
00249                IF ( .NOT.SYMB_ZERO )
00250      $              Y( IY ) = Y( IY ) + SIGN( SAFE1, Y( IY ) )
00251 
00252                IY = IY + INCY
00253             END DO
00254          END IF
00255       ELSE
00256          IF( TRANS.EQ.ILATRANS( 'N' ) )THEN
00257             DO I = 1, LENY
00258                IF ( BETA .EQ. ZERO ) THEN
00259                   SYMB_ZERO = .TRUE.
00260                   Y( IY ) = 0.0D+0
00261                ELSE IF ( Y( IY ) .EQ. ZERO ) THEN
00262                   SYMB_ZERO = .TRUE.
00263                ELSE
00264                   SYMB_ZERO = .FALSE.
00265                   Y( IY ) = BETA * ABS( Y( IY ) )
00266                END IF
00267                IF ( ALPHA .NE. ZERO ) THEN
00268                   JX = KX
00269                   DO J = 1, LENX
00270                      TEMP = ABS( A( I, J ) )
00271                      SYMB_ZERO = SYMB_ZERO .AND.
00272      $                    ( X( JX ) .EQ. ZERO .OR. TEMP .EQ. ZERO )
00273 
00274                      Y( IY ) = Y( IY ) + ALPHA*ABS( X( JX ) )*TEMP
00275                      JX = JX + INCX
00276                   END DO
00277                END IF
00278 
00279                IF (.NOT.SYMB_ZERO)
00280      $              Y( IY ) = Y( IY ) + SIGN( SAFE1, Y( IY ) )
00281 
00282                IY = IY + INCY
00283             END DO
00284          ELSE
00285             DO I = 1, LENY
00286                IF ( BETA .EQ. ZERO ) THEN
00287                   SYMB_ZERO = .TRUE.
00288                   Y( IY ) = 0.0D+0
00289                ELSE IF ( Y( IY ) .EQ. ZERO ) THEN
00290                   SYMB_ZERO = .TRUE.
00291                ELSE
00292                   SYMB_ZERO = .FALSE.
00293                   Y( IY ) = BETA * ABS( Y( IY ) )
00294                END IF
00295                IF ( ALPHA .NE. ZERO ) THEN
00296                   JX = KX
00297                   DO J = 1, LENX
00298                      TEMP = ABS( A( J, I ) )
00299                      SYMB_ZERO = SYMB_ZERO .AND.
00300      $                    ( X( JX ) .EQ. ZERO .OR. TEMP .EQ. ZERO )
00301 
00302                      Y( IY ) = Y( IY ) + ALPHA*ABS( X( JX ) )*TEMP
00303                      JX = JX + INCX
00304                   END DO
00305                END IF
00306 
00307                IF (.NOT.SYMB_ZERO)
00308      $              Y( IY ) = Y( IY ) + SIGN( SAFE1, Y( IY ) )
00309 
00310                IY = IY + INCY
00311             END DO
00312          END IF
00313 
00314       END IF
00315 *
00316       RETURN
00317 *
00318 *     End of DLA_GEAMV
00319 *
00320       END