LAPACK  3.10.1
LAPACK: Linear Algebra PACKage
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 *> \param[out] IHI
89 *> \verbatim
90 *> IHI is INTEGER
91 *> ILO and IHI are set to integers such that on exit
92 *> A(i,j) = 0 if i > j and j = 1,...,ILO-1 or I = IHI+1,...,N.
93 *> If JOB = 'N' or 'S', ILO = 1 and IHI = N.
94 *> \endverbatim
95 *>
96 *> \param[out] SCALE
97 *> \verbatim
98 *> SCALE is REAL array, dimension (N)
99 *> Details of the permutations and scaling factors applied to
100 *> A. If P(j) is the index of the row and column interchanged
101 *> with row and column j and D(j) is the scaling factor
102 *> applied to row and column j, then
103 *> SCALE(j) = P(j) for j = 1,...,ILO-1
104 *> = D(j) for j = ILO,...,IHI
105 *> = P(j) for j = IHI+1,...,N.
106 *> The order in which the interchanges are made is N to IHI+1,
107 *> then 1 to ILO-1.
108 *> \endverbatim
109 *>
110 *> \param[out] INFO
111 *> \verbatim
112 *> INFO is INTEGER
113 *> = 0: successful exit.
114 *> < 0: if INFO = -i, the i-th argument had an illegal value.
115 *> \endverbatim
116 *
117 * Authors:
118 * ========
119 *
120 *> \author Univ. of Tennessee
121 *> \author Univ. of California Berkeley
122 *> \author Univ. of Colorado Denver
123 *> \author NAG Ltd.
124 *
125 *> \ingroup complexGEcomputational
126 *
127 *> \par Further Details:
128 * =====================
129 *>
130 *> \verbatim
131 *>
132 *> The permutations consist of row and column interchanges which put
133 *> the matrix in the form
134 *>
135 *> ( T1 X Y )
136 *> P A P = ( 0 B Z )
137 *> ( 0 0 T2 )
138 *>
139 *> where T1 and T2 are upper triangular matrices whose eigenvalues lie
140 *> along the diagonal. The column indices ILO and IHI mark the starting
141 *> and ending columns of the submatrix B. Balancing consists of applying
142 *> a diagonal similarity transformation inv(D) * B * D to make the
143 *> 1-norms of each row of B and its corresponding column nearly equal.
144 *> The output matrix is
145 *>
146 *> ( T1 X*D Y )
147 *> ( 0 inv(D)*B*D inv(D)*Z ).
148 *> ( 0 0 T2 )
149 *>
150 *> Information about the permutations P and the diagonal matrix D is
151 *> returned in the vector SCALE.
152 *>
153 *> This subroutine is based on the EISPACK routine CBAL.
154 *>
155 *> Modified by Tzu-Yi Chen, Computer Science Division, University of
156 *> California at Berkeley, USA
157 *> \endverbatim
158 *>
159 * =====================================================================
160  SUBROUTINE cgebal( JOB, N, A, LDA, ILO, IHI, SCALE, INFO )
161 *
162 * -- LAPACK computational routine --
163 * -- LAPACK is a software package provided by Univ. of Tennessee, --
164 * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
165 *
166 * .. Scalar Arguments ..
167  CHARACTER JOB
168  INTEGER IHI, ILO, INFO, LDA, N
169 * ..
170 * .. Array Arguments ..
171  REAL SCALE( * )
172  COMPLEX A( LDA, * )
173 * ..
174 *
175 * =====================================================================
176 *
177 * .. Parameters ..
178  REAL ZERO, ONE
179  parameter( zero = 0.0e+0, one = 1.0e+0 )
180  REAL SCLFAC
181  parameter( sclfac = 2.0e+0 )
182  REAL FACTOR
183  parameter( factor = 0.95e+0 )
184 * ..
185 * .. Local Scalars ..
186  LOGICAL NOCONV
187  INTEGER I, ICA, IEXC, IRA, J, K, L, M
188  REAL C, CA, F, G, R, RA, S, SFMAX1, SFMAX2, SFMIN1,
189  $ SFMIN2
190 * ..
191 * .. External Functions ..
192  LOGICAL SISNAN, LSAME
193  INTEGER ICAMAX
194  REAL SLAMCH, SCNRM2
195  EXTERNAL sisnan, lsame, icamax, slamch, scnrm2
196 * ..
197 * .. External Subroutines ..
198  EXTERNAL csscal, cswap, xerbla
199 * ..
200 * .. Intrinsic Functions ..
201  INTRINSIC abs, aimag, max, min, real
202 *
203 * Test the input parameters
204 *
205  info = 0
206  IF( .NOT.lsame( job, 'N' ) .AND. .NOT.lsame( job, 'P' ) .AND.
207  $ .NOT.lsame( job, 'S' ) .AND. .NOT.lsame( job, 'B' ) ) THEN
208  info = -1
209  ELSE IF( n.LT.0 ) THEN
210  info = -2
211  ELSE IF( lda.LT.max( 1, n ) ) THEN
212  info = -4
213  END IF
214  IF( info.NE.0 ) THEN
215  CALL xerbla( 'CGEBAL', -info )
216  RETURN
217  END IF
218 *
219  k = 1
220  l = n
221 *
222  IF( n.EQ.0 )
223  $ GO TO 210
224 *
225  IF( lsame( job, 'N' ) ) THEN
226  DO 10 i = 1, n
227  scale( i ) = one
228  10 CONTINUE
229  GO TO 210
230  END IF
231 *
232  IF( lsame( job, 'S' ) )
233  $ GO TO 120
234 *
235 * Permutation to isolate eigenvalues if possible
236 *
237  GO TO 50
238 *
239 * Row and column exchange.
240 *
241  20 CONTINUE
242  scale( m ) = j
243  IF( j.EQ.m )
244  $ GO TO 30
245 *
246  CALL cswap( l, a( 1, j ), 1, a( 1, m ), 1 )
247  CALL cswap( n-k+1, a( j, k ), lda, a( m, k ), lda )
248 *
249  30 CONTINUE
250  GO TO ( 40, 80 )iexc
251 *
252 * Search for rows isolating an eigenvalue and push them down.
253 *
254  40 CONTINUE
255  IF( l.EQ.1 )
256  $ GO TO 210
257  l = l - 1
258 *
259  50 CONTINUE
260  DO 70 j = l, 1, -1
261 *
262  DO 60 i = 1, l
263  IF( i.EQ.j )
264  $ GO TO 60
265  IF( real( a( j, i ) ).NE.zero .OR. aimag( a( j, i ) ).NE.
266  $ zero )GO TO 70
267  60 CONTINUE
268 *
269  m = l
270  iexc = 1
271  GO TO 20
272  70 CONTINUE
273 *
274  GO TO 90
275 *
276 * Search for columns isolating an eigenvalue and push them left.
277 *
278  80 CONTINUE
279  k = k + 1
280 *
281  90 CONTINUE
282  DO 110 j = k, l
283 *
284  DO 100 i = k, l
285  IF( i.EQ.j )
286  $ GO TO 100
287  IF( real( a( i, j ) ).NE.zero .OR. aimag( a( i, j ) ).NE.
288  $ zero )GO TO 110
289  100 CONTINUE
290 *
291  m = k
292  iexc = 2
293  GO TO 20
294  110 CONTINUE
295 *
296  120 CONTINUE
297  DO 130 i = k, l
298  scale( i ) = one
299  130 CONTINUE
300 *
301  IF( lsame( job, 'P' ) )
302  $ GO TO 210
303 *
304 * Balance the submatrix in rows K to L.
305 *
306 * Iterative loop for norm reduction
307 *
308  sfmin1 = slamch( 'S' ) / slamch( 'P' )
309  sfmax1 = one / sfmin1
310  sfmin2 = sfmin1*sclfac
311  sfmax2 = one / sfmin2
312  140 CONTINUE
313  noconv = .false.
314 *
315  DO 200 i = k, l
316 *
317  c = scnrm2( l-k+1, a( k, i ), 1 )
318  r = scnrm2( l-k+1, a( i , k ), lda )
319  ica = icamax( l, a( 1, i ), 1 )
320  ca = abs( a( ica, i ) )
321  ira = icamax( n-k+1, a( i, k ), lda )
322  ra = abs( a( i, ira+k-1 ) )
323 *
324 * Guard against zero C or R due to underflow.
325 *
326  IF( c.EQ.zero .OR. r.EQ.zero )
327  $ GO TO 200
328  g = r / sclfac
329  f = one
330  s = c + r
331  160 CONTINUE
332  IF( c.GE.g .OR. max( f, c, ca ).GE.sfmax2 .OR.
333  $ min( r, g, ra ).LE.sfmin2 )GO TO 170
334  IF( sisnan( c+f+ca+r+g+ra ) ) THEN
335 *
336 * Exit if NaN to avoid infinite loop
337 *
338  info = -3
339  CALL xerbla( 'CGEBAL', -info )
340  RETURN
341  END IF
342  f = f*sclfac
343  c = c*sclfac
344  ca = ca*sclfac
345  r = r / sclfac
346  g = g / sclfac
347  ra = ra / sclfac
348  GO TO 160
349 *
350  170 CONTINUE
351  g = c / sclfac
352  180 CONTINUE
353  IF( g.LT.r .OR. max( r, ra ).GE.sfmax2 .OR.
354  $ min( f, c, g, ca ).LE.sfmin2 )GO TO 190
355  f = f / sclfac
356  c = c / sclfac
357  g = g / sclfac
358  ca = ca / sclfac
359  r = r*sclfac
360  ra = ra*sclfac
361  GO TO 180
362 *
363 * Now balance.
364 *
365  190 CONTINUE
366  IF( ( c+r ).GE.factor*s )
367  $ GO TO 200
368  IF( f.LT.one .AND. scale( i ).LT.one ) THEN
369  IF( f*scale( i ).LE.sfmin1 )
370  $ GO TO 200
371  END IF
372  IF( f.GT.one .AND. scale( i ).GT.one ) THEN
373  IF( scale( i ).GE.sfmax1 / f )
374  $ GO TO 200
375  END IF
376  g = one / f
377  scale( i ) = scale( i )*f
378  noconv = .true.
379 *
380  CALL csscal( n-k+1, g, a( i, k ), lda )
381  CALL csscal( l, f, a( 1, i ), 1 )
382 *
383  200 CONTINUE
384 *
385  IF( noconv )
386  $ GO TO 140
387 *
388  210 CONTINUE
389  ilo = k
390  ihi = l
391 *
392  RETURN
393 *
394 * End of CGEBAL
395 *
396  END
subroutine xerbla(SRNAME, INFO)
XERBLA
Definition: xerbla.f:60
subroutine csscal(N, SA, CX, INCX)
CSSCAL
Definition: csscal.f:78
subroutine cswap(N, CX, INCX, CY, INCY)
CSWAP
Definition: cswap.f:81
subroutine cgebal(JOB, N, A, LDA, ILO, IHI, SCALE, INFO)
CGEBAL
Definition: cgebal.f:161