--- rpl/lapack/lapack/zheequb.f 2016/08/27 15:34:49 1.11 +++ rpl/lapack/lapack/zheequb.f 2017/06/17 10:54:14 1.12 @@ -2,24 +2,24 @@ * * =========== DOCUMENTATION =========== * -* Online html documentation available at -* http://www.netlib.org/lapack/explore-html/ +* Online html documentation available at +* http://www.netlib.org/lapack/explore-html/ * *> \htmlonly -*> Download ZHEEQUB + dependencies -*> -*> [TGZ] -*> -*> [ZIP] -*> +*> Download ZHEEQUB + dependencies +*> +*> [TGZ] +*> +*> [ZIP] +*> *> [TXT] -*> \endhtmlonly +*> \endhtmlonly * * Definition: * =========== * * SUBROUTINE ZHEEQUB( UPLO, N, A, LDA, S, SCOND, AMAX, WORK, INFO ) -* +* * .. Scalar Arguments .. * INTEGER INFO, LDA, N * DOUBLE PRECISION AMAX, SCOND @@ -29,7 +29,7 @@ * COMPLEX*16 A( LDA, * ), WORK( * ) * DOUBLE PRECISION S( * ) * .. -* +* * *> \par Purpose: * ============= @@ -37,12 +37,11 @@ *> \verbatim *> *> ZHEEQUB computes row and column scalings intended to equilibrate a -*> Hermitian matrix A and reduce its condition number -*> (with respect to the two-norm). S contains the scale factors, -*> S(i) = 1/sqrt(A(i,i)), chosen so that the scaled matrix B with -*> elements B(i,j) = S(i)*A(i,j)*S(j) has ones on the diagonal. This -*> choice of S puts the condition number of B within a factor N of the -*> smallest possible condition number over all possible diagonal +*> Hermitian matrix A (with respect to the Euclidean norm) and reduce +*> its condition number. The scale factors S are computed by the BIN +*> algorithm (see references) so that the scaled matrix B with elements +*> B(i,j) = S(i)*A(i,j)*S(j) has a condition number within a factor N of +*> the smallest possible condition number over all possible diagonal *> scalings. *> \endverbatim * @@ -52,28 +51,27 @@ *> \param[in] UPLO *> \verbatim *> UPLO is CHARACTER*1 -*> = 'U': Upper triangles of A and B are stored; -*> = 'L': Lower triangles of A and B are stored. +*> = 'U': Upper triangle of A is stored; +*> = 'L': Lower triangle of A is stored. *> \endverbatim *> *> \param[in] N *> \verbatim *> N is INTEGER -*> The order of the matrix A. N >= 0. +*> The order of the matrix A. N >= 0. *> \endverbatim *> *> \param[in] A *> \verbatim *> A is COMPLEX*16 array, dimension (LDA,N) -*> The N-by-N Hermitian matrix whose scaling -*> factors are to be computed. Only the diagonal elements of A -*> are referenced. +*> The N-by-N Hermitian matrix whose scaling factors are to be +*> computed. *> \endverbatim *> *> \param[in] LDA *> \verbatim *> LDA is INTEGER -*> The leading dimension of the array A. LDA >= max(1,N). +*> The leading dimension of the array A. LDA >= max(1,N). *> \endverbatim *> *> \param[out] S @@ -86,21 +84,21 @@ *> \verbatim *> SCOND is DOUBLE PRECISION *> If INFO = 0, S contains the ratio of the smallest S(i) to -*> the largest S(i). If SCOND >= 0.1 and AMAX is neither too +*> the largest S(i). If SCOND >= 0.1 and AMAX is neither too *> large nor too small, it is not worth scaling by S. *> \endverbatim *> *> \param[out] AMAX *> \verbatim *> AMAX is DOUBLE PRECISION -*> Absolute value of largest matrix element. If AMAX is very -*> close to overflow or very close to underflow, the matrix -*> should be scaled. +*> Largest absolute value of any matrix element. If AMAX is +*> very close to overflow or very close to underflow, the +*> matrix should be scaled. *> \endverbatim *> *> \param[out] WORK *> \verbatim -*> WORK is COMPLEX*16 array, dimension (3*N) +*> WORK is COMPLEX*16 array, dimension (2*N) *> \endverbatim *> *> \param[out] INFO @@ -114,19 +112,27 @@ * Authors: * ======== * -*> \author Univ. of Tennessee -*> \author Univ. of California Berkeley -*> \author Univ. of Colorado Denver -*> \author NAG Ltd. +*> \author Univ. of Tennessee +*> \author Univ. of California Berkeley +*> \author Univ. of Colorado Denver +*> \author NAG Ltd. * *> \date April 2012 * *> \ingroup complex16HEcomputational * +*> \par References: +* ================ +*> +*> Livne, O.E. and Golub, G.H., "Scaling by Binormalization", \n +*> Numerical Algorithms, vol. 35, no. 1, pp. 97-120, January 2004. \n +*> DOI 10.1023/B:NUMA.0000016606.32820.69 \n +*> Tech report version: http://citeseerx.ist.psu.edu/viewdoc/summary?doi=10.1.1.3.1679 +*> * ===================================================================== SUBROUTINE ZHEEQUB( UPLO, N, A, LDA, S, SCOND, AMAX, WORK, INFO ) * -* -- LAPACK computational routine (version 3.4.1) -- +* -- LAPACK computational routine (version 3.7.0) -- * -- LAPACK is a software package provided by Univ. of Tennessee, -- * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..-- * April 2012 @@ -145,14 +151,14 @@ * * .. Parameters .. DOUBLE PRECISION ONE, ZERO - PARAMETER ( ONE = 1.0D+0, ZERO = 0.0D+0 ) + PARAMETER ( ONE = 1.0D0, ZERO = 0.0D0 ) INTEGER MAX_ITER PARAMETER ( MAX_ITER = 100 ) * .. * .. Local Scalars .. INTEGER I, J, ITER - DOUBLE PRECISION AVG, STD, TOL, C0, C1, C2, T, U, SI, D, - $ BASE, SMIN, SMAX, SMLNUM, BIGNUM, SCALE, SUMSQ + DOUBLE PRECISION AVG, STD, TOL, C0, C1, C2, T, U, SI, D, BASE, + $ SMIN, SMAX, SMLNUM, BIGNUM, SCALE, SUMSQ LOGICAL UP COMPLEX*16 ZDUM * .. @@ -172,20 +178,22 @@ * .. * .. Statement Function Definitions .. CABS1( ZDUM ) = ABS( DBLE( ZDUM ) ) + ABS( DIMAG( ZDUM ) ) +* .. +* .. Executable Statements .. * -* Test input parameters. +* Test the input parameters. * INFO = 0 - IF (.NOT. ( LSAME( UPLO, 'U' ) .OR. LSAME( UPLO, 'L' ) ) ) THEN - INFO = -1 + IF ( .NOT. ( LSAME( UPLO, 'U' ) .OR. LSAME( UPLO, 'L' ) ) ) THEN + INFO = -1 ELSE IF ( N .LT. 0 ) THEN - INFO = -2 + INFO = -2 ELSE IF ( LDA .LT. MAX( 1, N ) ) THEN - INFO = -4 + INFO = -4 END IF IF ( INFO .NE. 0 ) THEN - CALL XERBLA( 'ZHEEQUB', -INFO ) - RETURN + CALL XERBLA( 'ZHEEQUB', -INFO ) + RETURN END IF UP = LSAME( UPLO, 'U' ) @@ -194,12 +202,12 @@ * Quick return if possible. * IF ( N .EQ. 0 ) THEN - SCOND = ONE - RETURN + SCOND = ONE + RETURN END IF DO I = 1, N - S( I ) = ZERO + S( I ) = ZERO END DO AMAX = ZERO @@ -220,102 +228,100 @@ DO I = J+1, N S( I ) = MAX( S( I ), CABS1( A( I, J ) ) ) S( J ) = MAX( S( J ), CABS1( A( I, J ) ) ) - AMAX = MAX( AMAX, CABS1( A(I, J ) ) ) + AMAX = MAX( AMAX, CABS1( A( I, J ) ) ) END DO END DO END IF DO J = 1, N - S( J ) = 1.0D+0 / S( J ) + S( J ) = 1.0D0 / S( J ) END DO TOL = ONE / SQRT( 2.0D0 * N ) DO ITER = 1, MAX_ITER - SCALE = 0.0D+0 - SUMSQ = 0.0D+0 -* beta = |A|s - DO I = 1, N - WORK( I ) = ZERO - END DO - IF ( UP ) THEN - DO J = 1, N - DO I = 1, J-1 - T = CABS1( A( I, J ) ) - WORK( I ) = WORK( I ) + CABS1( A( I, J ) ) * S( J ) - WORK( J ) = WORK( J ) + CABS1( A( I, J ) ) * S( I ) - END DO - WORK( J ) = WORK( J ) + CABS1( A( J, J ) ) * S( J ) - END DO - ELSE - DO J = 1, N - WORK( J ) = WORK( J ) + CABS1( A( J, J ) ) * S( J ) - DO I = J+1, N - T = CABS1( A( I, J ) ) - WORK( I ) = WORK( I ) + CABS1( A( I, J ) ) * S( J ) - WORK( J ) = WORK( J ) + CABS1( A( I, J ) ) * S( I ) - END DO - END DO - END IF - -* avg = s^T beta / n - AVG = 0.0D+0 - DO I = 1, N - AVG = AVG + S( I )*WORK( I ) - END DO - AVG = AVG / N - - STD = 0.0D+0 - DO I = 2*N+1, 3*N - WORK( I ) = S( I-2*N ) * WORK( I-2*N ) - AVG - END DO - CALL ZLASSQ( N, WORK( 2*N+1 ), 1, SCALE, SUMSQ ) - STD = SCALE * SQRT( SUMSQ / N ) - - IF ( STD .LT. TOL * AVG ) GOTO 999 - - DO I = 1, N - T = CABS1( A( I, I ) ) - SI = S( I ) - C2 = ( N-1 ) * T - C1 = ( N-2 ) * ( WORK( I ) - T*SI ) - C0 = -(T*SI)*SI + 2*WORK( I )*SI - N*AVG - - D = C1*C1 - 4*C0*C2 - IF ( D .LE. 0 ) THEN - INFO = -1 - RETURN - END IF - SI = -2*C0 / ( C1 + SQRT( D ) ) - - D = SI - S(I) - U = ZERO - IF ( UP ) THEN - DO J = 1, I - T = CABS1( A( J, I ) ) - U = U + S( J )*T - WORK( J ) = WORK( J ) + D*T - END DO - DO J = I+1,N - T = CABS1( A( I, J ) ) - U = U + S( J )*T - WORK( J ) = WORK( J ) + D*T - END DO - ELSE - DO J = 1, I - T = CABS1( A( I, J ) ) - U = U + S( J )*T - WORK( J ) = WORK( J ) + D*T + SCALE = 0.0D0 + SUMSQ = 0.0D0 +* beta = |A|s + DO I = 1, N + WORK( I ) = ZERO + END DO + IF ( UP ) THEN + DO J = 1, N + DO I = 1, J-1 + WORK( I ) = WORK( I ) + CABS1( A( I, J ) ) * S( J ) + WORK( J ) = WORK( J ) + CABS1( A( I, J ) ) * S( I ) + END DO + WORK( J ) = WORK( J ) + CABS1( A( J, J ) ) * S( J ) END DO - DO J = I+1,N - T = CABS1( A( J, I ) ) - U = U + S( J )*T - WORK( J ) = WORK( J ) + D*T + ELSE + DO J = 1, N + WORK( J ) = WORK( J ) + CABS1( A( J, J ) ) * S( J ) + DO I = J+1, N + WORK( I ) = WORK( I ) + CABS1( A( I, J ) ) * S( J ) + WORK( J ) = WORK( J ) + CABS1( A( I, J ) ) * S( I ) + END DO END DO - END IF - AVG = AVG + ( U + WORK( I ) ) * D / N - S( I ) = SI - END DO + END IF +* avg = s^T beta / n + AVG = 0.0D0 + DO I = 1, N + AVG = AVG + S( I )*WORK( I ) + END DO + AVG = AVG / N + + STD = 0.0D0 + DO I = N+1, N + WORK( I ) = S( I-N ) * WORK( I-N ) - AVG + END DO + CALL ZLASSQ( N, WORK( N+1 ), 1, SCALE, SUMSQ ) + STD = SCALE * SQRT( SUMSQ / N ) + + IF ( STD .LT. TOL * AVG ) GOTO 999 + + DO I = 1, N + T = CABS1( A( I, I ) ) + SI = S( I ) + C2 = ( N-1 ) * T + C1 = ( N-2 ) * ( WORK( I ) - T*SI ) + C0 = -(T*SI)*SI + 2*WORK( I )*SI - N*AVG + D = C1*C1 - 4*C0*C2 + + IF ( D .LE. 0 ) THEN + INFO = -1 + RETURN + END IF + SI = -2*C0 / ( C1 + SQRT( D ) ) + + D = SI - S( I ) + U = ZERO + IF ( UP ) THEN + DO J = 1, I + T = CABS1( A( J, I ) ) + U = U + S( J )*T + WORK( J ) = WORK( J ) + D*T + END DO + DO J = I+1,N + T = CABS1( A( I, J ) ) + U = U + S( J )*T + WORK( J ) = WORK( J ) + D*T + END DO + ELSE + DO J = 1, I + T = CABS1( A( I, J ) ) + U = U + S( J )*T + WORK( J ) = WORK( J ) + D*T + END DO + DO J = I+1,N + T = CABS1( A( J, I ) ) + U = U + S( J )*T + WORK( J ) = WORK( J ) + D*T + END DO + END IF + + AVG = AVG + ( U + WORK( I ) ) * D / N + S( I ) = SI + END DO END DO 999 CONTINUE @@ -328,10 +334,10 @@ BASE = DLAMCH( 'B' ) U = ONE / LOG( BASE ) DO I = 1, N - S( I ) = BASE ** INT( U * LOG( S( I ) * T ) ) - SMIN = MIN( SMIN, S( I ) ) - SMAX = MAX( SMAX, S( I ) ) + S( I ) = BASE ** INT( U * LOG( S( I ) * T ) ) + SMIN = MIN( SMIN, S( I ) ) + SMAX = MAX( SMAX, S( I ) ) END DO SCOND = MAX( SMIN, SMLNUM ) / MIN( SMAX, BIGNUM ) - +* END