Previous Next Contents Generated Index Doc Set Home



Solution to a Linear System in a Symmetric Positive Definite Matrix (Simple Driver)

The subroutines described in this section solve a linear system AX = B for a real symmetric (or Hermitian) positive definite matrix A and general matrices B and X. Note that the expert driver xPOSVX is also available.

Calling Sequence

CALL DPOSV 
(UPLO, N, NRHS, DA, LDA, DB, LDB, INFO)
CALL SPOSV 
(UPLO, N, NRHS, SA, LDA, SB, LDB, INFO)
CALL ZPOSV 
(UPLO, N, NRHS, ZA, LDA, ZB, LDB, INFO)
CALL CPOSV 
(UPLO, N, NRHS, CA, LDA, CB, LDB, INFO)






void dposv 
(char uplo, int n, int nrhs, double *da, int lda, 
double *db, int ldb, int *info)
void sposv 
(char uplo, int n, int nrhs, float *sa, int lda, float 
*sb, int ldb, int *info)
void zposv 
(char uplo, int n, int nrhs, doublecomplex *za, int 
lda, doublecomplex *zb, int ldb, int *info)
void cposv 
(char uplo, int n, int nrhs, complex *ca, int lda, 
complex *cb, int ldb, 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.

NRHS

Number of right-hand sides, equal to the number of columns of the matrix B. NRHS 0.

xA

On entry, the upper or lower triangle of the matrix A.
On exit, the Cholesky factorization of the matrix A as computed by xPOTRF.

LDA

Leading dimension of the array A as specified in a dimension or type statement. LDA max(1, N).

xB

On entry, the N×NRHS right-hand side matrix B.
On exit, the N×NRHS solution matrix X.

LDB

Leading dimension of the array B as specified in a dimension or type statement. LDB max(1, N).

INFO

On exit:

INFO = 0

Subroutine completed normally.

INFO < 0

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

INFO > 0

The leading minor of order i of A, where i = INFO, is not positive definite. The factorization has not been completed, and the solution could not be computed.

Sample Program




      PROGRAM TEST
      IMPLICIT NONE
C
      INTEGER           LDA, LDB, N, NRHS
      PARAMETER        (N = 4)
      PARAMETER        (LDA = N)
      PARAMETER        (LDB = N)
      PARAMETER        (NRHS = 1)
C
      DOUBLE PRECISION  A(LDA,N), B(LDB,NRHS)
      INTEGER           ICOL, INFO, IROW
C
      EXTERNAL          DPOSV
C
C     Initialize the array A to store in symmetric form the 4x4
C     symmetric positive definite matrix A shown below. 
C     Initialize the array B to store the right hand side 
C     vector b shown below.
C
C          2  -1   0   0          6
C     A = -1   2  -1   0     b = 12
C          0  -1   2  -1         12
C          0   0  -1   2          6
C
      DATA A /  2.0D0, 3*8D8, -1.0D0, 2.0D0, 2*8D8, 0.0D0, -1.0D0,
     $          2.0D0, 8D8, 0.0D0, 0.0D0, -1.0D0, 2.0D0 /
      DATA B / 6.0D0, 1.2D1, 1.2D1, 6.0D0 /
C
C     Print the initial values of the arrays.
C
      PRINT 1000
      DO 100, IROW = 1, N
        PRINT 1010, (A(ICOL,IROW), ICOL = 1, IROW - 1),
     $              (A(IROW,ICOL), ICOL = IROW, N)
  100 CONTINUE
      PRINT 1020
      PRINT 1010, ((A(IROW,ICOL), ICOL = 1, N), IROW = 1, LDA)
      PRINT 1030
      PRINT 1040, B
C
C     Solve Ax=b and print the solution.
C
      CALL DPOSV ('UPPER TRIANGLE OF A STORED', N, NRHS, A, LDA,
     $            B, LDB, INFO)
      IF (INFO .NE. 0) THEN
        PRINT 1050, INFO
        STOP 1
      END IF
      PRINT 1060
      PRINT 1040, B
C
 1000 FORMAT (1X, 'A in full form:')
 1010 FORMAT (4(3X, F6.3))
 1020 FORMAT (/1X, 'A in symmetric form:  (* in unused elements)')
 1030 FORMAT (/1X, 'b:')
 1040 FORMAT (1X, F6.3)
 1050 FORMAT (1X, 'Error solving Ax=b, INFO = ', I5)
 1060 FORMAT (/1X, 'x:')
C
      END
 

Sample Output

 
 A in full form:
    2.000   -1.000    0.000    0.000
   -1.000    2.000   -1.000    0.000
    0.000   -1.000    2.000   -1.000
    0.000    0.000   -1.000    2.000



 A in symmetric form:  (* in unused elements)
    2.000   -1.000    0.000    0.000
   ******    2.000   -1.000    0.000
   ******   ******    2.000   -1.000
   ******   ******   ******    2.000



 b:
  6.000
 12.000
 12.000
  6.000



 x:
 18.000
 30.000
 30.000
 18.000






Previous Next Contents Generated Index Doc Set Home