LAPACK 3.12.1
LAPACK: Linear Algebra PACKage
Loading...
Searching...
No Matches

◆ zggbal()

subroutine zggbal ( character job,
integer n,
complex*16, dimension( lda, * ) a,
integer lda,
complex*16, dimension( ldb, * ) b,
integer ldb,
integer ilo,
integer ihi,
double precision, dimension( * ) lscale,
double precision, dimension( * ) rscale,
double precision, dimension( * ) work,
integer info )

ZGGBAL

Download ZGGBAL + dependencies [TGZ] [ZIP] [TXT]

Purpose:
!>
!> ZGGBAL balances a pair of general complex matrices (A,B).  This
!> involves, first, permuting A and B by similarity transformations to
!> isolate eigenvalues in the first 1 to ILO$-$1 and last IHI+1 to N
!> elements on the diagonal; and second, applying a diagonal similarity
!> transformation to rows and columns ILO to IHI to make the rows
!> and columns as close in norm as possible. Both steps are optional.
!>
!> Balancing may reduce the 1-norm of the matrices, and improve the
!> accuracy of the computed eigenvalues and/or eigenvectors in the
!> generalized eigenvalue problem A*x = lambda*B*x.
!>
Parameters
[in]JOB
!>          JOB is CHARACTER*1
!>          Specifies the operations to be performed on A and B:
!>          = 'N':  none:  simply set ILO = 1, IHI = N, LSCALE(I) = 1.0
!>                  and RSCALE(I) = 1.0 for i=1,...,N;
!>          = 'P':  permute only;
!>          = 'S':  scale only;
!>          = 'B':  both permute and scale.
!>
[in]N
!>          N is INTEGER
!>          The order of the matrices A and B.  N >= 0.
!>
[in,out]A
!>          A is COMPLEX*16 array, dimension (LDA,N)
!>          On entry, the input matrix A.
!>          On exit, A is overwritten by the balanced matrix.
!>          If JOB = 'N', A is not referenced.
!>
[in]LDA
!>          LDA is INTEGER
!>          The leading dimension of the array A. LDA >= max(1,N).
!>
[in,out]B
!>          B is COMPLEX*16 array, dimension (LDB,N)
!>          On entry, the input matrix B.
!>          On exit, B is overwritten by the balanced matrix.
!>          If JOB = 'N', B is not referenced.
!>
[in]LDB
!>          LDB is INTEGER
!>          The leading dimension of the array B. LDB >= max(1,N).
!>
[out]ILO
!>          ILO is INTEGER
!>
[out]IHI
!>          IHI is INTEGER
!>          ILO and IHI are set to integers such that on exit
!>          A(i,j) = 0 and B(i,j) = 0 if i > j and
!>          j = 1,...,ILO-1 or i = IHI+1,...,N.
!>          If JOB = 'N' or 'S', ILO = 1 and IHI = N.
!>
[out]LSCALE
!>          LSCALE is DOUBLE PRECISION array, dimension (N)
!>          Details of the permutations and scaling factors applied
!>          to the left side of A and B.  If P(j) is the index of the
!>          row interchanged with row j, and D(j) is the scaling factor
!>          applied to row j, then
!>            LSCALE(j) = P(j)    for J = 1,...,ILO-1
!>                      = D(j)    for J = ILO,...,IHI
!>                      = P(j)    for J = IHI+1,...,N.
!>          The order in which the interchanges are made is N to IHI+1,
!>          then 1 to ILO-1.
!>
[out]RSCALE
!>          RSCALE is DOUBLE PRECISION array, dimension (N)
!>          Details of the permutations and scaling factors applied
!>          to the right side of A and B.  If P(j) is the index of the
!>          column interchanged with column j, and D(j) is the scaling
!>          factor applied to column j, then
!>            RSCALE(j) = P(j)    for J = 1,...,ILO-1
!>                      = D(j)    for J = ILO,...,IHI
!>                      = P(j)    for J = IHI+1,...,N.
!>          The order in which the interchanges are made is N to IHI+1,
!>          then 1 to ILO-1.
!>
[out]WORK
!>          WORK is DOUBLE PRECISION array, dimension (lwork)
!>          lwork must be at least max(1,6*N) when JOB = 'S' or 'B', and
!>          at least 1 when JOB = 'N' or 'P'.
!>
[out]INFO
!>          INFO is INTEGER
!>          = 0:  successful exit
!>          < 0:  if INFO = -i, the i-th argument had an illegal value.
!>
Author
Univ. of Tennessee
Univ. of California Berkeley
Univ. of Colorado Denver
NAG Ltd.
Further Details:
!>
!>  See R.C. WARD, Balancing the generalized eigenvalue problem,
!>                 SIAM J. Sci. Stat. Comp. 2 (1981), 141-152.
!>

Definition at line 173 of file zggbal.f.

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 LSCALE( * ), RSCALE( * ), WORK( * )
186 COMPLEX*16 A( LDA, * ), B( LDB, * )
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 COMPLEX*16 CZERO
197 parameter( czero = ( 0.0d+0, 0.0d+0 ) )
198* ..
199* .. Local Scalars ..
200 INTEGER I, ICAB, IFLOW, IP1, IR, IRAB, IT, J, JC, JP1,
201 $ K, KOUNT, L, LCAB, LM1, LRAB, LSFMAX, LSFMIN,
202 $ M, NR, NRP2
203 DOUBLE PRECISION ALPHA, BASL, BETA, CAB, CMAX, COEF, COEF2,
204 $ COEF5, COR, EW, EWC, GAMMA, PGAMMA, RAB, SFMAX,
205 $ SFMIN, SUM, T, TA, TB, TC
206 COMPLEX*16 CDUM
207* ..
208* .. External Functions ..
209 LOGICAL LSAME
210 INTEGER IZAMAX
211 DOUBLE PRECISION DDOT, DLAMCH
212 EXTERNAL lsame, izamax, ddot, dlamch
213* ..
214* .. External Subroutines ..
215 EXTERNAL daxpy, dscal, xerbla, zdscal, zswap
216* ..
217* .. Intrinsic Functions ..
218 INTRINSIC abs, dble, dimag, int, log10, max, min, sign
219* ..
220* .. Statement Functions ..
221 DOUBLE PRECISION CABS1
222* ..
223* .. Statement Function definitions ..
224 cabs1( cdum ) = abs( dble( cdum ) ) + abs( dimag( cdum ) )
225* ..
226* .. Executable Statements ..
227*
228* Test the input parameters
229*
230 info = 0
231 IF( .NOT.lsame( job, 'N' ) .AND.
232 $ .NOT.lsame( job, 'P' ) .AND.
233 $ .NOT.lsame( job, 'S' ) .AND.
234 $ .NOT.lsame( job, 'B' ) ) THEN
235 info = -1
236 ELSE IF( n.LT.0 ) THEN
237 info = -2
238 ELSE IF( lda.LT.max( 1, n ) ) THEN
239 info = -4
240 ELSE IF( ldb.LT.max( 1, n ) ) THEN
241 info = -6
242 END IF
243 IF( info.NE.0 ) THEN
244 CALL xerbla( 'ZGGBAL', -info )
245 RETURN
246 END IF
247*
248* Quick return if possible
249*
250 IF( n.EQ.0 ) THEN
251 ilo = 1
252 ihi = n
253 RETURN
254 END IF
255*
256 IF( n.EQ.1 ) THEN
257 ilo = 1
258 ihi = n
259 lscale( 1 ) = one
260 rscale( 1 ) = one
261 RETURN
262 END IF
263*
264 IF( lsame( job, 'N' ) ) THEN
265 ilo = 1
266 ihi = n
267 DO 10 i = 1, n
268 lscale( i ) = one
269 rscale( i ) = one
270 10 CONTINUE
271 RETURN
272 END IF
273*
274 k = 1
275 l = n
276 IF( lsame( job, 'S' ) )
277 $ GO TO 190
278*
279 GO TO 30
280*
281* Permute the matrices A and B to isolate the eigenvalues.
282*
283* Find row with one nonzero in columns 1 through L
284*
285 20 CONTINUE
286 l = lm1
287 IF( l.NE.1 )
288 $ GO TO 30
289*
290 rscale( 1 ) = 1
291 lscale( 1 ) = 1
292 GO TO 190
293*
294 30 CONTINUE
295 lm1 = l - 1
296 DO 80 i = l, 1, -1
297 DO 40 j = 1, lm1
298 jp1 = j + 1
299 IF( a( i, j ).NE.czero .OR. b( i, j ).NE.czero )
300 $ GO TO 50
301 40 CONTINUE
302 j = l
303 GO TO 70
304*
305 50 CONTINUE
306 DO 60 j = jp1, l
307 IF( a( i, j ).NE.czero .OR. b( i, j ).NE.czero )
308 $ GO TO 80
309 60 CONTINUE
310 j = jp1 - 1
311*
312 70 CONTINUE
313 m = l
314 iflow = 1
315 GO TO 160
316 80 CONTINUE
317 GO TO 100
318*
319* Find column with one nonzero in rows K through N
320*
321 90 CONTINUE
322 k = k + 1
323*
324 100 CONTINUE
325 DO 150 j = k, l
326 DO 110 i = k, lm1
327 ip1 = i + 1
328 IF( a( i, j ).NE.czero .OR. b( i, j ).NE.czero )
329 $ GO TO 120
330 110 CONTINUE
331 i = l
332 GO TO 140
333 120 CONTINUE
334 DO 130 i = ip1, l
335 IF( a( i, j ).NE.czero .OR. b( i, j ).NE.czero )
336 $ GO TO 150
337 130 CONTINUE
338 i = ip1 - 1
339 140 CONTINUE
340 m = k
341 iflow = 2
342 GO TO 160
343 150 CONTINUE
344 GO TO 190
345*
346* Permute rows M and I
347*
348 160 CONTINUE
349 lscale( m ) = i
350 IF( i.EQ.m )
351 $ GO TO 170
352 CALL zswap( n-k+1, a( i, k ), lda, a( m, k ), lda )
353 CALL zswap( n-k+1, b( i, k ), ldb, b( m, k ), ldb )
354*
355* Permute columns M and J
356*
357 170 CONTINUE
358 rscale( m ) = j
359 IF( j.EQ.m )
360 $ GO TO 180
361 CALL zswap( l, a( 1, j ), 1, a( 1, m ), 1 )
362 CALL zswap( l, b( 1, j ), 1, b( 1, m ), 1 )
363*
364 180 CONTINUE
365 GO TO ( 20, 90 )iflow
366*
367 190 CONTINUE
368 ilo = k
369 ihi = l
370*
371 IF( lsame( job, 'P' ) ) THEN
372 DO 195 i = ilo, ihi
373 lscale( i ) = one
374 rscale( i ) = one
375 195 CONTINUE
376 RETURN
377 END IF
378*
379 IF( ilo.EQ.ihi )
380 $ RETURN
381*
382* Balance the submatrix in rows ILO to IHI.
383*
384 nr = ihi - ilo + 1
385 DO 200 i = ilo, ihi
386 rscale( i ) = zero
387 lscale( i ) = zero
388*
389 work( i ) = zero
390 work( i+n ) = zero
391 work( i+2*n ) = zero
392 work( i+3*n ) = zero
393 work( i+4*n ) = zero
394 work( i+5*n ) = zero
395 200 CONTINUE
396*
397* Compute right side vector in resulting linear equations
398*
399 basl = log10( sclfac )
400 DO 240 i = ilo, ihi
401 DO 230 j = ilo, ihi
402 IF( a( i, j ).EQ.czero ) THEN
403 ta = zero
404 GO TO 210
405 END IF
406 ta = log10( cabs1( a( i, j ) ) ) / basl
407*
408 210 CONTINUE
409 IF( b( i, j ).EQ.czero ) THEN
410 tb = zero
411 GO TO 220
412 END IF
413 tb = log10( cabs1( b( i, j ) ) ) / basl
414*
415 220 CONTINUE
416 work( i+4*n ) = work( i+4*n ) - ta - tb
417 work( j+5*n ) = work( j+5*n ) - ta - tb
418 230 CONTINUE
419 240 CONTINUE
420*
421 coef = one / dble( 2*nr )
422 coef2 = coef*coef
423 coef5 = half*coef2
424 nrp2 = nr + 2
425 beta = zero
426 it = 1
427*
428* Start generalized conjugate gradient iteration
429*
430 250 CONTINUE
431*
432 gamma = ddot( nr, work( ilo+4*n ), 1, work( ilo+4*n ), 1 ) +
433 $ ddot( nr, work( ilo+5*n ), 1, work( ilo+5*n ), 1 )
434*
435 ew = zero
436 ewc = zero
437 DO 260 i = ilo, ihi
438 ew = ew + work( i+4*n )
439 ewc = ewc + work( i+5*n )
440 260 CONTINUE
441*
442 gamma = coef*gamma - coef2*( ew**2+ewc**2 ) - coef5*( ew-ewc )**2
443 IF( gamma.EQ.zero )
444 $ GO TO 350
445 IF( it.NE.1 )
446 $ beta = gamma / pgamma
447 t = coef5*( ewc-three*ew )
448 tc = coef5*( ew-three*ewc )
449*
450 CALL dscal( nr, beta, work( ilo ), 1 )
451 CALL dscal( nr, beta, work( ilo+n ), 1 )
452*
453 CALL daxpy( nr, coef, work( ilo+4*n ), 1, work( ilo+n ), 1 )
454 CALL daxpy( nr, coef, work( ilo+5*n ), 1, work( ilo ), 1 )
455*
456 DO 270 i = ilo, ihi
457 work( i ) = work( i ) + tc
458 work( i+n ) = work( i+n ) + t
459 270 CONTINUE
460*
461* Apply matrix to vector
462*
463 DO 300 i = ilo, ihi
464 kount = 0
465 sum = zero
466 DO 290 j = ilo, ihi
467 IF( a( i, j ).EQ.czero )
468 $ GO TO 280
469 kount = kount + 1
470 sum = sum + work( j )
471 280 CONTINUE
472 IF( b( i, j ).EQ.czero )
473 $ GO TO 290
474 kount = kount + 1
475 sum = sum + work( j )
476 290 CONTINUE
477 work( i+2*n ) = dble( kount )*work( i+n ) + sum
478 300 CONTINUE
479*
480 DO 330 j = ilo, ihi
481 kount = 0
482 sum = zero
483 DO 320 i = ilo, ihi
484 IF( a( i, j ).EQ.czero )
485 $ GO TO 310
486 kount = kount + 1
487 sum = sum + work( i+n )
488 310 CONTINUE
489 IF( b( i, j ).EQ.czero )
490 $ GO TO 320
491 kount = kount + 1
492 sum = sum + work( i+n )
493 320 CONTINUE
494 work( j+3*n ) = dble( kount )*work( j ) + sum
495 330 CONTINUE
496*
497 sum = ddot( nr, work( ilo+n ), 1, work( ilo+2*n ), 1 ) +
498 $ ddot( nr, work( ilo ), 1, work( ilo+3*n ), 1 )
499 alpha = gamma / sum
500*
501* Determine correction to current iteration
502*
503 cmax = zero
504 DO 340 i = ilo, ihi
505 cor = alpha*work( i+n )
506 IF( abs( cor ).GT.cmax )
507 $ cmax = abs( cor )
508 lscale( i ) = lscale( i ) + cor
509 cor = alpha*work( i )
510 IF( abs( cor ).GT.cmax )
511 $ cmax = abs( cor )
512 rscale( i ) = rscale( i ) + cor
513 340 CONTINUE
514 IF( cmax.LT.half )
515 $ GO TO 350
516*
517 CALL daxpy( nr, -alpha, work( ilo+2*n ), 1, work( ilo+4*n ),
518 $ 1 )
519 CALL daxpy( nr, -alpha, work( ilo+3*n ), 1, work( ilo+5*n ),
520 $ 1 )
521*
522 pgamma = gamma
523 it = it + 1
524 IF( it.LE.nrp2 )
525 $ GO TO 250
526*
527* End generalized conjugate gradient iteration
528*
529 350 CONTINUE
530 sfmin = dlamch( 'S' )
531 sfmax = one / sfmin
532 lsfmin = int( log10( sfmin ) / basl+one )
533 lsfmax = int( log10( sfmax ) / basl )
534 DO 360 i = ilo, ihi
535 irab = izamax( n-ilo+1, a( i, ilo ), lda )
536 rab = abs( a( i, irab+ilo-1 ) )
537 irab = izamax( n-ilo+1, b( i, ilo ), ldb )
538 rab = max( rab, abs( b( i, irab+ilo-1 ) ) )
539 lrab = int( log10( rab+sfmin ) / basl+one )
540 ir = int(lscale( i ) + sign( half, lscale( i ) ))
541 ir = min( max( ir, lsfmin ), lsfmax, lsfmax-lrab )
542 lscale( i ) = sclfac**ir
543 icab = izamax( ihi, a( 1, i ), 1 )
544 cab = abs( a( icab, i ) )
545 icab = izamax( ihi, b( 1, i ), 1 )
546 cab = max( cab, abs( b( icab, i ) ) )
547 lcab = int( log10( cab+sfmin ) / basl+one )
548 jc = int(rscale( i ) + sign( half, rscale( i ) ))
549 jc = min( max( jc, lsfmin ), lsfmax, lsfmax-lcab )
550 rscale( i ) = sclfac**jc
551 360 CONTINUE
552*
553* Row scaling of matrices A and B
554*
555 DO 370 i = ilo, ihi
556 CALL zdscal( n-ilo+1, lscale( i ), a( i, ilo ), lda )
557 CALL zdscal( n-ilo+1, lscale( i ), b( i, ilo ), ldb )
558 370 CONTINUE
559*
560* Column scaling of matrices A and B
561*
562 DO 380 j = ilo, ihi
563 CALL zdscal( ihi, rscale( j ), a( 1, j ), 1 )
564 CALL zdscal( ihi, rscale( j ), b( 1, j ), 1 )
565 380 CONTINUE
566*
567 RETURN
568*
569* End of ZGGBAL
570*
xerbla
subroutine xerbla(srname, info)
Definition cblat2.f:3285
daxpy
subroutine daxpy(n, da, dx, incx, dy, incy)
DAXPY
Definition daxpy.f:89
ddot
double precision function ddot(n, dx, incx, dy, incy)
DDOT
Definition ddot.f:82
izamax
integer function izamax(n, zx, incx)
IZAMAX
Definition izamax.f:71
dlamch
double precision function dlamch(cmach)
DLAMCH
Definition dlamch.f:69
lsame
logical function lsame(ca, cb)
LSAME
Definition lsame.f:48
dscal
subroutine dscal(n, da, dx, incx)
DSCAL
Definition dscal.f:79
zdscal
subroutine zdscal(n, da, zx, incx)
ZDSCAL
Definition zdscal.f:78
zswap
subroutine zswap(n, zx, incx, zy, incy)
ZSWAP
Definition zswap.f:81
Here is the call graph for this function:
Here is the caller graph for this function: