001:       SUBROUTINE CSYSV( UPLO, N, NRHS, A, LDA, IPIV, B, LDB, WORK,
002:      $                  LWORK, INFO )
003: *
004: *  -- LAPACK driver routine (version 3.2) --
005: *  -- LAPACK is a software package provided by Univ. of Tennessee,    --
006: *  -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
007: *     November 2006
008: *
009: *     .. Scalar Arguments ..
010:       CHARACTER          UPLO
011:       INTEGER            INFO, LDA, LDB, LWORK, N, NRHS
012: *     ..
013: *     .. Array Arguments ..
014:       INTEGER            IPIV( * )
015:       COMPLEX            A( LDA, * ), B( LDB, * ), WORK( * )
016: *     ..
017: *
018: *  Purpose
019: *  =======
020: *
021: *  CSYSV computes the solution to a complex system of linear equations
022: *     A * X = B,
023: *  where A is an N-by-N symmetric matrix and X and B are N-by-NRHS
024: *  matrices.
025: *
026: *  The diagonal pivoting method is used to factor A as
027: *     A = U * D * U**T,  if UPLO = 'U', or
028: *     A = L * D * L**T,  if UPLO = 'L',
029: *  where U (or L) is a product of permutation and unit upper (lower)
030: *  triangular matrices, and D is symmetric and block diagonal with 
031: *  1-by-1 and 2-by-2 diagonal blocks.  The factored form of A is then
032: *  used to solve the system of equations A * X = B.
033: *
034: *  Arguments
035: *  =========
036: *
037: *  UPLO    (input) CHARACTER*1
038: *          = 'U':  Upper triangle of A is stored;
039: *          = 'L':  Lower triangle of A is stored.
040: *
041: *  N       (input) INTEGER
042: *          The number of linear equations, i.e., the order of the
043: *          matrix A.  N >= 0.
044: *
045: *  NRHS    (input) INTEGER
046: *          The number of right hand sides, i.e., the number of columns
047: *          of the matrix B.  NRHS >= 0.
048: *
049: *  A       (input/output) COMPLEX array, dimension (LDA,N)
050: *          On entry, the symmetric matrix A.  If UPLO = 'U', the leading
051: *          N-by-N upper triangular part of A contains the upper
052: *          triangular part of the matrix A, and the strictly lower
053: *          triangular part of A is not referenced.  If UPLO = 'L', the
054: *          leading N-by-N lower triangular part of A contains the lower
055: *          triangular part of the matrix A, and the strictly upper
056: *          triangular part of A is not referenced.
057: *
058: *          On exit, if INFO = 0, the block diagonal matrix D and the
059: *          multipliers used to obtain the factor U or L from the
060: *          factorization A = U*D*U**T or A = L*D*L**T as computed by
061: *          CSYTRF.
062: *
063: *  LDA     (input) INTEGER
064: *          The leading dimension of the array A.  LDA >= max(1,N).
065: *
066: *  IPIV    (output) INTEGER array, dimension (N)
067: *          Details of the interchanges and the block structure of D, as
068: *          determined by CSYTRF.  If IPIV(k) > 0, then rows and columns
069: *          k and IPIV(k) were interchanged, and D(k,k) is a 1-by-1
070: *          diagonal block.  If UPLO = 'U' and IPIV(k) = IPIV(k-1) < 0,
071: *          then rows and columns k-1 and -IPIV(k) were interchanged and
072: *          D(k-1:k,k-1:k) is a 2-by-2 diagonal block.  If UPLO = 'L' and
073: *          IPIV(k) = IPIV(k+1) < 0, then rows and columns k+1 and
074: *          -IPIV(k) were interchanged and D(k:k+1,k:k+1) is a 2-by-2
075: *          diagonal block.
076: *
077: *  B       (input/output) COMPLEX array, dimension (LDB,NRHS)
078: *          On entry, the N-by-NRHS right hand side matrix B.
079: *          On exit, if INFO = 0, the N-by-NRHS solution matrix X.
080: *
081: *  LDB     (input) INTEGER
082: *          The leading dimension of the array B.  LDB >= max(1,N).
083: *
084: *  WORK    (workspace/output) COMPLEX array, dimension (MAX(1,LWORK))
085: *          On exit, if INFO = 0, WORK(1) returns the optimal LWORK.
086: *
087: *  LWORK   (input) INTEGER
088: *          The length of WORK.  LWORK >= 1, and for best performance
089: *          LWORK >= max(1,N*NB), where NB is the optimal blocksize for
090: *          CSYTRF.
091: *
092: *          If LWORK = -1, then a workspace query is assumed; the routine
093: *          only calculates the optimal size of the WORK array, returns
094: *          this value as the first entry of the WORK array, and no error
095: *          message related to LWORK is issued by XERBLA.
096: *
097: *  INFO    (output) INTEGER
098: *          = 0: successful exit
099: *          < 0: if INFO = -i, the i-th argument had an illegal value
100: *          > 0: if INFO = i, D(i,i) is exactly zero.  The factorization
101: *               has been completed, but the block diagonal matrix D is
102: *               exactly singular, so the solution could not be computed.
103: *
104: *  =====================================================================
105: *
106: *     .. Local Scalars ..
107:       LOGICAL            LQUERY
108:       INTEGER            LWKOPT, NB
109: *     ..
110: *     .. External Functions ..
111:       LOGICAL            LSAME
112:       INTEGER            ILAENV
113:       EXTERNAL           ILAENV, LSAME
114: *     ..
115: *     .. External Subroutines ..
116:       EXTERNAL           CSYTRF, CSYTRS, XERBLA
117: *     ..
118: *     .. Intrinsic Functions ..
119:       INTRINSIC          MAX
120: *     ..
121: *     .. Executable Statements ..
122: *
123: *     Test the input parameters.
124: *
125:       INFO = 0
126:       LQUERY = ( LWORK.EQ.-1 )
127:       IF( .NOT.LSAME( UPLO, 'U' ) .AND. .NOT.LSAME( UPLO, 'L' ) ) THEN
128:          INFO = -1
129:       ELSE IF( N.LT.0 ) THEN
130:          INFO = -2
131:       ELSE IF( NRHS.LT.0 ) THEN
132:          INFO = -3
133:       ELSE IF( LDA.LT.MAX( 1, N ) ) THEN
134:          INFO = -5
135:       ELSE IF( LDB.LT.MAX( 1, N ) ) THEN
136:          INFO = -8
137:       ELSE IF( LWORK.LT.1 .AND. .NOT.LQUERY ) THEN
138:          INFO = -10
139:       END IF
140: *
141:       IF( INFO.EQ.0 ) THEN
142:          IF( N.EQ.0 ) THEN
143:             LWKOPT = 1
144:          ELSE
145:             NB = ILAENV( 1, 'CSYTRF', UPLO, N, -1, -1, -1 )
146:             LWKOPT = N*NB
147:          END IF
148:          WORK( 1 ) = LWKOPT
149:       END IF
150: *
151:       IF( INFO.NE.0 ) THEN
152:          CALL XERBLA( 'CSYSV ', -INFO )
153:          RETURN
154:       ELSE IF( LQUERY ) THEN
155:          RETURN
156:       END IF
157: *
158: *     Compute the factorization A = U*D*U' or A = L*D*L'.
159: *
160:       CALL CSYTRF( UPLO, N, A, LDA, IPIV, WORK, LWORK, INFO )
161:       IF( INFO.EQ.0 ) THEN
162: *
163: *        Solve the system A*X = B, overwriting B with X.
164: *
165:          CALL CSYTRS( UPLO, N, NRHS, A, LDA, IPIV, B, LDB, INFO )
166: *
167:       END IF
168: *
169:       WORK( 1 ) = LWKOPT
170: *
171:       RETURN
172: *
173: *     End of CSYSV
174: *
175:       END
176: