Skip to content

Commit 38f7703

Browse files
langoulangou
authored
Merge pull request #808 from eprovst/gebal
Refactor xGEBAL
2 parents 7ac31fa + 9c489fd commit 38f7703

File tree

4 files changed

+654
-574
lines changed

4 files changed

+654
-574
lines changed

‎SRC/cgebal.f

Lines changed: 167 additions & 145 deletions
Original file line numberDiff line numberDiff line change
@@ -85,6 +85,7 @@
8585
*> \verbatim
8686
*> ILO is INTEGER
8787
*> \endverbatim
88+
*>
8889
*> \param[out] IHI
8990
*> \verbatim
9091
*> IHI is INTEGER
@@ -154,6 +155,9 @@
154155
*>
155156
*> Modified by Tzu-Yi Chen, Computer Science Division, University of
156157
*> California at Berkeley, USA
158+
*>
159+
*> Refactored by Evert Provoost, Department of Computer Science,
160+
*> KU Leuven, Belgium
157161
*> \endverbatim
158162
*>
159163
* =====================================================================
@@ -183,8 +187,8 @@ SUBROUTINE CGEBAL( JOB, N, A, LDA, ILO, IHI, SCALE, INFO )
183187
PARAMETER ( FACTOR = 0.95E+0 )
184188
* ..
185189
* .. Local Scalars ..
186-
LOGICAL NOCONV
187-
INTEGER I, ICA, IEXC, IRA, J, K, L, M
190+
LOGICAL NOCONV, CANSWAP
191+
INTEGER I, ICA, IRA, J, K, L
188192
REAL C, CA, F, G, R, RA, S, SFMAX1, SFMAX2, SFMIN1,
189193
$ SFMIN2
190194
* ..
@@ -195,10 +199,10 @@ SUBROUTINE CGEBAL( JOB, N, A, LDA, ILO, IHI, SCALE, INFO )
195199
EXTERNAL SISNAN, LSAME, ICAMAX, SLAMCH, SCNRM2
196200
* ..
197201
* .. External Subroutines ..
198-
EXTERNAL CSSCAL, CSWAP, XERBLA
202+
EXTERNAL XERBLA, CSSCAL, CSWAP
199203
* ..
200204
* .. Intrinsic Functions ..
201-
INTRINSIC ABS, AIMAG, MAX, MIN, REAL
205+
INTRINSIC ABS, REAL, AIMAG, MAX, MIN
202206
*
203207
* Test the input parameters
204208
*
@@ -216,176 +220,194 @@ SUBROUTINE CGEBAL( JOB, N, A, LDA, ILO, IHI, SCALE, INFO )
216220
RETURN
217221
END IF
218222
*
219-
K = 1
220-
L = N
223+
* Quick returns.
221224
*
222-
IF( N.EQ.0 )
223-
$ GO TO 210
225+
IF( N.EQ.0 ) THEN
226+
ILO = 1
227+
IHI = 0
228+
RETURN
229+
END IF
224230
*
225231
IF( LSAME( JOB, 'N' ) ) THEN
226-
DO 10 I = 1, N
232+
DO I = 1, N
227233
SCALE( I ) = ONE
228-
10 CONTINUE
229-
GO TO 210
234+
END DO
235+
ILO = 1
236+
IHI = N
237+
RETURN
230238
END IF
231239
*
232-
IF( LSAME( JOB, 'S' ) )
233-
$ GO TO 120
234-
*
235-
* Permutation to isolate eigenvalues if possible
236-
*
237-
GO TO 50
238-
*
239-
* Row and column exchange.
240-
*
241-
20 CONTINUE
242-
SCALE( M ) = J
243-
IF( J.EQ.M )
244-
$ GO TO 30
245-
*
246-
CALL CSWAP( L, A( 1, J ), 1, A( 1, M ), 1 )
247-
CALL CSWAP( N-K+1, A( J, K ), LDA, A( M, K ), LDA )
248-
*
249-
30 CONTINUE
250-
GO TO ( 40, 80 )IEXC
251-
*
252-
* Search for rows isolating an eigenvalue and push them down.
253-
*
254-
40 CONTINUE
255-
IF( L.EQ.1 )
256-
$ GO TO 210
257-
L = L - 1
258-
*
259-
50 CONTINUE
260-
DO 70 J = L, 1, -1
240+
* Permutation to isolate eigenvalues if possible.
261241
*
262-
DO 60 I = 1, L
263-
IF( I.EQ.J )
264-
$ GO TO 60
265-
IF( REAL( A( J, I ) ).NE.ZERO .OR. AIMAG( A( J, I ) ).NE.
266-
$ ZERO )GO TO 70
267-
60 CONTINUE
268-
*
269-
M = L
270-
IEXC = 1
271-
GO TO 20
272-
70 CONTINUE
273-
*
274-
GO TO 90
242+
K = 1
243+
L = N
275244
*
276-
* Search for columns isolating an eigenvalue and push them left.
245+
IF( .NOT.LSAME( JOB, 'S' ) ) THEN
277246
*
278-
80 CONTINUE
279-
K = K + 1
247+
* Row and column exchange.
280248
*
281-
90 CONTINUE
282-
DO 110 J = K, L
249+
NOCONV = .TRUE.
250+
DO WHILE( NOCONV )
251+
*
252+
* Search for rows isolating an eigenvalue and push them down.
253+
*
254+
NOCONV = .FALSE.
255+
DO I = L, 1, -1
256+
CANSWAP = .TRUE.
257+
DO J = 1, L
258+
IF( I.NE.J .AND. ( REAL( A( I, J ) ).NE.ZERO .OR.
259+
$ AIMAG( A( I, J ) ).NE.ZERO ) ) THEN
260+
CANSWAP = .FALSE.
261+
EXIT
262+
END IF
263+
END DO
264+
*
265+
IF( CANSWAP ) THEN
266+
SCALE( L ) = I
267+
IF( I.NE.L ) THEN
268+
CALL CSWAP( L, A( 1, I ), 1, A( 1, L ), 1 )
269+
CALL CSWAP( N-K+1, A( I, K ), LDA, A( L, K ), LDA )
270+
END IF
271+
NOCONV = .TRUE.
272+
*
273+
IF( L.EQ.1 ) THEN
274+
ILO = 1
275+
IHI = 1
276+
RETURN
277+
END IF
278+
*
279+
L = L - 1
280+
END IF
281+
END DO
282+
*
283+
END DO
284+
285+
NOCONV = .TRUE.
286+
DO WHILE( NOCONV )
287+
*
288+
* Search for columns isolating an eigenvalue and push them left.
289+
*
290+
NOCONV = .FALSE.
291+
DO J = K, L
292+
CANSWAP = .TRUE.
293+
DO I = K, L
294+
IF( I.NE.J .AND. ( REAL( A( I, J ) ).NE.ZERO .OR.
295+
$ AIMAG( A( I, J ) ).NE.ZERO ) ) THEN
296+
CANSWAP = .FALSE.
297+
EXIT
298+
END IF
299+
END DO
300+
*
301+
IF( CANSWAP ) THEN
302+
SCALE( K ) = J
303+
IF( J.NE.K ) THEN
304+
CALL CSWAP( L, A( 1, J ), 1, A( 1, K ), 1 )
305+
CALL CSWAP( N-K+1, A( J, K ), LDA, A( K, K ), LDA )
306+
END IF
307+
NOCONV = .TRUE.
308+
*
309+
K = K + 1
310+
END IF
311+
END DO
312+
*
313+
END DO
283314
*
284-
DO 100 I = K, L
285-
IF( I.EQ.J )
286-
$ GO TO 100
287-
IF( REAL( A( I, J ) ).NE.ZERO .OR. AIMAG( A( I, J ) ).NE.
288-
$ ZERO )GO TO 110
289-
100 CONTINUE
315+
END IF
290316
*
291-
M = K
292-
IEXC = 2
293-
GO TO 20
294-
110 CONTINUE
317+
* Initialize SCALE for non-permuted submatrix.
295318
*
296-
120 CONTINUE
297-
DO 130 I = K, L
319+
DO I = K, L
298320
SCALE( I ) = ONE
299-
130 CONTINUE
321+
END DO
300322
*
301-
IF( LSAME( JOB, 'P' ) )
302-
$ GO TO 210
323+
* If we only had to permute, we are done.
324+
*
325+
IF( LSAME( JOB, 'P' ) ) THEN
326+
ILO = K
327+
IHI = L
328+
RETURN
329+
END IF
303330
*
304331
* Balance the submatrix in rows K to L.
305332
*
306-
* Iterative loop for norm reduction
333+
* Iterative loop for norm reduction.
307334
*
308335
SFMIN1 = SLAMCH( 'S' ) / SLAMCH( 'P' )
309336
SFMAX1 = ONE / SFMIN1
310337
SFMIN2 = SFMIN1*SCLFAC
311338
SFMAX2 = ONE / SFMIN2
312-
140 CONTINUE
313-
NOCONV = .FALSE.
314-
*
315-
DO 200 I = K, L
316-
*
317-
C = SCNRM2( L-K+1, A( K, I ), 1 )
318-
R = SCNRM2( L-K+1, A( I , K ), LDA )
319-
ICA = ICAMAX( L, A( 1, I ), 1 )
320-
CA = ABS( A( ICA, I ) )
321-
IRA = ICAMAX( N-K+1, A( I, K ), LDA )
322-
RA = ABS( A( I, IRA+K-1 ) )
323-
*
324-
* Guard against zero C or R due to underflow.
325-
*
326-
IF( C.EQ.ZERO .OR. R.EQ.ZERO )
327-
$ GO TO 200
328-
G = R / SCLFAC
329-
F = ONE
330-
S = C + R
331-
160 CONTINUE
332-
IF( C.GE.G .OR. MAX( F, C, CA ).GE.SFMAX2 .OR.
333-
$ MIN( R, G, RA ).LE.SFMIN2 )GO TO 170
334-
IF( SISNAN( C+F+CA+R+G+RA ) ) THEN
335339
*
336-
* Exit if NaN to avoid infinite loop
340+
NOCONV = .TRUE.
341+
DO WHILE( NOCONV )
342+
NOCONV = .FALSE.
337343
*
338-
INFO = -3
339-
CALL XERBLA( 'CGEBAL', -INFO )
340-
RETURN
341-
END IF
342-
F = F*SCLFAC
343-
C = C*SCLFAC
344-
CA = CA*SCLFAC
345-
R = R / SCLFAC
346-
G = G / SCLFAC
347-
RA = RA / SCLFAC
348-
GO TO 160
349-
*
350-
170 CONTINUE
351-
G = C / SCLFAC
352-
180 CONTINUE
353-
IF( G.LT.R .OR. MAX( R, RA ).GE.SFMAX2 .OR.
354-
$ MIN( F, C, G, CA ).LE.SFMIN2 )GO TO 190
355-
F = F / SCLFAC
356-
C = C / SCLFAC
357-
G = G / SCLFAC
358-
CA = CA / SCLFAC
359-
R = R*SCLFAC
360-
RA = RA*SCLFAC
361-
GO TO 180
362-
*
363-
* Now balance.
364-
*
365-
190 CONTINUE
366-
IF( ( C+R ).GE.FACTOR*S )
367-
$ GO TO 200
368-
IF( F.LT.ONE .AND. SCALE( I ).LT.ONE ) THEN
369-
IF( F*SCALE( I ).LE.SFMIN1 )
370-
$ GO TO 200
371-
END IF
372-
IF( F.GT.ONE .AND. SCALE( I ).GT.ONE ) THEN
373-
IF( SCALE( I ).GE.SFMAX1 / F )
374-
$ GO TO 200
375-
END IF
376-
G = ONE / F
377-
SCALE( I ) = SCALE( I )*F
378-
NOCONV = .TRUE.
344+
DO I = K, L
379345
*
380-
CALL CSSCAL( N-K+1, G, A( I, K ), LDA )
381-
CALL CSSCAL( L, F, A( 1, I ), 1 )
346+
C = SCNRM2( L-K+1, A( K, I ), 1 )
347+
R = SCNRM2( L-K+1, A( I, K ), LDA )
348+
ICA = ICAMAX( L, A( 1, I ), 1 )
349+
CA = ABS( A( ICA, I ) )
350+
IRA = ICAMAX( N-K+1, A( I, K ), LDA )
351+
RA = ABS( A( I, IRA+K-1 ) )
382352
*
383-
200 CONTINUE
353+
* Guard against zero C or R due to underflow.
354+
*
355+
IF( C.EQ.ZERO .OR. R.EQ.ZERO ) CYCLE
356+
*
357+
* Exit if NaN to avoid infinite loop
384358
*
385-
IF( NOCONV )
386-
$ GO TO 140
359+
IF( SISNAN( C+CA+R+RA ) ) THEN
360+
INFO = -3
361+
CALL XERBLA( 'CGEBAL', -INFO )
362+
RETURN
363+
END IF
364+
*
365+
G = R / SCLFAC
366+
F = ONE
367+
S = C + R
368+
*
369+
DO WHILE( C.LT.G .AND. MAX( F, C, CA ).LT.SFMAX2 .AND.
370+
$ MIN( R, G, RA ).GT.SFMIN2 )
371+
F = F*SCLFAC
372+
C = C*SCLFAC
373+
CA = CA*SCLFAC
374+
R = R / SCLFAC
375+
G = G / SCLFAC
376+
RA = RA / SCLFAC
377+
END DO
378+
*
379+
G = C / SCLFAC
380+
*
381+
DO WHILE( G.GE.R .AND. MAX( R, RA ).LT.SFMAX2 .AND.
382+
$ MIN( F, C, G, CA ).GT.SFMIN2 )
383+
F = F / SCLFAC
384+
C = C / SCLFAC
385+
G = G / SCLFAC
386+
CA = CA / SCLFAC
387+
R = R*SCLFAC
388+
RA = RA*SCLFAC
389+
END DO
390+
*
391+
* Now balance.
392+
*
393+
IF( ( C+R ).GE.FACTOR*S ) CYCLE
394+
IF( F.LT.ONE .AND. SCALE( I ).LT.ONE ) THEN
395+
IF( F*SCALE( I ).LE.SFMIN1 ) CYCLE
396+
END IF
397+
IF( F.GT.ONE .AND. SCALE( I ).GT.ONE ) THEN
398+
IF( SCALE( I ).GE.SFMAX1 / F ) CYCLE
399+
END IF
400+
G = ONE / F
401+
SCALE( I ) = SCALE( I )*F
402+
NOCONV = .TRUE.
403+
*
404+
CALL CSSCAL( N-K+1, G, A( I, K ), LDA )
405+
CALL CSSCAL( L, F, A( 1, I ), 1 )
406+
*
407+
END DO
408+
*
409+
END DO
387410
*
388-
210 CONTINUE
389411
ILO = K
390412
IHI = L
391413
*

0 commit comments

Comments
 (0)