← Back to branch summary

~ubuntu-branches/ubuntu/hoary/scilab/hoary

« back to all changes in this revision

Viewing changes to routines/lapack/zgebal.f

  • Committer: Bazaar Package Importer
  • Author(s): Torsten Werner
  • Date: 2005-01-09 22:58:21 UTC
  • mfrom: (1.1.1 upstream)
  • Revision ID: james.westby@ubuntu.com-20050109225821-473xr8vhgugxxx5j
Tags: 3.0-12
http://bugs.debian.org/255869
changed configure.in to build scilab's own malloc.o, closes: #255869

Show diffs side-by-side

added added

removed removed

Lines of Context:
routines/lapack/zgebal.f
 
1
      SUBROUTINE ZGEBAL( 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   SCALE( * )
 
14
      COMPLEX*16         A( LDA, * )
 
15
*     ..
 
16
*
 
17
*  Purpose
 
18
*  =======
 
19
*
 
20
*  ZGEBAL balances a general complex matrix A.  This involves, first,
 
21
*  permuting A by a similarity transformation to isolate eigenvalues
 
22
*  in the first 1 to ILO-1 and last IHI+1 to N elements on the
 
23
*  diagonal; and second, applying a diagonal similarity transformation
 
24
*  to rows and columns ILO to IHI to make the rows and columns as
 
25
*  close in norm as possible.  Both steps are optional.
 
26
*
 
27
*  Balancing may reduce the 1-norm of the matrix, and improve the
 
28
*  accuracy of the computed eigenvalues and/or eigenvectors.
 
29
*
 
30
*  Arguments
 
31
*  =========
 
32
*
 
33
*  JOB     (input) CHARACTER*1
 
34
*          Specifies the operations to be performed on A:
 
35
*          = 'N':  none:  simply set ILO = 1, IHI = N, SCALE(I) = 1.0
 
36
*                  for i = 1,...,N;
 
37
*          = 'P':  permute only;
 
38
*          = 'S':  scale only;
 
39
*          = 'B':  both permute and scale.
 
40
*
 
41
*  N       (input) INTEGER
 
42
*          The order of the matrix A.  N >= 0.
 
43
*
 
44
*  A       (input/output) COMPLEX*16 array, dimension (LDA,N)
 
45
*          On entry, the input matrix A.
 
46
*          On exit,  A is overwritten by the balanced matrix.
 
47
*          If JOB = 'N', A is not referenced.
 
48
*          See Further Details.
 
49
*
 
50
*  LDA     (input) INTEGER
 
51
*          The leading dimension of the array A.  LDA >= max(1,N).
 
52
*
 
53
*  ILO     (output) INTEGER
 
54
*  IHI     (output) INTEGER
 
55
*          ILO and IHI are set to integers such that on exit
 
56
*          A(i,j) = 0 if i > j and j = 1,...,ILO-1 or I = IHI+1,...,N.
 
57
*          If JOB = 'N' or 'S', ILO = 1 and IHI = N.
 
58
*
 
59
*  SCALE   (output) DOUBLE PRECISION array, dimension (N)
 
60
*          Details of the permutations and scaling factors applied to
 
61
*          A.  If P(j) is the index of the row and column interchanged
 
62
*          with row and column j and D(j) is the scaling factor
 
63
*          applied to row and column j, then
 
64
*          SCALE(j) = P(j)    for j = 1,...,ILO-1
 
65
*                   = D(j)    for j = ILO,...,IHI
 
66
*                   = P(j)    for j = IHI+1,...,N.
 
67
*          The order in which the interchanges are made is N to IHI+1,
 
68
*          then 1 to ILO-1.
 
69
*
 
70
*  INFO    (output) INTEGER
 
71
*          = 0:  successful exit.
 
72
*          < 0:  if INFO = -i, the i-th argument had an illegal value.
 
73
*
 
74
*  Further Details
 
75
*  ===============
 
76
*
 
77
*  The permutations consist of row and column interchanges which put
 
78
*  the matrix in the form
 
79
*
 
80
*             ( T1   X   Y  )
 
81
*     P A P = (  0   B   Z  )
 
82
*             (  0   0   T2 )
 
83
*
 
84
*  where T1 and T2 are upper triangular matrices whose eigenvalues lie
 
85
*  along the diagonal.  The column indices ILO and IHI mark the starting
 
86
*  and ending columns of the submatrix B. Balancing consists of applying
 
87
*  a diagonal similarity transformation inv(D) * B * D to make the
 
88
*  1-norms of each row of B and its corresponding column nearly equal.
 
89
*  The output matrix is
 
90
*
 
91
*     ( T1     X*D          Y    )
 
92
*     (  0  inv(D)*B*D  inv(D)*Z ).
 
93
*     (  0      0           T2   )
 
94
*
 
95
*  Information about the permutations P and the diagonal matrix D is
 
96
*  returned in the vector SCALE.
 
97
*
 
98
*  This subroutine is based on the EISPACK routine CBAL.
 
99
*
 
100
*  Modified by Tzu-Yi Chen, Computer Science Division, University of
 
101
*    California at Berkeley, USA
 
102
*
 
103
*  =====================================================================
 
104
*
 
105
*     .. Parameters ..
 
106
      DOUBLE PRECISION   ZERO, ONE
 
107
      PARAMETER          ( ZERO = 0.0D+0, ONE = 1.0D+0 )
 
108
      DOUBLE PRECISION   SCLFAC
 
109
      PARAMETER          ( SCLFAC = 0.8D+1 )
 
110
      DOUBLE PRECISION   FACTOR
 
111
      PARAMETER          ( FACTOR = 0.95D+0 )
 
112
*     ..
 
113
*     .. Local Scalars ..
 
114
      LOGICAL            NOCONV
 
115
      INTEGER            I, ICA, IEXC, IRA, J, K, L, M
 
116
      DOUBLE PRECISION   C, CA, F, G, R, RA, S, SFMAX1, SFMAX2, SFMIN1,
 
117
     $                   SFMIN2
 
118
      COMPLEX*16         CDUM
 
119
*     ..
 
120
*     .. External Functions ..
 
121
      LOGICAL            LSAME
 
122
      INTEGER            IZAMAX
 
123
      DOUBLE PRECISION   DLAMCH
 
124
      EXTERNAL           LSAME, IZAMAX, DLAMCH
 
125
*     ..
 
126
*     .. External Subroutines ..
 
127
      EXTERNAL           XERBLA, ZDSCAL, ZSWAP
 
128
*     ..
 
129
*     .. Intrinsic Functions ..
 
130
      INTRINSIC          ABS, DBLE, DIMAG, MAX, MIN
 
131
*     ..
 
132
*     .. Statement Functions ..
 
133
      DOUBLE PRECISION   CABS1
 
134
*     ..
 
135
*     .. Statement Function definitions ..
 
136
      CABS1( CDUM ) = ABS( DBLE( CDUM ) ) + ABS( DIMAG( CDUM ) )
 
137
*     ..
 
138
*     .. Executable Statements ..
 
139
*
 
140
*     Test the input parameters
 
141
*
 
142
      INFO = 0
 
143
      IF( .NOT.LSAME( JOB, 'N' ) .AND. .NOT.LSAME( JOB, 'P' ) .AND.
 
144
     $    .NOT.LSAME( JOB, 'S' ) .AND. .NOT.LSAME( JOB, 'B' ) ) THEN
 
145
         INFO = -1
 
146
      ELSE IF( N.LT.0 ) THEN
 
147
         INFO = -2
 
148
      ELSE IF( LDA.LT.MAX( 1, N ) ) THEN
 
149
         INFO = -4
 
150
      END IF
 
151
      IF( INFO.NE.0 ) THEN
 
152
         CALL XERBLA( 'ZGEBAL', -INFO )
 
153
         RETURN
 
154
      END IF
 
155
*
 
156
      K = 1
 
157
      L = N
 
158
*
 
159
      IF( N.EQ.0 )
 
160
     $   GO TO 210
 
161
*
 
162
      IF( LSAME( JOB, 'N' ) ) THEN
 
163
         DO 10 I = 1, N
 
164
            SCALE( I ) = ONE
 
165
   10    CONTINUE
 
166
         GO TO 210
 
167
      END IF
 
168
*
 
169
      IF( LSAME( JOB, 'S' ) )
 
170
     $   GO TO 120
 
171
*
 
172
*     Permutation to isolate eigenvalues if possible
 
173
*
 
174
      GO TO 50
 
175
*
 
176
*     Row and column exchange.
 
177
*
 
178
   20 CONTINUE
 
179
      SCALE( M ) = J
 
180
      IF( J.EQ.M )
 
181
     $   GO TO 30
 
182
*
 
183
      CALL ZSWAP( L, A( 1, J ), 1, A( 1, M ), 1 )
 
184
      CALL ZSWAP( N-K+1, A( J, K ), LDA, A( M, K ), LDA )
 
185
*
 
186
   30 CONTINUE
 
187
      GO TO ( 40, 80 )IEXC
 
188
*
 
189
*     Search for rows isolating an eigenvalue and push them down.
 
190
*
 
191
   40 CONTINUE
 
192
      IF( L.EQ.1 )
 
193
     $   GO TO 210
 
194
      L = L - 1
 
195
*
 
196
   50 CONTINUE
 
197
      DO 70 J = L, 1, -1
 
198
*
 
199
         DO 60 I = 1, L
 
200
            IF( I.EQ.J )
 
201
     $         GO TO 60
 
202
            IF( DBLE( A( J, I ) ).NE.ZERO .OR. DIMAG( A( J, I ) ).NE.
 
203
     $          ZERO )GO TO 70
 
204
   60    CONTINUE
 
205
*
 
206
         M = L
 
207
         IEXC = 1
 
208
         GO TO 20
 
209
   70 CONTINUE
 
210
*
 
211
      GO TO 90
 
212
*
 
213
*     Search for columns isolating an eigenvalue and push them left.
 
214
*
 
215
   80 CONTINUE
 
216
      K = K + 1
 
217
*
 
218
   90 CONTINUE
 
219
      DO 110 J = K, L
 
220
*
 
221
         DO 100 I = K, L
 
222
            IF( I.EQ.J )
 
223
     $         GO TO 100
 
224
            IF( DBLE( A( I, J ) ).NE.ZERO .OR. DIMAG( A( I, J ) ).NE.
 
225
     $          ZERO )GO TO 110
 
226
  100    CONTINUE
 
227
*
 
228
         M = K
 
229
         IEXC = 2
 
230
         GO TO 20
 
231
  110 CONTINUE
 
232
*
 
233
  120 CONTINUE
 
234
      DO 130 I = K, L
 
235
         SCALE( I ) = ONE
 
236
  130 CONTINUE
 
237
*
 
238
      IF( LSAME( JOB, 'P' ) )
 
239
     $   GO TO 210
 
240
*
 
241
*     Balance the submatrix in rows K to L.
 
242
*
 
243
*     Iterative loop for norm reduction
 
244
*
 
245
      SFMIN1 = DLAMCH( 'S' ) / DLAMCH( 'P' )
 
246
      SFMAX1 = ONE / SFMIN1
 
247
      SFMIN2 = SFMIN1*SCLFAC
 
248
      SFMAX2 = ONE / SFMIN2
 
249
  140 CONTINUE
 
250
      NOCONV = .FALSE.
 
251
*
 
252
      DO 200 I = K, L
 
253
         C = ZERO
 
254
         R = ZERO
 
255
*
 
256
         DO 150 J = K, L
 
257
            IF( J.EQ.I )
 
258
     $         GO TO 150
 
259
            C = C + CABS1( A( J, I ) )
 
260
            R = R + CABS1( A( I, J ) )
 
261
  150    CONTINUE
 
262
         ICA = IZAMAX( L, A( 1, I ), 1 )
 
263
         CA = ABS( A( ICA, I ) )
 
264
         IRA = IZAMAX( N-K+1, A( I, K ), LDA )
 
265
         RA = ABS( A( I, IRA+K-1 ) )
 
266
*
 
267
*        Guard against zero C or R due to underflow.
 
268
*
 
269
         IF( C.EQ.ZERO .OR. R.EQ.ZERO )
 
270
     $      GO TO 200
 
271
         G = R / SCLFAC
 
272
         F = ONE
 
273
         S = C + R
 
274
  160    CONTINUE
 
275
         IF( C.GE.G .OR. MAX( F, C, CA ).GE.SFMAX2 .OR.
 
276
     $       MIN( R, G, RA ).LE.SFMIN2 )GO TO 170
 
277
         F = F*SCLFAC
 
278
         C = C*SCLFAC
 
279
         CA = CA*SCLFAC
 
280
         R = R / SCLFAC
 
281
         G = G / SCLFAC
 
282
         RA = RA / SCLFAC
 
283
         GO TO 160
 
284
*
 
285
  170    CONTINUE
 
286
         G = C / SCLFAC
 
287
  180    CONTINUE
 
288
         IF( G.LT.R .OR. MAX( R, RA ).GE.SFMAX2 .OR.
 
289
     $       MIN( F, C, G, CA ).LE.SFMIN2 )GO TO 190
 
290
         F = F / SCLFAC
 
291
         C = C / SCLFAC
 
292
         G = G / SCLFAC
 
293
         CA = CA / SCLFAC
 
294
         R = R*SCLFAC
 
295
         RA = RA*SCLFAC
 
296
         GO TO 180
 
297
*
 
298
*        Now balance.
 
299
*
 
300
  190    CONTINUE
 
301
         IF( ( C+R ).GE.FACTOR*S )
 
302
     $      GO TO 200
 
303
         IF( F.LT.ONE .AND. SCALE( I ).LT.ONE ) THEN
 
304
            IF( F*SCALE( I ).LE.SFMIN1 )
 
305
     $         GO TO 200
 
306
         END IF
 
307
         IF( F.GT.ONE .AND. SCALE( I ).GT.ONE ) THEN
 
308
            IF( SCALE( I ).GE.SFMAX1 / F )
 
309
     $         GO TO 200
 
310
         END IF
 
311
         G = ONE / F
 
312
         SCALE( I ) = SCALE( I )*F
 
313
         NOCONV = .TRUE.
 
314
*
 
315
         CALL ZDSCAL( N-K+1, G, A( I, K ), LDA )
 
316
         CALL ZDSCAL( L, F, A( 1, I ), 1 )
 
317
*
 
318
  200 CONTINUE
 
319
*
 
320
      IF( NOCONV )
 
321
     $   GO TO 140
 
322
*
 
323
  210 CONTINUE
 
324
      ILO = K
 
325
      IHI = L
 
326
*
 
327
      RETURN
 
328
*
 
329
*     End of ZGEBAL
 
330
*
 
331
      END