NAME
zpbco - compute a Cholesky factorization and condition
number of a symmetric positive definite matrix A in banded
storage. If the condition number is not needed then xPBFA
is slightly faster. It is typical to follow a call to xPBCO
with a call to xPBSL to solve Ax = b or to xPBDI to compute
the determinant of A.
SYNOPSIS
SUBROUTINE DPBCO (DA, LDA, N, NDIAG, DRCOND, DWORK, INFO)
SUBROUTINE SPBCO (SA, LDA, N, NDIAG, SRCOND, SWORK, INFO)
SUBROUTINE ZPBCO (ZA, LDA, N, NDIAG, DRCOND, ZWORK, INFO)
SUBROUTINE CPBCO (CA, LDA, N, NDIAG, SRCOND, CWORK, INFO)
#include <sunperf.h>
void dpbco(double *abd, int lda, int n, int m, double
*drcond, int *info) ;
void spbco(float *abd, int lda, int n, int m, float *srcond,
int *info) ;
void zpbco(doublecomplex *abd, int lda, int n, int m, double
*drcond, int *info) ;
void cpbco(complex *abd, int lda, int n, int m, float
*srcond, int *info) ;
ARGUMENTS
xA On entry, the upper triangle of the matrix A. On
exit, a Cholesky factorization of the matrix A.
LDA Leading dimension of the array A as specified in a
dimension or type statement. LDA >= NDIAG + 1.
N Order of the matrix A. N * 0.
NDIAG Number of diagonals. N-1 >= NDIAG >= 0 but if N =
0 then NDIAG = 0.
xRCOND On exit, an estimate of the reciprocal condition
number of A. 0.0 <= RCOND <= 1.0. As the value
of RCOND gets smaller, operations with A such as
solving Ax = b may become less stable. If RCOND
satisfies RCOND + 1.0 = 1.0 then A may be singular
to working precision.
xWORK Scratch array with a dimension of N.
INFO On exit:
INFO = 0 Subroutine completed normally.
INFO * 0 Returns a value k if the leading minor
of order k is not positive definite.
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), B(N), RCOND, WORK(N)
INTEGER ICOL, INFO, IROW
C
EXTERNAL DPBCO, DPBSL
C
C Initialize the array A to store in banded storage mode
C the matrix A shown below. Initialize the array B to
C store the vector B shown below.
C
C 2 -1 0 0 60
C A = -1 2 -1 0 b = 60
C 0 -1 2 -1 60
C 0 0 -1 2 60
C
DATA A / 8D8, 2.0D0, -1.0D0, 2.0D0, -1.0D0, 2.0D0, -1.0D0, 2.0D0 /
DATA B / 4*6.0D1 /
C
PRINT 1000
PRINT 1010, A(2,1), A(1,2)
PRINT 1010, A(3,1), A(2,2), A(1,3)
PRINT 1020, A(3,2), A(2,3), A(1,4)
PRINT 1030, A(3,3), A(2,4)
PRINT 1040
PRINT 1010, ((A(IROW,ICOL), ICOL = 1, N), IROW = 1, LDA)
PRINT 1050
PRINT 1060, B
CALL DPBCO (A, LDA, N, NDIAG, RCOND, WORK, INFO)
IF (INFO .EQ. 0) THEN
IF ((RCOND + 1.0D0) .EQ. 1.0D0) THEN
PRINT 1100
END IF
CALL DPBSL (A, LDA, N, NDIAG, B)
PRINT 1070, RCOND
PRINT 1080
PRINT 1060, B
ELSE
PRINT 1090
END IF
C
1000 FORMAT (1X, 'A in full form:')
1010 FORMAT (4(3X, F5.1))
1020 FORMAT (8X, 3(3X, F5.1))
1030 FORMAT (16X, 3(3X, F5.1))
1040 FORMAT (/1X, 'A in banded form: (* in unused entries)')
1050 FORMAT (/1X, 'b:')
1060 FORMAT (3X, F5.1)
1070 FORMAT (/1X, 'Reciprocal condition number of A:', F5.1)
1080 FORMAT (/1X, 'A**(-1) * b:')
1090 FORMAT (/1X, 'A is not positive definite.')
1100 FORMAT (1X, 'A may be singular to working precision.')
C
END
SAMPLE OUTPUT
A in full form:
2.0 -1.0
-1.0 2.0 -1.0
-1.0 2.0 -1.0
-1.0 2.0
A in banded form: (* in unused entries)
***** -1.0 -1.0 -1.0
2.0 2.0 2.0 2.0
b:
60.0
60.0
60.0
60.0
Reciprocal condition number of A: 0.1
A**(-1) * b:
120.0
180.0
180.0
120.0
|
Закладки на сайте Проследить за страницей |
Created 1996-2025 by Maxim Chirkov Добавить, Поддержать, Вебмастеру |