001:       SUBROUTINE ZHPSVX( FACT, UPLO, N, NRHS, AP, AFP, IPIV, B, LDB, X,
002:      $                   LDX, RCOND, FERR, BERR, WORK, RWORK, 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          FACT, UPLO
011:       INTEGER            INFO, LDB, LDX, N, NRHS
012:       DOUBLE PRECISION   RCOND
013: *     ..
014: *     .. Array Arguments ..
015:       INTEGER            IPIV( * )
016:       DOUBLE PRECISION   BERR( * ), FERR( * ), RWORK( * )
017:       COMPLEX*16         AFP( * ), AP( * ), B( LDB, * ), WORK( * ),
018:      $                   X( LDX, * )
019: *     ..
020: *
021: *  Purpose
022: *  =======
023: *
024: *  ZHPSVX uses the diagonal pivoting factorization A = U*D*U**H or
025: *  A = L*D*L**H to compute the solution to a complex system of linear
026: *  equations A * X = B, where A is an N-by-N Hermitian matrix stored
027: *  in packed format and X and B are N-by-NRHS matrices.
028: *
029: *  Error bounds on the solution and a condition estimate are also
030: *  provided.
031: *
032: *  Description
033: *  ===========
034: *
035: *  The following steps are performed:
036: *
037: *  1. If FACT = 'N', the diagonal pivoting method is used to factor A as
038: *        A = U * D * U**H,  if UPLO = 'U', or
039: *        A = L * D * L**H,  if UPLO = 'L',
040: *     where U (or L) is a product of permutation and unit upper (lower)
041: *     triangular matrices and D is Hermitian and block diagonal with
042: *     1-by-1 and 2-by-2 diagonal blocks.
043: *
044: *  2. If some D(i,i)=0, so that D is exactly singular, then the routine
045: *     returns with INFO = i. Otherwise, the factored form of A is used
046: *     to estimate the condition number of the matrix A.  If the
047: *     reciprocal of the condition number is less than machine precision,
048: *     INFO = N+1 is returned as a warning, but the routine still goes on
049: *     to solve for X and compute error bounds as described below.
050: *
051: *  3. The system of equations is solved for X using the factored form
052: *     of A.
053: *
054: *  4. Iterative refinement is applied to improve the computed solution
055: *     matrix and calculate error bounds and backward error estimates
056: *     for it.
057: *
058: *  Arguments
059: *  =========
060: *
061: *  FACT    (input) CHARACTER*1
062: *          Specifies whether or not the factored form of A has been
063: *          supplied on entry.
064: *          = 'F':  On entry, AFP and IPIV contain the factored form of
065: *                  A.  AFP and IPIV will not be modified.
066: *          = 'N':  The matrix A will be copied to AFP and factored.
067: *
068: *  UPLO    (input) CHARACTER*1
069: *          = 'U':  Upper triangle of A is stored;
070: *          = 'L':  Lower triangle of A is stored.
071: *
072: *  N       (input) INTEGER
073: *          The number of linear equations, i.e., the order of the
074: *          matrix A.  N >= 0.
075: *
076: *  NRHS    (input) INTEGER
077: *          The number of right hand sides, i.e., the number of columns
078: *          of the matrices B and X.  NRHS >= 0.
079: *
080: *  AP      (input) COMPLEX*16 array, dimension (N*(N+1)/2)
081: *          The upper or lower triangle of the Hermitian matrix A, packed
082: *          columnwise in a linear array.  The j-th column of A is stored
083: *          in the array AP as follows:
084: *          if UPLO = 'U', AP(i + (j-1)*j/2) = A(i,j) for 1<=i<=j;
085: *          if UPLO = 'L', AP(i + (j-1)*(2*n-j)/2) = A(i,j) for j<=i<=n.
086: *          See below for further details.
087: *
088: *  AFP     (input or output) COMPLEX*16 array, dimension (N*(N+1)/2)
089: *          If FACT = 'F', then AFP is an input argument and on entry
090: *          contains the block diagonal matrix D and the multipliers used
091: *          to obtain the factor U or L from the factorization
092: *          A = U*D*U**H or A = L*D*L**H as computed by ZHPTRF, stored as
093: *          a packed triangular matrix in the same storage format as A.
094: *
095: *          If FACT = 'N', then AFP is an output argument and on exit
096: *          contains the block diagonal matrix D and the multipliers used
097: *          to obtain the factor U or L from the factorization
098: *          A = U*D*U**H or A = L*D*L**H as computed by ZHPTRF, stored as
099: *          a packed triangular matrix in the same storage format as A.
100: *
101: *  IPIV    (input or output) INTEGER array, dimension (N)
102: *          If FACT = 'F', then IPIV is an input argument and on entry
103: *          contains details of the interchanges and the block structure
104: *          of D, as determined by ZHPTRF.
105: *          If IPIV(k) > 0, then rows and columns k and IPIV(k) were
106: *          interchanged and D(k,k) is a 1-by-1 diagonal block.
107: *          If UPLO = 'U' and IPIV(k) = IPIV(k-1) < 0, then rows and
108: *          columns k-1 and -IPIV(k) were interchanged and D(k-1:k,k-1:k)
109: *          is a 2-by-2 diagonal block.  If UPLO = 'L' and IPIV(k) =
110: *          IPIV(k+1) < 0, then rows and columns k+1 and -IPIV(k) were
111: *          interchanged and D(k:k+1,k:k+1) is a 2-by-2 diagonal block.
112: *
113: *          If FACT = 'N', then IPIV is an output argument and on exit
114: *          contains details of the interchanges and the block structure
115: *          of D, as determined by ZHPTRF.
116: *
117: *  B       (input) COMPLEX*16 array, dimension (LDB,NRHS)
118: *          The N-by-NRHS right hand side matrix B.
119: *
120: *  LDB     (input) INTEGER
121: *          The leading dimension of the array B.  LDB >= max(1,N).
122: *
123: *  X       (output) COMPLEX*16 array, dimension (LDX,NRHS)
124: *          If INFO = 0 or INFO = N+1, the N-by-NRHS solution matrix X.
125: *
126: *  LDX     (input) INTEGER
127: *          The leading dimension of the array X.  LDX >= max(1,N).
128: *
129: *  RCOND   (output) DOUBLE PRECISION
130: *          The estimate of the reciprocal condition number of the matrix
131: *          A.  If RCOND is less than the machine precision (in
132: *          particular, if RCOND = 0), the matrix is singular to working
133: *          precision.  This condition is indicated by a return code of
134: *          INFO > 0.
135: *
136: *  FERR    (output) DOUBLE PRECISION array, dimension (NRHS)
137: *          The estimated forward error bound for each solution vector
138: *          X(j) (the j-th column of the solution matrix X).
139: *          If XTRUE is the true solution corresponding to X(j), FERR(j)
140: *          is an estimated upper bound for the magnitude of the largest
141: *          element in (X(j) - XTRUE) divided by the magnitude of the
142: *          largest element in X(j).  The estimate is as reliable as
143: *          the estimate for RCOND, and is almost always a slight
144: *          overestimate of the true error.
145: *
146: *  BERR    (output) DOUBLE PRECISION array, dimension (NRHS)
147: *          The componentwise relative backward error of each solution
148: *          vector X(j) (i.e., the smallest relative change in
149: *          any element of A or B that makes X(j) an exact solution).
150: *
151: *  WORK    (workspace) COMPLEX*16 array, dimension (2*N)
152: *
153: *  RWORK   (workspace) DOUBLE PRECISION array, dimension (N)
154: *
155: *  INFO    (output) INTEGER
156: *          = 0: successful exit
157: *          < 0: if INFO = -i, the i-th argument had an illegal value
158: *          > 0:  if INFO = i, and i is
159: *                <= N:  D(i,i) is exactly zero.  The factorization
160: *                       has been completed but the factor D is exactly
161: *                       singular, so the solution and error bounds could
162: *                       not be computed. RCOND = 0 is returned.
163: *                = N+1: D is nonsingular, but RCOND is less than machine
164: *                       precision, meaning that the matrix is singular
165: *                       to working precision.  Nevertheless, the
166: *                       solution and error bounds are computed because
167: *                       there are a number of situations where the
168: *                       computed solution can be more accurate than the
169: *                       value of RCOND would suggest.
170: *
171: *  Further Details
172: *  ===============
173: *
174: *  The packed storage scheme is illustrated by the following example
175: *  when N = 4, UPLO = 'U':
176: *
177: *  Two-dimensional storage of the Hermitian matrix A:
178: *
179: *     a11 a12 a13 a14
180: *         a22 a23 a24
181: *             a33 a34     (aij = conjg(aji))
182: *                 a44
183: *
184: *  Packed storage of the upper triangle of A:
185: *
186: *  AP = [ a11, a12, a22, a13, a23, a33, a14, a24, a34, a44 ]
187: *
188: *  =====================================================================
189: *
190: *     .. Parameters ..
191:       DOUBLE PRECISION   ZERO
192:       PARAMETER          ( ZERO = 0.0D+0 )
193: *     ..
194: *     .. Local Scalars ..
195:       LOGICAL            NOFACT
196:       DOUBLE PRECISION   ANORM
197: *     ..
198: *     .. External Functions ..
199:       LOGICAL            LSAME
200:       DOUBLE PRECISION   DLAMCH, ZLANHP
201:       EXTERNAL           LSAME, DLAMCH, ZLANHP
202: *     ..
203: *     .. External Subroutines ..
204:       EXTERNAL           XERBLA, ZCOPY, ZHPCON, ZHPRFS, ZHPTRF, ZHPTRS,
205:      $                   ZLACPY
206: *     ..
207: *     .. Intrinsic Functions ..
208:       INTRINSIC          MAX
209: *     ..
210: *     .. Executable Statements ..
211: *
212: *     Test the input parameters.
213: *
214:       INFO = 0
215:       NOFACT = LSAME( FACT, 'N' )
216:       IF( .NOT.NOFACT .AND. .NOT.LSAME( FACT, 'F' ) ) THEN
217:          INFO = -1
218:       ELSE IF( .NOT.LSAME( UPLO, 'U' ) .AND. .NOT.LSAME( UPLO, 'L' ) )
219:      $          THEN
220:          INFO = -2
221:       ELSE IF( N.LT.0 ) THEN
222:          INFO = -3
223:       ELSE IF( NRHS.LT.0 ) THEN
224:          INFO = -4
225:       ELSE IF( LDB.LT.MAX( 1, N ) ) THEN
226:          INFO = -9
227:       ELSE IF( LDX.LT.MAX( 1, N ) ) THEN
228:          INFO = -11
229:       END IF
230:       IF( INFO.NE.0 ) THEN
231:          CALL XERBLA( 'ZHPSVX', -INFO )
232:          RETURN
233:       END IF
234: *
235:       IF( NOFACT ) THEN
236: *
237: *        Compute the factorization A = U*D*U' or A = L*D*L'.
238: *
239:          CALL ZCOPY( N*( N+1 ) / 2, AP, 1, AFP, 1 )
240:          CALL ZHPTRF( UPLO, N, AFP, IPIV, INFO )
241: *
242: *        Return if INFO is non-zero.
243: *
244:          IF( INFO.GT.0 )THEN
245:             RCOND = ZERO
246:             RETURN
247:          END IF
248:       END IF
249: *
250: *     Compute the norm of the matrix A.
251: *
252:       ANORM = ZLANHP( 'I', UPLO, N, AP, RWORK )
253: *
254: *     Compute the reciprocal of the condition number of A.
255: *
256:       CALL ZHPCON( UPLO, N, AFP, IPIV, ANORM, RCOND, WORK, INFO )
257: *
258: *     Compute the solution vectors X.
259: *
260:       CALL ZLACPY( 'Full', N, NRHS, B, LDB, X, LDX )
261:       CALL ZHPTRS( UPLO, N, NRHS, AFP, IPIV, X, LDX, INFO )
262: *
263: *     Use iterative refinement to improve the computed solutions and
264: *     compute error bounds and backward error estimates for them.
265: *
266:       CALL ZHPRFS( UPLO, N, NRHS, AP, AFP, IPIV, B, LDB, X, LDX, FERR,
267:      $             BERR, WORK, RWORK, INFO )
268: *
269: *     Set INFO = N+1 if the matrix is singular to working precision.
270: *
271:       IF( RCOND.LT.DLAMCH( 'Epsilon' ) )
272:      $   INFO = N + 1
273: *
274:       RETURN
275: *
276: *     End of ZHPSVX
277: *
278:       END
279: