LAPACK 3.11.0
LAPACK: Linear Algebra PACKage
Loading...
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
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