LAPACK  3.4.2
LAPACK: Linear Algebra PACKage
 All Files Functions Groups
zgebal.f
Go to the documentation of this file.
1 *> \brief \b ZGEBAL
2 *
3 * =========== DOCUMENTATION ===========
4 *
5 * Online html documentation available at
6 * http://www.netlib.org/lapack/explore-html/
7 *
8 *> \htmlonly
9 *> Download ZGEBAL + dependencies
10 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/zgebal.f">
11 *> [TGZ]</a>
12 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/zgebal.f">
13 *> [ZIP]</a>
14 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/zgebal.f">
15 *> [TXT]</a>
16 *> \endhtmlonly
17 *
18 * Definition:
19 * ===========
20 *
21 * SUBROUTINE ZGEBAL( 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 * DOUBLE PRECISION SCALE( * )
29 * COMPLEX*16 A( LDA, * )
30 * ..
31 *
32 *
33 *> \par Purpose:
34 * =============
35 *>
36 *> \verbatim
37 *>
38 *> ZGEBAL 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*16 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 *> \endverbatim
87 *>
88 *> \param[out] IHI
89 *> \verbatim
90 *> ILO and IHI are set to INTEGER such that on exit
91 *> A(i,j) = 0 if i > j and j = 1,...,ILO-1 or I = IHI+1,...,N.
92 *> If JOB = 'N' or 'S', ILO = 1 and IHI = N.
93 *> \endverbatim
94 *>
95 *> \param[out] SCALE
96 *> \verbatim
97 *> SCALE is DOUBLE PRECISION array, dimension (N)
98 *> Details of the permutations and scaling factors applied to
99 *> A. If P(j) is the index of the row and column interchanged
100 *> with row and column j and D(j) is the scaling factor
101 *> applied to row and column j, then
102 *> SCALE(j) = P(j) for j = 1,...,ILO-1
103 *> = D(j) for j = ILO,...,IHI
104 *> = P(j) for j = IHI+1,...,N.
105 *> The order in which the interchanges are made is N to IHI+1,
106 *> then 1 to ILO-1.
107 *> \endverbatim
108 *>
109 *> \param[out] INFO
110 *> \verbatim
111 *> INFO is INTEGER
112 *> = 0: successful exit.
113 *> < 0: if INFO = -i, the i-th argument had an illegal value.
114 *> \endverbatim
115 *
116 * Authors:
117 * ========
118 *
119 *> \author Univ. of Tennessee
120 *> \author Univ. of California Berkeley
121 *> \author Univ. of Colorado Denver
122 *> \author NAG Ltd.
123 *
124 *> \date November 2011
125 *
126 *> \ingroup complex16GEcomputational
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 *> \endverbatim
159 *>
160 * =====================================================================
161  SUBROUTINE zgebal( JOB, N, A, LDA, ILO, IHI, SCALE, INFO )
162 *
163 * -- LAPACK computational routine (version 3.4.0) --
164 * -- LAPACK is a software package provided by Univ. of Tennessee, --
165 * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
166 * November 2011
167 *
168 * .. Scalar Arguments ..
169  CHARACTER job
170  INTEGER ihi, ilo, info, lda, n
171 * ..
172 * .. Array Arguments ..
173  DOUBLE PRECISION scale( * )
174  COMPLEX*16 a( lda, * )
175 * ..
176 *
177 * =====================================================================
178 *
179 * .. Parameters ..
180  DOUBLE PRECISION zero, one
181  parameter( zero = 0.0d+0, one = 1.0d+0 )
182  DOUBLE PRECISION sclfac
183  parameter( sclfac = 2.0d+0 )
184  DOUBLE PRECISION factor
185  parameter( factor = 0.95d+0 )
186 * ..
187 * .. Local Scalars ..
188  LOGICAL noconv
189  INTEGER i, ica, iexc, ira, j, k, l, m
190  DOUBLE PRECISION c, ca, f, g, r, ra, s, sfmax1, sfmax2, sfmin1,
191  $ sfmin2
192  COMPLEX*16 cdum
193 * ..
194 * .. External Functions ..
195  LOGICAL disnan, lsame
196  INTEGER izamax
197  DOUBLE PRECISION dlamch
198  EXTERNAL disnan, lsame, izamax, dlamch
199 * ..
200 * .. External Subroutines ..
201  EXTERNAL xerbla, zdscal, zswap
202 * ..
203 * .. Intrinsic Functions ..
204  INTRINSIC abs, dble, dimag, max, min
205 * ..
206 * .. Statement Functions ..
207  DOUBLE PRECISION cabs1
208 * ..
209 * .. Statement Function definitions ..
210  cabs1( cdum ) = abs( dble( cdum ) ) + abs( dimag( cdum ) )
211 * ..
212 * .. Executable Statements ..
213 *
214 * Test the input parameters
215 *
216  info = 0
217  IF( .NOT.lsame( job, 'N' ) .AND. .NOT.lsame( job, 'P' ) .AND.
218  $ .NOT.lsame( job, 'S' ) .AND. .NOT.lsame( job, 'B' ) ) THEN
219  info = -1
220  ELSE IF( n.LT.0 ) THEN
221  info = -2
222  ELSE IF( lda.LT.max( 1, n ) ) THEN
223  info = -4
224  END IF
225  IF( info.NE.0 ) THEN
226  CALL xerbla( 'ZGEBAL', -info )
227  return
228  END IF
229 *
230  k = 1
231  l = n
232 *
233  IF( n.EQ.0 )
234  $ go to 210
235 *
236  IF( lsame( job, 'N' ) ) THEN
237  DO 10 i = 1, n
238  scale( i ) = one
239  10 continue
240  go to 210
241  END IF
242 *
243  IF( lsame( job, 'S' ) )
244  $ go to 120
245 *
246 * Permutation to isolate eigenvalues if possible
247 *
248  go to 50
249 *
250 * Row and column exchange.
251 *
252  20 continue
253  scale( m ) = j
254  IF( j.EQ.m )
255  $ go to 30
256 *
257  CALL zswap( l, a( 1, j ), 1, a( 1, m ), 1 )
258  CALL zswap( n-k+1, a( j, k ), lda, a( m, k ), lda )
259 *
260  30 continue
261  go to( 40, 80 )iexc
262 *
263 * Search for rows isolating an eigenvalue and push them down.
264 *
265  40 continue
266  IF( l.EQ.1 )
267  $ go to 210
268  l = l - 1
269 *
270  50 continue
271  DO 70 j = l, 1, -1
272 *
273  DO 60 i = 1, l
274  IF( i.EQ.j )
275  $ go to 60
276  IF( dble( a( j, i ) ).NE.zero .OR. dimag( a( j, i ) ).NE.
277  $ zero )go to 70
278  60 continue
279 *
280  m = l
281  iexc = 1
282  go to 20
283  70 continue
284 *
285  go to 90
286 *
287 * Search for columns isolating an eigenvalue and push them left.
288 *
289  80 continue
290  k = k + 1
291 *
292  90 continue
293  DO 110 j = k, l
294 *
295  DO 100 i = k, l
296  IF( i.EQ.j )
297  $ go to 100
298  IF( dble( a( i, j ) ).NE.zero .OR. dimag( a( i, j ) ).NE.
299  $ zero )go to 110
300  100 continue
301 *
302  m = k
303  iexc = 2
304  go to 20
305  110 continue
306 *
307  120 continue
308  DO 130 i = k, l
309  scale( i ) = one
310  130 continue
311 *
312  IF( lsame( job, 'P' ) )
313  $ go to 210
314 *
315 * Balance the submatrix in rows K to L.
316 *
317 * Iterative loop for norm reduction
318 *
319  sfmin1 = dlamch( 'S' ) / dlamch( 'P' )
320  sfmax1 = one / sfmin1
321  sfmin2 = sfmin1*sclfac
322  sfmax2 = one / sfmin2
323  140 continue
324  noconv = .false.
325 *
326  DO 200 i = k, l
327  c = zero
328  r = zero
329 *
330  DO 150 j = k, l
331  IF( j.EQ.i )
332  $ go to 150
333  c = c + cabs1( a( j, i ) )
334  r = r + cabs1( a( i, j ) )
335  150 continue
336  ica = izamax( l, a( 1, i ), 1 )
337  ca = abs( a( ica, i ) )
338  ira = izamax( n-k+1, a( i, k ), lda )
339  ra = abs( a( i, ira+k-1 ) )
340 *
341 * Guard against zero C or R due to underflow.
342 *
343  IF( c.EQ.zero .OR. r.EQ.zero )
344  $ go to 200
345  g = r / sclfac
346  f = one
347  s = c + r
348  160 continue
349  IF( c.GE.g .OR. max( f, c, ca ).GE.sfmax2 .OR.
350  $ min( r, g, ra ).LE.sfmin2 )go to 170
351  IF( disnan( c+f+ca+r+g+ra ) ) THEN
352 *
353 * Exit if NaN to avoid infinite loop
354 *
355  info = -3
356  CALL xerbla( 'ZGEBAL', -info )
357  return
358  END IF
359  f = f*sclfac
360  c = c*sclfac
361  ca = ca*sclfac
362  r = r / sclfac
363  g = g / sclfac
364  ra = ra / sclfac
365  go to 160
366 *
367  170 continue
368  g = c / sclfac
369  180 continue
370  IF( g.LT.r .OR. max( r, ra ).GE.sfmax2 .OR.
371  $ min( f, c, g, ca ).LE.sfmin2 )go to 190
372  f = f / sclfac
373  c = c / sclfac
374  g = g / sclfac
375  ca = ca / sclfac
376  r = r*sclfac
377  ra = ra*sclfac
378  go to 180
379 *
380 * Now balance.
381 *
382  190 continue
383  IF( ( c+r ).GE.factor*s )
384  $ go to 200
385  IF( f.LT.one .AND. scale( i ).LT.one ) THEN
386  IF( f*scale( i ).LE.sfmin1 )
387  $ go to 200
388  END IF
389  IF( f.GT.one .AND. scale( i ).GT.one ) THEN
390  IF( scale( i ).GE.sfmax1 / f )
391  $ go to 200
392  END IF
393  g = one / f
394  scale( i ) = scale( i )*f
395  noconv = .true.
396 *
397  CALL zdscal( n-k+1, g, a( i, k ), lda )
398  CALL zdscal( l, f, a( 1, i ), 1 )
399 *
400  200 continue
401 *
402  IF( noconv )
403  $ go to 140
404 *
405  210 continue
406  ilo = k
407  ihi = l
408 *
409  return
410 *
411 * End of ZGEBAL
412 *
413  END