SUBROUTINE DLAQSP( UPLO, N, AP, S, SCOND, AMAX, EQUED )
*
* -- LAPACK auxiliary routine (version 3.1) --
* Univ. of Tennessee, Univ. of California Berkeley and NAG Ltd..
* November 2006
*
* .. Scalar Arguments ..
CHARACTER EQUED, UPLO
INTEGER N
DOUBLE PRECISION AMAX, SCOND
* ..
* .. Array Arguments ..
DOUBLE PRECISION AP( * ), S( * )
* ..
*
* Purpose
* =======
*
* DLAQSP equilibrates a symmetric matrix A using the scaling factors
* in the vector S.
*
* Arguments
* =========
*
* UPLO (input) CHARACTER*1
* Specifies whether the upper or lower triangular part of the
* symmetric matrix A is stored.
* = 'U': Upper triangular
* = 'L': Lower triangular
*
* N (input) INTEGER
* The order of the matrix A. N >= 0.
*
* AP (input/output) DOUBLE PRECISION array, dimension (N*(N+1)/2)
* On entry, the upper or lower triangle of the symmetric matrix
* A, packed columnwise in a linear array. The j-th column of A
* is stored in the array AP as follows:
* if UPLO = 'U', AP(i + (j-1)*j/2) = A(i,j) for 1<=i<=j;
* if UPLO = 'L', AP(i + (j-1)*(2n-j)/2) = A(i,j) for j<=i<=n.
*
* On exit, the equilibrated matrix: diag(S) * A * diag(S), in
* the same storage format as A.
*
* S (input) DOUBLE PRECISION array, dimension (N)
* The scale factors for A.
*
* SCOND (input) DOUBLE PRECISION
* Ratio of the smallest S(i) to the largest S(i).
*
* AMAX (input) DOUBLE PRECISION
* Absolute value of largest matrix entry.
*
* EQUED (output) CHARACTER*1
* Specifies whether or not equilibration was done.
* = 'N': No equilibration.
* = 'Y': Equilibration was done, i.e., A has been replaced by
* diag(S) * A * diag(S).
*
* Internal Parameters
* ===================
*
* THRESH is a threshold value used to decide if scaling should be done
* based on the ratio of the scaling factors. If SCOND < THRESH,
* scaling is done.
*
* LARGE and SMALL are threshold values used to decide if scaling should
* be done based on the absolute size of the largest matrix element.
* If AMAX > LARGE or AMAX < SMALL, scaling is done.
*
* =====================================================================
*
* .. Parameters ..
DOUBLE PRECISION ONE, THRESH
PARAMETER ( ONE = 1.0D+0, THRESH = 0.1D+0 )
* ..
* .. Local Scalars ..
INTEGER I, J, JC
DOUBLE PRECISION CJ, LARGE, SMALL
* ..
* .. External Functions ..
LOGICAL LSAME
DOUBLE PRECISION DLAMCH
EXTERNAL LSAME, DLAMCH
* ..
* .. Executable Statements ..
*
* Quick return if possible
*
IF( N.LE.0 ) THEN
EQUED = 'N'
RETURN
END IF
*
* Initialize LARGE and SMALL.
*
SMALL = DLAMCH( 'Safe minimum' ) / DLAMCH( 'Precision' )
LARGE = ONE / SMALL
*
IF( SCOND.GE.THRESH .AND. AMAX.GE.SMALL .AND. AMAX.LE.LARGE ) THEN
*
* No equilibration
*
EQUED = 'N'
ELSE
*
* Replace A by diag(S) * A * diag(S).
*
IF( LSAME( UPLO, 'U' ) ) THEN
*
* Upper triangle of A is stored.
*
JC = 1
DO 20 J = 1, N
CJ = S( J )
DO 10 I = 1, J
AP( JC+I-1 ) = CJ*S( I )*AP( JC+I-1 )
10 CONTINUE
JC = JC + J
20 CONTINUE
ELSE
*
* Lower triangle of A is stored.
*
JC = 1
DO 40 J = 1, N
CJ = S( J )
DO 30 I = J, N
AP( JC+I-J ) = CJ*S( I )*AP( JC+I-J )
30 CONTINUE
JC = JC + N - J + 1
40 CONTINUE
END IF
EQUED = 'Y'
END IF
*
RETURN
*
* End of DLAQSP
*
END