LAPACK  3.10.0 LAPACK: Linear Algebra PACKage

## ◆ sggbal()

 subroutine sggbal ( character JOB, integer N, real, dimension( lda, * ) A, integer LDA, real, dimension( ldb, * ) B, integer LDB, integer ILO, integer IHI, real, dimension( * ) LSCALE, real, dimension( * ) RSCALE, real, dimension( * ) WORK, integer INFO )

SGGBAL

Purpose:
``` SGGBAL balances a pair of general real 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 REAL 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 REAL 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 REAL 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 REAL 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 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] WORK ``` WORK is REAL 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.```
Further Details:
```  See R.C. WARD, Balancing the generalized eigenvalue problem,
SIAM J. Sci. Stat. Comp. 2 (1981), 141-152.```

Definition at line 175 of file sggbal.f.

177 *
178 * -- LAPACK computational routine --
179 * -- LAPACK is a software package provided by Univ. of Tennessee, --
180 * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
181 *
182 * .. Scalar Arguments ..
183  CHARACTER JOB
184  INTEGER IHI, ILO, INFO, LDA, LDB, N
185 * ..
186 * .. Array Arguments ..
187  REAL A( LDA, * ), B( LDB, * ), LSCALE( * ),
188  \$ RSCALE( * ), WORK( * )
189 * ..
190 *
191 * =====================================================================
192 *
193 * .. Parameters ..
194  REAL ZERO, HALF, ONE
195  parameter( zero = 0.0e+0, half = 0.5e+0, one = 1.0e+0 )
196  REAL THREE, SCLFAC
197  parameter( three = 3.0e+0, sclfac = 1.0e+1 )
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  REAL ALPHA, BASL, BETA, CAB, CMAX, COEF, COEF2,
204  \$ COEF5, COR, EW, EWC, GAMMA, PGAMMA, RAB, SFMAX,
205  \$ SFMIN, SUM, T, TA, TB, TC
206 * ..
207 * .. External Functions ..
208  LOGICAL LSAME
209  INTEGER ISAMAX
210  REAL SDOT, SLAMCH
211  EXTERNAL lsame, isamax, sdot, slamch
212 * ..
213 * .. External Subroutines ..
214  EXTERNAL saxpy, sscal, sswap, xerbla
215 * ..
216 * .. Intrinsic Functions ..
217  INTRINSIC abs, int, log10, max, min, real, sign
218 * ..
219 * .. Executable Statements ..
220 *
221 * Test the input parameters
222 *
223  info = 0
224  IF( .NOT.lsame( job, 'N' ) .AND. .NOT.lsame( job, 'P' ) .AND.
225  \$ .NOT.lsame( job, 'S' ) .AND. .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( 'SGGBAL', -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 sswap( n-k+1, a( i, k ), lda, a( m, k ), lda )
344  CALL sswap( 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 sswap( l, a( 1, j ), 1, a( 1, m ), 1 )
353  CALL sswap( 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 / real( 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 = sdot( nr, work( ilo+4*n ), 1, work( ilo+4*n ), 1 ) +
420  \$ sdot( 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 sscal( nr, beta, work( ilo ), 1 )
438  CALL sscal( nr, beta, work( ilo+n ), 1 )
439 *
440  CALL saxpy( nr, coef, work( ilo+4*n ), 1, work( ilo+n ), 1 )
441  CALL saxpy( 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 ) = real( 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 ) = real( kount )*work( j ) + sum
482  330 CONTINUE
483 *
484  sum = sdot( nr, work( ilo+n ), 1, work( ilo+2*n ), 1 ) +
485  \$ sdot( 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 saxpy( nr, -alpha, work( ilo+2*n ), 1, work( ilo+4*n ), 1 )
505  CALL saxpy( nr, -alpha, work( ilo+3*n ), 1, work( ilo+5*n ), 1 )
506 *
507  pgamma = gamma
508  it = it + 1
509  IF( it.LE.nrp2 )
510  \$ GO TO 250
511 *
512 * End generalized conjugate gradient iteration
513 *
514  350 CONTINUE
515  sfmin = slamch( 'S' )
516  sfmax = one / sfmin
517  lsfmin = int( log10( sfmin ) / basl+one )
518  lsfmax = int( log10( sfmax ) / basl )
519  DO 360 i = ilo, ihi
520  irab = isamax( n-ilo+1, a( i, ilo ), lda )
521  rab = abs( a( i, irab+ilo-1 ) )
522  irab = isamax( n-ilo+1, b( i, ilo ), ldb )
523  rab = max( rab, abs( b( i, irab+ilo-1 ) ) )
524  lrab = int( log10( rab+sfmin ) / basl+one )
525  ir = lscale( i ) + sign( half, lscale( i ) )
526  ir = min( max( ir, lsfmin ), lsfmax, lsfmax-lrab )
527  lscale( i ) = sclfac**ir
528  icab = isamax( ihi, a( 1, i ), 1 )
529  cab = abs( a( icab, i ) )
530  icab = isamax( ihi, b( 1, i ), 1 )
531  cab = max( cab, abs( b( icab, i ) ) )
532  lcab = int( log10( cab+sfmin ) / basl+one )
533  jc = rscale( i ) + sign( half, rscale( i ) )
534  jc = min( max( jc, lsfmin ), lsfmax, lsfmax-lcab )
535  rscale( i ) = sclfac**jc
536  360 CONTINUE
537 *
538 * Row scaling of matrices A and B
539 *
540  DO 370 i = ilo, ihi
541  CALL sscal( n-ilo+1, lscale( i ), a( i, ilo ), lda )
542  CALL sscal( n-ilo+1, lscale( i ), b( i, ilo ), ldb )
543  370 CONTINUE
544 *
545 * Column scaling of matrices A and B
546 *
547  DO 380 j = ilo, ihi
548  CALL sscal( ihi, rscale( j ), a( 1, j ), 1 )
549  CALL sscal( ihi, rscale( j ), b( 1, j ), 1 )
550  380 CONTINUE
551 *
552  RETURN
553 *
554 * End of SGGBAL
555 *
integer function isamax(N, SX, INCX)
ISAMAX
Definition: isamax.f:71
subroutine xerbla(SRNAME, INFO)
XERBLA
Definition: xerbla.f:60
logical function lsame(CA, CB)
LSAME
Definition: lsame.f:53
subroutine sswap(N, SX, INCX, SY, INCY)
SSWAP
Definition: sswap.f:82
subroutine sscal(N, SA, SX, INCX)
SSCAL
Definition: sscal.f:79
real function sdot(N, SX, INCX, SY, INCY)
SDOT
Definition: sdot.f:82
subroutine saxpy(N, SA, SX, INCX, SY, INCY)
SAXPY
Definition: saxpy.f:89
real function slamch(CMACH)
SLAMCH
Definition: slamch.f:68
Here is the call graph for this function:
Here is the caller graph for this function: