1       SUBROUTINE DLAQSP( UPLO, N, AP, S, SCOND, AMAX, EQUED )
  2 *
  3 *  -- LAPACK auxiliary routine (version 3.2) --
  4 *  -- LAPACK is a software package provided by Univ. of Tennessee,    --
  5 *  -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
  6 *     November 2006
  7 *
  8 *     .. Scalar Arguments ..
  9       CHARACTER          EQUED, UPLO
 10       INTEGER            N
 11       DOUBLE PRECISION   AMAX, SCOND
 12 *     ..
 13 *     .. Array Arguments ..
 14       DOUBLE PRECISION   AP( * ), S( * )
 15 *     ..
 16 *
 17 *  Purpose
 18 *  =======
 19 *
 20 *  DLAQSP equilibrates a symmetric matrix A using the scaling factors
 21 *  in the vector S.
 22 *
 23 *  Arguments
 24 *  =========
 25 *
 26 *  UPLO    (input) CHARACTER*1
 27 *          Specifies whether the upper or lower triangular part of the
 28 *          symmetric matrix A is stored.
 29 *          = 'U':  Upper triangular
 30 *          = 'L':  Lower triangular
 31 *
 32 *  N       (input) INTEGER
 33 *          The order of the matrix A.  N >= 0.
 34 *
 35 *  AP      (input/output) DOUBLE PRECISION array, dimension (N*(N+1)/2)
 36 *          On entry, the upper or lower triangle of the symmetric matrix
 37 *          A, packed columnwise in a linear array.  The j-th column of A
 38 *          is stored in the array AP as follows:
 39 *          if UPLO = 'U', AP(i + (j-1)*j/2) = A(i,j) for 1<=i<=j;
 40 *          if UPLO = 'L', AP(i + (j-1)*(2n-j)/2) = A(i,j) for j<=i<=n.
 41 *
 42 *          On exit, the equilibrated matrix:  diag(S) * A * diag(S), in
 43 *          the same storage format as A.
 44 *
 45 *  S       (input) DOUBLE PRECISION array, dimension (N)
 46 *          The scale factors for A.
 47 *
 48 *  SCOND   (input) DOUBLE PRECISION
 49 *          Ratio of the smallest S(i) to the largest S(i).
 50 *
 51 *  AMAX    (input) DOUBLE PRECISION
 52 *          Absolute value of largest matrix entry.
 53 *
 54 *  EQUED   (output) CHARACTER*1
 55 *          Specifies whether or not equilibration was done.
 56 *          = 'N':  No equilibration.
 57 *          = 'Y':  Equilibration was done, i.e., A has been replaced by
 58 *                  diag(S) * A * diag(S).
 59 *
 60 *  Internal Parameters
 61 *  ===================
 62 *
 63 *  THRESH is a threshold value used to decide if scaling should be done
 64 *  based on the ratio of the scaling factors.  If SCOND < THRESH,
 65 *  scaling is done.
 66 *
 67 *  LARGE and SMALL are threshold values used to decide if scaling should
 68 *  be done based on the absolute size of the largest matrix element.
 69 *  If AMAX > LARGE or AMAX < SMALL, scaling is done.
 70 *
 71 *  =====================================================================
 72 *
 73 *     .. Parameters ..
 74       DOUBLE PRECISION   ONE, THRESH
 75       PARAMETER          ( ONE = 1.0D+0, THRESH = 0.1D+0 )
 76 *     ..
 77 *     .. Local Scalars ..
 78       INTEGER            I, J, JC
 79       DOUBLE PRECISION   CJ, LARGE, SMALL
 80 *     ..
 81 *     .. External Functions ..
 82       LOGICAL            LSAME
 83       DOUBLE PRECISION   DLAMCH
 84       EXTERNAL           LSAME, DLAMCH
 85 *     ..
 86 *     .. Executable Statements ..
 87 *
 88 *     Quick return if possible
 89 *
 90       IF( N.LE.0 ) THEN
 91          EQUED = 'N'
 92          RETURN
 93       END IF
 94 *
 95 *     Initialize LARGE and SMALL.
 96 *
 97       SMALL = DLAMCH( 'Safe minimum' ) / DLAMCH( 'Precision' )
 98       LARGE = ONE / SMALL
 99 *
100       IF( SCOND.GE.THRESH .AND. AMAX.GE.SMALL .AND. AMAX.LE.LARGE ) THEN
101 *
102 *        No equilibration
103 *
104          EQUED = 'N'
105       ELSE
106 *
107 *        Replace A by diag(S) * A * diag(S).
108 *
109          IF( LSAME( UPLO, 'U' ) ) THEN
110 *
111 *           Upper triangle of A is stored.
112 *
113             JC = 1
114             DO 20 J = 1, N
115                CJ = S( J )
116                DO 10 I = 1, J
117                   AP( JC+I-1 ) = CJ*S( I )*AP( JC+I-1 )
118    10          CONTINUE
119                JC = JC + J
120    20       CONTINUE
121          ELSE
122 *
123 *           Lower triangle of A is stored.
124 *
125             JC = 1
126             DO 40 J = 1, N
127                CJ = S( J )
128                DO 30 I = J, N
129                   AP( JC+I-J ) = CJ*S( I )*AP( JC+I-J )
130    30          CONTINUE
131                JC = JC + N - J + 1
132    40       CONTINUE
133          END IF
134          EQUED = 'Y'
135       END IF
136 *
137       RETURN
138 *
139 *     End of DLAQSP
140 *
141       END