LAPACK 3.3.0

dgbequ.f

Go to the documentation of this file.
00001       SUBROUTINE DGBEQU( M, N, KL, KU, AB, LDAB, R, C, ROWCND, COLCND,
00002      $                   AMAX, INFO )
00003 *
00004 *  -- LAPACK routine (version 3.2) --
00005 *  -- LAPACK is a software package provided by Univ. of Tennessee,    --
00006 *  -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
00007 *     November 2006
00008 *
00009 *     .. Scalar Arguments ..
00010       INTEGER            INFO, KL, KU, LDAB, M, N
00011       DOUBLE PRECISION   AMAX, COLCND, ROWCND
00012 *     ..
00013 *     .. Array Arguments ..
00014       DOUBLE PRECISION   AB( LDAB, * ), C( * ), R( * )
00015 *     ..
00016 *
00017 *  Purpose
00018 *  =======
00019 *
00020 *  DGBEQU computes row and column scalings intended to equilibrate an
00021 *  M-by-N band matrix A and reduce its condition number.  R returns the
00022 *  row scale factors and C the column scale factors, chosen to try to
00023 *  make the largest element in each row and column of the matrix B with
00024 *  elements B(i,j)=R(i)*A(i,j)*C(j) have absolute value 1.
00025 *
00026 *  R(i) and C(j) are restricted to be between SMLNUM = smallest safe
00027 *  number and BIGNUM = largest safe number.  Use of these scaling
00028 *  factors is not guaranteed to reduce the condition number of A but
00029 *  works well in practice.
00030 *
00031 *  Arguments
00032 *  =========
00033 *
00034 *  M       (input) INTEGER
00035 *          The number of rows of the matrix A.  M >= 0.
00036 *
00037 *  N       (input) INTEGER
00038 *          The number of columns of the matrix A.  N >= 0.
00039 *
00040 *  KL      (input) INTEGER
00041 *          The number of subdiagonals within the band of A.  KL >= 0.
00042 *
00043 *  KU      (input) INTEGER
00044 *          The number of superdiagonals within the band of A.  KU >= 0.
00045 *
00046 *  AB      (input) DOUBLE PRECISION array, dimension (LDAB,N)
00047 *          The band matrix A, stored in rows 1 to KL+KU+1.  The j-th
00048 *          column of A is stored in the j-th column of the array AB as
00049 *          follows:
00050 *          AB(ku+1+i-j,j) = A(i,j) for max(1,j-ku)<=i<=min(m,j+kl).
00051 *
00052 *  LDAB    (input) INTEGER
00053 *          The leading dimension of the array AB.  LDAB >= KL+KU+1.
00054 *
00055 *  R       (output) DOUBLE PRECISION array, dimension (M)
00056 *          If INFO = 0, or INFO > M, R contains the row scale factors
00057 *          for A.
00058 *
00059 *  C       (output) DOUBLE PRECISION array, dimension (N)
00060 *          If INFO = 0, C contains the column scale factors for A.
00061 *
00062 *  ROWCND  (output) DOUBLE PRECISION
00063 *          If INFO = 0 or INFO > M, ROWCND contains the ratio of the
00064 *          smallest R(i) to the largest R(i).  If ROWCND >= 0.1 and
00065 *          AMAX is neither too large nor too small, it is not worth
00066 *          scaling by R.
00067 *
00068 *  COLCND  (output) DOUBLE PRECISION
00069 *          If INFO = 0, COLCND contains the ratio of the smallest
00070 *          C(i) to the largest C(i).  If COLCND >= 0.1, it is not
00071 *          worth scaling by C.
00072 *
00073 *  AMAX    (output) DOUBLE PRECISION
00074 *          Absolute value of largest matrix element.  If AMAX is very
00075 *          close to overflow or very close to underflow, the matrix
00076 *          should be scaled.
00077 *
00078 *  INFO    (output) INTEGER
00079 *          = 0:  successful exit
00080 *          < 0:  if INFO = -i, the i-th argument had an illegal value
00081 *          > 0:  if INFO = i, and i is
00082 *                <= M:  the i-th row of A is exactly zero
00083 *                >  M:  the (i-M)-th column of A is exactly zero
00084 *
00085 *  =====================================================================
00086 *
00087 *     .. Parameters ..
00088       DOUBLE PRECISION   ONE, ZERO
00089       PARAMETER          ( ONE = 1.0D+0, ZERO = 0.0D+0 )
00090 *     ..
00091 *     .. Local Scalars ..
00092       INTEGER            I, J, KD
00093       DOUBLE PRECISION   BIGNUM, RCMAX, RCMIN, SMLNUM
00094 *     ..
00095 *     .. External Functions ..
00096       DOUBLE PRECISION   DLAMCH
00097       EXTERNAL           DLAMCH
00098 *     ..
00099 *     .. External Subroutines ..
00100       EXTERNAL           XERBLA
00101 *     ..
00102 *     .. Intrinsic Functions ..
00103       INTRINSIC          ABS, MAX, MIN
00104 *     ..
00105 *     .. Executable Statements ..
00106 *
00107 *     Test the input parameters
00108 *
00109       INFO = 0
00110       IF( M.LT.0 ) THEN
00111          INFO = -1
00112       ELSE IF( N.LT.0 ) THEN
00113          INFO = -2
00114       ELSE IF( KL.LT.0 ) THEN
00115          INFO = -3
00116       ELSE IF( KU.LT.0 ) THEN
00117          INFO = -4
00118       ELSE IF( LDAB.LT.KL+KU+1 ) THEN
00119          INFO = -6
00120       END IF
00121       IF( INFO.NE.0 ) THEN
00122          CALL XERBLA( 'DGBEQU', -INFO )
00123          RETURN
00124       END IF
00125 *
00126 *     Quick return if possible
00127 *
00128       IF( M.EQ.0 .OR. N.EQ.0 ) THEN
00129          ROWCND = ONE
00130          COLCND = ONE
00131          AMAX = ZERO
00132          RETURN
00133       END IF
00134 *
00135 *     Get machine constants.
00136 *
00137       SMLNUM = DLAMCH( 'S' )
00138       BIGNUM = ONE / SMLNUM
00139 *
00140 *     Compute row scale factors.
00141 *
00142       DO 10 I = 1, M
00143          R( I ) = ZERO
00144    10 CONTINUE
00145 *
00146 *     Find the maximum element in each row.
00147 *
00148       KD = KU + 1
00149       DO 30 J = 1, N
00150          DO 20 I = MAX( J-KU, 1 ), MIN( J+KL, M )
00151             R( I ) = MAX( R( I ), ABS( AB( KD+I-J, J ) ) )
00152    20    CONTINUE
00153    30 CONTINUE
00154 *
00155 *     Find the maximum and minimum scale factors.
00156 *
00157       RCMIN = BIGNUM
00158       RCMAX = ZERO
00159       DO 40 I = 1, M
00160          RCMAX = MAX( RCMAX, R( I ) )
00161          RCMIN = MIN( RCMIN, R( I ) )
00162    40 CONTINUE
00163       AMAX = RCMAX
00164 *
00165       IF( RCMIN.EQ.ZERO ) THEN
00166 *
00167 *        Find the first zero scale factor and return an error code.
00168 *
00169          DO 50 I = 1, M
00170             IF( R( I ).EQ.ZERO ) THEN
00171                INFO = I
00172                RETURN
00173             END IF
00174    50    CONTINUE
00175       ELSE
00176 *
00177 *        Invert the scale factors.
00178 *
00179          DO 60 I = 1, M
00180             R( I ) = ONE / MIN( MAX( R( I ), SMLNUM ), BIGNUM )
00181    60    CONTINUE
00182 *
00183 *        Compute ROWCND = min(R(I)) / max(R(I))
00184 *
00185          ROWCND = MAX( RCMIN, SMLNUM ) / MIN( RCMAX, BIGNUM )
00186       END IF
00187 *
00188 *     Compute column scale factors
00189 *
00190       DO 70 J = 1, N
00191          C( J ) = ZERO
00192    70 CONTINUE
00193 *
00194 *     Find the maximum element in each column,
00195 *     assuming the row scaling computed above.
00196 *
00197       KD = KU + 1
00198       DO 90 J = 1, N
00199          DO 80 I = MAX( J-KU, 1 ), MIN( J+KL, M )
00200             C( J ) = MAX( C( J ), ABS( AB( KD+I-J, J ) )*R( I ) )
00201    80    CONTINUE
00202    90 CONTINUE
00203 *
00204 *     Find the maximum and minimum scale factors.
00205 *
00206       RCMIN = BIGNUM
00207       RCMAX = ZERO
00208       DO 100 J = 1, N
00209          RCMIN = MIN( RCMIN, C( J ) )
00210          RCMAX = MAX( RCMAX, C( J ) )
00211   100 CONTINUE
00212 *
00213       IF( RCMIN.EQ.ZERO ) THEN
00214 *
00215 *        Find the first zero scale factor and return an error code.
00216 *
00217          DO 110 J = 1, N
00218             IF( C( J ).EQ.ZERO ) THEN
00219                INFO = M + J
00220                RETURN
00221             END IF
00222   110    CONTINUE
00223       ELSE
00224 *
00225 *        Invert the scale factors.
00226 *
00227          DO 120 J = 1, N
00228             C( J ) = ONE / MIN( MAX( C( J ), SMLNUM ), BIGNUM )
00229   120    CONTINUE
00230 *
00231 *        Compute COLCND = min(C(J)) / max(C(J))
00232 *
00233          COLCND = MAX( RCMIN, SMLNUM ) / MIN( RCMAX, BIGNUM )
00234       END IF
00235 *
00236       RETURN
00237 *
00238 *     End of DGBEQU
00239 *
00240       END
 All Files Functions