NAME
zsico - compute the UDU factorization and condition number
of a symmetric matrix A. If the condition number is not
needed then xSIFA is slightly faster. It is typical to fol-
low a call to xSICO with a call to xSISL to solve Ax = b or
to xSIDI to compute the determinant, inverse, and inertia of
A.
SYNOPSIS
SUBROUTINE DSICO (DA, LDA, N, IPIVOT, DRCOND, DWORK)
SUBROUTINE SSICO (SA, LDA, N, IPIVOT, SRCOND, SWORK)
SUBROUTINE ZSICO (ZA, LDA, N, IPIVOT, DRCOND, ZWORK)
SUBROUTINE CSICO (CA, LDA, N, IPIVOT, SRCOND, CWORK)
#include <sunperf.h>
void dsico(double *da, int lda, int n, int *kpvt, double
*rcond) ;
void ssico(float *sa, int lda, int n, int *kpvt, float
*rcond) ;
void zsico(doublecomplex *za, int lda, int n, int *kpvt,
double *rcond) ;
void csico(complex *ca, int lda, int n, int *kpvt, float
*rcond) ;
ARGUMENTS
xA On entry, the upper triangle of the matrix A. On
exit, a UDU factorization of the matrix A. The
strict lower triangle of A is not referenced.
LDA Leading dimension of the array A as specified in a
dimension or type statement. LDA >= max(1,N).
N Order of the matrix A. N >= 0.
IPIVOT On exit, a vector of pivot indices.
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.
SAMPLE PROGRAM
PROGRAM TEST
IMPLICIT NONE
C
INTEGER LDA, N
PARAMETER (N = 4)
PARAMETER (LDA = 5)
C
DOUBLE PRECISION A(LDA,N), B(N), RCOND, WORK(N)
INTEGER ICOL, IPIVOT(N), IROW
C
EXTERNAL DSICO, DSISL
C
C Initialize the array A to store the matrix A shown below.
C Initialize the array B to store the vector b shown below.
C
C -.5 -.5 -.5 -.5 12
C A = -.5 -1.5 -1.5 -1.5 b = 6
C -.5 -1.5 -2.5 -2.5 6
C -.5 -1.5 -2.5 -3.5 12
C
DATA A / -5.0D-1, 4*8D8, -5.0D-1, -1.5D0, 3*8D8, -5.0D-1,
$ -1.5D0, -2.5D0, 2*8D8, -5.0D-1, -1.5D0, -2.5D0,
$ -3.5D0, 8D8 /
DATA B / 1.2D1, 6.0D0, 6.0D0, 1.2D1 /
C
PRINT 1000
DO 100, IROW = 1, N
PRINT 1010, (A(ICOL,IROW), ICOL = 1, IROW),
$ (A(IROW,ICOL), ICOL = IROW + 1, N)
100 CONTINUE
PRINT 1020
PRINT 1010, ((A(IROW,ICOL), ICOL = 1, N), IROW = 1, N)
PRINT 1030
PRINT 1040, B
CALL DSICO (A, LDA, N, IPIVOT, RCOND, WORK)
PRINT 1050, RCOND
IF ((RCOND + 1.0D0) .EQ. RCOND) THEN
PRINT 1060
END IF
PRINT 1070, 1.0D0 / RCOND
CALL DSISL (A, LDA, N, IPIVOT, B)
PRINT 1080
PRINT 1040, B
C
1000 FORMAT (1X, 'A in full form:')
1010 FORMAT (4(3X, F5.1))
1020 FORMAT (/1X, 'A in symmetric form: (* in unused elements)')
1030 FORMAT (/1X, 'b:')
1040 FORMAT (3X, F5.1)
1050 FORMAT (/1X, 'Reciprocal condition number of A:', F6.3)
1060 FORMAT (1X, 'A may be singular to working precision.')
1070 FORMAT (1X, 'Condition number of A: ', F6.3)
1080 FORMAT (/1X, 'A**(-1) * b:')
C
END
SAMPLE OUTPUT
A in full form:
-0.5 -0.5 -0.5 -0.5
-0.5 -1.5 -1.5 -1.5
-0.5 -1.5 -2.5 -2.5
-0.5 -1.5 -2.5 -3.5
A in symmetric form: (* in unused elements)
-0.5 -0.5 -0.5 -0.5
***** -1.5 -1.5 -1.5
***** ***** -2.5 -2.5
***** ***** ***** -3.5
b:
12.0
6.0
6.0
12.0
Reciprocal condition number of A: 0.031
Condition number of A: 32.000
A**(-1) * b:
-30.0
6.0
6.0
-6.0
|
Закладки на сайте Проследить за страницей |
Created 1996-2025 by Maxim Chirkov Добавить, Поддержать, Вебмастеру |