001:       SUBROUTINE DGEBAL( JOB, N, A, LDA, ILO, IHI, SCALE, INFO )
002: *
003: *  -- LAPACK routine (version 3.2) --
004: *  -- LAPACK is a software package provided by Univ. of Tennessee,    --
005: *  -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
006: *     November 2006
007: *
008: *     .. Scalar Arguments ..
009:       CHARACTER          JOB
010:       INTEGER            IHI, ILO, INFO, LDA, N
011: *     ..
012: *     .. Array Arguments ..
013:       DOUBLE PRECISION   A( LDA, * ), SCALE( * )
014: *     ..
015: *
016: *  Purpose
017: *  =======
018: *
019: *  DGEBAL balances a general real matrix A.  This involves, first,
020: *  permuting A by a similarity transformation to isolate eigenvalues
021: *  in the first 1 to ILO-1 and last IHI+1 to N elements on the
022: *  diagonal; and second, applying a diagonal similarity transformation
023: *  to rows and columns ILO to IHI to make the rows and columns as
024: *  close in norm as possible.  Both steps are optional.
025: *
026: *  Balancing may reduce the 1-norm of the matrix, and improve the
027: *  accuracy of the computed eigenvalues and/or eigenvectors.
028: *
029: *  Arguments
030: *  =========
031: *
032: *  JOB     (input) CHARACTER*1
033: *          Specifies the operations to be performed on A:
034: *          = 'N':  none:  simply set ILO = 1, IHI = N, SCALE(I) = 1.0
035: *                  for i = 1,...,N;
036: *          = 'P':  permute only;
037: *          = 'S':  scale only;
038: *          = 'B':  both permute and scale.
039: *
040: *  N       (input) INTEGER
041: *          The order of the matrix A.  N >= 0.
042: *
043: *  A       (input/output) DOUBLE PRECISION array, dimension (LDA,N)
044: *          On entry, the input matrix A.
045: *          On exit,  A is overwritten by the balanced matrix.
046: *          If JOB = 'N', A is not referenced.
047: *          See Further Details.
048: *
049: *  LDA     (input) INTEGER
050: *          The leading dimension of the array A.  LDA >= max(1,N).
051: *
052: *  ILO     (output) INTEGER
053: *  IHI     (output) INTEGER
054: *          ILO and IHI are set to integers such that on exit
055: *          A(i,j) = 0 if i > j and j = 1,...,ILO-1 or I = IHI+1,...,N.
056: *          If JOB = 'N' or 'S', ILO = 1 and IHI = N.
057: *
058: *  SCALE   (output) DOUBLE PRECISION array, dimension (N)
059: *          Details of the permutations and scaling factors applied to
060: *          A.  If P(j) is the index of the row and column interchanged
061: *          with row and column j and D(j) is the scaling factor
062: *          applied to row and column j, then
063: *          SCALE(j) = P(j)    for j = 1,...,ILO-1
064: *                   = D(j)    for j = ILO,...,IHI
065: *                   = P(j)    for j = IHI+1,...,N.
066: *          The order in which the interchanges are made is N to IHI+1,
067: *          then 1 to ILO-1.
068: *
069: *  INFO    (output) INTEGER
070: *          = 0:  successful exit.
071: *          < 0:  if INFO = -i, the i-th argument had an illegal value.
072: *
073: *  Further Details
074: *  ===============
075: *
076: *  The permutations consist of row and column interchanges which put
077: *  the matrix in the form
078: *
079: *             ( T1   X   Y  )
080: *     P A P = (  0   B   Z  )
081: *             (  0   0   T2 )
082: *
083: *  where T1 and T2 are upper triangular matrices whose eigenvalues lie
084: *  along the diagonal.  The column indices ILO and IHI mark the starting
085: *  and ending columns of the submatrix B. Balancing consists of applying
086: *  a diagonal similarity transformation inv(D) * B * D to make the
087: *  1-norms of each row of B and its corresponding column nearly equal.
088: *  The output matrix is
089: *
090: *     ( T1     X*D          Y    )
091: *     (  0  inv(D)*B*D  inv(D)*Z ).
092: *     (  0      0           T2   )
093: *
094: *  Information about the permutations P and the diagonal matrix D is
095: *  returned in the vector SCALE.
096: *
097: *  This subroutine is based on the EISPACK routine BALANC.
098: *
099: *  Modified by Tzu-Yi Chen, Computer Science Division, University of
100: *    California at Berkeley, USA
101: *
102: *  =====================================================================
103: *
104: *     .. Parameters ..
105:       DOUBLE PRECISION   ZERO, ONE
106:       PARAMETER          ( ZERO = 0.0D+0, ONE = 1.0D+0 )
107:       DOUBLE PRECISION   SCLFAC
108:       PARAMETER          ( SCLFAC = 2.0D+0 )
109:       DOUBLE PRECISION   FACTOR
110:       PARAMETER          ( FACTOR = 0.95D+0 )
111: *     ..
112: *     .. Local Scalars ..
113:       LOGICAL            NOCONV
114:       INTEGER            I, ICA, IEXC, IRA, J, K, L, M
115:       DOUBLE PRECISION   C, CA, F, G, R, RA, S, SFMAX1, SFMAX2, SFMIN1,
116:      $                   SFMIN2
117: *     ..
118: *     .. External Functions ..
119:       LOGICAL            LSAME
120:       INTEGER            IDAMAX
121:       DOUBLE PRECISION   DLAMCH
122:       EXTERNAL           LSAME, IDAMAX, DLAMCH
123: *     ..
124: *     .. External Subroutines ..
125:       EXTERNAL           DSCAL, DSWAP, XERBLA
126: *     ..
127: *     .. Intrinsic Functions ..
128:       INTRINSIC          ABS, MAX, MIN
129: *     ..
130: *     .. Executable Statements ..
131: *
132: *     Test the input parameters
133: *
134:       INFO = 0
135:       IF( .NOT.LSAME( JOB, 'N' ) .AND. .NOT.LSAME( JOB, 'P' ) .AND.
136:      $    .NOT.LSAME( JOB, 'S' ) .AND. .NOT.LSAME( JOB, 'B' ) ) THEN
137:          INFO = -1
138:       ELSE IF( N.LT.0 ) THEN
139:          INFO = -2
140:       ELSE IF( LDA.LT.MAX( 1, N ) ) THEN
141:          INFO = -4
142:       END IF
143:       IF( INFO.NE.0 ) THEN
144:          CALL XERBLA( 'DGEBAL', -INFO )
145:          RETURN
146:       END IF
147: *
148:       K = 1
149:       L = N
150: *
151:       IF( N.EQ.0 )
152:      $   GO TO 210
153: *
154:       IF( LSAME( JOB, 'N' ) ) THEN
155:          DO 10 I = 1, N
156:             SCALE( I ) = ONE
157:    10    CONTINUE
158:          GO TO 210
159:       END IF
160: *
161:       IF( LSAME( JOB, 'S' ) )
162:      $   GO TO 120
163: *
164: *     Permutation to isolate eigenvalues if possible
165: *
166:       GO TO 50
167: *
168: *     Row and column exchange.
169: *
170:    20 CONTINUE
171:       SCALE( M ) = J
172:       IF( J.EQ.M )
173:      $   GO TO 30
174: *
175:       CALL DSWAP( L, A( 1, J ), 1, A( 1, M ), 1 )
176:       CALL DSWAP( N-K+1, A( J, K ), LDA, A( M, K ), LDA )
177: *
178:    30 CONTINUE
179:       GO TO ( 40, 80 )IEXC
180: *
181: *     Search for rows isolating an eigenvalue and push them down.
182: *
183:    40 CONTINUE
184:       IF( L.EQ.1 )
185:      $   GO TO 210
186:       L = L - 1
187: *
188:    50 CONTINUE
189:       DO 70 J = L, 1, -1
190: *
191:          DO 60 I = 1, L
192:             IF( I.EQ.J )
193:      $         GO TO 60
194:             IF( A( J, I ).NE.ZERO )
195:      $         GO TO 70
196:    60    CONTINUE
197: *
198:          M = L
199:          IEXC = 1
200:          GO TO 20
201:    70 CONTINUE
202: *
203:       GO TO 90
204: *
205: *     Search for columns isolating an eigenvalue and push them left.
206: *
207:    80 CONTINUE
208:       K = K + 1
209: *
210:    90 CONTINUE
211:       DO 110 J = K, L
212: *
213:          DO 100 I = K, L
214:             IF( I.EQ.J )
215:      $         GO TO 100
216:             IF( A( I, J ).NE.ZERO )
217:      $         GO TO 110
218:   100    CONTINUE
219: *
220:          M = K
221:          IEXC = 2
222:          GO TO 20
223:   110 CONTINUE
224: *
225:   120 CONTINUE
226:       DO 130 I = K, L
227:          SCALE( I ) = ONE
228:   130 CONTINUE
229: *
230:       IF( LSAME( JOB, 'P' ) )
231:      $   GO TO 210
232: *
233: *     Balance the submatrix in rows K to L.
234: *
235: *     Iterative loop for norm reduction
236: *
237:       SFMIN1 = DLAMCH( 'S' ) / DLAMCH( 'P' )
238:       SFMAX1 = ONE / SFMIN1
239:       SFMIN2 = SFMIN1*SCLFAC
240:       SFMAX2 = ONE / SFMIN2
241:   140 CONTINUE
242:       NOCONV = .FALSE.
243: *
244:       DO 200 I = K, L
245:          C = ZERO
246:          R = ZERO
247: *
248:          DO 150 J = K, L
249:             IF( J.EQ.I )
250:      $         GO TO 150
251:             C = C + ABS( A( J, I ) )
252:             R = R + ABS( A( I, J ) )
253:   150    CONTINUE
254:          ICA = IDAMAX( L, A( 1, I ), 1 )
255:          CA = ABS( A( ICA, I ) )
256:          IRA = IDAMAX( N-K+1, A( I, K ), LDA )
257:          RA = ABS( A( I, IRA+K-1 ) )
258: *
259: *        Guard against zero C or R due to underflow.
260: *
261:          IF( C.EQ.ZERO .OR. R.EQ.ZERO )
262:      $      GO TO 200
263:          G = R / SCLFAC
264:          F = ONE
265:          S = C + R
266:   160    CONTINUE
267:          IF( C.GE.G .OR. MAX( F, C, CA ).GE.SFMAX2 .OR.
268:      $       MIN( R, G, RA ).LE.SFMIN2 )GO TO 170
269:          F = F*SCLFAC
270:          C = C*SCLFAC
271:          CA = CA*SCLFAC
272:          R = R / SCLFAC
273:          G = G / SCLFAC
274:          RA = RA / SCLFAC
275:          GO TO 160
276: *
277:   170    CONTINUE
278:          G = C / SCLFAC
279:   180    CONTINUE
280:          IF( G.LT.R .OR. MAX( R, RA ).GE.SFMAX2 .OR.
281:      $       MIN( F, C, G, CA ).LE.SFMIN2 )GO TO 190
282:          F = F / SCLFAC
283:          C = C / SCLFAC
284:          G = G / SCLFAC
285:          CA = CA / SCLFAC
286:          R = R*SCLFAC
287:          RA = RA*SCLFAC
288:          GO TO 180
289: *
290: *        Now balance.
291: *
292:   190    CONTINUE
293:          IF( ( C+R ).GE.FACTOR*S )
294:      $      GO TO 200
295:          IF( F.LT.ONE .AND. SCALE( I ).LT.ONE ) THEN
296:             IF( F*SCALE( I ).LE.SFMIN1 )
297:      $         GO TO 200
298:          END IF
299:          IF( F.GT.ONE .AND. SCALE( I ).GT.ONE ) THEN
300:             IF( SCALE( I ).GE.SFMAX1 / F )
301:      $         GO TO 200
302:          END IF
303:          G = ONE / F
304:          SCALE( I ) = SCALE( I )*F
305:          NOCONV = .TRUE.
306: *
307:          CALL DSCAL( N-K+1, G, A( I, K ), LDA )
308:          CALL DSCAL( L, F, A( 1, I ), 1 )
309: *
310:   200 CONTINUE
311: *
312:       IF( NOCONV )
313:      $   GO TO 140
314: *
315:   210 CONTINUE
316:       ILO = K
317:       IHI = L
318: *
319:       RETURN
320: *
321: *     End of DGEBAL
322: *
323:       END
324: