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
395