1       SUBROUTINE ZPBEQU( UPLO, N, KD, AB, LDAB, S, SCOND, AMAX, INFO )
  2 *
  3 *  -- LAPACK 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          UPLO
 10       INTEGER            INFO, KD, LDAB, N
 11       DOUBLE PRECISION   AMAX, SCOND
 12 *     ..
 13 *     .. Array Arguments ..
 14       DOUBLE PRECISION   S( * )
 15       COMPLEX*16         AB( LDAB, * )
 16 *     ..
 17 *
 18 *  Purpose
 19 *  =======
 20 *
 21 *  ZPBEQU computes row and column scalings intended to equilibrate a
 22 *  Hermitian positive definite band matrix A and reduce its condition
 23 *  number (with respect to the two-norm).  S contains the scale factors,
 24 *  S(i) = 1/sqrt(A(i,i)), chosen so that the scaled matrix B with
 25 *  elements B(i,j) = S(i)*A(i,j)*S(j) has ones on the diagonal.  This
 26 *  choice of S puts the condition number of B within a factor N of the
 27 *  smallest possible condition number over all possible diagonal
 28 *  scalings.
 29 *
 30 *  Arguments
 31 *  =========
 32 *
 33 *  UPLO    (input) CHARACTER*1
 34 *          = 'U':  Upper triangular of A is stored;
 35 *          = 'L':  Lower triangular of A is stored.
 36 *
 37 *  N       (input) INTEGER
 38 *          The order of the matrix A.  N >= 0.
 39 *
 40 *  KD      (input) INTEGER
 41 *          The number of superdiagonals of the matrix A if UPLO = 'U',
 42 *          or the number of subdiagonals if UPLO = 'L'.  KD >= 0.
 43 *
 44 *  AB      (input) COMPLEX*16 array, dimension (LDAB,N)
 45 *          The upper or lower triangle of the Hermitian band matrix A,
 46 *          stored in the first KD+1 rows of the array.  The j-th column
 47 *          of A is stored in the j-th column of the array AB as follows:
 48 *          if UPLO = 'U', AB(kd+1+i-j,j) = A(i,j) for max(1,j-kd)<=i<=j;
 49 *          if UPLO = 'L', AB(1+i-j,j)    = A(i,j) for j<=i<=min(n,j+kd).
 50 *
 51 *  LDAB     (input) INTEGER
 52 *          The leading dimension of the array A.  LDAB >= KD+1.
 53 *
 54 *  S       (output) DOUBLE PRECISION array, dimension (N)
 55 *          If INFO = 0, S contains the scale factors for A.
 56 *
 57 *  SCOND   (output) DOUBLE PRECISION
 58 *          If INFO = 0, S contains the ratio of the smallest S(i) to
 59 *          the largest S(i).  If SCOND >= 0.1 and AMAX is neither too
 60 *          large nor too small, it is not worth scaling by S.
 61 *
 62 *  AMAX    (output) DOUBLE PRECISION
 63 *          Absolute value of largest matrix element.  If AMAX is very
 64 *          close to overflow or very close to underflow, the matrix
 65 *          should be scaled.
 66 *
 67 *  INFO    (output) INTEGER
 68 *          = 0:  successful exit
 69 *          < 0:  if INFO = -i, the i-th argument had an illegal value.
 70 *          > 0:  if INFO = i, the i-th diagonal element is nonpositive.
 71 *
 72 *  =====================================================================
 73 *
 74 *     .. Parameters ..
 75       DOUBLE PRECISION   ZERO, ONE
 76       PARAMETER          ( ZERO = 0.0D+0, ONE = 1.0D+0 )
 77 *     ..
 78 *     .. Local Scalars ..
 79       LOGICAL            UPPER
 80       INTEGER            I, J
 81       DOUBLE PRECISION   SMIN
 82 *     ..
 83 *     .. External Functions ..
 84       LOGICAL            LSAME
 85       EXTERNAL           LSAME
 86 *     ..
 87 *     .. External Subroutines ..
 88       EXTERNAL           XERBLA
 89 *     ..
 90 *     .. Intrinsic Functions ..
 91       INTRINSIC          DBLEMAXMINSQRT
 92 *     ..
 93 *     .. Executable Statements ..
 94 *
 95 *     Test the input parameters.
 96 *
 97       INFO = 0
 98       UPPER = LSAME( UPLO, 'U' )
 99       IF.NOT.UPPER .AND. .NOT.LSAME( UPLO, 'L' ) ) THEN
100          INFO = -1
101       ELSE IF( N.LT.0 ) THEN
102          INFO = -2
103       ELSE IF( KD.LT.0 ) THEN
104          INFO = -3
105       ELSE IF( LDAB.LT.KD+1 ) THEN
106          INFO = -5
107       END IF
108       IF( INFO.NE.0 ) THEN
109          CALL XERBLA( 'ZPBEQU'-INFO )
110          RETURN
111       END IF
112 *
113 *     Quick return if possible
114 *
115       IF( N.EQ.0 ) THEN
116          SCOND = ONE
117          AMAX = ZERO
118          RETURN
119       END IF
120 *
121       IF( UPPER ) THEN
122          J = KD + 1
123       ELSE
124          J = 1
125       END IF
126 *
127 *     Initialize SMIN and AMAX.
128 *
129       S( 1 ) = DBLE( AB( J, 1 ) )
130       SMIN = S( 1 )
131       AMAX = S( 1 )
132 *
133 *     Find the minimum and maximum diagonal elements.
134 *
135       DO 10 I = 2, N
136          S( I ) = DBLE( AB( J, I ) )
137          SMIN = MIN( SMIN, S( I ) )
138          AMAX = MAX( AMAX, S( I ) )
139    10 CONTINUE
140 *
141       IF( SMIN.LE.ZERO ) THEN
142 *
143 *        Find the first non-positive diagonal element and return.
144 *
145          DO 20 I = 1, N
146             IF( S( I ).LE.ZERO ) THEN
147                INFO = I
148                RETURN
149             END IF
150    20    CONTINUE
151       ELSE
152 *
153 *        Set the scale factors to the reciprocals
154 *        of the diagonal elements.
155 *
156          DO 30 I = 1, N
157             S( I ) = ONE / SQRT( S( I ) )
158    30    CONTINUE
159 *
160 *        Compute SCOND = min(S(I)) / max(S(I))
161 *
162          SCOND = SQRT( SMIN ) / SQRT( AMAX )
163       END IF
164       RETURN
165 *
166 *     End of ZPBEQU
167 *
168       END