LAPACK 3.3.0

cgbequ.f

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