--- rpl/lapack/lapack/zsyequb.f 2011/11/21 20:43:21 1.5 +++ rpl/lapack/lapack/zsyequb.f 2018/05/29 07:18:36 1.14 @@ -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 ZSYEQUB + dependencies -*> -*> [TGZ] -*> -*> [ZIP] -*> +*> Download ZSYEQUB + dependencies +*> +*> [TGZ] +*> +*> [ZIP] +*> *> [TXT] -*> \endhtmlonly +*> \endhtmlonly * * Definition: * =========== * * SUBROUTINE ZSYEQUB( 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 *> *> ZSYEQUB computes row and column scalings intended to equilibrate a -*> symmetric 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 +*> symmetric 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,30 +51,27 @@ *> \param[in] UPLO *> \verbatim *> UPLO is CHARACTER*1 -*> Specifies whether the details of the factorization are stored -*> as an upper or lower triangular matrix. -*> = 'U': Upper triangular, form is A = U*D*U**T; -*> = 'L': Lower triangular, form is A = L*D*L**T. +*> = '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 symmetric matrix whose scaling -*> factors are to be computed. Only the diagonal elements of A -*> are referenced. +*> The N-by-N symmetric 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 @@ -88,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 @@ -116,12 +112,12 @@ * 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 November 2011 +*> \date November 2017 * *> \ingroup complex16SYcomputational * @@ -130,16 +126,16 @@ *> *> 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://ruready.utah.edu/archive/papers/bin.pdf +*> 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 ZSYEQUB( UPLO, N, A, LDA, S, SCOND, AMAX, WORK, INFO ) * -* -- LAPACK computational routine (version 3.4.0) -- +* -- LAPACK computational routine (version 3.8.0) -- * -- LAPACK is a software package provided by Univ. of Tennessee, -- * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..-- -* November 2011 +* November 2017 * * .. Scalar Arguments .. INTEGER INFO, LDA, N @@ -172,7 +168,7 @@ EXTERNAL DLAMCH, LSAME * .. * .. External Subroutines .. - EXTERNAL ZLASSQ + EXTERNAL ZLASSQ, XERBLA * .. * .. Intrinsic Functions .. INTRINSIC ABS, DBLE, DIMAG, INT, LOG, MAX, MIN, SQRT @@ -180,7 +176,7 @@ * .. Statement Functions .. DOUBLE PRECISION CABS1 * .. -* Statement Function Definitions +* .. Statement Function Definitions .. CABS1( ZDUM ) = ABS( DBLE( ZDUM ) ) + ABS( DIMAG( ZDUM ) ) * .. * .. Executable Statements .. @@ -189,15 +185,15 @@ * INFO = 0 IF ( .NOT. ( LSAME( UPLO, 'U' ) .OR. LSAME( UPLO, 'L' ) ) ) THEN - INFO = -1 + 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( 'ZSYEQUB', -INFO ) - RETURN + CALL XERBLA( 'ZSYEQUB', -INFO ) + RETURN END IF UP = LSAME( UPLO, 'U' ) @@ -206,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 @@ -222,7 +218,7 @@ S( J ) = MAX( S( J ), CABS1( A( I, J ) ) ) AMAX = MAX( AMAX, CABS1( A( I, J ) ) ) END DO - S( J ) = MAX( S( J ), CABS1( A( J, J) ) ) + S( J ) = MAX( S( J ), CABS1( A( J, J ) ) ) AMAX = MAX( AMAX, CABS1( A( J, J ) ) ) END DO ELSE @@ -231,102 +227,101 @@ AMAX = MAX( AMAX, CABS1( A( J, J ) ) ) DO I = J+1, N S( I ) = MAX( S( I ), CABS1( A( I, J ) ) ) - S( J ) = MAX( S( J ), CABS1 (A( I, J ) ) ) + S( J ) = MAX( S( J ), 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 = N+1, 2*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 + 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, 2*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 @@ -339,9 +334,9 @@ 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 ) *