LAPACK 3.12.0
LAPACK: Linear Algebra PACKage
Loading...
Searching...
No Matches
cgebal.f
Go to the documentation of this file.
1*> \brief \b CGEBAL
2*
3* =========== DOCUMENTATION ===========
4*
5* Online html documentation available at
6* http://www.netlib.org/lapack/explore-html/
7*
8*> \htmlonly
9*> Download CGEBAL + dependencies
10*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/cgebal.f">
11*> [TGZ]</a>
12*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/cgebal.f">
13*> [ZIP]</a>
14*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/cgebal.f">
15*> [TXT]</a>
16*> \endhtmlonly
17*
18* Definition:
19* ===========
20*
21* SUBROUTINE CGEBAL( JOB, N, A, LDA, ILO, IHI, SCALE, INFO )
22*
23* .. Scalar Arguments ..
24* CHARACTER JOB
25* INTEGER IHI, ILO, INFO, LDA, N
26* ..
27* .. Array Arguments ..
28* REAL SCALE( * )
29* COMPLEX A( LDA, * )
30* ..
31*
32*
33*> \par Purpose:
34* =============
35*>
36*> \verbatim
37*>
38*> CGEBAL balances a general complex matrix A. This involves, first,
39*> permuting A by a similarity transformation to isolate eigenvalues
40*> in the first 1 to ILO-1 and last IHI+1 to N elements on the
41*> diagonal; and second, applying a diagonal similarity transformation
42*> to rows and columns ILO to IHI to make the rows and columns as
43*> close in norm as possible. Both steps are optional.
44*>
45*> Balancing may reduce the 1-norm of the matrix, and improve the
46*> accuracy of the computed eigenvalues and/or eigenvectors.
47*> \endverbatim
48*
49* Arguments:
50* ==========
51*
52*> \param[in] JOB
53*> \verbatim
54*> JOB is CHARACTER*1
55*> Specifies the operations to be performed on A:
56*> = 'N': none: simply set ILO = 1, IHI = N, SCALE(I) = 1.0
57*> for i = 1,...,N;
58*> = 'P': permute only;
59*> = 'S': scale only;
60*> = 'B': both permute and scale.
61*> \endverbatim
62*>
63*> \param[in] N
64*> \verbatim
65*> N is INTEGER
66*> The order of the matrix A. N >= 0.
67*> \endverbatim
68*>
69*> \param[in,out] A
70*> \verbatim
71*> A is COMPLEX array, dimension (LDA,N)
72*> On entry, the input matrix A.
73*> On exit, A is overwritten by the balanced matrix.
74*> If JOB = 'N', A is not referenced.
75*> See Further Details.
76*> \endverbatim
77*>
78*> \param[in] LDA
79*> \verbatim
80*> LDA is INTEGER
81*> The leading dimension of the array A. LDA >= max(1,N).
82*> \endverbatim
83*>
84*> \param[out] ILO
85*> \verbatim
86*> ILO is INTEGER
87*> \endverbatim
88*>
89*> \param[out] IHI
90*> \verbatim
91*> IHI is INTEGER
92*> ILO and IHI are set to integers such that on exit
93*> A(i,j) = 0 if i > j and j = 1,...,ILO-1 or I = IHI+1,...,N.
94*> If JOB = 'N' or 'S', ILO = 1 and IHI = N.
95*> \endverbatim
96*>
97*> \param[out] SCALE
98*> \verbatim
99*> SCALE is REAL array, dimension (N)
100*> Details of the permutations and scaling factors applied to
101*> A. If P(j) is the index of the row and column interchanged
102*> with row and column j and D(j) is the scaling factor
103*> applied to row and column j, then
104*> SCALE(j) = P(j) for j = 1,...,ILO-1
105*> = D(j) for j = ILO,...,IHI
106*> = P(j) for j = IHI+1,...,N.
107*> The order in which the interchanges are made is N to IHI+1,
108*> then 1 to ILO-1.
109*> \endverbatim
110*>
111*> \param[out] INFO
112*> \verbatim
113*> INFO is INTEGER
114*> = 0: successful exit.
115*> < 0: if INFO = -i, the i-th argument had an illegal value.
116*> \endverbatim
117*
118* Authors:
119* ========
120*
121*> \author Univ. of Tennessee
122*> \author Univ. of California Berkeley
123*> \author Univ. of Colorado Denver
124*> \author NAG Ltd.
125*
126*> \ingroup gebal
127*
128*> \par Further Details:
129* =====================
130*>
131*> \verbatim
132*>
133*> The permutations consist of row and column interchanges which put
134*> the matrix in the form
135*>
136*> ( T1 X Y )
137*> P A P = ( 0 B Z )
138*> ( 0 0 T2 )
139*>
140*> where T1 and T2 are upper triangular matrices whose eigenvalues lie
141*> along the diagonal. The column indices ILO and IHI mark the starting
142*> and ending columns of the submatrix B. Balancing consists of applying
143*> a diagonal similarity transformation inv(D) * B * D to make the
144*> 1-norms of each row of B and its corresponding column nearly equal.
145*> The output matrix is
146*>
147*> ( T1 X*D Y )
148*> ( 0 inv(D)*B*D inv(D)*Z ).
149*> ( 0 0 T2 )
150*>
151*> Information about the permutations P and the diagonal matrix D is
152*> returned in the vector SCALE.
153*>
154*> This subroutine is based on the EISPACK routine CBAL.
155*>
156*> Modified by Tzu-Yi Chen, Computer Science Division, University of
157*> California at Berkeley, USA
158*>
159*> Refactored by Evert Provoost, Department of Computer Science,
160*> KU Leuven, Belgium
161*> \endverbatim
162*>
163* =====================================================================
164 SUBROUTINE cgebal( JOB, N, A, LDA, ILO, IHI, SCALE, INFO )
165*
166* -- LAPACK computational routine --
167* -- LAPACK is a software package provided by Univ. of Tennessee, --
168* -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
169*
170* .. Scalar Arguments ..
171 CHARACTER JOB
172 INTEGER IHI, ILO, INFO, LDA, N
173* ..
174* .. Array Arguments ..
175 REAL SCALE( * )
176 COMPLEX A( LDA, * )
177* ..
178*
179* =====================================================================
180*
181* .. Parameters ..
182 REAL ZERO, ONE
183 parameter( zero = 0.0e+0, one = 1.0e+0 )
184 REAL SCLFAC
185 parameter( sclfac = 2.0e+0 )
186 REAL FACTOR
187 parameter( factor = 0.95e+0 )
188* ..
189* .. Local Scalars ..
190 LOGICAL NOCONV, CANSWAP
191 INTEGER I, ICA, IRA, J, K, L
192 REAL C, CA, F, G, R, RA, S, SFMAX1, SFMAX2, SFMIN1,
193 $ SFMIN2
194* ..
195* .. External Functions ..
196 LOGICAL SISNAN, LSAME
197 INTEGER ICAMAX
198 REAL SLAMCH, SCNRM2
199 EXTERNAL sisnan, lsame, icamax, slamch, scnrm2
200* ..
201* .. External Subroutines ..
202 EXTERNAL xerbla, csscal, cswap
203* ..
204* .. Intrinsic Functions ..
205 INTRINSIC abs, real, aimag, max, min
206*
207* Test the input parameters
208*
209 info = 0
210 IF( .NOT.lsame( job, 'N' ) .AND. .NOT.lsame( job, 'P' ) .AND.
211 $ .NOT.lsame( job, 'S' ) .AND. .NOT.lsame( job, 'B' ) ) THEN
212 info = -1
213 ELSE IF( n.LT.0 ) THEN
214 info = -2
215 ELSE IF( lda.LT.max( 1, n ) ) THEN
216 info = -4
217 END IF
218 IF( info.NE.0 ) THEN
219 CALL xerbla( 'CGEBAL', -info )
220 RETURN
221 END IF
222*
223* Quick returns.
224*
225 IF( n.EQ.0 ) THEN
226 ilo = 1
227 ihi = 0
228 RETURN
229 END IF
230*
231 IF( lsame( job, 'N' ) ) THEN
232 DO i = 1, n
233 scale( i ) = one
234 END DO
235 ilo = 1
236 ihi = n
237 RETURN
238 END IF
239*
240* Permutation to isolate eigenvalues if possible.
241*
242 k = 1
243 l = n
244*
245 IF( .NOT.lsame( job, 'S' ) ) THEN
246*
247* Row and column exchange.
248*
249 noconv = .true.
250 DO WHILE( noconv )
251*
252* Search for rows isolating an eigenvalue and push them down.
253*
254 noconv = .false.
255 DO i = l, 1, -1
256 canswap = .true.
257 DO j = 1, l
258 IF( i.NE.j .AND. ( real( a( i, j ) ).NE.zero .OR.
259 $ aimag( a( i, j ) ).NE.zero ) ) THEN
260 canswap = .false.
261 EXIT
262 END IF
263 END DO
264*
265 IF( canswap ) THEN
266 scale( l ) = i
267 IF( i.NE.l ) THEN
268 CALL cswap( l, a( 1, i ), 1, a( 1, l ), 1 )
269 CALL cswap( n-k+1, a( i, k ), lda, a( l, k ), lda )
270 END IF
271 noconv = .true.
272*
273 IF( l.EQ.1 ) THEN
274 ilo = 1
275 ihi = 1
276 RETURN
277 END IF
278*
279 l = l - 1
280 END IF
281 END DO
282*
283 END DO
284
285 noconv = .true.
286 DO WHILE( noconv )
287*
288* Search for columns isolating an eigenvalue and push them left.
289*
290 noconv = .false.
291 DO j = k, l
292 canswap = .true.
293 DO i = k, l
294 IF( i.NE.j .AND. ( real( a( i, j ) ).NE.zero .OR.
295 $ aimag( a( i, j ) ).NE.zero ) ) THEN
296 canswap = .false.
297 EXIT
298 END IF
299 END DO
300*
301 IF( canswap ) THEN
302 scale( k ) = j
303 IF( j.NE.k ) THEN
304 CALL cswap( l, a( 1, j ), 1, a( 1, k ), 1 )
305 CALL cswap( n-k+1, a( j, k ), lda, a( k, k ), lda )
306 END IF
307 noconv = .true.
308*
309 k = k + 1
310 END IF
311 END DO
312*
313 END DO
314*
315 END IF
316*
317* Initialize SCALE for non-permuted submatrix.
318*
319 DO i = k, l
320 scale( i ) = one
321 END DO
322*
323* If we only had to permute, we are done.
324*
325 IF( lsame( job, 'P' ) ) THEN
326 ilo = k
327 ihi = l
328 RETURN
329 END IF
330*
331* Balance the submatrix in rows K to L.
332*
333* Iterative loop for norm reduction.
334*
335 sfmin1 = slamch( 'S' ) / slamch( 'P' )
336 sfmax1 = one / sfmin1
337 sfmin2 = sfmin1*sclfac
338 sfmax2 = one / sfmin2
339*
340 noconv = .true.
341 DO WHILE( noconv )
342 noconv = .false.
343*
344 DO i = k, l
345*
346 c = scnrm2( l-k+1, a( k, i ), 1 )
347 r = scnrm2( l-k+1, a( i, k ), lda )
348 ica = icamax( l, a( 1, i ), 1 )
349 ca = abs( a( ica, i ) )
350 ira = icamax( n-k+1, a( i, k ), lda )
351 ra = abs( a( i, ira+k-1 ) )
352*
353* Guard against zero C or R due to underflow.
354*
355 IF( c.EQ.zero .OR. r.EQ.zero ) cycle
356*
357* Exit if NaN to avoid infinite loop
358*
359 IF( sisnan( c+ca+r+ra ) ) THEN
360 info = -3
361 CALL xerbla( 'CGEBAL', -info )
362 RETURN
363 END IF
364*
365 g = r / sclfac
366 f = one
367 s = c + r
368*
369 DO WHILE( c.LT.g .AND. max( f, c, ca ).LT.sfmax2 .AND.
370 $ min( r, g, ra ).GT.sfmin2 )
371 f = f*sclfac
372 c = c*sclfac
373 ca = ca*sclfac
374 r = r / sclfac
375 g = g / sclfac
376 ra = ra / sclfac
377 END DO
378*
379 g = c / sclfac
380*
381 DO WHILE( g.GE.r .AND. max( r, ra ).LT.sfmax2 .AND.
382 $ min( f, c, g, ca ).GT.sfmin2 )
383 f = f / sclfac
384 c = c / sclfac
385 g = g / sclfac
386 ca = ca / sclfac
387 r = r*sclfac
388 ra = ra*sclfac
389 END DO
390*
391* Now balance.
392*
393 IF( ( c+r ).GE.factor*s ) cycle
394 IF( f.LT.one .AND. scale( i ).LT.one ) THEN
395 IF( f*scale( i ).LE.sfmin1 ) cycle
396 END IF
397 IF( f.GT.one .AND. scale( i ).GT.one ) THEN
398 IF( scale( i ).GE.sfmax1 / f ) cycle
399 END IF
400 g = one / f
401 scale( i ) = scale( i )*f
402 noconv = .true.
403*
404 CALL csscal( n-k+1, g, a( i, k ), lda )
405 CALL csscal( l, f, a( 1, i ), 1 )
406*
407 END DO
408*
409 END DO
410*
411 ilo = k
412 ihi = l
413*
414 RETURN
415*
416* End of CGEBAL
417*
418 END
subroutine xerbla(srname, info)
Definition cblat2.f:3285
subroutine cgebal(job, n, a, lda, ilo, ihi, scale, info)
CGEBAL
Definition cgebal.f:165
subroutine csscal(n, sa, cx, incx)
CSSCAL
Definition csscal.f:78
subroutine cswap(n, cx, incx, cy, incy)
CSWAP
Definition cswap.f:81