LAPACK  3.4.1
LAPACK: Linear Algebra PACKage
csyequb.f
Go to the documentation of this file.
00001 *> \brief \b CSYEQUB
00002 *
00003 *  =========== DOCUMENTATION ===========
00004 *
00005 * Online html documentation available at 
00006 *            http://www.netlib.org/lapack/explore-html/ 
00007 *
00008 *> \htmlonly
00009 *> Download CSYEQUB + dependencies 
00010 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/csyequb.f"> 
00011 *> [TGZ]</a> 
00012 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/csyequb.f"> 
00013 *> [ZIP]</a> 
00014 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/csyequb.f"> 
00015 *> [TXT]</a>
00016 *> \endhtmlonly 
00017 *
00018 *  Definition:
00019 *  ===========
00020 *
00021 *       SUBROUTINE CSYEQUB( UPLO, N, A, LDA, S, SCOND, AMAX, WORK, INFO )
00022 * 
00023 *       .. Scalar Arguments ..
00024 *       INTEGER            INFO, LDA, N
00025 *       REAL               AMAX, SCOND
00026 *       CHARACTER          UPLO
00027 *       ..
00028 *       .. Array Arguments ..
00029 *       COMPLEX            A( LDA, * ), WORK( * )
00030 *       REAL               S( * )
00031 *       ..
00032 *  
00033 *
00034 *> \par Purpose:
00035 *  =============
00036 *>
00037 *> \verbatim
00038 *>
00039 *> CSYEQUB computes row and column scalings intended to equilibrate a
00040 *> symmetric matrix A and reduce its condition number
00041 *> (with respect to the two-norm).  S contains the scale factors,
00042 *> S(i) = 1/sqrt(A(i,i)), chosen so that the scaled matrix B with
00043 *> elements B(i,j) = S(i)*A(i,j)*S(j) has ones on the diagonal.  This
00044 *> choice of S puts the condition number of B within a factor N of the
00045 *> smallest possible condition number over all possible diagonal
00046 *> scalings.
00047 *> \endverbatim
00048 *
00049 *  Arguments:
00050 *  ==========
00051 *
00052 *> \param[in] UPLO
00053 *> \verbatim
00054 *>          UPLO is CHARACTER*1
00055 *>          Specifies whether the details of the factorization are stored
00056 *>          as an upper or lower triangular matrix.
00057 *>          = 'U':  Upper triangular, form is A = U*D*U**T;
00058 *>          = 'L':  Lower triangular, form is A = L*D*L**T.
00059 *> \endverbatim
00060 *>
00061 *> \param[in] N
00062 *> \verbatim
00063 *>          N is INTEGER
00064 *>          The order of the matrix A.  N >= 0.
00065 *> \endverbatim
00066 *>
00067 *> \param[in] A
00068 *> \verbatim
00069 *>          A is COMPLEX array, dimension (LDA,N)
00070 *>          The N-by-N symmetric matrix whose scaling
00071 *>          factors are to be computed.  Only the diagonal elements of A
00072 *>          are referenced.
00073 *> \endverbatim
00074 *>
00075 *> \param[in] LDA
00076 *> \verbatim
00077 *>          LDA is INTEGER
00078 *>          The leading dimension of the array A.  LDA >= max(1,N).
00079 *> \endverbatim
00080 *>
00081 *> \param[out] S
00082 *> \verbatim
00083 *>          S is REAL array, dimension (N)
00084 *>          If INFO = 0, S contains the scale factors for A.
00085 *> \endverbatim
00086 *>
00087 *> \param[out] SCOND
00088 *> \verbatim
00089 *>          SCOND is REAL
00090 *>          If INFO = 0, S contains the ratio of the smallest S(i) to
00091 *>          the largest S(i).  If SCOND >= 0.1 and AMAX is neither too
00092 *>          large nor too small, it is not worth scaling by S.
00093 *> \endverbatim
00094 *>
00095 *> \param[out] AMAX
00096 *> \verbatim
00097 *>          AMAX is REAL
00098 *>          Absolute value of largest matrix element.  If AMAX is very
00099 *>          close to overflow or very close to underflow, the matrix
00100 *>          should be scaled.
00101 *> \endverbatim
00102 *>
00103 *> \param[out] WORK
00104 *> \verbatim
00105 *>          WORK is COMPLEX array, dimension (3*N)
00106 *> \endverbatim
00107 *>
00108 *> \param[out] INFO
00109 *> \verbatim
00110 *>          INFO is INTEGER
00111 *>          = 0:  successful exit
00112 *>          < 0:  if INFO = -i, the i-th argument had an illegal value
00113 *>          > 0:  if INFO = i, the i-th diagonal element is nonpositive.
00114 *> \endverbatim
00115 *
00116 *  Authors:
00117 *  ========
00118 *
00119 *> \author Univ. of Tennessee 
00120 *> \author Univ. of California Berkeley 
00121 *> \author Univ. of Colorado Denver 
00122 *> \author NAG Ltd. 
00123 *
00124 *> \date November 2011
00125 *
00126 *> \ingroup complexSYcomputational
00127 *
00128 *> \par References:
00129 *  ================
00130 *>
00131 *>  Livne, O.E. and Golub, G.H., "Scaling by Binormalization", \n
00132 *>  Numerical Algorithms, vol. 35, no. 1, pp. 97-120, January 2004. \n
00133 *>  DOI 10.1023/B:NUMA.0000016606.32820.69 \n
00134 *>  Tech report version: http://ruready.utah.edu/archive/papers/bin.pdf
00135 *>
00136 *  =====================================================================
00137       SUBROUTINE CSYEQUB( UPLO, N, A, LDA, S, SCOND, AMAX, WORK, INFO )
00138 *
00139 *  -- LAPACK computational routine (version 3.4.0) --
00140 *  -- LAPACK is a software package provided by Univ. of Tennessee,    --
00141 *  -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
00142 *     November 2011
00143 *
00144 *     .. Scalar Arguments ..
00145       INTEGER            INFO, LDA, N
00146       REAL               AMAX, SCOND
00147       CHARACTER          UPLO
00148 *     ..
00149 *     .. Array Arguments ..
00150       COMPLEX            A( LDA, * ), WORK( * )
00151       REAL               S( * )
00152 *     ..
00153 *
00154 *  =====================================================================
00155 *
00156 *     .. Parameters ..
00157       REAL               ONE, ZERO
00158       PARAMETER          ( ONE = 1.0E0, ZERO = 0.0E0 )
00159       INTEGER            MAX_ITER
00160       PARAMETER          ( MAX_ITER = 100 )
00161 *     ..
00162 *     .. Local Scalars ..
00163       INTEGER            I, J, ITER
00164       REAL               AVG, STD, TOL, C0, C1, C2, T, U, SI, D, BASE,
00165      $                   SMIN, SMAX, SMLNUM, BIGNUM, SCALE, SUMSQ
00166       LOGICAL            UP
00167       COMPLEX            ZDUM
00168 *     ..
00169 *     .. External Functions ..
00170       REAL               SLAMCH
00171       LOGICAL            LSAME
00172       EXTERNAL           LSAME, SLAMCH
00173 *     ..
00174 *     .. External Subroutines ..
00175       EXTERNAL           CLASSQ
00176 *     ..
00177 *     .. Intrinsic Functions ..
00178       INTRINSIC          ABS, AIMAG, INT, LOG, MAX, MIN, REAL, SQRT
00179 *     ..
00180 *     .. Statement Functions ..
00181       REAL               CABS1
00182 *     ..
00183 *     Statement Function Definitions
00184       CABS1( ZDUM ) = ABS( REAL( ZDUM ) ) + ABS( AIMAG( ZDUM ) )
00185 *     ..
00186 *     .. Executable Statements ..
00187 *
00188 *     Test the input parameters.
00189 *
00190       INFO = 0
00191       IF ( .NOT. ( LSAME( UPLO, 'U' ) .OR. LSAME( UPLO, 'L' ) ) ) THEN
00192         INFO = -1
00193       ELSE IF ( N .LT. 0 ) THEN
00194         INFO = -2
00195       ELSE IF ( LDA .LT. MAX( 1, N ) ) THEN
00196         INFO = -4
00197       END IF
00198       IF ( INFO .NE. 0 ) THEN
00199         CALL XERBLA( 'CSYEQUB', -INFO )
00200         RETURN
00201       END IF
00202 
00203       UP = LSAME( UPLO, 'U' )
00204       AMAX = ZERO
00205 *
00206 *     Quick return if possible.
00207 *
00208       IF ( N .EQ. 0 ) THEN
00209         SCOND = ONE
00210         RETURN
00211       END IF
00212 
00213       DO I = 1, N
00214         S( I ) = ZERO
00215       END DO
00216 
00217       AMAX = ZERO
00218       IF ( UP ) THEN
00219          DO J = 1, N
00220             DO I = 1, J-1
00221                S( I ) = MAX( S( I ), CABS1( A( I, J ) ) )
00222                S( J ) = MAX( S( J ), CABS1( A( I, J ) ) )
00223                AMAX = MAX( AMAX, CABS1( A( I, J ) ) )
00224             END DO
00225             S( J ) = MAX( S( J ), CABS1( A( J, J) ) )
00226             AMAX = MAX( AMAX, CABS1( A( J, J ) ) )
00227          END DO
00228       ELSE
00229          DO J = 1, N
00230             S( J ) = MAX( S( J ), CABS1( A( J, J ) ) )
00231             AMAX = MAX( AMAX, CABS1( A( J, J ) ) )
00232             DO I = J+1, N
00233                S( I ) = MAX( S( I ), CABS1( A( I, J ) ) )
00234                S( J ) = MAX( S( J ), CABS1 (A( I, J ) ) )
00235                AMAX = MAX( AMAX, CABS1( A( I, J ) ) )
00236             END DO
00237          END DO
00238       END IF
00239       DO J = 1, N
00240          S( J ) = 1.0 / S( J )
00241       END DO
00242 
00243       TOL = ONE / SQRT( 2.0E0 * N )
00244 
00245       DO ITER = 1, MAX_ITER
00246          SCALE = 0.0
00247          SUMSQ = 0.0
00248 *       beta = |A|s
00249         DO I = 1, N
00250            WORK( I ) = ZERO
00251         END DO
00252         IF ( UP ) THEN
00253            DO J = 1, N
00254               DO I = 1, J-1
00255                  T = CABS1( A( I, J ) )
00256                  WORK( I ) = WORK( I ) + CABS1( A( I, J ) ) * S( J )
00257                  WORK( J ) = WORK( J ) + CABS1( A( I, J ) ) * S( I )
00258               END DO
00259               WORK( J ) = WORK( J ) + CABS1( A( J, J ) ) * S( J )
00260            END DO
00261         ELSE
00262            DO J = 1, N
00263               WORK( J ) = WORK( J ) + CABS1( A( J, J ) ) * S( J )
00264               DO I = J+1, N
00265                  T = CABS1( A( I, J ) )
00266                  WORK( I ) = WORK( I ) + CABS1( A( I, J ) ) * S( J )
00267                  WORK( J ) = WORK( J ) + CABS1( A( I, J ) ) * S( I )
00268               END DO
00269            END DO
00270         END IF
00271 
00272 *       avg = s^T beta / n
00273         AVG = 0.0
00274         DO I = 1, N
00275           AVG = AVG + S( I )*WORK( I )
00276         END DO
00277         AVG = AVG / N
00278 
00279         STD = 0.0
00280         DO I = N+1, 2*N
00281            WORK( I ) = S( I-N ) * WORK( I-N ) - AVG
00282         END DO
00283         CALL CLASSQ( N, WORK( N+1 ), 1, SCALE, SUMSQ )
00284         STD = SCALE * SQRT( SUMSQ / N )
00285 
00286         IF ( STD .LT. TOL * AVG ) GOTO 999
00287 
00288         DO I = 1, N
00289           T = CABS1( A( I, I ) )
00290           SI = S( I )
00291           C2 = ( N-1 ) * T
00292           C1 = ( N-2 ) * ( WORK( I ) - T*SI )
00293           C0 = -(T*SI)*SI + 2*WORK( I )*SI - N*AVG
00294           D = C1*C1 - 4*C0*C2
00295 
00296           IF ( D .LE. 0 ) THEN
00297             INFO = -1
00298             RETURN
00299           END IF
00300           SI = -2*C0 / ( C1 + SQRT( D ) )
00301 
00302           D = SI - S( I )
00303           U = ZERO
00304           IF ( UP ) THEN
00305             DO J = 1, I
00306               T = CABS1( A( J, I ) )
00307               U = U + S( J )*T
00308               WORK( J ) = WORK( J ) + D*T
00309             END DO
00310             DO J = I+1,N
00311               T = CABS1( A( I, J ) )
00312               U = U + S( J )*T
00313               WORK( J ) = WORK( J ) + D*T
00314             END DO
00315           ELSE
00316             DO J = 1, I
00317               T = CABS1( A( I, J ) )
00318               U = U + S( J )*T
00319               WORK( J ) = WORK( J ) + D*T
00320             END DO
00321             DO J = I+1,N
00322               T = CABS1( A( J, I ) )
00323               U = U + S( J )*T
00324               WORK( J ) = WORK( J ) + D*T
00325             END DO
00326           END IF
00327           AVG = AVG + ( U + WORK( I ) ) * D / N
00328           S( I ) = SI
00329         END DO
00330       END DO
00331 
00332  999  CONTINUE
00333 
00334       SMLNUM = SLAMCH( 'SAFEMIN' )
00335       BIGNUM = ONE / SMLNUM
00336       SMIN = BIGNUM
00337       SMAX = ZERO
00338       T = ONE / SQRT( AVG )
00339       BASE = SLAMCH( 'B' )
00340       U = ONE / LOG( BASE )
00341       DO I = 1, N
00342         S( I ) = BASE ** INT( U * LOG( S( I ) * T ) )
00343         SMIN = MIN( SMIN, S( I ) )
00344         SMAX = MAX( SMAX, S( I ) )
00345       END DO
00346       SCOND = MAX( SMIN, SMLNUM ) / MIN( SMAX, BIGNUM )
00347 *
00348       END
 All Files Functions