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
CVSweb interface <joel.bertrand@systella.fr>