LAPACK: SRC/csyequb.f Source File

Go to the documentation of this file.
00001       SUBROUTINE CSYEQUB( UPLO, N, A, LDA, S, SCOND, AMAX, WORK, INFO )
00002 *
00003 *     -- LAPACK routine (version 3.2.2)                                 --
00004 *     -- Contributed by James Demmel, Deaglan Halligan, Yozo Hida and --
00005 *     -- Jason Riedy of Univ. of California Berkeley.                 --
00006 *     -- June 2010                                                    --
00007 *
00008 *     -- LAPACK is a software package provided by Univ. of Tennessee, --
00009 *     -- Univ. of California Berkeley and NAG Ltd.                    --
00010 *
00011       IMPLICIT NONE
00012 *     ..
00013 *     .. Scalar Arguments ..
00014       INTEGER            INFO, LDA, N
00015       REAL               AMAX, SCOND
00016       CHARACTER          UPLO
00017 *     ..
00018 *     .. Array Arguments ..
00019       COMPLEX            A( LDA, * ), WORK( * )
00020       REAL               S( * )
00021 *     ..
00022 *
00023 *  Purpose
00024 *  =======
00025 *
00026 *  CSYEQUB computes row and column scalings intended to equilibrate a
00027 *  symmetric matrix A and reduce its condition number
00028 *  (with respect to the two-norm).  S contains the scale factors,
00029 *  S(i) = 1/sqrt(A(i,i)), chosen so that the scaled matrix B with
00030 *  elements B(i,j) = S(i)*A(i,j)*S(j) has ones on the diagonal.  This
00031 *  choice of S puts the condition number of B within a factor N of the
00032 *  smallest possible condition number over all possible diagonal
00033 *  scalings.
00034 *
00035 *  Arguments
00036 *  =========
00037 *
00038 *  UPLO    (input) CHARACTER*1
00039 *          Specifies whether the details of the factorization are stored
00040 *          as an upper or lower triangular matrix.
00041 *          = 'U':  Upper triangular, form is A = U*D*U**T;
00042 *          = 'L':  Lower triangular, form is A = L*D*L**T.
00043 *
00044 *  N       (input) INTEGER
00045 *          The order of the matrix A.  N >= 0.
00046 *
00047 *  A       (input) COMPLEX array, dimension (LDA,N)
00048 *          The N-by-N symmetric matrix whose scaling
00049 *          factors are to be computed.  Only the diagonal elements of A
00050 *          are referenced.
00051 *
00052 *  LDA     (input) INTEGER
00053 *          The leading dimension of the array A.  LDA >= max(1,N).
00054 *
00055 *  S       (output) REAL array, dimension (N)
00056 *          If INFO = 0, S contains the scale factors for A.
00057 *
00058 *  SCOND   (output) REAL
00059 *          If INFO = 0, S contains the ratio of the smallest S(i) to
00060 *          the largest S(i).  If SCOND >= 0.1 and AMAX is neither too
00061 *          large nor too small, it is not worth scaling by S.
00062 *
00063 *  AMAX    (output) REAL
00064 *          Absolute value of largest matrix element.  If AMAX is very
00065 *          close to overflow or very close to underflow, the matrix
00066 *          should be scaled.
00067 *
00068 *  WORK    (workspace) COMPLEX array, dimension (3*N)
00069 *
00070 *  INFO    (output) INTEGER
00071 *          = 0:  successful exit
00072 *          < 0:  if INFO = -i, the i-th argument had an illegal value
00073 *          > 0:  if INFO = i, the i-th diagonal element is nonpositive.
00074 *
00075 *  Further Details
00076 *  ======= =======
00077 *
00078 *  Reference: Livne, O.E. and Golub, G.H., "Scaling by Binormalization",
00079 *  Numerical Algorithms, vol. 35, no. 1, pp. 97-120, January 2004.
00080 *  DOI 10.1023/B:NUMA.0000016606.32820.69
00081 *  Tech report version: http://ruready.utah.edu/archive/papers/bin.pdf
00082 *
00083 *  =====================================================================
00084 *
00085 *     .. Parameters ..
00086       REAL               ONE, ZERO
00087       PARAMETER          ( ONE = 1.0E0, ZERO = 0.0E0 )
00088       INTEGER            MAX_ITER
00089       PARAMETER          ( MAX_ITER = 100 )
00090 *     ..
00091 *     .. Local Scalars ..
00092       INTEGER            I, J, ITER
00093       REAL               AVG, STD, TOL, C0, C1, C2, T, U, SI, D, BASE,
00094      $                   SMIN, SMAX, SMLNUM, BIGNUM, SCALE, SUMSQ
00095       LOGICAL            UP
00096       COMPLEX            ZDUM
00097 *     ..
00098 *     .. External Functions ..
00099       REAL               SLAMCH
00100       LOGICAL            LSAME
00101       EXTERNAL           LSAME, SLAMCH
00102 *     ..
00103 *     .. External Subroutines ..
00104       EXTERNAL           CLASSQ
00105 *     ..
00106 *     .. Intrinsic Functions ..
00107       INTRINSIC          ABS, AIMAG, INT, LOG, MAX, MIN, REAL, SQRT
00108 *     ..
00109 *     .. Statement Functions ..
00110       REAL               CABS1
00111 *     ..
00112 *     Statement Function Definitions
00113       CABS1( ZDUM ) = ABS( REAL( ZDUM ) ) + ABS( AIMAG( ZDUM ) )
00114 *     ..
00115 *     .. Executable Statements ..
00116 *
00117 *     Test the input parameters.
00118 *
00119       INFO = 0
00120       IF ( .NOT. ( LSAME( UPLO, 'U' ) .OR. LSAME( UPLO, 'L' ) ) ) THEN
00121         INFO = -1
00122       ELSE IF ( N .LT. 0 ) THEN
00123         INFO = -2
00124       ELSE IF ( LDA .LT. MAX( 1, N ) ) THEN
00125         INFO = -4
00126       END IF
00127       IF ( INFO .NE. 0 ) THEN
00128         CALL XERBLA( 'CSYEQUB', -INFO )
00129         RETURN
00130       END IF
00131 
00132       UP = LSAME( UPLO, 'U' )
00133       AMAX = ZERO
00134 *
00135 *     Quick return if possible.
00136 *
00137       IF ( N .EQ. 0 ) THEN
00138         SCOND = ONE
00139         RETURN
00140       END IF
00141 
00142       DO I = 1, N
00143         S( I ) = ZERO
00144       END DO
00145 
00146       AMAX = ZERO
00147       IF ( UP ) THEN
00148          DO J = 1, N
00149             DO I = 1, J-1
00150                S( I ) = MAX( S( I ), CABS1( A( I, J ) ) )
00151                S( J ) = MAX( S( J ), CABS1( A( I, J ) ) )
00152                AMAX = MAX( AMAX, CABS1( A( I, J ) ) )
00153             END DO
00154             S( J ) = MAX( S( J ), CABS1( A( J, J) ) )
00155             AMAX = MAX( AMAX, CABS1( A( J, J ) ) )
00156          END DO
00157       ELSE
00158          DO J = 1, N
00159             S( J ) = MAX( S( J ), CABS1( A( J, J ) ) )
00160             AMAX = MAX( AMAX, CABS1( A( J, J ) ) )
00161             DO I = J+1, N
00162                S( I ) = MAX( S( I ), CABS1( A( I, J ) ) )
00163                S( J ) = MAX( S( J ), CABS1 (A( I, J ) ) )
00164                AMAX = MAX( AMAX, CABS1( A( I, J ) ) )
00165             END DO
00166          END DO
00167       END IF
00168       DO J = 1, N
00169          S( J ) = 1.0 / S( J )
00170       END DO
00171 
00172       TOL = ONE / SQRT( 2.0E0 * N )
00173 
00174       DO ITER = 1, MAX_ITER
00175          SCALE = 0.0
00176          SUMSQ = 0.0
00177 *       beta = |A|s
00178         DO I = 1, N
00179            WORK( I ) = ZERO
00180         END DO
00181         IF ( UP ) THEN
00182            DO J = 1, N
00183               DO I = 1, J-1
00184                  T = CABS1( A( I, J ) )
00185                  WORK( I ) = WORK( I ) + CABS1( A( I, J ) ) * S( J )
00186                  WORK( J ) = WORK( J ) + CABS1( A( I, J ) ) * S( I )
00187               END DO
00188               WORK( J ) = WORK( J ) + CABS1( A( J, J ) ) * S( J )
00189            END DO
00190         ELSE
00191            DO J = 1, N
00192               WORK( J ) = WORK( J ) + CABS1( A( J, J ) ) * S( J )
00193               DO I = J+1, N
00194                  T = CABS1( A( I, J ) )
00195                  WORK( I ) = WORK( I ) + CABS1( A( I, J ) ) * S( J )
00196                  WORK( J ) = WORK( J ) + CABS1( A( I, J ) ) * S( I )
00197               END DO
00198            END DO
00199         END IF
00200 
00201 *       avg = s^T beta / n
00202         AVG = 0.0
00203         DO I = 1, N
00204           AVG = AVG + S( I )*WORK( I )
00205         END DO
00206         AVG = AVG / N
00207 
00208         STD = 0.0
00209         DO I = N+1, 2*N
00210            WORK( I ) = S( I-N ) * WORK( I-N ) - AVG
00211         END DO
00212         CALL CLASSQ( N, WORK( N+1 ), 1, SCALE, SUMSQ )
00213         STD = SCALE * SQRT( SUMSQ / N )
00214 
00215         IF ( STD .LT. TOL * AVG ) GOTO 999
00216 
00217         DO I = 1, N
00218           T = CABS1( A( I, I ) )
00219           SI = S( I )
00220           C2 = ( N-1 ) * T
00221           C1 = ( N-2 ) * ( WORK( I ) - T*SI )
00222           C0 = -(T*SI)*SI + 2*WORK( I )*SI - N*AVG
00223           D = C1*C1 - 4*C0*C2
00224 
00225           IF ( D .LE. 0 ) THEN
00226             INFO = -1
00227             RETURN
00228           END IF
00229           SI = -2*C0 / ( C1 + SQRT( D ) )
00230 
00231           D = SI - S( I )
00232           U = ZERO
00233           IF ( UP ) THEN
00234             DO J = 1, I
00235               T = CABS1( A( J, I ) )
00236               U = U + S( J )*T
00237               WORK( J ) = WORK( J ) + D*T
00238             END DO
00239             DO J = I+1,N
00240               T = CABS1( A( I, J ) )
00241               U = U + S( J )*T
00242               WORK( J ) = WORK( J ) + D*T
00243             END DO
00244           ELSE
00245             DO J = 1, I
00246               T = CABS1( A( I, J ) )
00247               U = U + S( J )*T
00248               WORK( J ) = WORK( J ) + D*T
00249             END DO
00250             DO J = I+1,N
00251               T = CABS1( A( J, I ) )
00252               U = U + S( J )*T
00253               WORK( J ) = WORK( J ) + D*T
00254             END DO
00255           END IF
00256           AVG = AVG + ( U + WORK( I ) ) * D / N
00257           S( I ) = SI
00258         END DO
00259       END DO
00260 
00261  999  CONTINUE
00262 
00263       SMLNUM = SLAMCH( 'SAFEMIN' )
00264       BIGNUM = ONE / SMLNUM
00265       SMIN = BIGNUM
00266       SMAX = ZERO
00267       T = ONE / SQRT( AVG )
00268       BASE = SLAMCH( 'B' )
00269       U = ONE / LOG( BASE )
00270       DO I = 1, N
00271         S( I ) = BASE ** INT( U * LOG( S( I ) * T ) )
00272         SMIN = MIN( SMIN, S( I ) )
00273         SMAX = MAX( SMAX, S( I ) )
00274       END DO
00275       SCOND = MAX( SMIN, SMLNUM ) / MIN( SMAX, BIGNUM )
00276 *
00277       END