KTH framework for Nek5000 toolboxes; testing version  0.0.1
dgebal.f
Go to the documentation of this file.
1  SUBROUTINE dgebal( JOB, N, A, LDA, ILO, IHI, SCALE, INFO )
2 *
3 * -- LAPACK routine (version 3.0) --
4 * Univ. of Tennessee, Univ. of California Berkeley, NAG Ltd.,
5 * Courant Institute, Argonne National Lab, and Rice University
6 * June 30, 1999
7 *
8 * .. Scalar Arguments ..
9  CHARACTER JOB
10  INTEGER IHI, ILO, INFO, LDA, N
11 * ..
12 * .. Array Arguments ..
13  DOUBLE PRECISION A( LDA, * ), SCALE( * )
14 * ..
15 *
16 * Purpose
17 * =======
18 *
19 * DGEBAL balances a general real matrix A. This involves, first,
20 * permuting A by a similarity transformation to isolate eigenvalues
21 * in the first 1 to ILO-1 and last IHI+1 to N elements on the
22 * diagonal; and second, applying a diagonal similarity transformation
23 * to rows and columns ILO to IHI to make the rows and columns as
24 * close in norm as possible. Both steps are optional.
25 *
26 * Balancing may reduce the 1-norm of the matrix, and improve the
27 * accuracy of the computed eigenvalues and/or eigenvectors.
28 *
29 * Arguments
30 * =========
31 *
32 * JOB (input) CHARACTER*1
33 * Specifies the operations to be performed on A:
34 * = 'N': none: simply set ILO = 1, IHI = N, SCALE(I) = 1.0
35 * for i = 1,...,N;
36 * = 'P': permute only;
37 * = 'S': scale only;
38 * = 'B': both permute and scale.
39 *
40 * N (input) INTEGER
41 * The order of the matrix A. N >= 0.
42 *
43 * A (input/output) DOUBLE PRECISION array, dimension (LDA,N)
44 * On entry, the input matrix A.
45 * On exit, A is overwritten by the balanced matrix.
46 * If JOB = 'N', A is not referenced.
47 * See Further Details.
48 *
49 * LDA (input) INTEGER
50 * The leading dimension of the array A. LDA >= max(1,N).
51 *
52 * ILO (output) INTEGER
53 * IHI (output) INTEGER
54 * ILO and IHI are set to integers such that on exit
55 * A(i,j) = 0 if i > j and j = 1,...,ILO-1 or I = IHI+1,...,N.
56 * If JOB = 'N' or 'S', ILO = 1 and IHI = N.
57 *
58 * SCALE (output) DOUBLE PRECISION array, dimension (N)
59 * Details of the permutations and scaling factors applied to
60 * A. If P(j) is the index of the row and column interchanged
61 * with row and column j and D(j) is the scaling factor
62 * applied to row and column j, then
63 * SCALE(j) = P(j) for j = 1,...,ILO-1
64 * = D(j) for j = ILO,...,IHI
65 * = P(j) for j = IHI+1,...,N.
66 * The order in which the interchanges are made is N to IHI+1,
67 * then 1 to ILO-1.
68 *
69 * INFO (output) INTEGER
70 * = 0: successful exit.
71 * < 0: if INFO = -i, the i-th argument had an illegal value.
72 *
73 * Further Details
74 * ===============
75 *
76 * The permutations consist of row and column interchanges which put
77 * the matrix in the form
78 *
79 * ( T1 X Y )
80 * P A P = ( 0 B Z )
81 * ( 0 0 T2 )
82 *
83 * where T1 and T2 are upper triangular matrices whose eigenvalues lie
84 * along the diagonal. The column indices ILO and IHI mark the starting
85 * and ending columns of the submatrix B. Balancing consists of applying
86 * a diagonal similarity transformation inv(D) * B * D to make the
87 * 1-norms of each row of B and its corresponding column nearly equal.
88 * The output matrix is
89 *
90 * ( T1 X*D Y )
91 * ( 0 inv(D)*B*D inv(D)*Z ).
92 * ( 0 0 T2 )
93 *
94 * Information about the permutations P and the diagonal matrix D is
95 * returned in the vector SCALE.
96 *
97 * This subroutine is based on the EISPACK routine BALANC.
98 *
99 * 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 = 0.8d+1 )
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
subroutine dgebal(JOB, N, A, LDA, ILO, IHI, SCALE, INFO)
Definition: dgebal.f:2
subroutine dscal(n, da, dx, incx)
Definition: dscal.f:2
subroutine dswap(n, dx, incx, dy, incy)
Definition: dswap.f:2
subroutine xerbla(SRNAME, INFO)
Definition: xerbla.f:2