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