LAPACK: SRC/zla_gbamv.f Source File

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