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