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