Book a Demo!
CoCalc Logo Icon
StoreFeaturesDocsShareSupportNewsAboutPoliciesSign UpSign In
ElmerCSC
GitHub Repository: ElmerCSC/elmerfem
Path: blob/devel/mathlibs/src/lapack/cgebal.f
5212 views
1
SUBROUTINE CGEBAL( JOB, N, A, LDA, ILO, IHI, SCALE, INFO )
2
*
3
* -- LAPACK routine (version 3.0) --
4
* Univ. of Tennessee, Univ. of California Berkeley, NAG Ltd.,
5
* Courant Institute, Argonne National Lab, and Rice University
6
* June 30, 1999
7
*
8
* .. Scalar Arguments ..
9
CHARACTER JOB
10
INTEGER IHI, ILO, INFO, LDA, N
11
* ..
12
* .. Array Arguments ..
13
REAL SCALE( * )
14
COMPLEX A( LDA, * )
15
* ..
16
*
17
* Purpose
18
* =======
19
*
20
* CGEBAL balances a general complex matrix A. This involves, first,
21
* permuting A by a similarity transformation to isolate eigenvalues
22
* in the first 1 to ILO-1 and last IHI+1 to N elements on the
23
* diagonal; and second, applying a diagonal similarity transformation
24
* to rows and columns ILO to IHI to make the rows and columns as
25
* close in norm as possible. Both steps are optional.
26
*
27
* Balancing may reduce the 1-norm of the matrix, and improve the
28
* accuracy of the computed eigenvalues and/or eigenvectors.
29
*
30
* Arguments
31
* =========
32
*
33
* JOB (input) CHARACTER*1
34
* Specifies the operations to be performed on A:
35
* = 'N': none: simply set ILO = 1, IHI = N, SCALE(I) = 1.0
36
* for i = 1,...,N;
37
* = 'P': permute only;
38
* = 'S': scale only;
39
* = 'B': both permute and scale.
40
*
41
* N (input) INTEGER
42
* The order of the matrix A. N >= 0.
43
*
44
* A (input/output) COMPLEX array, dimension (LDA,N)
45
* On entry, the input matrix A.
46
* On exit, A is overwritten by the balanced matrix.
47
* If JOB = 'N', A is not referenced.
48
* See Further Details.
49
*
50
* LDA (input) INTEGER
51
* The leading dimension of the array A. LDA >= max(1,N).
52
*
53
* ILO (output) INTEGER
54
* IHI (output) INTEGER
55
* ILO and IHI are set to integers such that on exit
56
* A(i,j) = 0 if i > j and j = 1,...,ILO-1 or I = IHI+1,...,N.
57
* If JOB = 'N' or 'S', ILO = 1 and IHI = N.
58
*
59
* SCALE (output) REAL array, dimension (N)
60
* Details of the permutations and scaling factors applied to
61
* A. If P(j) is the index of the row and column interchanged
62
* with row and column j and D(j) is the scaling factor
63
* applied to row and column j, then
64
* SCALE(j) = P(j) for j = 1,...,ILO-1
65
* = D(j) for j = ILO,...,IHI
66
* = P(j) for j = IHI+1,...,N.
67
* The order in which the interchanges are made is N to IHI+1,
68
* then 1 to ILO-1.
69
*
70
* INFO (output) INTEGER
71
* = 0: successful exit.
72
* < 0: if INFO = -i, the i-th argument had an illegal value.
73
*
74
* Further Details
75
* ===============
76
*
77
* The permutations consist of row and column interchanges which put
78
* the matrix in the form
79
*
80
* ( T1 X Y )
81
* P A P = ( 0 B Z )
82
* ( 0 0 T2 )
83
*
84
* where T1 and T2 are upper triangular matrices whose eigenvalues lie
85
* along the diagonal. The column indices ILO and IHI mark the starting
86
* and ending columns of the submatrix B. Balancing consists of applying
87
* a diagonal similarity transformation inv(D) * B * D to make the
88
* 1-norms of each row of B and its corresponding column nearly equal.
89
* The output matrix is
90
*
91
* ( T1 X*D Y )
92
* ( 0 inv(D)*B*D inv(D)*Z ).
93
* ( 0 0 T2 )
94
*
95
* Information about the permutations P and the diagonal matrix D is
96
* returned in the vector SCALE.
97
*
98
* This subroutine is based on the EISPACK routine CBAL.
99
*
100
* Modified by Tzu-Yi Chen, Computer Science Division, University of
101
* California at Berkeley, USA
102
*
103
* =====================================================================
104
*
105
* .. Parameters ..
106
REAL ZERO, ONE
107
PARAMETER ( ZERO = 0.0E+0, ONE = 1.0E+0 )
108
REAL SCLFAC
109
PARAMETER ( SCLFAC = 0.8E+1 )
110
REAL FACTOR
111
PARAMETER ( FACTOR = 0.95E+0 )
112
* ..
113
* .. Local Scalars ..
114
LOGICAL NOCONV
115
INTEGER I, ICA, IEXC, IRA, J, K, L, M
116
REAL C, CA, F, G, R, RA, S, SFMAX1, SFMAX2, SFMIN1,
117
$ SFMIN2
118
COMPLEX CDUM
119
* ..
120
* .. External Functions ..
121
LOGICAL LSAME
122
INTEGER ICAMAX
123
REAL SLAMCH
124
EXTERNAL LSAME, ICAMAX, SLAMCH
125
* ..
126
* .. External Subroutines ..
127
EXTERNAL CSSCAL, CSWAP, XERBLA
128
* ..
129
* .. Intrinsic Functions ..
130
INTRINSIC ABS, AIMAG, MAX, MIN, REAL
131
* ..
132
* .. Statement Functions ..
133
REAL CABS1
134
* ..
135
* .. Statement Function definitions ..
136
CABS1( CDUM ) = ABS( REAL( CDUM ) ) + ABS( AIMAG( CDUM ) )
137
* ..
138
* .. Executable Statements ..
139
*
140
* Test the input parameters
141
*
142
INFO = 0
143
IF( .NOT.LSAME( JOB, 'N' ) .AND. .NOT.LSAME( JOB, 'P' ) .AND.
144
$ .NOT.LSAME( JOB, 'S' ) .AND. .NOT.LSAME( JOB, 'B' ) ) THEN
145
INFO = -1
146
ELSE IF( N.LT.0 ) THEN
147
INFO = -2
148
ELSE IF( LDA.LT.MAX( 1, N ) ) THEN
149
INFO = -4
150
END IF
151
IF( INFO.NE.0 ) THEN
152
CALL XERBLA( 'CGEBAL', -INFO )
153
RETURN
154
END IF
155
*
156
K = 1
157
L = N
158
*
159
IF( N.EQ.0 )
160
$ GO TO 210
161
*
162
IF( LSAME( JOB, 'N' ) ) THEN
163
DO 10 I = 1, N
164
SCALE( I ) = ONE
165
10 CONTINUE
166
GO TO 210
167
END IF
168
*
169
IF( LSAME( JOB, 'S' ) )
170
$ GO TO 120
171
*
172
* Permutation to isolate eigenvalues if possible
173
*
174
GO TO 50
175
*
176
* Row and column exchange.
177
*
178
20 CONTINUE
179
SCALE( M ) = J
180
IF( J.EQ.M )
181
$ GO TO 30
182
*
183
CALL CSWAP( L, A( 1, J ), 1, A( 1, M ), 1 )
184
CALL CSWAP( N-K+1, A( J, K ), LDA, A( M, K ), LDA )
185
*
186
30 CONTINUE
187
GO TO ( 40, 80 )IEXC
188
*
189
* Search for rows isolating an eigenvalue and push them down.
190
*
191
40 CONTINUE
192
IF( L.EQ.1 )
193
$ GO TO 210
194
L = L - 1
195
*
196
50 CONTINUE
197
DO 70 J = L, 1, -1
198
*
199
DO 60 I = 1, L
200
IF( I.EQ.J )
201
$ GO TO 60
202
IF( REAL( A( J, I ) ).NE.ZERO .OR. AIMAG( A( J, I ) ).NE.
203
$ ZERO )GO TO 70
204
60 CONTINUE
205
*
206
M = L
207
IEXC = 1
208
GO TO 20
209
70 CONTINUE
210
*
211
GO TO 90
212
*
213
* Search for columns isolating an eigenvalue and push them left.
214
*
215
80 CONTINUE
216
K = K + 1
217
*
218
90 CONTINUE
219
DO 110 J = K, L
220
*
221
DO 100 I = K, L
222
IF( I.EQ.J )
223
$ GO TO 100
224
IF( REAL( A( I, J ) ).NE.ZERO .OR. AIMAG( A( I, J ) ).NE.
225
$ ZERO )GO TO 110
226
100 CONTINUE
227
*
228
M = K
229
IEXC = 2
230
GO TO 20
231
110 CONTINUE
232
*
233
120 CONTINUE
234
DO 130 I = K, L
235
SCALE( I ) = ONE
236
130 CONTINUE
237
*
238
IF( LSAME( JOB, 'P' ) )
239
$ GO TO 210
240
*
241
* Balance the submatrix in rows K to L.
242
*
243
* Iterative loop for norm reduction
244
*
245
SFMIN1 = SLAMCH( 'S' ) / SLAMCH( 'P' )
246
SFMAX1 = ONE / SFMIN1
247
SFMIN2 = SFMIN1*SCLFAC
248
SFMAX2 = ONE / SFMIN2
249
140 CONTINUE
250
NOCONV = .FALSE.
251
*
252
DO 200 I = K, L
253
C = ZERO
254
R = ZERO
255
*
256
DO 150 J = K, L
257
IF( J.EQ.I )
258
$ GO TO 150
259
C = C + CABS1( A( J, I ) )
260
R = R + CABS1( A( I, J ) )
261
150 CONTINUE
262
ICA = ICAMAX( L, A( 1, I ), 1 )
263
CA = ABS( A( ICA, I ) )
264
IRA = ICAMAX( N-K+1, A( I, K ), LDA )
265
RA = ABS( A( I, IRA+K-1 ) )
266
*
267
* Guard against zero C or R due to underflow.
268
*
269
IF( C.EQ.ZERO .OR. R.EQ.ZERO )
270
$ GO TO 200
271
G = R / SCLFAC
272
F = ONE
273
S = C + R
274
160 CONTINUE
275
IF( C.GE.G .OR. MAX( F, C, CA ).GE.SFMAX2 .OR.
276
$ MIN( R, G, RA ).LE.SFMIN2 )GO TO 170
277
F = F*SCLFAC
278
C = C*SCLFAC
279
CA = CA*SCLFAC
280
R = R / SCLFAC
281
G = G / SCLFAC
282
RA = RA / SCLFAC
283
GO TO 160
284
*
285
170 CONTINUE
286
G = C / SCLFAC
287
180 CONTINUE
288
IF( G.LT.R .OR. MAX( R, RA ).GE.SFMAX2 .OR.
289
$ MIN( F, C, G, CA ).LE.SFMIN2 )GO TO 190
290
F = F / SCLFAC
291
C = C / SCLFAC
292
G = G / SCLFAC
293
CA = CA / SCLFAC
294
R = R*SCLFAC
295
RA = RA*SCLFAC
296
GO TO 180
297
*
298
* Now balance.
299
*
300
190 CONTINUE
301
IF( ( C+R ).GE.FACTOR*S )
302
$ GO TO 200
303
IF( F.LT.ONE .AND. SCALE( I ).LT.ONE ) THEN
304
IF( F*SCALE( I ).LE.SFMIN1 )
305
$ GO TO 200
306
END IF
307
IF( F.GT.ONE .AND. SCALE( I ).GT.ONE ) THEN
308
IF( SCALE( I ).GE.SFMAX1 / F )
309
$ GO TO 200
310
END IF
311
G = ONE / F
312
SCALE( I ) = SCALE( I )*F
313
NOCONV = .TRUE.
314
*
315
CALL CSSCAL( N-K+1, G, A( I, K ), LDA )
316
CALL CSSCAL( L, F, A( 1, I ), 1 )
317
*
318
200 CONTINUE
319
*
320
IF( NOCONV )
321
$ GO TO 140
322
*
323
210 CONTINUE
324
ILO = K
325
IHI = L
326
*
327
RETURN
328
*
329
* End of CGEBAL
330
*
331
END
332
333