LAPACK 3.12.0 LAPACK: Linear Algebra PACKage
Searching...
No Matches
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
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 gebal
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*>
157*> Refactored by Evert Provoost, Department of Computer Science,
158*> KU Leuven, Belgium
159*> \endverbatim
160*>
161* =====================================================================
162 SUBROUTINE sgebal( JOB, N, A, LDA, ILO, IHI, SCALE, INFO )
163*
164* -- LAPACK computational routine --
165* -- LAPACK is a software package provided by Univ. of Tennessee, --
166* -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
167*
168* .. Scalar Arguments ..
169 CHARACTER JOB
170 INTEGER IHI, ILO, INFO, LDA, N
171* ..
172* .. Array Arguments ..
173 REAL A( LDA, * ), SCALE( * )
174* ..
175*
176* =====================================================================
177*
178* .. Parameters ..
179 REAL ZERO, ONE
180 parameter( zero = 0.0e+0, one = 1.0e+0 )
181 REAL SCLFAC
182 parameter( sclfac = 2.0e+0 )
183 REAL FACTOR
184 parameter( factor = 0.95e+0 )
185* ..
186* .. Local Scalars ..
187 LOGICAL NOCONV, CANSWAP
188 INTEGER I, ICA, IRA, J, K, L
189 REAL C, CA, F, G, R, RA, S, SFMAX1, SFMAX2, SFMIN1,
190 \$ SFMIN2
191* ..
192* .. External Functions ..
193 LOGICAL SISNAN, LSAME
194 INTEGER ISAMAX
195 REAL SLAMCH, SNRM2
196 EXTERNAL sisnan, lsame, isamax, slamch, snrm2
197* ..
198* .. External Subroutines ..
199 EXTERNAL sscal, sswap, xerbla
200* ..
201* .. Intrinsic Functions ..
202 INTRINSIC abs, max, min
203* ..
204* Test the input parameters
205*
206 info = 0
207 IF( .NOT.lsame( job, 'N' ) .AND. .NOT.lsame( job, 'P' ) .AND.
208 \$ .NOT.lsame( job, 'S' ) .AND. .NOT.lsame( job, 'B' ) ) THEN
209 info = -1
210 ELSE IF( n.LT.0 ) THEN
211 info = -2
212 ELSE IF( lda.LT.max( 1, n ) ) THEN
213 info = -4
214 END IF
215 IF( info.NE.0 ) THEN
216 CALL xerbla( 'SGEBAL', -info )
217 RETURN
218 END IF
219*
220* Quick returns.
221*
222 IF( n.EQ.0 ) THEN
223 ilo = 1
224 ihi = 0
225 RETURN
226 END IF
227*
228 IF( lsame( job, 'N' ) ) THEN
229 DO i = 1, n
230 scale( i ) = one
231 END DO
232 ilo = 1
233 ihi = n
234 RETURN
235 END IF
236*
237* Permutation to isolate eigenvalues if possible.
238*
239 k = 1
240 l = n
241*
242 IF( .NOT.lsame( job, 'S' ) ) THEN
243*
244* Row and column exchange.
245*
246 noconv = .true.
247 DO WHILE( noconv )
248*
249* Search for rows isolating an eigenvalue and push them down.
250*
251 noconv = .false.
252 DO i = l, 1, -1
253 canswap = .true.
254 DO j = 1, l
255 IF( i.NE.j .AND. a( i, j ).NE.zero ) THEN
256 canswap = .false.
257 EXIT
258 END IF
259 END DO
260*
261 IF( canswap ) THEN
262 scale( l ) = i
263 IF( i.NE.l ) THEN
264 CALL sswap( l, a( 1, i ), 1, a( 1, l ), 1 )
265 CALL sswap( n-k+1, a( i, k ), lda, a( l, k ), lda )
266 END IF
267 noconv = .true.
268*
269 IF( l.EQ.1 ) THEN
270 ilo = 1
271 ihi = 1
272 RETURN
273 END IF
274*
275 l = l - 1
276 END IF
277 END DO
278*
279 END DO
280
281 noconv = .true.
282 DO WHILE( noconv )
283*
284* Search for columns isolating an eigenvalue and push them left.
285*
286 noconv = .false.
287 DO j = k, l
288 canswap = .true.
289 DO i = k, l
290 IF( i.NE.j .AND. a( i, j ).NE.zero ) THEN
291 canswap = .false.
292 EXIT
293 END IF
294 END DO
295*
296 IF( canswap ) THEN
297 scale( k ) = j
298 IF( j.NE.k ) THEN
299 CALL sswap( l, a( 1, j ), 1, a( 1, k ), 1 )
300 CALL sswap( n-k+1, a( j, k ), lda, a( k, k ), lda )
301 END IF
302 noconv = .true.
303*
304 k = k + 1
305 END IF
306 END DO
307*
308 END DO
309*
310 END IF
311*
312* Initialize SCALE for non-permuted submatrix.
313*
314 DO i = k, l
315 scale( i ) = one
316 END DO
317*
318* If we only had to permute, we are done.
319*
320 IF( lsame( job, 'P' ) ) THEN
321 ilo = k
322 ihi = l
323 RETURN
324 END IF
325*
326* Balance the submatrix in rows K to L.
327*
328* Iterative loop for norm reduction.
329*
330 sfmin1 = slamch( 'S' ) / slamch( 'P' )
331 sfmax1 = one / sfmin1
332 sfmin2 = sfmin1*sclfac
333 sfmax2 = one / sfmin2
334*
335 noconv = .true.
336 DO WHILE( noconv )
337 noconv = .false.
338*
339 DO i = k, l
340*
341 c = snrm2( l-k+1, a( k, i ), 1 )
342 r = snrm2( l-k+1, a( i, k ), lda )
343 ica = isamax( l, a( 1, i ), 1 )
344 ca = abs( a( ica, i ) )
345 ira = isamax( n-k+1, a( i, k ), lda )
346 ra = abs( a( i, ira+k-1 ) )
347*
348* Guard against zero C or R due to underflow.
349*
350 IF( c.EQ.zero .OR. r.EQ.zero ) cycle
351*
352* Exit if NaN to avoid infinite loop
353*
354 IF( sisnan( c+ca+r+ra ) ) THEN
355 info = -3
356 CALL xerbla( 'SGEBAL', -info )
357 RETURN
358 END IF
359*
360 g = r / sclfac
361 f = one
362 s = c + r
363*
364 DO WHILE( c.LT.g .AND. max( f, c, ca ).LT.sfmax2 .AND.
365 \$ min( r, g, ra ).GT.sfmin2 )
366 f = f*sclfac
367 c = c*sclfac
368 ca = ca*sclfac
369 r = r / sclfac
370 g = g / sclfac
371 ra = ra / sclfac
372 END DO
373*
374 g = c / sclfac
375*
376 DO WHILE( g.GE.r .AND. max( r, ra ).LT.sfmax2 .AND.
377 \$ min( f, c, g, ca ).GT.sfmin2 )
378 f = f / sclfac
379 c = c / sclfac
380 g = g / sclfac
381 ca = ca / sclfac
382 r = r*sclfac
383 ra = ra*sclfac
384 END DO
385*
386* Now balance.
387*
388 IF( ( c+r ).GE.factor*s ) cycle
389 IF( f.LT.one .AND. scale( i ).LT.one ) THEN
390 IF( f*scale( i ).LE.sfmin1 ) cycle
391 END IF
392 IF( f.GT.one .AND. scale( i ).GT.one ) THEN
393 IF( scale( i ).GE.sfmax1 / f ) cycle
394 END IF
395 g = one / f
396 scale( i ) = scale( i )*f
397 noconv = .true.
398*
399 CALL sscal( n-k+1, g, a( i, k ), lda )
400 CALL sscal( l, f, a( 1, i ), 1 )
401*
402 END DO
403*
404 END DO
405*
406 ilo = k
407 ihi = l
408*
409 RETURN
410*
411* End of SGEBAL
412*
413 END
subroutine xerbla(srname, info)
Definition cblat2.f:3285
subroutine sgebal(job, n, a, lda, ilo, ihi, scale, info)
SGEBAL
Definition sgebal.f:163
subroutine sscal(n, sa, sx, incx)
SSCAL
Definition sscal.f:79
subroutine sswap(n, sx, incx, sy, incy)
SSWAP
Definition sswap.f:82