Sun Performance Library Reference Manual: 1

Equilibration Scale Factors for a Symmetric Positive Definite Matrix in Banded Storage

The subroutines described in this section compute equilibration scale factors for a real symmetric (or Hermitian) positive definite matrix A in banded storage. Equilibration is intended to reduce the condition number of A with respect to the 2-norm. SCALE contains scale factors chosen so that the scaled matrix B with elements B(i,j) = SCALE(i) × A(i,j) × SCALE(j) has ones on the diagonal. This choice of SCALE puts the condition number of B within a factor N of the smallest possible condition number over all possible diagonal scalings. Note that these subroutines only compute the scale factors for A. Additional code must be written to actually equilibrate A.

Calling Sequence


CALL DPBEQU	(UPLO, N, NDIAG, DA, LDA, DSCALE, DSCOND, DAMAX, INFO)
CALL SPBEQU	(UPLO, N, NDIAG, SA, LDA, SSCALE, SSCOND, SAMAX, INFO)
CALL ZPBEQU	(UPLO, N, NDIAG, ZA, LDA, DSCALE, DSCOND, DAMAX, INFO)
CALL CPBEQU	(UPLO, N, NDIAG, CA, LDA, SSCALE, SSCOND, SAMAX, INFO)

void dpbequ	(char uplo, int n, int ndiag, double da, int lda, double dscale, double dscond, double damax, int *info)
void spbequ	(char uplo, int n, int ndiag, float sa, int lda, float sscale, float sscond, float samax, int *info)
void zpbequ	(char uplo, int n, int ndiag, doublecomplex za, int lda, double dscale, double dscond, double damax, int *info)
void cpbequ	(char uplo, int n, int ndiag, complex ca, int lda, float sscale, float sscond, float samax, int *info)

Arguments

UPLO

Indicates whether xA contains the upper or lower triangle of the matrix. The legal values for UPLO are listed below. Any values not listed below are illegal.

'U' or 'u'

xA contains the upper triangle.

'L' or 'l'

xA contains the lower triangle.

N

Order of the matrix A. N 0.

NDIAG

Number of superdiagonals or subdiagonals of the matrix A. N-1 NDIAG 0 but if N = 0 then NDIAG = 0.

If UPLO = 'U' or 'u', NDIAG is the number of superdiagonals.

If UPLO = 'L' or 'l', NDIAG is the number of subdiagonals.

xA

Upper or lower triangle of the matrix A.

LDA

Leading dimension of the array A as specified in a dimension or type statement. LDA NDIAG + 1.

xSCALE

On exit, the N scale factors for A.

xSCOND

On exit, the ratio of the smallest SCALE(i) to the largest SCALE(i). If SCOND 0.1 and AMAX is not close to overflow or underflow, it is not worth scaling by SCALE.

xAMAX

On exit, the absolute value of the largest element of A. If AMAX is very close to overflow or underflow, A should be scaled regardless of the value of SCOND.

INFO

On exit:

INFO = 0

Subroutine completed normally.

INFO < 0

The ith argument, where i = |INFO|, had an illegal value.

INFO > 0

The ith diagonal entry, where i = INFO, is non-positive.

UPLO	Indicates whether xA contains the upper or lower triangle of the matrix. The legal values for UPLO are listed below. Any values not listed below are illegal.
	'U' or 'u'	xA contains the upper triangle.
	'L' or 'l'	xA contains the lower triangle.
N	Order of the matrix A. N 0.
NDIAG	Number of superdiagonals or subdiagonals of the matrix A. N-1 NDIAG 0 but if N = 0 then NDIAG = 0.
	If UPLO = 'U' or 'u', NDIAG is the number of superdiagonals.
	If UPLO = 'L' or 'l', NDIAG is the number of subdiagonals.
xA	Upper or lower triangle of the matrix A.
LDA	Leading dimension of the array A as specified in a dimension or type statement. LDA NDIAG + 1.
xSCALE	On exit, the N scale factors for A.
xSCOND	On exit, the ratio of the smallest SCALE(i) to the largest SCALE(i). If SCOND 0.1 and AMAX is not close to overflow or underflow, it is not worth scaling by SCALE.
xAMAX	On exit, the absolute value of the largest element of A. If AMAX is very close to overflow or underflow, A should be scaled regardless of the value of SCOND.
INFO	On exit:
	INFO = 0	Subroutine completed normally.
	INFO < 0	The ith argument, where i = \|INFO\|, had an illegal value.
	INFO > 0	The ith diagonal entry, where i = INFO, is non-positive.

Sample Program



PROGRAM TEST
IMPLICIT NONE
C
INTEGER LDA, N, NDIAG
PARAMETER (N = 4)
PARAMETER (NDIAG = 1)
PARAMETER (LDA = NDIAG + 1)
C
DOUBLE PRECISION A(LDA,N), AMAX, SCALE(N), SCOND
INTEGER ICOL, INFO, IROW
C
EXTERNAL DPBEQU
INTRINSIC ABS
C
C Initialize the array A to store in symmetric banded form
C the 4x4 symmetric positive definite matrix A with one
C subdiagonal and one superdiagonal shown below.
C
C 1024 -2 0 0
C A = -2 128 -2 0
C 0 -2 16 -2
C 0 0 -2 1
C
DATA A / 8D8, 1.024D3, -2.0D0, 1.28D2,
$ -2.0D0, 1.6D1, -2.0D0, 1.0D0 /
C
C Print the initial values of the arrays.
C
PRINT 1000
PRINT 1010, A(2,1), A(1,2), 0.0D0, 0.0D0
PRINT 1010, A(1,2), A(2,2), A(1,3), 0.0D0
PRINT 1010, 0.0D0, A(1,3), A(2,3), A(1,4)
PRINT 1010, 0.0D0, 0.0D0, A(1,4), A(2,4)
PRINT 1020
PRINT 1010, ((A(IROW,ICOL), ICOL = 1, N), IROW = 1, LDA)
C
C Compute the scale factors to use to equilibrate A.
C
CALL DPBEQU ('UPPER TRIANGLE OF A STORED', N, NDIAG, A,
$ LDA, SCALE, SCOND, AMAX, INFO)
IF (INFO .NE. 0) THEN
PRINT 1030, INFO
STOP 1
END IF
C
C Apply the scale factors to A then print the result.
C
DO 100, ICOL = 1, N
A(NDIAG+1,ICOL) = A(NDIAG+1,ICOL) * (SCALE(ICOL) ** 2)
100 CONTINUE
DO 110, ICOL = 2, N
A(1,ICOL) = A(1,ICOL) * SCALE(ICOL-1) * SCALE(ICOL)
110 CONTINUE
PRINT 1040
PRINT 1050, SCALE
PRINT 1060
PRINT 1010, A(2,1), A(1,2), 0.0D0, 0.0D0
PRINT 1010, A(1,2), A(2,2), A(1,3), 0.0D0
PRINT 1010, 0.0D0, A(1,3), A(2,3), A(1,4)
PRINT 1010, 0.0D0, 0.0D0, A(1,4), A(2,4)
C
1000 FORMAT (1X, 'A in full form:')
1010 FORMAT (4(3X, F8.3))
1020 FORMAT (/1X, 'A in banded form: (* in unused elements)')
1030 FORMAT (1X, 'Error computing scale factors of A, INFO =',
$ 1X, I5)
1040 FORMAT (/1X, 'Scale factors for A:')
1050 FORMAT (4X, F8.5)
1060 FORMAT (/1X, 'Scaled A:')
C
END

Sample Output



A in full form:
1024.000 -2.000 0.000 0.000
-2.000 128.000 -2.000 0.000
0.000 -2.000 16.000 -2.000
0.000 0.000 -2.000 1.000

A in banded form: (* in unused elements)
******** -2.000 -2.000 -2.000
1024.000 128.000 16.000 1.000

Scale factors for A:
0.03125
0.08839
0.25000
1.00000

Scaled A:
1.000 -0.006 0.000 0.000
-0.006 1.000 -0.044 0.000
0.000 -0.044 1.000 -0.500
0.000 0.000 -0.500 1.000