LAPACK: SRC/dla_gbamv.f Source File

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