LAPACK 3.12.1
LAPACK: Linear Algebra PACKage
Loading...
Searching...
No Matches
dggbal.f
Go to the documentation of this file.
1*> \brief \b DGGBAL
2*
3* =========== DOCUMENTATION ===========
4*
5* Online html documentation available at
6* http://www.netlib.org/lapack/explore-html/
7*
8*> Download DGGBAL + dependencies
9*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/dggbal.f">
10*> [TGZ]</a>
11*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/dggbal.f">
12*> [ZIP]</a>
13*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/dggbal.f">
14*> [TXT]</a>
15*
16* Definition:
17* ===========
18*
19* SUBROUTINE DGGBAL( JOB, N, A, LDA, B, LDB, ILO, IHI, LSCALE,
20* RSCALE, WORK, INFO )
21*
22* .. Scalar Arguments ..
23* CHARACTER JOB
24* INTEGER IHI, ILO, INFO, LDA, LDB, N
25* ..
26* .. Array Arguments ..
27* DOUBLE PRECISION A( LDA, * ), B( LDB, * ), LSCALE( * ),
28* $ RSCALE( * ), WORK( * )
29* ..
30*
31*
32*> \par Purpose:
33* =============
34*>
35*> \verbatim
36*>
37*> DGGBAL balances a pair of general real matrices (A,B). This
38*> involves, first, permuting A and B by similarity transformations to
39*> isolate eigenvalues in the first 1 to ILO$-$1 and last IHI+1 to N
40*> elements on the diagonal; and second, applying a diagonal similarity
41*> transformation to rows and columns ILO to IHI to make the rows
42*> and columns as close in norm as possible. Both steps are optional.
43*>
44*> Balancing may reduce the 1-norm of the matrices, and improve the
45*> accuracy of the computed eigenvalues and/or eigenvectors in the
46*> generalized eigenvalue problem A*x = lambda*B*x.
47*> \endverbatim
48*
49* Arguments:
50* ==========
51*
52*> \param[in] JOB
53*> \verbatim
54*> JOB is CHARACTER*1
55*> Specifies the operations to be performed on A and B:
56*> = 'N': none: simply set ILO = 1, IHI = N, LSCALE(I) = 1.0
57*> and RSCALE(I) = 1.0 for i = 1,...,N.
58*> = 'P': permute only;
59*> = 'S': scale only;
60*> = 'B': both permute and scale.
61*> \endverbatim
62*>
63*> \param[in] N
64*> \verbatim
65*> N is INTEGER
66*> The order of the matrices A and B. N >= 0.
67*> \endverbatim
68*>
69*> \param[in,out] A
70*> \verbatim
71*> A is DOUBLE PRECISION array, dimension (LDA,N)
72*> On entry, the input matrix A.
73*> On exit, A is overwritten by the balanced matrix.
74*> If JOB = 'N', A is not referenced.
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[in,out] B
84*> \verbatim
85*> B is DOUBLE PRECISION array, dimension (LDB,N)
86*> On entry, the input matrix B.
87*> On exit, B is overwritten by the balanced matrix.
88*> If JOB = 'N', B is not referenced.
89*> \endverbatim
90*>
91*> \param[in] LDB
92*> \verbatim
93*> LDB is INTEGER
94*> The leading dimension of the array B. LDB >= max(1,N).
95*> \endverbatim
96*>
97*> \param[out] ILO
98*> \verbatim
99*> ILO is INTEGER
100*> \endverbatim
101*>
102*> \param[out] IHI
103*> \verbatim
104*> IHI is INTEGER
105*> ILO and IHI are set to integers such that on exit
106*> A(i,j) = 0 and B(i,j) = 0 if i > j and
107*> j = 1,...,ILO-1 or i = IHI+1,...,N.
108*> If JOB = 'N' or 'S', ILO = 1 and IHI = N.
109*> \endverbatim
110*>
111*> \param[out] LSCALE
112*> \verbatim
113*> LSCALE is DOUBLE PRECISION array, dimension (N)
114*> Details of the permutations and scaling factors applied
115*> to the left side of A and B. If P(j) is the index of the
116*> row interchanged with row j, and D(j)
117*> is the scaling factor applied to row j, then
118*> LSCALE(j) = P(j) for J = 1,...,ILO-1
119*> = D(j) for J = ILO,...,IHI
120*> = P(j) for J = IHI+1,...,N.
121*> The order in which the interchanges are made is N to IHI+1,
122*> then 1 to ILO-1.
123*> \endverbatim
124*>
125*> \param[out] RSCALE
126*> \verbatim
127*> RSCALE is DOUBLE PRECISION array, dimension (N)
128*> Details of the permutations and scaling factors applied
129*> to the right side of A and B. If P(j) is the index of the
130*> column interchanged with column j, and D(j)
131*> is the scaling factor applied to column j, then
132*> LSCALE(j) = P(j) for J = 1,...,ILO-1
133*> = D(j) for J = ILO,...,IHI
134*> = P(j) for J = IHI+1,...,N.
135*> The order in which the interchanges are made is N to IHI+1,
136*> then 1 to ILO-1.
137*> \endverbatim
138*>
139*> \param[out] WORK
140*> \verbatim
141*> WORK is DOUBLE PRECISION array, dimension (lwork)
142*> lwork must be at least max(1,6*N) when JOB = 'S' or 'B', and
143*> at least 1 when JOB = 'N' or 'P'.
144*> \endverbatim
145*>
146*> \param[out] INFO
147*> \verbatim
148*> INFO is INTEGER
149*> = 0: successful exit
150*> < 0: if INFO = -i, the i-th argument had an illegal value.
151*> \endverbatim
152*
153* Authors:
154* ========
155*
156*> \author Univ. of Tennessee
157*> \author Univ. of California Berkeley
158*> \author Univ. of Colorado Denver
159*> \author NAG Ltd.
160*
161*> \ingroup ggbal
162*
163*> \par Further Details:
164* =====================
165*>
166*> \verbatim
167*>
168*> See R.C. WARD, Balancing the generalized eigenvalue problem,
169*> SIAM J. Sci. Stat. Comp. 2 (1981), 141-152.
170*> \endverbatim
171*>
172* =====================================================================
173 SUBROUTINE dggbal( JOB, N, A, LDA, B, LDB, ILO, IHI, LSCALE,
174 $ RSCALE, WORK, INFO )
175*
176* -- LAPACK computational routine --
177* -- LAPACK is a software package provided by Univ. of Tennessee, --
178* -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
179*
180* .. Scalar Arguments ..
181 CHARACTER JOB
182 INTEGER IHI, ILO, INFO, LDA, LDB, N
183* ..
184* .. Array Arguments ..
185 DOUBLE PRECISION A( LDA, * ), B( LDB, * ), LSCALE( * ),
186 $ rscale( * ), work( * )
187* ..
188*
189* =====================================================================
190*
191* .. Parameters ..
192 DOUBLE PRECISION ZERO, HALF, ONE
193 parameter( zero = 0.0d+0, half = 0.5d+0, one = 1.0d+0 )
194 DOUBLE PRECISION THREE, SCLFAC
195 parameter( three = 3.0d+0, sclfac = 1.0d+1 )
196* ..
197* .. Local Scalars ..
198 INTEGER I, ICAB, IFLOW, IP1, IR, IRAB, IT, J, JC, JP1,
199 $ k, kount, l, lcab, lm1, lrab, lsfmax, lsfmin,
200 $ m, nr, nrp2
201 DOUBLE PRECISION ALPHA, BASL, BETA, CAB, CMAX, COEF, COEF2,
202 $ coef5, cor, ew, ewc, gamma, pgamma, rab, sfmax,
203 $ sfmin, sum, t, ta, tb, tc
204* ..
205* .. External Functions ..
206 LOGICAL LSAME
207 INTEGER IDAMAX
208 DOUBLE PRECISION DDOT, DLAMCH
209 EXTERNAL lsame, idamax, ddot, dlamch
210* ..
211* .. External Subroutines ..
212 EXTERNAL daxpy, dscal, dswap, xerbla
213* ..
214* .. Intrinsic Functions ..
215 INTRINSIC abs, dble, int, log10, max, min, sign
216* ..
217* .. Executable Statements ..
218*
219* Test the input parameters
220*
221 info = 0
222 IF( .NOT.lsame( job, 'N' ) .AND.
223 $ .NOT.lsame( job, 'P' ) .AND.
224 $ .NOT.lsame( job, 'S' ) .AND.
225 $ .NOT.lsame( job, 'B' ) ) THEN
226 info = -1
227 ELSE IF( n.LT.0 ) THEN
228 info = -2
229 ELSE IF( lda.LT.max( 1, n ) ) THEN
230 info = -4
231 ELSE IF( ldb.LT.max( 1, n ) ) THEN
232 info = -6
233 END IF
234 IF( info.NE.0 ) THEN
235 CALL xerbla( 'DGGBAL', -info )
236 RETURN
237 END IF
238*
239* Quick return if possible
240*
241 IF( n.EQ.0 ) THEN
242 ilo = 1
243 ihi = n
244 RETURN
245 END IF
246*
247 IF( n.EQ.1 ) THEN
248 ilo = 1
249 ihi = n
250 lscale( 1 ) = one
251 rscale( 1 ) = one
252 RETURN
253 END IF
254*
255 IF( lsame( job, 'N' ) ) THEN
256 ilo = 1
257 ihi = n
258 DO 10 i = 1, n
259 lscale( i ) = one
260 rscale( i ) = one
261 10 CONTINUE
262 RETURN
263 END IF
264*
265 k = 1
266 l = n
267 IF( lsame( job, 'S' ) )
268 $ GO TO 190
269*
270 GO TO 30
271*
272* Permute the matrices A and B to isolate the eigenvalues.
273*
274* Find row with one nonzero in columns 1 through L
275*
276 20 CONTINUE
277 l = lm1
278 IF( l.NE.1 )
279 $ GO TO 30
280*
281 rscale( 1 ) = one
282 lscale( 1 ) = one
283 GO TO 190
284*
285 30 CONTINUE
286 lm1 = l - 1
287 DO 80 i = l, 1, -1
288 DO 40 j = 1, lm1
289 jp1 = j + 1
290 IF( a( i, j ).NE.zero .OR. b( i, j ).NE.zero )
291 $ GO TO 50
292 40 CONTINUE
293 j = l
294 GO TO 70
295*
296 50 CONTINUE
297 DO 60 j = jp1, l
298 IF( a( i, j ).NE.zero .OR. b( i, j ).NE.zero )
299 $ GO TO 80
300 60 CONTINUE
301 j = jp1 - 1
302*
303 70 CONTINUE
304 m = l
305 iflow = 1
306 GO TO 160
307 80 CONTINUE
308 GO TO 100
309*
310* Find column with one nonzero in rows K through N
311*
312 90 CONTINUE
313 k = k + 1
314*
315 100 CONTINUE
316 DO 150 j = k, l
317 DO 110 i = k, lm1
318 ip1 = i + 1
319 IF( a( i, j ).NE.zero .OR. b( i, j ).NE.zero )
320 $ GO TO 120
321 110 CONTINUE
322 i = l
323 GO TO 140
324 120 CONTINUE
325 DO 130 i = ip1, l
326 IF( a( i, j ).NE.zero .OR. b( i, j ).NE.zero )
327 $ GO TO 150
328 130 CONTINUE
329 i = ip1 - 1
330 140 CONTINUE
331 m = k
332 iflow = 2
333 GO TO 160
334 150 CONTINUE
335 GO TO 190
336*
337* Permute rows M and I
338*
339 160 CONTINUE
340 lscale( m ) = i
341 IF( i.EQ.m )
342 $ GO TO 170
343 CALL dswap( n-k+1, a( i, k ), lda, a( m, k ), lda )
344 CALL dswap( n-k+1, b( i, k ), ldb, b( m, k ), ldb )
345*
346* Permute columns M and J
347*
348 170 CONTINUE
349 rscale( m ) = j
350 IF( j.EQ.m )
351 $ GO TO 180
352 CALL dswap( l, a( 1, j ), 1, a( 1, m ), 1 )
353 CALL dswap( l, b( 1, j ), 1, b( 1, m ), 1 )
354*
355 180 CONTINUE
356 GO TO ( 20, 90 )iflow
357*
358 190 CONTINUE
359 ilo = k
360 ihi = l
361*
362 IF( lsame( job, 'P' ) ) THEN
363 DO 195 i = ilo, ihi
364 lscale( i ) = one
365 rscale( i ) = one
366 195 CONTINUE
367 RETURN
368 END IF
369*
370 IF( ilo.EQ.ihi )
371 $ RETURN
372*
373* Balance the submatrix in rows ILO to IHI.
374*
375 nr = ihi - ilo + 1
376 DO 200 i = ilo, ihi
377 rscale( i ) = zero
378 lscale( i ) = zero
379*
380 work( i ) = zero
381 work( i+n ) = zero
382 work( i+2*n ) = zero
383 work( i+3*n ) = zero
384 work( i+4*n ) = zero
385 work( i+5*n ) = zero
386 200 CONTINUE
387*
388* Compute right side vector in resulting linear equations
389*
390 basl = log10( sclfac )
391 DO 240 i = ilo, ihi
392 DO 230 j = ilo, ihi
393 tb = b( i, j )
394 ta = a( i, j )
395 IF( ta.EQ.zero )
396 $ GO TO 210
397 ta = log10( abs( ta ) ) / basl
398 210 CONTINUE
399 IF( tb.EQ.zero )
400 $ GO TO 220
401 tb = log10( abs( tb ) ) / basl
402 220 CONTINUE
403 work( i+4*n ) = work( i+4*n ) - ta - tb
404 work( j+5*n ) = work( j+5*n ) - ta - tb
405 230 CONTINUE
406 240 CONTINUE
407*
408 coef = one / dble( 2*nr )
409 coef2 = coef*coef
410 coef5 = half*coef2
411 nrp2 = nr + 2
412 beta = zero
413 it = 1
414*
415* Start generalized conjugate gradient iteration
416*
417 250 CONTINUE
418*
419 gamma = ddot( nr, work( ilo+4*n ), 1, work( ilo+4*n ), 1 ) +
420 $ ddot( nr, work( ilo+5*n ), 1, work( ilo+5*n ), 1 )
421*
422 ew = zero
423 ewc = zero
424 DO 260 i = ilo, ihi
425 ew = ew + work( i+4*n )
426 ewc = ewc + work( i+5*n )
427 260 CONTINUE
428*
429 gamma = coef*gamma - coef2*( ew**2+ewc**2 ) - coef5*( ew-ewc )**2
430 IF( gamma.EQ.zero )
431 $ GO TO 350
432 IF( it.NE.1 )
433 $ beta = gamma / pgamma
434 t = coef5*( ewc-three*ew )
435 tc = coef5*( ew-three*ewc )
436*
437 CALL dscal( nr, beta, work( ilo ), 1 )
438 CALL dscal( nr, beta, work( ilo+n ), 1 )
439*
440 CALL daxpy( nr, coef, work( ilo+4*n ), 1, work( ilo+n ), 1 )
441 CALL daxpy( nr, coef, work( ilo+5*n ), 1, work( ilo ), 1 )
442*
443 DO 270 i = ilo, ihi
444 work( i ) = work( i ) + tc
445 work( i+n ) = work( i+n ) + t
446 270 CONTINUE
447*
448* Apply matrix to vector
449*
450 DO 300 i = ilo, ihi
451 kount = 0
452 sum = zero
453 DO 290 j = ilo, ihi
454 IF( a( i, j ).EQ.zero )
455 $ GO TO 280
456 kount = kount + 1
457 sum = sum + work( j )
458 280 CONTINUE
459 IF( b( i, j ).EQ.zero )
460 $ GO TO 290
461 kount = kount + 1
462 sum = sum + work( j )
463 290 CONTINUE
464 work( i+2*n ) = dble( kount )*work( i+n ) + sum
465 300 CONTINUE
466*
467 DO 330 j = ilo, ihi
468 kount = 0
469 sum = zero
470 DO 320 i = ilo, ihi
471 IF( a( i, j ).EQ.zero )
472 $ GO TO 310
473 kount = kount + 1
474 sum = sum + work( i+n )
475 310 CONTINUE
476 IF( b( i, j ).EQ.zero )
477 $ GO TO 320
478 kount = kount + 1
479 sum = sum + work( i+n )
480 320 CONTINUE
481 work( j+3*n ) = dble( kount )*work( j ) + sum
482 330 CONTINUE
483*
484 sum = ddot( nr, work( ilo+n ), 1, work( ilo+2*n ), 1 ) +
485 $ ddot( nr, work( ilo ), 1, work( ilo+3*n ), 1 )
486 alpha = gamma / sum
487*
488* Determine correction to current iteration
489*
490 cmax = zero
491 DO 340 i = ilo, ihi
492 cor = alpha*work( i+n )
493 IF( abs( cor ).GT.cmax )
494 $ cmax = abs( cor )
495 lscale( i ) = lscale( i ) + cor
496 cor = alpha*work( i )
497 IF( abs( cor ).GT.cmax )
498 $ cmax = abs( cor )
499 rscale( i ) = rscale( i ) + cor
500 340 CONTINUE
501 IF( cmax.LT.half )
502 $ GO TO 350
503*
504 CALL daxpy( nr, -alpha, work( ilo+2*n ), 1, work( ilo+4*n ),
505 $ 1 )
506 CALL daxpy( nr, -alpha, work( ilo+3*n ), 1, work( ilo+5*n ),
507 $ 1 )
508*
509 pgamma = gamma
510 it = it + 1
511 IF( it.LE.nrp2 )
512 $ GO TO 250
513*
514* End generalized conjugate gradient iteration
515*
516 350 CONTINUE
517 sfmin = dlamch( 'S' )
518 sfmax = one / sfmin
519 lsfmin = int( log10( sfmin ) / basl+one )
520 lsfmax = int( log10( sfmax ) / basl )
521 DO 360 i = ilo, ihi
522 irab = idamax( n-ilo+1, a( i, ilo ), lda )
523 rab = abs( a( i, irab+ilo-1 ) )
524 irab = idamax( n-ilo+1, b( i, ilo ), ldb )
525 rab = max( rab, abs( b( i, irab+ilo-1 ) ) )
526 lrab = int( log10( rab+sfmin ) / basl+one )
527 ir = int(lscale( i ) + sign( half, lscale( i ) ))
528 ir = min( max( ir, lsfmin ), lsfmax, lsfmax-lrab )
529 lscale( i ) = sclfac**ir
530 icab = idamax( ihi, a( 1, i ), 1 )
531 cab = abs( a( icab, i ) )
532 icab = idamax( ihi, b( 1, i ), 1 )
533 cab = max( cab, abs( b( icab, i ) ) )
534 lcab = int( log10( cab+sfmin ) / basl+one )
535 jc = int(rscale( i ) + sign( half, rscale( i ) ))
536 jc = min( max( jc, lsfmin ), lsfmax, lsfmax-lcab )
537 rscale( i ) = sclfac**jc
538 360 CONTINUE
539*
540* Row scaling of matrices A and B
541*
542 DO 370 i = ilo, ihi
543 CALL dscal( n-ilo+1, lscale( i ), a( i, ilo ), lda )
544 CALL dscal( n-ilo+1, lscale( i ), b( i, ilo ), ldb )
545 370 CONTINUE
546*
547* Column scaling of matrices A and B
548*
549 DO 380 j = ilo, ihi
550 CALL dscal( ihi, rscale( j ), a( 1, j ), 1 )
551 CALL dscal( ihi, rscale( j ), b( 1, j ), 1 )
552 380 CONTINUE
553*
554 RETURN
555*
556* End of DGGBAL
557*
558 END
xerbla
subroutine xerbla(srname, info)
Definition cblat2.f:3285
daxpy
subroutine daxpy(n, da, dx, incx, dy, incy)
DAXPY
Definition daxpy.f:89
dggbal
subroutine dggbal(job, n, a, lda, b, ldb, ilo, ihi, lscale, rscale, work, info)
DGGBAL
Definition dggbal.f:175
dscal
subroutine dscal(n, da, dx, incx)
DSCAL
Definition dscal.f:79
dswap
subroutine dswap(n, dx, incx, dy, incy)
DSWAP
Definition dswap.f:82